diff -r 495448cca603 -r 6aa955ad8122 src/fe_model/p2_local_model.rs
--- a/src/fe_model/p2_local_model.rs	Thu May 01 08:40:33 2025 -0500
+++ b/src/fe_model/p2_local_model.rs	Thu May 01 13:06:58 2025 -0500
@@ -2,21 +2,21 @@
 Second order polynomical (P2) models on real intervals and planar 2D simplices.
 */
 
-use crate::types::*;
-use crate::loc::Loc;
-use crate::sets::{Set,NPolygon,SpannedHalfspace};
-use crate::linsolve::*;
+use super::base::{LocalModel, RealLocalModel};
 use crate::euclidean::Euclidean;
 use crate::instance::Instance;
-use super::base::{LocalModel,RealLocalModel};
+use crate::linsolve::*;
+use crate::loc::Loc;
 use crate::sets::Cube;
+use crate::sets::{NPolygon, Set, SpannedHalfspace};
+use crate::types::*;
 use numeric_literals::replace_float_literals;
 
 /// Type for simplices of arbitrary dimension `N`.
 ///
 /// The type parameter `D` indicates the number of nodes. (Rust's const generics do not currently
 /// allow its automatic calculation from `N`.)
-pub struct Simplex<F : Float, const N : usize, const D : usize>(pub [Loc<F, N>; D]);
+pub struct Simplex<F: Float, const N: usize, const D: usize>(pub [Loc<N, F>; D]);
 /// A two-dimensional planar simplex
 pub type PlanarSimplex<F> = Simplex<F, 2, 3>;
 /// A real interval
@@ -25,26 +25,29 @@
 /// Calculates (a+b)/2
 #[inline]
 #[replace_float_literals(F::cast_from(literal))]
-pub(crate) fn midpoint<F : Float, const N : usize>(a : &Loc<F,N>, b : &Loc<F,N>) -> Loc<F, N> {
-    (a+b)/2.0
+pub(crate) fn midpoint<F: Float, const N: usize>(a: &Loc<N, F>, b: &Loc<N, F>) -> Loc<N, F> {
+    (a + b) / 2.0
 }
 
-impl<'a, F : Float> Set<Loc<F,1>> for RealInterval<F> {
+impl<'a, F: Float> Set<Loc<1, F>> for RealInterval<F> {
     #[inline]
-    fn contains<I : Instance<Loc<F, 1>>>(&self, z : I) -> bool {
+    fn contains<I: Instance<Loc<1, F>>>(&self, z: I) -> bool {
         let &Loc([x]) = z.ref_instance();
         let &[Loc([x0]), Loc([x1])] = &self.0;
         (x0 < x && x < x1) || (x1 < x && x < x0)
     }
 }
 
-impl<'a, F : Float> Set<Loc<F,2>> for PlanarSimplex<F> {
+impl<'a, F: Float> Set<Loc<2, F>> for PlanarSimplex<F> {
     #[inline]
-    fn contains<I : Instance<Loc<F, 2>>>(&self, z : I) -> bool {
+    fn contains<I: Instance<Loc<2, F>>>(&self, z: I) -> bool {
         let &[x0, x1, x2] = &self.0;
-        NPolygon([[x0, x1].spanned_halfspace(),
-                  [x1, x2].spanned_halfspace(),
-                  [x2, x0].spanned_halfspace()]).contains(z)
+        NPolygon([
+            [x0, x1].spanned_halfspace(),
+            [x1, x2].spanned_halfspace(),
+            [x2, x0].spanned_halfspace(),
+        ])
+        .contains(z)
     }
 }
 
@@ -58,121 +61,118 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> P2Powers for Loc<F, 1> {
-    type Output = Loc<F, 1>;
-    type Full = Loc<F, 3>;
-    type Diff = Loc<Loc<F, 1>, 1>;
+impl<F: Float> P2Powers for Loc<1, F> {
+    type Output = Loc<1, F>;
+    type Full = Loc<3, F>;
+    type Diff = Loc<1, Loc<1, F>>;
 
     #[inline]
     fn p2powers(&self) -> Self::Output {
         let &Loc([x0]) = self;
-        [x0*x0].into()
+        [x0 * x0].into()
     }
 
     #[inline]
     fn p2powers_full(&self) -> Self::Full {
         let &Loc([x0]) = self;
-        [1.0, x0, x0*x0].into()
+        [1.0, x0, x0 * x0].into()
     }
 
     #[inline]
     fn p2powers_diff(&self) -> Self::Diff {
         let &Loc([x0]) = self;
-        [[x0+x0].into()].into()
+        [[x0 + x0].into()].into()
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> P2Powers for Loc<F, 2> {
-    type Output = Loc<F, 3>;
-    type Full = Loc<F, 6>;
-    type Diff = Loc<Loc<F, 3>, 2>;
+impl<F: Float> P2Powers for Loc<2, F> {
+    type Output = Loc<3, F>;
+    type Full = Loc<6, F>;
+    type Diff = Loc<2, Loc<3, F>>;
 
     #[inline]
     fn p2powers(&self) -> Self::Output {
         let &Loc([x0, x1]) = self;
-        [x0*x0, x0*x1, x1*x1].into()
+        [x0 * x0, x0 * x1, x1 * x1].into()
     }
 
     #[inline]
     fn p2powers_full(&self) -> Self::Full {
         let &Loc([x0, x1]) = self;
-        [1.0, x0, x1, x0*x0, x0*x1, x1*x1].into()
+        [1.0, x0, x1, x0 * x0, x0 * x1, x1 * x1].into()
     }
 
     #[inline]
     fn p2powers_diff(&self) -> Self::Diff {
         let &Loc([x0, x1]) = self;
-        [[x0+x0, x1, 0.0].into(), [0.0, x0, x1+x1].into()].into()
+        [[x0 + x0, x1, 0.0].into(), [0.0, x0, x1 + x1].into()].into()
     }
 }
 
 /// A trait for generating second order polynomial model of dimension `N` on `Self´.
 ///
 /// `Self` should present a subset aset of elements of the type [`Loc`]`<F, N>`.
-pub trait P2Model<F : Num, const N : usize> {
+pub trait P2Model<F: Num, const N: usize> {
     /// Implementation type of the second order polynomical model.
     /// Typically a [`P2LocalModel`].
-    type Model : LocalModel<Loc<F,N>,F>;
+    type Model: LocalModel<Loc<N, F>, F>;
     /// Generates a second order polynomial model of the function `g` on `Self`.
-    fn p2_model<G : Fn(&Loc<F, N>) -> F>(&self, g : G) -> Self::Model;
+    fn p2_model<G: Fn(&Loc<N, F>) -> F>(&self, g: G) -> Self::Model;
 }
 
 /// A local second order polynomical model of dimension `N` with `E` edges
-pub struct P2LocalModel<F : Num, const N : usize, const E : usize/*, const V : usize, const Q : usize*/> {
-    a0 : F,
-    a1 : Loc<F, N>,
-    a2 :  Loc<F, E>,
-    //node_values : Loc<F, V>,
-    //edge_values : Loc<F, Q>,
+pub struct P2LocalModel<
+    F: Num,
+    const N: usize,
+    const E: usize, /*, const V : usize, const Q : usize*/
+> {
+    a0: F,
+    a1: Loc<N, F>,
+    a2: Loc<E, F>,
+    //node_values : Loc<V, F>,
+    //edge_values : Loc<Q, F>,
 }
 
 //
 // 1D planar model construction
 //
 
-impl<F : Float> RealInterval<F> {
+impl<F: Float> RealInterval<F> {
     #[inline]
-    fn midpoints(&self) -> [Loc<F, 1>; 1] {
+    fn midpoints(&self) -> [Loc<1, F>; 1] {
         let [ref n0, ref n1] = &self.0;
         let n01 = midpoint(n0, n1);
         [n01]
     }
 }
 
-impl<F : Float> P2LocalModel<F, 1, 1/*, 2, 0*/> {
+impl<F: Float> P2LocalModel<F, 1, 1 /*, 2, 0*/> {
     /// Creates a new 1D second order polynomical model based on three nodal coordinates and
     /// corresponding function values.
     #[inline]
-    pub fn new(
-        &[n0, n1, n01] : &[Loc<F, 1>; 3],
-        &[v0, v1, v01] : &[F; 3],
-    ) -> Self {
-        let p = move |x : &Loc<F, 1>, v : F| {
+    pub fn new(&[n0, n1, n01]: &[Loc<1, F>; 3], &[v0, v1, v01]: &[F; 3]) -> Self {
+        let p = move |x: &Loc<1, F>, v: F| {
             let Loc([c, d, e]) = x.p2powers_full();
             [c, d, e, v]
         };
-        let [a0, a1, a11] = linsolve([
-            p(&n0, v0),
-            p(&n1, v1),
-            p(&n01, v01)
-        ]);
+        let [a0, a1, a11] = linsolve([p(&n0, v0), p(&n1, v1), p(&n01, v01)]);
         P2LocalModel {
-            a0 : a0,
-            a1 : [a1].into(),
-            a2 : [a11].into(),
+            a0: a0,
+            a1: [a1].into(),
+            a2: [a11].into(),
             //node_values : [v0, v1].into(),
             //edge_values: [].into(),
         }
     }
 }
 
-impl<F : Float> P2Model<F,1> for RealInterval<F> {
-    type Model = P2LocalModel<F, 1, 1/*, 2, 0*/>;
+impl<F: Float> P2Model<F, 1> for RealInterval<F> {
+    type Model = P2LocalModel<F, 1, 1 /*, 2, 0*/>;
 
     #[inline]
-    fn p2_model<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> Self::Model {
-        let [n01] =  self.midpoints();
+    fn p2_model<G: Fn(&Loc<1, F>) -> F>(&self, g: G) -> Self::Model {
+        let [n01] = self.midpoints();
         let [n0, n1] = self.0;
         let vals = [g(&n0), g(&n1), g(&n01)];
         let nodes = [n0, n1, n01];
@@ -184,10 +184,10 @@
 // 2D planar model construction
 //
 
-impl<F : Float> PlanarSimplex<F> {
+impl<F: Float> PlanarSimplex<F> {
     #[inline]
     /// Returns the midpoints of all the edges of the simplex
-    fn midpoints(&self) -> [Loc<F, 2>; 3] {
+    fn midpoints(&self) -> [Loc<2, F>; 3] {
         let [ref n0, ref n1, ref n2] = &self.0;
         let n01 = midpoint(n0, n1);
         let n12 = midpoint(n1, n2);
@@ -196,15 +196,15 @@
     }
 }
 
-impl<F : Float> P2LocalModel<F, 2, 3/*, 3, 3*/> {
+impl<F: Float> P2LocalModel<F, 2, 3 /*, 3, 3*/> {
     /// Creates a new 2D second order polynomical model based on six nodal coordinates and
     /// corresponding function values.
     #[inline]
     pub fn new(
-        &[n0, n1, n2, n01, n12, n20] : &[Loc<F, 2>; 6],
-        &[v0, v1, v2, v01, v12, v20] : &[F; 6],
+        &[n0, n1, n2, n01, n12, n20]: &[Loc<2, F>; 6],
+        &[v0, v1, v2, v01, v12, v20]: &[F; 6],
     ) -> Self {
-        let p = move |x : &Loc<F,2>, v :F| {
+        let p = move |x: &Loc<2, F>, v: F| {
             let Loc([c, d, e, f, g, h]) = x.p2powers_full();
             [c, d, e, f, g, h, v]
         };
@@ -217,20 +217,20 @@
             p(&n20, v20),
         ]);
         P2LocalModel {
-            a0 : a0,
-            a1 : [a1, a2].into(),
-            a2 : [a11, a12, a22].into(),
+            a0: a0,
+            a1: [a1, a2].into(),
+            a2: [a11, a12, a22].into(),
             //node_values : [v0, v1, v2].into(),
             //edge_values: [v01, v12, v20].into(),
         }
     }
 }
 
-impl<F : Float> P2Model<F,2> for PlanarSimplex<F> {
-    type Model = P2LocalModel<F, 2, 3/*, 3, 3*/>;
+impl<F: Float> P2Model<F, 2> for PlanarSimplex<F> {
+    type Model = P2LocalModel<F, 2, 3 /*, 3, 3*/>;
 
     #[inline]
-    fn p2_model<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> Self::Model {
+    fn p2_model<G: Fn(&Loc<2, F>) -> F>(&self, g: G) -> Self::Model {
         let midpoints = self.midpoints();
         let [ref n0, ref n1, ref n2] = self.0;
         let [ref n01, ref n12, ref n20] = midpoints;
@@ -242,125 +242,137 @@
 
 macro_rules! impl_local_model {
     ($n:literal, $e:literal, $v:literal, $q:literal) => {
-        impl<F : Float> LocalModel<Loc<F, $n>, F> for P2LocalModel<F, $n, $e/*, $v, $q*/> {
+        impl<F: Float> LocalModel<Loc<$n, F>, F> for P2LocalModel<F, $n, /*, $v, $q*/ $e> {
             #[inline]
-            fn value(&self, x : &Loc<F,$n>) -> F {
+            fn value(&self, x: &Loc<$n, F>) -> F {
                 self.a0 + x.dot(&self.a1) + x.p2powers().dot(&self.a2)
             }
 
             #[inline]
-            fn differential(&self, x : &Loc<F,$n>) -> Loc<F,$n> {
+            fn differential(&self, x: &Loc<$n, F>) -> Loc<$n, F> {
                 self.a1 + x.p2powers_diff().map(|di| di.dot(&self.a2))
             }
         }
-    }
+    };
 }
 
 impl_local_model!(1, 1, 2, 0);
 impl_local_model!(2, 3, 3, 3);
 
-
 //
 // Minimisation
 //
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> P2LocalModel<F, 1, 1/*, 2, 0*/> {
+impl<F: Float> P2LocalModel<F, 1, 1 /*, 2, 0*/> {
     /// Minimises the model along the edge `[x0, x1]`.
     #[inline]
-    fn minimise_edge(&self, x0 : Loc<F, 1>, x1 : Loc<F,1>) -> (Loc<F,1>, F) {
-        let &P2LocalModel{
-            a1 : Loc([a1]),
-            a2 : Loc([a11]),
+    fn minimise_edge(&self, x0: Loc<1, F>, x1: Loc<1, F>) -> (Loc<1, F>, F) {
+        let &P2LocalModel {
+            a1: Loc([a1]),
+            a2: Loc([a11]),
             //node_values : Loc([v0, v1]),
             ..
-         } = self;
+        } = self;
         // We do this in cases, first trying for an interior solution, then edges.
         // For interior solution, first check determinant; no point trying if non-positive
         if a11 > 0.0 {
             // An interior solution x[1] has to satisfy
             //   2a₁₁*x[1] + a₁ =0
             // This gives
-            let t = -a1/(2.0*a11);
+            let t = -a1 / (2.0 * a11);
             let (Loc([t0]), Loc([t1])) = (x0, x1);
             if (t0 <= t && t <= t1) || (t1 <= t && t <= t0) {
                 let x = [t].into();
                 let v = self.value(&x);
-                return (x, v)
+                return (x, v);
             }
         }
 
         let v0 = self.value(&x0);
         let v1 = self.value(&x1);
-        if v0 < v1 { (x0, v0) } else { (x1, v1) }
+        if v0 < v1 {
+            (x0, v0)
+        } else {
+            (x1, v1)
+        }
     }
 }
 
-impl<'a, F : Float> RealLocalModel<RealInterval<F>,Loc<F,1>,F>
-for P2LocalModel<F, 1, 1/*, 2, 0*/> {
+impl<'a, F: Float> RealLocalModel<RealInterval<F>, Loc<1, F>, F>
+    for P2LocalModel<F, 1, 1 /*, 2, 0*/>
+{
     #[inline]
-    fn minimise(&self, &Simplex([x0, x1]) : &RealInterval<F>) -> (Loc<F,1>, F) {
+    fn minimise(&self, &Simplex([x0, x1]): &RealInterval<F>) -> (Loc<1, F>, F) {
         self.minimise_edge(x0, x1)
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> P2LocalModel<F, 2, 3/*, 3, 3*/> {
+impl<F: Float> P2LocalModel<F, 2, 3 /*, 3, 3*/> {
     /// Minimise the 2D model along the edge `[x0, x1] = {x0 + t(x1 - x0) | t ∈ [0, 1] }`.
     #[inline]
-    fn minimise_edge(&self, x0 : &Loc<F,2>, x1 : &Loc<F,2>/*, v0 : F, v1 : F*/) -> (Loc<F,2>, F) {
+    fn minimise_edge(
+        &self,
+        x0: &Loc<2, F>,
+        x1: &Loc<2, F>, /*, v0 : F, v1 : F*/
+    ) -> (Loc<2, F>, F) {
         let &P2LocalModel {
             a0,
-            a1 : Loc([a1, a2]),
-            a2 : Loc([a11, a12, a22]),
+            a1: Loc([a1, a2]),
+            a2: Loc([a11, a12, a22]),
             ..
-         } = self;
+        } = self;
         let &Loc([x00, x01]) = x0;
-        let d@Loc([d0, d1]) = x1 - x0;
-        let b0 = a0 + a1*x00 + a2*x01 + a11*x00*x00 + a12*x00*x01 + a22*x01*x01;
-        let b1 = a1*d0 + a2*d1 + 2.0*a11*d0*x00 + a12*(d0*x01 + d1*x00) + 2.0*a22*d1*x01;
-        let b11 = a11*d0*d0 + a12*d0*d1 + a22*d1*d1;
+        let d @ Loc([d0, d1]) = x1 - x0;
+        let b0 = a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + a12 * x00 * x01 + a22 * x01 * x01;
+        let b1 = a1 * d0
+            + a2 * d1
+            + 2.0 * a11 * d0 * x00
+            + a12 * (d0 * x01 + d1 * x00)
+            + 2.0 * a22 * d1 * x01;
+        let b11 = a11 * d0 * d0 + a12 * d0 * d1 + a22 * d1 * d1;
         let edge_1d_model = P2LocalModel {
-            a0 : b0,
-            a1 : Loc([b1]),
-            a2 : Loc([b11]),
+            a0: b0,
+            a1: Loc([b1]),
+            a2: Loc([b11]),
             //node_values : Loc([v0, v1]),
         };
         let (Loc([t]), v) = edge_1d_model.minimise_edge(0.0.into(), 1.0.into());
-        (x0 + d*t, v)
+        (x0 + d * t, v)
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float> RealLocalModel<PlanarSimplex<F>,Loc<F,2>,F>
-for P2LocalModel<F, 2, 3/*, 3, 3*/> {
+impl<'a, F: Float> RealLocalModel<PlanarSimplex<F>, Loc<2, F>, F>
+    for P2LocalModel<F, 2, 3 /*, 3, 3*/>
+{
     #[inline]
-    fn minimise(&self, el : &PlanarSimplex<F>) -> (Loc<F,2>, F) {
+    fn minimise(&self, el: &PlanarSimplex<F>) -> (Loc<2, F>, F) {
         let &P2LocalModel {
-            a1 : Loc([a1, a2]),
-            a2 : Loc([a11, a12, a22]),
+            a1: Loc([a1, a2]),
+            a2: Loc([a11, a12, a22]),
             //node_values : Loc([v0, v1, v2]),
             ..
-         } = self;
+        } = self;
 
         // We do this in cases, first trying for an interior solution, then edges.
         // For interior solution, first check determinant; no point trying if non-positive
-        let r = 2.0*(a11*a22-a12*a12);
+        let r = 2.0 * (a11 * a22 - a12 * a12);
         if r > 0.0 {
             // An interior solution (x[1], x[2]) has to satisfy
             //   2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0
             // This gives
-            let x = [(a22*a1-a12*a2)/r, (a12*a1-a11*a2)/r].into();
+            let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into();
             if el.contains(&x) {
-                return (x, self.value(&x))
+                return (x, self.value(&x));
             }
         }
 
         let &[ref x0, ref x1, ref x2] = &el.0;
         let mut min_edge = self.minimise_edge(x0, x1);
-        let more_edge = [self.minimise_edge(x1, x2),
-                         self.minimise_edge(x2, x0)];
-        
+        let more_edge = [self.minimise_edge(x1, x2), self.minimise_edge(x2, x0)];
+
         for edge in more_edge {
             if edge.1 < min_edge.1 {
                 min_edge = edge;
@@ -372,35 +384,38 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float> RealLocalModel<Cube<F, 2>,Loc<F,2>,F>
-for P2LocalModel<F, 2, 3/*, 3, 3*/> {
+impl<'a, F: Float> RealLocalModel<Cube<2, F>, Loc<2, F>, F>
+    for P2LocalModel<F, 2, 3 /*, 3, 3*/>
+{
     #[inline]
-    fn minimise(&self, el : &Cube<F, 2>) -> (Loc<F,2>, F) {
+    fn minimise(&self, el: &Cube<2, F>) -> (Loc<2, F>, F) {
         let &P2LocalModel {
-            a1 : Loc([a1, a2]),
-            a2 : Loc([a11, a12, a22]),
+            a1: Loc([a1, a2]),
+            a2: Loc([a11, a12, a22]),
             //node_values : Loc([v0, v1, v2]),
             ..
-         } = self;
+        } = self;
 
         // We do this in cases, first trying for an interior solution, then edges.
         // For interior solution, first check determinant; no point trying if non-positive
-        let r = 2.0*(a11*a22-a12*a12);
+        let r = 2.0 * (a11 * a22 - a12 * a12);
         if r > 0.0 {
             // An interior solution (x[1], x[2]) has to satisfy
             //   2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0
             // This gives
-            let x = [(a22*a1-a12*a2)/r, (a12*a1-a11*a2)/r].into();
+            let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into();
             if el.contains(&x) {
-                return (x, self.value(&x))
+                return (x, self.value(&x));
             }
         }
 
         let [x0, x1, x2, x3] = el.corners();
         let mut min_edge = self.minimise_edge(&x0, &x1);
-        let more_edge = [self.minimise_edge(&x1, &x2),
-                         self.minimise_edge(&x2, &x3),
-                         self.minimise_edge(&x3, &x0)];
+        let more_edge = [
+            self.minimise_edge(&x1, &x2),
+            self.minimise_edge(&x2, &x3),
+            self.minimise_edge(&x3, &x0),
+        ];
 
         for edge in more_edge {
             if edge.1 < min_edge.1 {
@@ -422,7 +437,7 @@
         let domain = Simplex(vertices);
         // A simple quadratic function for which the approximation is exact on reals,
         // and appears exact on f64 as well.
-        let f = |&Loc([x]) : &Loc<f64, 1>| x*x + x + 1.0;
+        let f = |&Loc([x]): &Loc<1, f64>| x * x + x + 1.0;
         let model = domain.p2_model(f);
         let xs = [Loc([0.5]), Loc([0.25])];
 
@@ -439,7 +454,7 @@
         let domain = Simplex(vertices);
         // A simple quadratic function for which the approximation is exact on reals,
         // and appears exact on f64 as well.
-        let f = |&Loc([x, y]) : &Loc<f64, 2>| - (x*x + x*y + x - 2.0 * y + 1.0);
+        let f = |&Loc([x, y]): &Loc<2, f64>| -(x * x + x * y + x - 2.0 * y + 1.0);
         let model = domain.p2_model(f);
         let xs = [Loc([0.5, 0.5]), Loc([0.25, 0.25])];