--- a/src/prox_penalty/radon_squared.rs Tue Apr 08 13:31:39 2025 -0500
+++ b/src/prox_penalty/radon_squared.rs Fri May 08 16:47:58 2026 -0500
@@ -4,148 +4,78 @@
Instead of the $𝒟$-norm of `fb.rs`, this uses a standard Radon norm for the proximal map.
*/
+use super::{InsertionConfig, ProxPenalty, ProxTerm, StepLengthBound, StepLengthBoundPD};
+use crate::dataterm::QuadraticDataTerm;
+use crate::forward_model::ForwardModel;
+use crate::measures::merging::SpikeMerging;
+use crate::measures::{DiscreteMeasure, Radon};
+use crate::regularisation::RadonSquaredRegTerm;
+use crate::types::*;
+use alg_tools::bounds::MinMaxMapping;
+use alg_tools::error::DynResult;
+use alg_tools::instance::{Instance, Space};
+use alg_tools::iterate::{AlgIterator, AlgIteratorIteration};
+use alg_tools::linops::BoundedLinear;
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::{Norm, L2};
use numeric_literals::replace_float_literals;
-use serde::{Serialize, Deserialize};
-use nalgebra::DVector;
-
-use alg_tools::iterate::{
- AlgIteratorIteration,
- AlgIterator
-};
-use alg_tools::norms::{L2, Norm};
-use alg_tools::linops::Mapping;
-use alg_tools::bisection_tree::{
- BTFN,
- Bounds,
- BTSearch,
- SupportGenerator,
- LocalAnalysis,
-};
-use alg_tools::mapping::RealMapping;
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-
-use crate::types::*;
-use crate::measures::{
- RNDM,
- DeltaMeasure,
- Radon,
-};
-use crate::measures::merging::SpikeMerging;
-use crate::regularisation::RegTerm;
-use crate::forward_model::{
- ForwardModel,
- AdjointProductBoundedBy
-};
-use super::{
- FBGenericConfig,
- ProxPenalty,
-};
+use serde::{Deserialize, Serialize};
/// Radon-norm squared proximal penalty
-#[derive(Copy,Clone,Serialize,Deserialize)]
+#[derive(Copy, Clone, Serialize, Deserialize)]
pub struct RadonSquared;
#[replace_float_literals(F::cast_from(literal))]
-impl<F, GA, BTA, S, Reg, const N : usize>
-ProxPenalty<F, BTFN<F, GA, BTA, N>, Reg, N> for RadonSquared
+impl<Domain, F, M, Reg> ProxPenalty<Domain, M, Reg, F> for RadonSquared
where
- F : Float + ToNalgebraRealField,
- GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
- BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
- S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
- Reg : RegTerm<F, N>,
- RNDM<F, N> : SpikeMerging<F>,
+ Domain: Space + Clone + PartialEq + 'static,
+ for<'a> &'a Domain: Instance<Domain>,
+ F: Float + ToNalgebraRealField,
+ M: MinMaxMapping<Domain, F>,
+ Reg: RadonSquaredRegTerm<Domain, F>,
+ DiscreteMeasure<Domain, F>: SpikeMerging<F>,
{
- type ReturnMapping = BTFN<F, GA, BTA, N>;
+ type ReturnMapping = M;
+
+ fn prox_type() -> ProxTerm {
+ ProxTerm::RadonSquared
+ }
fn insert_and_reweigh<I>(
&self,
- μ : &mut RNDM<F, N>,
- τv : &mut BTFN<F, GA, BTA, N>,
- μ_base : &RNDM<F, N>,
- ν_delta: Option<&RNDM<F, N>>,
- τ : F,
- ε : F,
- config : &FBGenericConfig<F>,
- reg : &Reg,
- _state : &AlgIteratorIteration<I>,
- stats : &mut IterInfo<F, N>,
- ) -> (Option<Self::ReturnMapping>, bool)
+ μ: &mut DiscreteMeasure<Domain, F>,
+ τv: &mut M,
+ τ: F,
+ ε: F,
+ config: &InsertionConfig<F>,
+ reg: &Reg,
+ _state: &AlgIteratorIteration<I>,
+ stats: &mut IterInfo<F>,
+ ) -> DynResult<(Option<Self::ReturnMapping>, bool)>
where
- I : AlgIterator
+ I: AlgIterator,
{
- let mut y = μ_base.masses_dvector();
+ let violation = reg.find_tolerance_violation(τv, τ, ε, true, config);
+ reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats);
- assert!(μ_base.len() <= μ.len());
-
- 'i_and_w: for i in 0..=1 {
- // Optimise weights
- if μ.len() > 0 {
- // Form finite-dimensional subproblem. The subproblem references to the original μ^k
- // from the beginning of the iteration are all contained in the immutable c and g.
- // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
- // problems have not yet been updated to sign change.
- let g̃ = DVector::from_iterator(μ.len(),
- μ.iter_locations()
- .map(|ζ| - F::to_nalgebra_mixed(τv.apply(ζ))));
- let mut x = μ.masses_dvector();
- y.extend(std::iter::repeat(0.0.to_nalgebra_mixed()).take(0.max(x.len()-y.len())));
- assert_eq!(y.len(), x.len());
- // Solve finite-dimensional subproblem.
- // TODO: This assumes that ν_delta has no common locations with μ-μ_base, to
- // ignore it.
- stats.inner_iters += reg.solve_findim_l1squared(&y, &g̃, τ, &mut x, ε, config);
-
- // Update masses of μ based on solution of finite-dimensional subproblem.
- μ.set_masses_dvector(&x);
- }
-
- if i>0 {
- // Simple debugging test to see if more inserts would be needed. Doesn't seem so.
- //let n = μ.dist_matching(μ_base);
- //println!("{:?}", reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n));
- break 'i_and_w
- }
-
- // Calculate ‖μ - μ_base‖_ℳ
- // TODO: This assumes that ν_delta has no common locations with μ-μ_base.
- let n = μ.dist_matching(μ_base) + ν_delta.map_or(0.0, |ν| ν.norm(Radon));
-
- // Find a spike to insert, if needed.
- // This only check the overall tolerances, not tolerances on support of μ-μ_base or μ,
- // which are supposed to have been guaranteed by the finite-dimensional weight optimisation.
- match reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n) {
- None => { break 'i_and_w },
- Some((ξ, _v_ξ, _in_bounds)) => {
- // Weight is found out by running the finite-dimensional optimisation algorithm
- // above
- *μ += DeltaMeasure { x : ξ, α : 0.0 };
- stats.inserted += 1;
- }
- };
- }
-
- (None, true)
+ Ok((None, true))
}
fn merge_spikes(
&self,
- μ : &mut RNDM<F, N>,
- τv : &mut BTFN<F, GA, BTA, N>,
- μ_base : &RNDM<F, N>,
- ν_delta: Option<&RNDM<F, N>>,
- τ : F,
- ε : F,
- config : &FBGenericConfig<F>,
- reg : &Reg,
- fitness : Option<impl Fn(&RNDM<F, N>) -> F>,
- ) -> usize
- {
+ μ: &mut DiscreteMeasure<Domain, F>,
+ τv: &mut M,
+ μ_base: &DiscreteMeasure<Domain, F>,
+ τ: F,
+ ε: F,
+ config: &InsertionConfig<F>,
+ reg: &Reg,
+ fitness: Option<impl Fn(&DiscreteMeasure<Domain, F>) -> F>,
+ ) -> usize {
if config.fitness_merging {
if let Some(f) = fitness {
- return μ.merge_spikes_fitness(config.merging, f, |&v| v)
- .1
+ return μ.merge_spikes_fitness(config.merging, f, |&v| v).1;
}
}
μ.merge_spikes(config.merging, |μ_candidate| {
@@ -156,10 +86,7 @@
// after μ_candidate's extra points.
// TODO: doesn't seem to work, maybe need to merge μ_base as well?
// Although that doesn't seem to make sense.
- let μ_radon = match ν_delta {
- None => μ_candidate.sub_matching(μ_base),
- Some(ν) => μ_candidate.sub_matching(μ_base) - ν,
- };
+ let μ_radon = μ_candidate.sub_matching(μ_base);
reg.verify_merge_candidate_radonsq(τv, μ_candidate, τ, ε, &config, &μ_radon)
//let n = μ_candidate.dist_matching(μ_base);
//reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n).is_none()
@@ -167,16 +94,27 @@
}
}
-
-impl<F, A, const N : usize> AdjointProductBoundedBy<RNDM<F, N>, RadonSquared>
-for A
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, A, Domain> StepLengthBound<F, QuadraticDataTerm<F, Domain, A>> for RadonSquared
where
- F : Float,
- A : ForwardModel<RNDM<F, N>, F>
+ F: Float + ToNalgebraRealField,
+ Domain: Space + Norm<Radon, F>,
+ A: ForwardModel<Domain, F> + BoundedLinear<Domain, Radon, L2, F>,
{
- type FloatType = F;
-
- fn adjoint_product_bound(&self, _ : &RadonSquared) -> Option<Self::FloatType> {
- self.opnorm_bound(Radon, L2).powi(2).into()
+ fn step_length_bound(&self, f: &QuadraticDataTerm<F, Domain, A>) -> DynResult<F> {
+ // TODO: direct squared calculation
+ Ok(f.operator().opnorm_bound(Radon, L2)?.powi(2))
}
}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, A, Domain> StepLengthBoundPD<F, A, DiscreteMeasure<Domain, F>> for RadonSquared
+where
+ Domain: Space + Clone + PartialEq + 'static,
+ F: Float + ToNalgebraRealField,
+ A: BoundedLinear<DiscreteMeasure<Domain, F>, Radon, L2, F>,
+{
+ fn step_length_bound_pd(&self, opA: &A) -> DynResult<F> {
+ opA.opnorm_bound(Radon, L2)
+ }
+}