Mercurial > repos > pointsource_algs / file diff

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dataterm.rs	Mon Feb 17 13:54:53 2025 -0500
@@ -0,0 +1,94 @@
+/*!
+Basid definitions for data terms
+*/
+
+use numeric_literals::replace_float_literals;
+
+use alg_tools::euclidean::Euclidean;
+use alg_tools::linops::GEMV;
+pub use alg_tools::norms::L1;
+use alg_tools::norms::Norm;
+use alg_tools::instance::{Instance, Space};
+
+use crate::types::*;
+pub use crate::types::L2Squared;
+use crate::measures::RNDM;
+
+/// Calculates the residual $Aμ-b$.
+#[replace_float_literals(F::cast_from(literal))]
+pub(crate) fn calculate_residual<
+    X : Space,
+    I : Instance<X>,
+    F : Float,
+    V : Euclidean<F> + Clone,
+    A : GEMV<F, X, Codomain = V>,
+>(
+    μ : I,
+    opA : &A,
+    b : &V
+) -> V {
+    let mut r = b.clone();
+    opA.gemv(&mut r, 1.0, μ, -1.0);
+    r
+}
+
+/// Calculates the residual $A(μ+μ_delta)-b$.
+#[replace_float_literals(F::cast_from(literal))]
+pub(crate) fn calculate_residual2<
+    F : Float,
+    X : Space,
+    I : Instance<X>,
+    J : Instance<X>,
+    V : Euclidean<F> + Clone,
+    A : GEMV<F, X, Codomain = V>,
+>(
+    μ : I,
+    μ_delta : J,
+    opA : &A,
+    b : &V
+) -> V {
+    let mut r = b.clone();
+    opA.gemv(&mut r, 1.0, μ, -1.0);
+    opA.gemv(&mut r, 1.0, μ_delta, 1.0);
+    r
+}
+
+
+/// Trait for data terms
+#[replace_float_literals(F::cast_from(literal))]
+pub trait DataTerm<F : Float, V, const N : usize> {
+    /// Calculates $F(y)$, where $F$ is the data fidelity.
+    fn calculate_fit(&self, _residual : &V) -> F;
+
+    /// Calculates $F(Aμ-b)$, where $F$ is the data fidelity.
+    fn calculate_fit_op<I, A : GEMV<F, RNDM<F, N>, Codomain = V>>(
+        &self,
+        μ : I,
+        opA : &A,
+        b : &V
+    ) -> F
+    where
+        V : Euclidean<F> + Clone,
+        I : Instance<RNDM<F, N>>,
+    {
+        let r = calculate_residual(μ, opA, b);
+        self.calculate_fit(&r)
+    }
+}
+
+impl<F : Float, V : Euclidean<F>, const N : usize>
+DataTerm<F, V, N>
+for L2Squared {
+    fn calculate_fit(&self, residual : &V) -> F {
+        residual.norm2_squared_div2()
+    }
+}
+
+
+impl<F : Float, V : Euclidean<F> + Norm<F, L1>, const N : usize>
+DataTerm<F, V, N>
+for L1 {
+    fn calculate_fit(&self, residual : &V) -> F {
+        residual.norm(L1)
+    }
+}