Merge pull request #42 from probsys/20251217-fsaad-sop-times

20251217 fsaad sop times

Author: fsaad

docs/src/tutorials/decomposition.ipynb

@@ -60,7 +60,7 @@ ax.scatter(df_test.ds, df_test.y, marker="o", color="w", edgecolor="k", label="T
 
 
 
-    Python: <matplotlib.collections.PathCollection object at 0x7367187bbe30>
+    Python: <matplotlib.collections.PathCollection object at 0x7ddc363cdb80>
 
 
 
@@ -109,7 +109,7 @@ ax.scatter(df_test.ds, df_test.y, marker="o", color="w", edgecolor="k", label="T
 
 
 
-    Python: <matplotlib.collections.PathCollection object at 0x7367183731a0>
+    Python: <matplotlib.collections.PathCollection object at 0x7ddc35e42880>
 
 
 
@@ -628,7 +628,7 @@ for (ax, m, f) in zip(axes, [model_lin, model_per, model_ge], [forecasts_lin, fo
 end
 axes[0].set_title("STRUCTURE: LIN")
 axes[1].set_title("STRUCTURE: PER")
-axes[2].set_title("STRUCTURE: GE")
+axes[2].set_title("STRUCTURE: GE");
 ```
 
 
@@ -637,12 +637,6 @@ axes[2].set_title("STRUCTURE: GE")
 
 
 
-
-
-
-
-
-
 ## Sum-of-Products Decomposition
 
 An third approach to decomposing kernels is using [`AutoGP.split_kernel_sop`](@ref), which is based on a sum-of-products decomposition of kernels.
@@ -724,18 +718,20 @@ nrows = AutoGP.num_particles(model)
 fig, axes = PythonPlot.subplots(ncols=3, nrows=nrows, figsize=(26,5*nrows))
 for particle=1:AutoGP.num_particles(model)
     f = forecasts_sum[forecasts_sum.particle .== particle,:]
-    for (component, ax) in zip(0:2, axes[particle-1])
+    for (component, ax) in zip([0,2,1], axes[particle-1])
         subdf = f[f.component .== component,:]
         ax.plot(subdf[!,"ds"], subdf[!,"y_mean"], color="k", linewidth=1)
         ax.fill_between(subdf.ds, subdf[!,"y_0.025"], subdf[!,"y_0.975"]; color="tab:blue", alpha=0.5)
-        ax.scatter(df_train.ds, df_train.y, marker="o", color="k", label="Observed Data")
-        ax.scatter(df_test.ds, df_test.y, marker="o", color="w", edgecolor="k", label="Test Data")
+        if component in [0,2]
+            ax.scatter(df_train.ds, df_train.y, marker="o", color="k", label="Observed Data")
+            ax.scatter(df_test.ds, df_test.y, marker="o", color="w", edgecolor="k", label="Test Data")
+            ax.set_ylim([5000, 7500])
+        end
         w = @sprintf("%1.2f", maximum(subdf.weight))
         ax.set_title("Particle $(particle) - Component $(title[component+1]) - Weight $(w)")
-        ax.set_ylim([5000, 7500])
     end
 end
-fig.savefig("forecasts-joint.png", dpi=150)
+fig.savefig("forecasts-joint.png", dpi=150);
 ```
 
 
@@ -744,27 +740,24 @@ fig.savefig("forecasts-joint.png", dpi=150)
 
 
 
-
-
-
-
-
-
 Using [`AutoGP.GP.predict_mvn_sum`](@ref) gives the overall mixture of Gaussians.
 
 
 ```julia
 ds_probe = vcat(df_train.ds, df_test.ds)
 mvn, indexes = AutoGP.predict_mvn_sum(model, ds_probe, AutoGP.GP.Periodic);
 
-overall_predictions = mean(mvn)[indexes.Y]
-seasonal_predictions = mean(mvn)[indexes.F[1]]
-nonseasonal_predictions = mean(mvn)[indexes.F[2]]
-fig, ax = PythonPlot.subplots(figsize=(10,6))
-ax.plot(ds_probe, seasonal_predictions, marker=".", color="k")
-ax.plot(ds_probe, nonseasonal_predictions, marker=".", color="k")
-ax.scatter(df_train.ds, df_train.y, marker=".")
-ax.scatter(df_test.ds, df_test.y, marker="o", edgecolor="k", facecolor="w")
+overall_predictions = mean(mvn.components[1])[indexes.Y]
+seasonal_predictions = mean(mvn.components[1])[indexes.F[1]]
+nonseasonal_predictions = mean(mvn.components[1])[indexes.F[2]]
+
+fig, ax = PythonPlot.subplots(figsize=(18,6), ncols=2)
+ax[0].plot(ds_probe, nonseasonal_predictions, marker=".", color="k")
+ax[0].plot(ds_probe, nonseasonal_predictions .+ seasonal_predictions , marker=".", color="k")
+ax[0].scatter(df_train.ds, df_train.y, marker=".")
+ax[0].scatter(df_test.ds, df_test.y, marker="o", edgecolor="k", facecolor="w");
+
+ax[1].plot(ds_probe, seasonal_predictions, marker=".", color="k")
 ```
 
 
@@ -776,5 +769,11 @@ ax.scatter(df_test.ds, df_test.y, marker="o", edgecolor="k", facecolor="w")
 
 
 
+    Python: [<matplotlib.lines.Line2D object at 0x7ddb977b57c0>]
+
 
 
+
+```julia
+
+```

docs/src/tutorials/decomposition.md

@@ -596,6 +596,8 @@ julia> split_kernel_sop(l*p + l*c, Periodic)
 (l*p, l*c)
 julia> split_kernel_sop(p*p, Periodic)
 (p*p, Constant(0))
+julia> split_kernel_sop((l+p)*(l+p), Periodic)
+(p*l+p*p+l*p, l*l)
 ```
 """
 function split_kernel_sop(node::Node, ::Type{T}) where {T<:LeafNode}
@@ -613,7 +615,19 @@ split_kernel_sop_helper(node::T, ::Type{T}) where {T<:LeafNode} = (node, nothing
 split_kernel_sop_helper(node::LeafNode, ::Type{T}) where {T<:LeafNode} = (nothing, node)
 
 function split_kernel_sop_helper(node::Times, ::Type{T}) where {T<:LeafNode}
-    return (has_leaf(node.left, T) || has_leaf(node.right, T)) ? (node, nothing) : (nothing, node)
+    (left_a, left_b) = split_kernel_sop_helper(node.left, T)
+    (right_a, right_b) = split_kernel_sop_helper(node.right, T)
+    mult(a, b) = (isnothing(a) || isnothing(b)) ? nothing : a * b
+    terms = [
+        mult(left_a, right_a),
+        mult(left_a, right_b),
+        mult(left_b, right_a),
+        mult(left_b, right_b),
+    ]
+    l_sop = merge_split_operand(node.left + node.right, terms[1], terms[2])
+    l_sop = merge_split_operand(node.left + node.right, l_sop, terms[3])
+    r_sop = terms[4]
+    return (l_sop, r_sop)
 end
 
 function split_kernel_sop_helper(node::B, ::Type{T}) where {B<:BinaryOpNode, T <: LeafNode}
@@ -625,7 +639,7 @@ function split_kernel_sop_helper(node::B, ::Type{T}) where {B<:BinaryOpNode, T <
 end
 
 # Helper function for split_kernel_sop_helper
-merge_split_operand(node::Plus, node_a, node_b) = @match (node_a, node_b) begin
+merge_split_operand(::Plus, node_a, node_b) = @match (node_a, node_b) begin
     (::Nothing, ::Nothing) => nothing
     (::Node, ::Nothing)    => node_a
     (::Nothing, ::Node)    => node_b

docs/src/tutorials/decomposition_files/decomposition_12_0.png

@@ -1006,11 +1006,17 @@ function predict_mvn_sum(
             noise_pred=noise_pred)
         mu = Distributions.mean(mvn.mvn)
         cov = Distributions.cov(mvn.mvn)
-        # TODO: We must only include the mean in at most ONE of the two
-        # components. For now, the mean is included in all three.
+        # Transform to data space.
         mu, cov = Transforms.unapply_mean_var(model.y_transform, mu, cov)
-        # Remove the mean from the first latent function.
-        # mu[mvn.indexes.F[1]] .+= model.y_transform.intercept / model.y_transform.slope
+        # Avoid double counting the linear offset by removing it from the first
+        # component. Derivation, where Y is data space and X is model space:
+        #   X = aY + b
+        #   X = F1 + F2
+        #   Y = (X-b)/a = F1/a + F2/a -b/a
+        # However, using unapply_mean gives
+        #   Y' = (F1-b)/a + (F2-b)/a = F1/a - b/a + F2/a - b/a
+        # Therefore Y' is double counting by b/a, which we add to F1.
+        mu[mvn.indexes.F[1]] .+= model.y_transform.intercept / model.y_transform.slope
         distributions[particle] = MvNormal(mu, cov)
         # Set the indexes.
         if !isnothing(indexes[1])