comparison: src/fe_model/p2_local_model.rs
src/fe_model/p2_local_model.rs
- changeset 208
- 6be69f736c79
- parent 206
- ce37ff3ce507
- child 209
- 060891c3f537
equal
deleted
inserted
replaced
| 207:032957f8d6de | 208:6be69f736c79 |
|---|---|
| 6 use crate::euclidean::Euclidean; | 6 use crate::euclidean::Euclidean; |
| 7 use crate::instance::Instance; | 7 use crate::instance::Instance; |
| 8 use crate::linsolve::*; | 8 use crate::linsolve::*; |
| 9 use crate::loc::Loc; | 9 use crate::loc::Loc; |
| 10 use crate::sets::Cube; | 10 use crate::sets::Cube; |
| 11 use crate::sets::{NPolygon, Set, SpannedHalfspace}; | 11 use crate::sets::Set; |
| 12 use crate::types::*; | 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 /// |
| 100 type Diff = Loc<2, Loc<3, F>>; | 100 type Diff = Loc<2, Loc<3, F>>; |
| 101 | 101 |
| 102 #[inline] | 102 #[inline] |
| 103 fn p2powers(&self) -> Self::Output { | 103 fn p2powers(&self) -> Self::Output { |
| 104 let &Loc([x0, x1]) = self; | 104 let &Loc([x0, x1]) = self; |
| 105 [x0 * x0, x0 * x1, x1 * x1].into() | 105 [x0 * x0, 2.0 * x0 * x1, x1 * x1].into() |
| 106 } | 106 } |
| 107 | 107 |
| 108 #[inline] | 108 #[inline] |
| 109 fn p2powers_full(&self) -> Self::Full { | 109 fn p2powers_full(&self) -> Self::Full { |
| 110 let &Loc([x0, x1]) = self; | 110 let &Loc([x0, x1]) = self; |
| 111 [1.0, x0, x1, x0 * x0, x0 * x1, x1 * x1].into() | 111 [1.0, x0, x1, x0 * x0, 2.0 * x0 * x1, x1 * x1].into() |
| 112 } | 112 } |
| 113 | 113 |
| 114 #[inline] | 114 #[inline] |
| 115 fn p2powers_diff(&self) -> Self::Diff { | 115 fn p2powers_diff(&self) -> Self::Diff { |
| 116 let &Loc([x0, x1]) = self; | 116 let &Loc([x0, x1]) = self; |
| 326 x1: &Loc<2, F>, /*, v0 : F, v1 : F*/ | 326 x1: &Loc<2, F>, /*, v0 : F, v1 : F*/ |
| 327 ) -> (Loc<2, F>, F) { | 327 ) -> (Loc<2, F>, F) { |
| 328 let &P2LocalModel { a0, a1: Loc([a1, a2]), a2: Loc([a11, a12, a22]), .. } = self; | 328 let &P2LocalModel { a0, a1: Loc([a1, a2]), a2: Loc([a11, a12, a22]), .. } = self; |
| 329 let &Loc([x00, x01]) = x0; | 329 let &Loc([x00, x01]) = x0; |
| 330 let d @ Loc([d0, d1]) = x1 - x0; | 330 let d @ Loc([d0, d1]) = x1 - x0; |
| 331 let b0 = a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + a12 * x00 * x01 + a22 * x01 * x01; | 331 let b0 = |
| 332 a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + 2.0 * a12 * x00 * x01 + a22 * x01 * x01; | |
| 332 let b1 = a1 * d0 | 333 let b1 = a1 * d0 |
| 333 + a2 * d1 | 334 + a2 * d1 |
| 334 + 2.0 * a11 * d0 * x00 | 335 + 2.0 * a11 * d0 * x00 |
| 335 + a12 * (d0 * x01 + d1 * x00) | 336 + 2.0 * a12 * (d0 * x01 + d1 * x00) |
| 336 + 2.0 * a22 * d1 * x01; | 337 + 2.0 * a22 * d1 * x01; |
| 337 let b11 = a11 * d0 * d0 + a12 * d0 * d1 + a22 * d1 * d1; | 338 let b11 = a11 * d0 * d0 + 2.0 * a12 * d0 * d1 + a22 * d1 * d1; |
| 338 let edge_1d_model = P2LocalModel { | 339 let edge_1d_model = P2LocalModel { |
| 339 a0: b0, | 340 a0: b0, |
| 340 a1: Loc([b1]), | 341 a1: Loc([b1]), |
| 341 a2: Loc([b11]), | 342 a2: Loc([b11]), |
| 342 //node_values : Loc([v0, v1]), | 343 //node_values : Loc([v0, v1]), |
| 350 impl<'a, F: Float> RealLocalModel<PlanarSimplex<F>, Loc<2, F>, F> | 351 impl<'a, F: Float> RealLocalModel<PlanarSimplex<F>, Loc<2, F>, F> |
| 351 for P2LocalModel<F, 2, 3 /*, 3, 3*/> | 352 for P2LocalModel<F, 2, 3 /*, 3, 3*/> |
| 352 { | 353 { |
| 353 #[inline] | 354 #[inline] |
| 354 fn minimise(&self, el: &PlanarSimplex<F>) -> (Loc<2, F>, F) { | 355 fn minimise(&self, el: &PlanarSimplex<F>) -> (Loc<2, F>, F) { |
| 356 let &P2LocalModel { | |
| 357 a1: Loc([a1, a2]), | |
| 358 a2: Loc([a11, a12, a22]), | |
| 359 //node_values : Loc([v0, v1, v2]), | |
| 360 .. | |
| 361 } = self; | |
| 362 | |
| 363 // We do this in cases, first trying for an interior solution, then edges. | |
| 364 // For interior solution, first check determinant; no point trying if non-positive | |
| 365 let r = 2.0 * (a11 * a22 - a12 * a12); | |
| 366 if r > 0.0 { | |
| 367 // An interior solution (x[1], x[2]) has to satisfy | |
| 368 // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[2] + 2a₁₂*x[1]+a₂=0 | |
| 369 // This gives | |
| 370 let x = [(a22 * a1 - a12 * a2) / r, (a11 * a2 - a12 * a1) / r].into(); | |
| 371 if el.contains(&x) { | |
| 372 return (x, self.value(&x)); | |
| 373 } | |
| 374 } | |
| 375 | |
| 376 let &[ref x0, ref x1, ref x2] = &el.0; | |
| 377 let mut min_edge = self.minimise_edge(x0, x1); | |
| 378 let more_edge = [self.minimise_edge(x1, x2), self.minimise_edge(x2, x0)]; | |
| 379 | |
| 380 for edge in more_edge { | |
| 381 if edge.1 < min_edge.1 { | |
| 382 min_edge = edge; | |
| 383 } | |
| 384 } | |
| 385 | |
| 386 min_edge | |
| 387 } | |
| 388 } | |
| 389 | |
| 390 #[replace_float_literals(F::cast_from(literal))] | |
| 391 impl<'a, F: Float> RealLocalModel<Cube<2, F>, Loc<2, F>, F> | |
| 392 for P2LocalModel<F, 2, 3 /*, 3, 3*/> | |
| 393 { | |
| 394 #[inline] | |
| 395 fn minimise(&self, el: &Cube<2, F>) -> (Loc<2, F>, F) { | |
| 355 let &P2LocalModel { | 396 let &P2LocalModel { |
| 356 a1: Loc([a1, a2]), | 397 a1: Loc([a1, a2]), |
| 357 a2: Loc([a11, a12, a22]), | 398 a2: Loc([a11, a12, a22]), |
| 358 //node_values : Loc([v0, v1, v2]), | 399 //node_values : Loc([v0, v1, v2]), |
| 359 .. | 400 .. |
| 370 if el.contains(&x) { | 411 if el.contains(&x) { |
| 371 return (x, self.value(&x)); | 412 return (x, self.value(&x)); |
| 372 } | 413 } |
| 373 } | 414 } |
| 374 | 415 |
| 375 let &[ref x0, ref x1, ref x2] = &el.0; | |
| 376 let mut min_edge = self.minimise_edge(x0, x1); | |
| 377 let more_edge = [self.minimise_edge(x1, x2), self.minimise_edge(x2, 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 #[replace_float_literals(F::cast_from(literal))] | |
| 390 impl<'a, F: Float> RealLocalModel<Cube<2, F>, Loc<2, F>, F> | |
| 391 for P2LocalModel<F, 2, 3 /*, 3, 3*/> | |
| 392 { | |
| 393 #[inline] | |
| 394 fn minimise(&self, el: &Cube<2, F>) -> (Loc<2, F>, F) { | |
| 395 let &P2LocalModel { | |
| 396 a1: Loc([a1, a2]), | |
| 397 a2: Loc([a11, a12, a22]), | |
| 398 //node_values : Loc([v0, v1, v2]), | |
| 399 .. | |
| 400 } = self; | |
| 401 | |
| 402 // We do this in cases, first trying for an interior solution, then edges. | |
| 403 // For interior solution, first check determinant; no point trying if non-positive | |
| 404 let r = 2.0 * (a11 * a22 - a12 * a12); | |
| 405 if r > 0.0 { | |
| 406 // An interior solution (x[1], x[2]) has to satisfy | |
| 407 // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 | |
| 408 // This gives | |
| 409 let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into(); | |
| 410 if el.contains(&x) { | |
| 411 return (x, self.value(&x)); | |
| 412 } | |
| 413 } | |
| 414 | |
| 415 let [x0, x1, x2, x3] = el.corners(); | 416 let [x0, x1, x2, x3] = el.corners(); |
| 416 let mut min_edge = self.minimise_edge(&x0, &x1); | 417 let mut min_edge = self.minimise_edge(&x0, &x1); |
| 417 let more_edge = [ | 418 let more_edge = [ |
| 418 self.minimise_edge(&x1, &x2), | 419 self.minimise_edge(&x1, &x2), |
| 419 self.minimise_edge(&x2, &x3), | 420 self.minimise_edge(&x2, &x3), |
