bug fix: fe_ls_parabolic edf computation correctly handles n_dofs_

AlePalu · AlePalu · commit 9c6dc1b9af74 · 2025-08-04T16:48:47.000+02:00
diff --git a/fdaPDE/core b/fdaPDE/core
@@ -1 +1 @@
-Subproject commit 7ed083b65046793d4ce58f238f0630011343c401
+Subproject commit a2a9c88b203b6f1b74bbc71d698f131f839f25ba
diff --git a/fdaPDE/src/solvers/fe_ls_parabolic.h b/fdaPDE/src/solvers/fe_ls_parabolic.h
@@ -422,16 +422,16 @@ class fe_ls_parabolic_mono {
                 for (int j = 0; j < r; ++j) { Us_->operator()(i, j) = rademacher(rng); }
             }
             Ys_ = Us_->transpose() * Psi_;
-            Bs_ = matrix_t::Zero(2 * n_dofs_, r);   // implicitly enforce homogeneous forcing
+            Bs_ = matrix_t::Zero(2 * m_ * n_dofs_, r);   // implicitly enforce homogeneous forcing
         }
         if (n_covs_ == 0) {
-            Bs_->topRows(n_dofs_) = -PsiNA().transpose() * D_ * W_ * (*Us_);
+            Bs_->topRows(m_ * n_dofs_) = -PsiNA().transpose() * D_ * W_ * (*Us_);
         } else {
-            Bs_->topRows(n_dofs_) = -PsiNA().transpose() * D_ * internals::lmbQ(W_, X_, invXtWX_, *Us_);
+            Bs_->topRows(m_ * n_dofs_) = -PsiNA().transpose() * D_ * internals::lmbQ(W_, X_, invXtWX_, *Us_);
         }
         matrix_t x = n_covs_ == 0 ? invA_.solve(*Bs_) : woodbury_system_solve(invA_, U_, XtWX_, V_, *Bs_);
         double trS = 0;   // monte carlo Tr[S] approximation
-        for (int i = 0; i < r; ++i) { trS += Ys_->row(i).dot(x.col(i).head(n_dofs_)); }
+        for (int i = 0; i < r; ++i) { trS += Ys_->row(i).dot(x.col(i).head(m_ * n_dofs_)); }
         return trS / r;
     }
     template <typename... LambdaT>

Original file line number	Diff line number	Diff line change
`@@ -422,16 +422,16 @@ class fe_ls_parabolic_mono {`
`422`	`422`	`for (int j = 0; j < r; ++j) { Us_->operator()(i, j) = rademacher(rng); }`
`423`	`423`	`}`
`424`	`424`	`Ys_ = Us_->transpose() * Psi_;`
`425`		`- Bs_ = matrix_t::Zero(2 * n_dofs_, r); // implicitly enforce homogeneous forcing`
	`425`	`+ Bs_ = matrix_t::Zero(2 * m_ * n_dofs_, r); // implicitly enforce homogeneous forcing`
`426`	`426`	`}`
`427`	`427`	`if (n_covs_ == 0) {`
`428`		`- Bs_->topRows(n_dofs_) = -PsiNA().transpose() * D_ * W_ * (*Us_);`
	`428`	`+ Bs_->topRows(m_ * n_dofs_) = -PsiNA().transpose() * D_ * W_ * (*Us_);`
`429`	`429`	`} else {`
`430`		`- Bs_->topRows(n_dofs_) = -PsiNA().transpose() * D_ * internals::lmbQ(W_, X_, invXtWX_, *Us_);`
	`430`	`+ Bs_->topRows(m_ * n_dofs_) = -PsiNA().transpose() * D_ * internals::lmbQ(W_, X_, invXtWX_, *Us_);`
`431`	`431`	`}`
`432`	`432`	`matrix_t x = n_covs_ == 0 ? invA_.solve(Bs_) : woodbury_system_solve(invA_, U_, XtWX_, V_, Bs_);`
`433`	`433`	`double trS = 0; // monte carlo Tr[S] approximation`
`434`		`- for (int i = 0; i < r; ++i) { trS += Ys_->row(i).dot(x.col(i).head(n_dofs_)); }`
	`434`	`+ for (int i = 0; i < r; ++i) { trS += Ys_->row(i).dot(x.col(i).head(m_ * n_dofs_)); }`
`435`	`435`	`return trS / r;`
`436`	`436`	`}`
`437`	`437`	`template <typename... LambdaT>`