Mercurial > repos > pointsource_algs / file comparison

comparison: src/prox_penalty/wave.rs

src/prox_penalty/wave.rs

changeset 75
677a5fd1b014
parent 63
7a8a55fd41c0
child 76
b921ed0ab99b
equal deleted inserted replaced
74:df92e78cc3f4 75:677a5fd1b014
50 ε: F, 50 ε: F,
51 config: &InsertionConfig<F>, 51 config: &InsertionConfig<F>,
52 reg: &Reg, 52 reg: &Reg,
53 state: &AlgIteratorIteration<I>, 53 state: &AlgIteratorIteration<I>,
54 stats: &mut IterInfo<F>, 54 stats: &mut IterInfo<F>,
55 ) -> DynResult<(Option<Self::ReturnMapping>, bool)> 55 ) -> DynResult<bool>
56 where 56 where
57 I: AlgIterator, 57 I: AlgIterator,
58 { 58 {
59 let op𝒟norm = self.opnorm_bound(Radon, Linfinity)?; 59 let op𝒟norm = self.opnorm_bound(Radon, Linfinity)?;
60 60
68 let μ_base = μ.clone(); 68 let μ_base = μ.clone();
69 let ω0 = self.apply(&μ_base); 69 let ω0 = self.apply(&μ_base);
70 70
71 // Add points to support until within error tolerance or maximum insertion count reached. 71 // Add points to support until within error tolerance or maximum insertion count reached.
72 let mut count = 0; 72 let mut count = 0;
73 let (within_tolerances, d) = 'insertion: loop { 73 let within_tolerances = 'insertion: loop {
74 if μ.len() > 0 { 74 if μ.len() > 0 {
75 // Form finite-dimensional subproblem. The subproblem references to the original μ^k 75 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
76 // from the beginning of the iteration are all contained in the immutable c and g. 76 // from the beginning of the iteration are all contained in the immutable c and g.
77 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional 77 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
78 // problems have not yet been updated to sign change. 78 // problems have not yet been updated to sign change.
113 }; 113 };
114 114
115 // Find a spike to insert, if needed 115 // Find a spike to insert, if needed
116 let (ξ, _v_ξ, in_bounds) = 116 let (ξ, _v_ξ, in_bounds) =
117 match reg.find_tolerance_violation(&mut d, τ, ε, skip_by_rough_check, config) { 117 match reg.find_tolerance_violation(&mut d, τ, ε, skip_by_rough_check, config) {
118 None => break 'insertion (true, d), 118 None => break 'insertion true,
119 Some(res) => res, 119 Some(res) => res,
120 }; 120 };
121 121
122 // Break if maximum insertion count reached 122 // Break if maximum insertion count reached
123 if count >= max_insertions { 123 if count >= max_insertions {
124 break 'insertion (in_bounds, d); 124 break 'insertion in_bounds;
125 } 125 }
126 126
127 // No point in optimising the weight here; the finite-dimensional algorithm is fast. 127 // No point in optimising the weight here; the finite-dimensional algorithm is fast.
128 *μ += DeltaMeasure { x: ξ, α: 0.0 }; 128 *μ += DeltaMeasure { x: ξ, α: 0.0 };
129 count += 1; 129 count += 1;
137 subproblem solution tolerance" 137 subproblem solution tolerance"
138 ); 138 );
139 println!("{}", err.red()); 139 println!("{}", err.red());
140 } 140 }
141 141
142 Ok((Some(d), within_tolerances)) 142 Ok(within_tolerances)
143 }
144
145 fn reweigh<I>(
146 &self,
147 μ: &mut DiscreteMeasure<Domain, F>,
148 τv: &mut M,
149 τ: F,
150 ε: F,
151 config: &InsertionConfig<F>,
152 reg: &Reg,
153 _state: &AlgIteratorIteration<I>,
154 stats: &mut IterInfo<F>,
155 ) -> DynResult<()>
156 where
157 I: AlgIterator,
158 {
159 if μ.len() > 0 {
160 let op𝒟norm = self.opnorm_bound(Radon, Linfinity)?;
161 let ω0 = self.apply(&*μ);
162
163 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
164 // from the beginning of the iteration are all contained in the immutable c and g.
165 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
166 // problems have not yet been updated to sign change.
167 let à = self.findim_matrix(μ.iter_locations());
168 let g̃ = DVector::from_iterator(
169 μ.len(),
170 μ.iter_locations()
171 .map(|ζ| ω0.apply(ζ) - τv.apply(ζ))
172 .map(F::to_nalgebra_mixed),
173 );
174 let mut x = μ.masses_dvector();
175
176 // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
177 // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
178 // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ
179 // ≤ sup_{|z|_2 ≤ 1} |z|_1 |𝒟| |x|_1 ≤ sup_{|z|_2 ≤ 1} n |z|_2 |𝒟| |x|_2
180 // = n |𝒟| |x|_2, where n is the number of points. Therefore
181 let Ã_normest = op𝒟norm * F::cast_from(μ.len());
182
183 // Solve finite-dimensional subproblem.
184 stats.inner_iters += reg.solve_findim(&Ã, &g̃, τ, &mut x, Ã_normest, ε, config);
185
186 // Update masses of μ based on solution of finite-dimensional subproblem.
187 μ.set_masses_dvector(&x);
188 }
189
190 Ok(())
143 } 191 }
144 192
145 fn merge_spikes( 193 fn merge_spikes(
146 &self, 194 &self,
147 μ: &mut DiscreteMeasure<Domain, F>, 195 μ: &mut DiscreteMeasure<Domain, F>,

Mercurial Repository: pointsource_algs

mercurial