# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1664131556 -10800
# Node ID 5aa5c279e341bbdc8f7bc9becb521b62bd990973
# Parent  eb3c7813b67acf6b6b89c71acd0b24594b0fccd6
Attempt at directly referring to measure masses as SliceStorageMut

diff -r eb3c7813b67a -r 5aa5c279e341 src/fb.rs
--- a/src/fb.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/fb.rs	Sun Sep 25 21:45:56 2022 +0300
@@ -457,7 +457,7 @@
     iterator : I,
     plotter : SeqPlotter<F, N>
 ) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
+where F : Float + ToNalgebraRealField<MixedType = F>,
       I : AlgIteratorFactory<IterInfo<F, N>>,
       for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
                                   //+ std::ops::Mul<F, Output=A::Observable>,  <-- FIXME: compiler overflow
@@ -521,7 +521,7 @@
     mut residual : A::Observable,
     mut specialisation : Spec,
 ) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
+where F : Float + ToNalgebraRealField<MixedType=F>,
       I : AlgIteratorFactory<IterInfo<F, N>>,
       Spec : FBSpecialisation<F, A::Observable, N>,
       A::Observable : std::ops::MulAssign<F>,
@@ -657,7 +657,7 @@
                                                μ.iter_locations()
                                                 .map(|ζ| minus_τv.apply(ζ) + ω0.apply(ζ))
                                                 .map(F::to_nalgebra_mixed));
-                let mut x = μ.masses_dvector();
+                let mut x = μ.masses_mut();
 
                 // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
                 // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
@@ -672,7 +672,7 @@
                                                 inner_τ, inner_it);
 
                 // Update masses of μ based on solution of finite-dimensional subproblem.
-                μ.set_masses_dvector(&x);
+                //μ.set_masses_dvector(&x);
             }
 
             // Form d = ω0 - τv - 𝒟μ = -𝒟(μ - μ^k) - τv for checking the proximate optimality
diff -r eb3c7813b67a -r 5aa5c279e341 src/frank_wolfe.rs
--- a/src/frank_wolfe.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/frank_wolfe.rs	Sun Sep 25 21:45:56 2022 +0300
@@ -149,13 +149,13 @@
     inner : &InnerSettings<F>,
     iterator : I
 ) -> usize
-where F : Float + ToNalgebraRealField,
+where F : Float + ToNalgebraRealField<MixedType=F>,
       I : AlgIteratorFactory<F>,
       A : ForwardModel<Loc<F, N>, F>
 {
     // Form and solve finite-dimensional subproblem.
     let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
-    let mut x = μ.masses_dvector();
+    let mut x = μ.masses_mut();
 
     // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
     // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
@@ -168,7 +168,7 @@
     let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
     let iters = quadratic_nonneg(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator);
     // Update masses of μ based on solution of finite-dimensional subproblem.
-    μ.set_masses_dvector(&x);
+    //μ.set_masses_dvector(&x);
 
     iters
 }
@@ -193,7 +193,7 @@
     iterator : I,
     mut plotter : SeqPlotter<F, N>,
 ) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
+where F : Float + ToNalgebraRealField<MixedType=F>,
       I : AlgIteratorFactory<IterInfo<F, N>>,
       for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
                                   //+ std::ops::Mul<F, Output=A::Observable>,  <-- FIXME: compiler overflow
diff -r eb3c7813b67a -r 5aa5c279e341 src/measures/discrete.rs
--- a/src/measures/discrete.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/measures/discrete.rs	Sun Sep 25 21:45:56 2022 +0300
@@ -7,7 +7,14 @@
 };
 use std::iter::Sum;
 use serde::ser::{Serializer, Serialize, SerializeSeq};
-use nalgebra::DVector;
+use nalgebra::{
+    DVector,
+    DVectorSliceMut,
+    Dynamic,
+    Const,
+    SliceStorageMut,
+    Matrix,
+};
 
 use alg_tools::norms::Norm;
 use alg_tools::tabledump::TableDump;
@@ -176,6 +183,25 @@
     }
 }
 
+impl<Domain, F : Float + ToNalgebraRealField<MixedType=F>> DiscreteMeasure<Domain, F> {
+    pub fn masses_mut(&mut self) -> DVectorSliceMut<'_, F::MixedType, Dynamic> {
+        let n = self.spikes.len();
+        unsafe {
+            let start = self.spikes.as_mut_ptr();
+            let next = start.add(1);
+            let ptr = &mut (*start).α as *mut F;
+            let ptrnext = &mut (*next).α as *mut F;
+            let stride = ptrnext.offset_from(ptr);
+            assert_eq!(start.add(stride as usize), next);
+            Matrix::from_data(SliceStorageMut::from_raw_parts(
+                ptr,
+                (Dynamic::new(n), Const),
+                (Dynamic::new(stride as usize), Dynamic::new(0))
+            ))
+        }
+    }
+}
+
 impl<Domain, F :Num> Index<usize> for DiscreteMeasure<Domain, F> {
     type Output = DeltaMeasure<Domain, F>;
     #[inline]
diff -r eb3c7813b67a -r 5aa5c279e341 src/pdps.rs
--- a/src/pdps.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/pdps.rs	Sun Sep 25 21:45:56 2022 +0300
@@ -311,7 +311,7 @@
     plotter : SeqPlotter<F, N>,
     dataterm : D,
 ) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
+where F : Float + ToNalgebraRealField<MixedType=F>,
       I : AlgIteratorFactory<IterInfo<F, N>>,
       for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>
                                   + std::ops::Add<A::Observable, Output=A::Observable>,
diff -r eb3c7813b67a -r 5aa5c279e341 src/subproblem.rs
--- a/src/subproblem.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/subproblem.rs	Sun Sep 25 21:45:56 2022 +0300
@@ -1,7 +1,19 @@
 //! Iterative algorithms for solving finite-dimensional subproblems.
 
 use serde::{Serialize, Deserialize};
-use nalgebra::{DVector, DMatrix};
+use nalgebra::{
+    Matrix,
+    Vector,
+    Storage,
+    StorageMut,
+    Const,
+    Dim,
+    DVector,
+    DMatrix,
+    DefaultAllocator,
+    Dynamic,
+};
+use nalgebra::base::allocator::Allocator;
 use numeric_literals::replace_float_literals;
 use itertools::{izip, Itertools};
 use colored::Colorize;
@@ -80,25 +92,33 @@
 /// The `λ` component of the model is handled in the proximal step instead of the gradient step
 /// for potential performance improvements.
 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn quadratic_nonneg_fb<F, I>(
-    mA : &DMatrix<F::MixedType>,
-    g : &DVector<F::MixedType>,
+pub fn quadratic_nonneg_fb<F, I, S, SA, SG, D>(
+    mA : &Matrix<F::MixedType, D, D, SA>,
+    g : &Vector<F::MixedType, D, SG>,
     //c_ : F,
     λ_ : F,
-    x : &mut DVector<F::MixedType>,
+    x : &mut Vector<F::MixedType, D, S>,
     τ_ : F,
     iterator : I
 ) -> usize
 where F : Float + ToNalgebraRealField,
-      I : AlgIteratorFactory<F>
+      I : AlgIteratorFactory<F>,
+      D : Dim,
+      S : StorageMut<F::MixedType, D, Const::<1>>,
+      SA : Storage<F::MixedType, D, D>,
+      SG : Storage<F::MixedType, D, Const::<1>>,
+      DefaultAllocator : Allocator<F::MixedType, D>
 {
-    let mut xprev = x.clone();
+    let mut xprev = x.clone_owned();
     //let c = c_.to_nalgebra_mixed();
     let λ = λ_.to_nalgebra_mixed();
     let τ = τ_.to_nalgebra_mixed();
     let τλ = τ * λ;
-    let mut v = DVector::zeros(x.len());
     let mut iters = 0;
+    let mut v = {
+        let (r, c) = x.shape_generic();
+        Vector::zeros_generic(r, c)
+    };
 
     iterator.iterate(|state| {
         // Replace `x` with $x - τ[Ax-g]= [x + τg]- τAx$
@@ -195,26 +215,33 @@
 /// forward-backward step.
 /// </p>
 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn quadratic_nonneg_ssn<F, I>(
-    mA : &DMatrix<F::MixedType>,
-    g : &DVector<F::MixedType>,
+pub fn quadratic_nonneg_ssn<F, I, S, SA, SG>(
+    mA : &Matrix<F::MixedType, Dynamic, Dynamic, SA>,
+    g : &Vector<F::MixedType, Dynamic, SG>,
     //c_ : F,
     λ_ : F,
-    x : &mut DVector<F::MixedType>,
+    x : &mut Vector<F::MixedType, Dynamic, S>,
     τ_ : F,
     iterator : I
 ) -> Result<usize, NumericalError>
-where F : Float + ToNalgebraRealField,
-      I : AlgIteratorFactory<F>
+where F : Float + ToNalgebraRealField<MixedType=F>,
+      I : AlgIteratorFactory<F>,
+      S : StorageMut<F::MixedType, Dynamic, Const::<1>>,
+      SA : Storage<F::MixedType, Dynamic, Dynamic>,
+      SG : Storage<F::MixedType, Dynamic, Const::<1>>,
 {
     let n = x.len();
-    let mut xprev = x.clone();
+    let mut xprev = x.clone_owned();
     let mut v = DVector::zeros(n);
     //let c = c_.to_nalgebra_mixed();
     let λ = λ_.to_nalgebra_mixed();
     let τ = τ_.to_nalgebra_mixed();
     let τλ = τ * λ;
     let mut inact : Vec<bool> = Vec::from_iter(std::iter::repeat(false).take(n));
+    let mut v = {
+        let (r, c) = x.shape_generic();
+        Vector::zeros_generic(r, c)
+    };
     let mut s = DVector::zeros(0);
     let mut decomp = nalgebra::linalg::LU::new(DMatrix::zeros(0, 0));
     let mut iters = 0;
@@ -346,18 +373,21 @@
 ///    cone_, <https://doi.org/10.1080/02331934.2014.929680>.
 ///
 /// This function returns the number of iterations taken.
-pub fn quadratic_nonneg<F, I>(
+pub fn quadratic_nonneg<F, I, S, SA, SG>(
     method : InnerMethod,
-    mA : &DMatrix<F::MixedType>,
-    g : &DVector<F::MixedType>,
+    mA : &Matrix<F::MixedType, Dynamic, Dynamic, SA>,
+    g : &Vector<F::MixedType, Dynamic, SG>,
     //c_ : F,
     λ : F,
-    x : &mut DVector<F::MixedType>,
+    x : &mut Vector<F::MixedType, Dynamic, S>,
     τ : F,
     iterator : I
 ) -> usize
-where F : Float + ToNalgebraRealField,
-      I : AlgIteratorFactory<F>
+where F : Float + ToNalgebraRealField<MixedType=F>,
+      I : AlgIteratorFactory<F>,
+      S : StorageMut<F::MixedType, Dynamic, Const::<1>>,
+      SA : Storage<F::MixedType, Dynamic, Dynamic>,
+      SG : Storage<F::MixedType, Dynamic, Const::<1>>,
 {
     
     match method {