Mercurial > repos > alg_tools / file comparison

comparison: src/fe_model/p2_local_model.rs

src/fe_model/p2_local_model.rs

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

Mercurial Repository: alg_tools

mercurial