Mercurial > repos > pointsource_algs / file diff

--- a/src/forward_model.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/forward_model.rs	Tue Nov 29 15:36:12 2022 +0200
@@ -53,14 +53,22 @@
     type Observable : Euclidean<F, Output=Self::Observable>
                       + AXPY<F>
                       + Clone;
+    /// Finite-dimensional model for weight optimisation
+    type FindimModel : Linear<DVector<F::MixedType>, Codomain=Self::Observable>;
 
-    /// Return A_*A and A_* b
+    /// Return A_*A and A_* b instantiated at the locations of the spikes of `μ`.
     fn findim_quadratic_model(
         &self,
         μ : &DiscreteMeasure<Domain, F>,
         b : &Self::Observable
     ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>);
 
+    /// Return A instantiated at the locations of the spikes of `μ`
+    fn findim_model(
+        &self,
+        μ : &DiscreteMeasure<Domain, F>,
+    ) -> Self::FindimModel;
+
     /// Write an observable into a file.
     fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError;
 
@@ -517,6 +525,7 @@
       ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
       Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
     type Observable = DVector<F>;
+    type FindimModel = DMatrix<F>;
 
     fn findim_quadratic_model(
         &self,
@@ -524,6 +533,15 @@
         b : &Self::Observable
     ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>) {
         assert_eq!(b.len(), self.n_sensors());
+        let mA = self.findim_model(μ);
+        let mAt = mA.transpose();
+        (&mAt * mA, &mAt * b)
+    }
+
+    fn findim_model(
+        &self,
+        μ : &DiscreteMeasure<Loc<F, N>, F>,
+    ) -> DMatrix<F::MixedType> {
         let mut mA = DMatrix::zeros(self.n_sensors(), μ.len());
         let grid = self.grid();
         for (mut mAcol, δ) in mA.column_iter_mut().zip(μ.iter_spikes()) {
@@ -532,8 +550,7 @@
                 mAcol[d] += sensor.apply(&δ.x);
             }
         }
-        let mAt = mA.transpose();
-        (&mAt * mA, &mAt * b)
+        mA
     }
 
     fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError {