Mercurial > repos > alg_tools / file comparison

comparison: src/fe_model/p2_local_model.rs

src/fe_model/p2_local_model.rs

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

Mercurial Repository: alg_tools

mercurial