Mon, 24 Oct 2022 09:41:43 +0300
Allow step closure of AlgIterators to indicate succesfull termination or failure.
0 | 1 | use crate::types::*; |
2 | use crate::loc::Loc; | |
3 | use crate::sets::{Set,NPolygon,SpannedHalfspace}; | |
4 | use crate::linsolve::*; | |
5 | use crate::norms::Dot; | |
6 | use super::base::{LocalModel,RealLocalModel}; | |
7 | use crate::sets::Cube; | |
8 | use numeric_literals::replace_float_literals; | |
9 | ||
10 | pub struct Simplex<F : Float, const N : usize, const D : usize>(pub [Loc<F, N>; D]); | |
11 | pub type PlanarSimplex<F> = Simplex<F,2,3>; | |
12 | pub type RealInterval<F> = Simplex<F,1,2>; | |
13 | ||
14 | #[inline] | |
15 | #[replace_float_literals(F::cast_from(literal))] | |
16 | pub(crate) fn midpoint<F : Float, const N : usize>(a : &Loc<F,N>, b : &Loc<F,N>) -> Loc<F, N> { | |
17 | (a+b)/2.0 | |
18 | } | |
19 | ||
20 | impl<'a, F : Float> Set<Loc<F,1>> for RealInterval<F> { | |
21 | #[inline] | |
22 | fn contains(&self, &Loc([x]) : &Loc<F,1>) -> bool { | |
23 | let &[Loc([x0]), Loc([x1])] = &self.0; | |
24 | (x0 < x && x < x1) || (x1 < x && x < x0) | |
25 | } | |
26 | } | |
27 | ||
28 | impl<'a, F : Float> Set<Loc<F,2>> for PlanarSimplex<F> { | |
29 | #[inline] | |
30 | fn contains(&self, x : &Loc<F,2>) -> bool { | |
31 | let &[x0, x1, x2] = &self.0; | |
32 | NPolygon([[x0, x1].spanned_halfspace(), | |
33 | [x1, x2].spanned_halfspace(), | |
34 | [x2, x0].spanned_halfspace()]).contains(x) | |
35 | } | |
36 | } | |
37 | ||
38 | trait P2Powers { | |
39 | type Output; | |
40 | type Diff; | |
41 | type Full; | |
42 | fn p2powers(&self) -> Self::Output; | |
43 | fn p2powers_full(&self) -> Self::Full; | |
44 | fn p2powers_diff(&self) -> Self::Diff; | |
45 | } | |
46 | ||
47 | #[replace_float_literals(F::cast_from(literal))] | |
48 | impl<F : Float> P2Powers for Loc<F, 1> { | |
49 | type Output = Loc<F, 1>; | |
50 | type Full = Loc<F, 3>; | |
51 | type Diff = Loc<Loc<F, 1>, 1>; | |
52 | ||
53 | #[inline] | |
54 | fn p2powers(&self) -> Self::Output { | |
55 | let &Loc([x0]) = self; | |
56 | [x0*x0].into() | |
57 | } | |
58 | ||
59 | #[inline] | |
60 | fn p2powers_full(&self) -> Self::Full { | |
61 | let &Loc([x0]) = self; | |
62 | [1.0, x0, x0*x0].into() | |
63 | } | |
64 | ||
65 | #[inline] | |
66 | fn p2powers_diff(&self) -> Self::Diff { | |
67 | let &Loc([x0]) = self; | |
68 | [[x0+x0].into()].into() | |
69 | } | |
70 | } | |
71 | ||
72 | #[replace_float_literals(F::cast_from(literal))] | |
73 | impl<F : Float> P2Powers for Loc<F, 2> { | |
74 | type Output = Loc<F, 3>; | |
75 | type Full = Loc<F, 6>; | |
76 | type Diff = Loc<Loc<F, 3>, 2>; | |
77 | ||
78 | #[inline] | |
79 | fn p2powers(&self) -> Self::Output { | |
80 | let &Loc([x0, x1]) = self; | |
81 | [x0*x0, x0*x1, x1*x1].into() | |
82 | } | |
83 | ||
84 | #[inline] | |
85 | fn p2powers_full(&self) -> Self::Full { | |
86 | let &Loc([x0, x1]) = self; | |
87 | [1.0, x0, x1, x0*x0, x0*x1, x1*x1].into() | |
88 | } | |
89 | ||
90 | #[inline] | |
91 | fn p2powers_diff(&self) -> Self::Diff { | |
92 | let &Loc([x0, x1]) = self; | |
93 | [[x0+x0, x1, 0.0].into(), [0.0, x0, x1+x1].into()].into() | |
94 | } | |
95 | } | |
96 | ||
97 | pub trait P2Model<F : Num, const N : usize> { | |
98 | type Model : LocalModel<Loc<F,N>,F>; | |
99 | fn p2_model<G : Fn(&Loc<F, N>) -> F>(&self, g : G) -> Self::Model; | |
100 | } | |
101 | ||
102 | pub struct P2LocalModel<F : Num, const N : usize, const E : usize, const V : usize, const Q : usize> { | |
103 | a0 : F, | |
104 | a1 : Loc<F, N>, | |
105 | a2 : Loc<F, E>, | |
106 | //node_values : Loc<F, V>, | |
107 | //edge_values : Loc<F, Q>, | |
108 | } | |
109 | ||
110 | // | |
111 | // 1D planar model construction | |
112 | // | |
113 | ||
114 | impl<F : Float> RealInterval<F> { | |
115 | #[inline] | |
116 | fn midpoints(&self) -> [Loc<F, 1>; 1] { | |
117 | let [ref n0, ref n1] = &self.0; | |
118 | let n01 = midpoint(n0, n1); | |
119 | [n01] | |
120 | } | |
121 | } | |
122 | ||
123 | impl<F : Float> P2LocalModel<F, 1, 1, 2, 0> { | |
124 | #[inline] | |
125 | pub fn new( | |
126 | &[n0, n1, n01] : &[Loc<F, 1>; 3], | |
127 | &[v0, v1, v01] : &[F; 3], | |
128 | ) -> Self { | |
129 | let p = move |x : &Loc<F, 1>, v : F| { | |
130 | let Loc([c, d, e]) = x.p2powers_full(); | |
131 | [c, d, e, v] | |
132 | }; | |
133 | let [a0, a1, a11] = linsolve([ | |
134 | p(&n0, v0), | |
135 | p(&n1, v1), | |
136 | p(&n01, v01) | |
137 | ]); | |
138 | P2LocalModel { | |
139 | a0 : a0, | |
140 | a1 : [a1].into(), | |
141 | a2 : [a11].into(), | |
142 | //node_values : [v0, v1].into(), | |
143 | //edge_values: [].into(), | |
144 | } | |
145 | } | |
146 | } | |
147 | ||
148 | impl<F : Float> P2Model<F,1> for RealInterval<F> { | |
149 | type Model = P2LocalModel<F, 1, 1, 2, 0>; | |
150 | ||
151 | #[inline] | |
152 | fn p2_model<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> Self::Model { | |
153 | let [n01] = self.midpoints(); | |
154 | let [n0, n1] = self.0; | |
155 | let vals = [g(&n0), g(&n1), g(&n01)]; | |
156 | let nodes = [n0, n1, n01]; | |
157 | Self::Model::new(&nodes, &vals) | |
158 | } | |
159 | } | |
160 | ||
161 | // | |
162 | // 2D planar model construction | |
163 | // | |
164 | ||
165 | impl<F : Float> PlanarSimplex<F> { | |
166 | #[inline] | |
167 | fn midpoints(&self) -> [Loc<F, 2>; 3] { | |
168 | let [ref n0, ref n1, ref n2] = &self.0; | |
169 | let n01 = midpoint(n0, n1); | |
170 | let n12 = midpoint(n1, n2); | |
171 | let n20 = midpoint(n2, n0); | |
172 | [n01, n12, n20] | |
173 | } | |
174 | } | |
175 | ||
176 | impl<F : Float> P2LocalModel<F, 2, 3, 3, 3> { | |
177 | #[inline] | |
178 | pub fn new( | |
179 | &[n0, n1, n2, n01, n12, n20] : &[Loc<F, 2>; 6], | |
180 | &[v0, v1, v2, v01, v12, v20] : &[F; 6], | |
181 | ) -> Self { | |
182 | let p = move |x : &Loc<F,2>, v :F| { | |
183 | let Loc([c, d, e, f, g, h]) = x.p2powers_full(); | |
184 | [c, d, e, f, g, h, v] | |
185 | }; | |
186 | let [a0, a1, a2, a11, a12, a22] = linsolve([ | |
187 | p(&n0, v0), | |
188 | p(&n1, v1), | |
189 | p(&n2, v2), | |
190 | p(&n01, v01), | |
191 | p(&n12, v12), | |
192 | p(&n20, v20), | |
193 | ]); | |
194 | P2LocalModel { | |
195 | a0 : a0, | |
196 | a1 : [a1, a2].into(), | |
197 | a2 : [a11, a12, a22].into(), | |
198 | //node_values : [v0, v1, v2].into(), | |
199 | //edge_values: [v01, v12, v20].into(), | |
200 | } | |
201 | } | |
202 | } | |
203 | ||
204 | impl<F : Float> P2Model<F,2> for PlanarSimplex<F> { | |
205 | type Model = P2LocalModel<F, 2, 3, 3, 3>; | |
206 | ||
207 | #[inline] | |
208 | fn p2_model<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> Self::Model { | |
209 | let midpoints = self.midpoints(); | |
210 | let [ref n0, ref n1, ref n2] = self.0; | |
211 | let [ref n01, ref n12, ref n20] = midpoints; | |
212 | let vals = [g(n0), g(n1), g(n2), g(n01), g(n12), g(n20)]; | |
213 | let nodes = [*n0, *n1, *n2, *n01, *n12, *n20]; | |
214 | Self::Model::new(&nodes, &vals) | |
215 | } | |
216 | } | |
217 | ||
218 | macro_rules! impl_local_model { | |
219 | ($n:literal, $e:literal, $v:literal, $q:literal) => { | |
220 | impl<F : Float> LocalModel<Loc<F, $n>, F> for P2LocalModel<F, $n, $e, $v, $q> { | |
221 | #[inline] | |
222 | fn value(&self, x : &Loc<F,$n>) -> F { | |
223 | self.a0 + x.dot(&self.a1) + x.p2powers().dot(&self.a2) | |
224 | } | |
225 | ||
226 | #[inline] | |
227 | fn differential(&self, x : &Loc<F,$n>) -> Loc<F,$n> { | |
228 | self.a1 + x.p2powers_diff().map(|di| di.dot(&self.a2)) | |
229 | } | |
230 | } | |
231 | } | |
232 | } | |
233 | ||
234 | impl_local_model!(1, 1, 2, 0); | |
235 | impl_local_model!(2, 3, 3, 3); | |
236 | ||
237 | ||
238 | // | |
239 | // Minimisation | |
240 | // | |
241 | ||
242 | #[replace_float_literals(F::cast_from(literal))] | |
243 | impl<F : Float> P2LocalModel<F, 1, 1, 2, 0> { | |
244 | #[inline] | |
245 | fn minimise_edge(&self, x0 : Loc<F, 1>, x1 : Loc<F,1>) -> (Loc<F,1>, F) { | |
246 | let &P2LocalModel{ | |
247 | a1 : Loc([a1]), | |
248 | a2 : Loc([a11]), | |
249 | //node_values : Loc([v0, v1]), | |
250 | .. | |
251 | } = self; | |
252 | // We do this in cases, first trying for an interior solution, then edges. | |
253 | // For interior solution, first check determinant; no point trying if non-positive | |
254 | if a11 > 0.0 { | |
255 | // An interior solution x[1] has to satisfy | |
256 | // 2a₁₁*x[1] + a₁ =0 | |
257 | // This gives | |
258 | let t = -a1/(2.0*a11); | |
259 | let (Loc([t0]), Loc([t1])) = (x0, x1); | |
260 | if (t0 <= t && t <= t1) || (t1 <= t && t <= t0) { | |
261 | let x = [t].into(); | |
262 | let v = self.value(&x); | |
263 | return (x, v) | |
264 | } | |
265 | } | |
266 | ||
267 | let v0 = self.value(&x0); | |
268 | let v1 = self.value(&x1); | |
269 | if v0 < v1 { (x0, v0) } else { (x1, v1) } | |
270 | } | |
271 | } | |
272 | ||
273 | impl<'a, F : Float> RealLocalModel<RealInterval<F>,Loc<F,1>,F> | |
274 | for P2LocalModel<F, 1, 1, 2, 0> { | |
275 | #[inline] | |
276 | fn minimise(&self, &Simplex([x0, x1]) : &RealInterval<F>) -> (Loc<F,1>, F) { | |
277 | self.minimise_edge(x0, x1) | |
278 | } | |
279 | } | |
280 | ||
281 | #[replace_float_literals(F::cast_from(literal))] | |
282 | impl<F : Float> P2LocalModel<F, 2, 3, 3, 3> { | |
283 | /// Minimise the 2D model on the edge {x0 + t(x1 - x0) | t ∈ [0, 1] } | |
284 | #[inline] | |
285 | fn minimise_edge(&self, x0 : &Loc<F,2>, x1 : &Loc<F,2>/*, v0 : F, v1 : F*/) -> (Loc<F,2>, F) { | |
286 | let &P2LocalModel { | |
287 | a0, | |
288 | a1 : Loc([a1, a2]), | |
289 | a2 : Loc([a11, a12, a22]), | |
290 | .. | |
291 | } = self; | |
292 | let &Loc([x00, x01]) = x0; | |
293 | let d@Loc([d0, d1]) = x1 - x0; | |
294 | let b0 = a0 + a1*x00 + a2*x01 + a11*x00*x00 + a12*x00*x01 + a22*x01*x01; | |
295 | let b1 = a1*d0 + a2*d1 + 2.0*a11*d0*x00 + a12*(d0*x01 + d1*x00) + 2.0*a22*d1*x01; | |
296 | let b11 = a11*d0*d0 + a12*d0*d1 + a22*d1*d1; | |
297 | let edge_1d_model = P2LocalModel { | |
298 | a0 : b0, | |
299 | a1 : Loc([b1]), | |
300 | a2 : Loc([b11]), | |
301 | //node_values : Loc([v0, v1]), | |
302 | }; | |
303 | let (Loc([t]), v) = edge_1d_model.minimise_edge(0.0.into(), 1.0.into()); | |
304 | (x0 + d*t, v) | |
305 | } | |
306 | } | |
307 | ||
308 | #[replace_float_literals(F::cast_from(literal))] | |
309 | impl<'a, F : Float> RealLocalModel<PlanarSimplex<F>,Loc<F,2>,F> | |
310 | for P2LocalModel<F, 2, 3, 3, 3> { | |
311 | #[inline] | |
312 | fn minimise(&self, el : &PlanarSimplex<F>) -> (Loc<F,2>, F) { | |
313 | let &P2LocalModel { | |
314 | a1 : Loc([a1, a2]), | |
315 | a2 : Loc([a11, a12, a22]), | |
316 | //node_values : Loc([v0, v1, v2]), | |
317 | .. | |
318 | } = self; | |
319 | ||
320 | // We do this in cases, first trying for an interior solution, then edges. | |
321 | // For interior solution, first check determinant; no point trying if non-positive | |
322 | let r = 2.0*(a11*a22-a12*a12); | |
323 | if r > 0.0 { | |
324 | // An interior solution (x[1], x[2]) has to satisfy | |
325 | // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 | |
326 | // This gives | |
327 | let x = [(a22*a1-a12*a2)/r, (a12*a1-a11*a2)/r].into(); | |
328 | if el.contains(&x) { | |
329 | return (x, self.value(&x)) | |
330 | } | |
331 | } | |
332 | ||
333 | let &[ref x0, ref x1, ref x2] = &el.0; | |
334 | let mut min_edge = self.minimise_edge(x0, x1); | |
335 | let more_edge = [self.minimise_edge(x1, x2), | |
336 | self.minimise_edge(x2, x0)]; | |
337 | ||
338 | for edge in more_edge { | |
339 | if edge.1 < min_edge.1 { | |
340 | min_edge = edge; | |
341 | } | |
342 | } | |
343 | ||
344 | min_edge | |
345 | } | |
346 | } | |
347 | ||
348 | #[replace_float_literals(F::cast_from(literal))] | |
349 | impl<'a, F : Float> RealLocalModel<Cube<F, 2>,Loc<F,2>,F> | |
350 | for P2LocalModel<F, 2, 3, 3, 3> { | |
351 | #[inline] | |
352 | fn minimise(&self, el : &Cube<F, 2>) -> (Loc<F,2>, F) { | |
353 | let &P2LocalModel { | |
354 | a1 : Loc([a1, a2]), | |
355 | a2 : Loc([a11, a12, a22]), | |
356 | //node_values : Loc([v0, v1, v2]), | |
357 | .. | |
358 | } = self; | |
359 | ||
360 | // We do this in cases, first trying for an interior solution, then edges. | |
361 | // For interior solution, first check determinant; no point trying if non-positive | |
362 | let r = 2.0*(a11*a22-a12*a12); | |
363 | if r > 0.0 { | |
364 | // An interior solution (x[1], x[2]) has to satisfy | |
365 | // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 | |
366 | // This gives | |
367 | let x = [(a22*a1-a12*a2)/r, (a12*a1-a11*a2)/r].into(); | |
368 | if el.contains(&x) { | |
369 | return (x, self.value(&x)) | |
370 | } | |
371 | } | |
372 | ||
373 | let [x0, x1, x2, x3] = el.corners(); | |
374 | let mut min_edge = self.minimise_edge(&x0, &x1); | |
375 | let more_edge = [self.minimise_edge(&x1, &x2), | |
376 | self.minimise_edge(&x2, &x3), | |
377 | self.minimise_edge(&x3, &x0)]; | |
378 | ||
379 | for edge in more_edge { | |
380 | if edge.1 < min_edge.1 { | |
381 | min_edge = edge; | |
382 | } | |
383 | } | |
384 | ||
385 | min_edge | |
386 | } | |
387 | } | |
388 | ||
389 | // impl<F : Float> P2Model<F, 1> for Cube<F, 1> { | |
390 | // type Model = CubeP2LocalModel<F,1,1,1>; | |
391 | ||
392 | // fn p2_model<G : FnMut(&Loc<F, 1>) -> F>(&self, fun : &G) -> Self::Model { | |
393 | // CubeP2LocalModel([Simplex(self.corners()).p2_model(fun)]) | |
394 | // } | |
395 | // } | |
396 | ||
397 | // impl<F : Float> P2Model<F, 2> for Cube<F, 2> { | |
398 | // type Model = CubeP2LocalModel<F,2,3,2>; | |
399 | ||
400 | // fn p2_model<G : FnMut(&Loc<F, 2>) -> F>(&self, fun : &G) -> Self::Model { | |
401 | // let [a, b, c, d] = self.corners(); | |
402 | // CubeP2LocalModel([Simplex([a, b, c]).p2_model(fun), | |
403 | // Simplex([b, c, d]).p2_model(fun)]) | |
404 | // } | |
405 | // } | |
406 | ||
407 | ||
408 | #[cfg(test)] | |
409 | mod tests { | |
410 | use super::*; | |
411 | ||
412 | #[test] | |
413 | fn p2_model_1d_test() { | |
414 | let vertices = [Loc([0.0]), Loc([1.0])]; | |
415 | let domain = Simplex(vertices); | |
416 | // A simple quadratic function for which the approximation is exact on reals, | |
417 | // and appears exact on f64 as well. | |
418 | let f = |&Loc([x]) : &Loc<f64, 1>| x*x + x + 1.0; | |
419 | let model = domain.p2_model(f); | |
420 | let xs = [Loc([0.5]), Loc([0.25])]; | |
421 | ||
422 | for x in vertices.iter().chain(xs.iter()) { | |
423 | assert_eq!(model.value(&x), f(&x)); | |
424 | } | |
425 | ||
426 | assert_eq!(model.minimise(&domain), (Loc([0.0]), 1.0)); | |
427 | } | |
428 | ||
429 | #[test] | |
430 | fn p2_model_2d_test() { | |
431 | let vertices = [Loc([0.0, 0.0]), Loc([1.0, 0.0]), Loc([1.0, 1.0])]; | |
432 | let domain = Simplex(vertices); | |
433 | // A simple quadratic function for which the approximation is exact on reals, | |
434 | // and appears exact on f64 as well. | |
1
df3901ec2f5d
Fix some unit tests after fundamental changes that made them invalid
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
435 | let f = |&Loc([x, y]) : &Loc<f64, 2>| - (x*x + x*y + x - 2.0 * y + 1.0); |
0 | 436 | let model = domain.p2_model(f); |
437 | let xs = [Loc([0.5, 0.5]), Loc([0.25, 0.25])]; | |
438 | ||
439 | for x in vertices.iter().chain(xs.iter()) { | |
440 | assert_eq!(model.value(&x), f(&x)); | |
441 | } | |
442 | ||
1
df3901ec2f5d
Fix some unit tests after fundamental changes that made them invalid
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
443 | assert_eq!(model.minimise(&domain), (Loc([1.0, 0.0]), -3.0)); |
0 | 444 | } |
445 | } |