Mercurial > repos > non-riemannian-opt / changeset

--- a/Cargo.toml	Thu Nov 07 16:52:05 2024 -0500
+++ b/Cargo.toml	Wed Dec 04 23:19:46 2024 -0500
@@ -23,3 +23,5 @@
 colored = "~2.0.0"
 image = "~0.24.3"
 serde_repr = "0.1"
+
+

--- a/src/cube.rs	Thu Nov 07 16:52:05 2024 -0500
+++ b/src/cube.rs	Wed Dec 04 23:19:46 2024 -0500
@@ -6,7 +6,7 @@
 use serde::Serialize;
 use alg_tools::loc::Loc;
 use alg_tools::norms::{Norm, L2};
-use crate::manifold::{ManifoldPoint, EmbeddedManifoldPoint};
+use crate::manifold::{EmbeddedManifoldPoint, FacedManifoldPoint, ManifoldPoint};
 
 /// All the difference faces of a [`OnCube`].
 #[derive(Copy, Clone, Debug, Eq, PartialEq, Serialize_repr, Deserialize_repr)]
@@ -292,9 +292,12 @@
         }
         (best_tan, best_len)
     }
+}
 
+impl FacedManifoldPoint for OnCube {
+    type Face = Face;
     /// Returns the face of this point.
-    pub fn face(&self) -> Face {
+    fn face(&self) -> Face {
         self.face
     }
 }

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/cylinder.rs	Wed Dec 04 23:19:46 2024 -0500
@@ -0,0 +1,932 @@
+/*!
+Implementation of the surface of a 3D cylinder as a [`ManifoldPoint`].
+*/
+
+use alg_tools::euclidean::Euclidean;
+use serde_repr::*;
+use serde::{Serialize, Deserialize};
+use alg_tools::loc::Loc;
+use alg_tools::norms::{Norm, L2};
+use alg_tools::types::Float;
+use crate::manifold::{ManifoldPoint, EmbeddedManifoldPoint, FacedManifoldPoint};
+use crate::newton::{newton_sym1x1, newton_sym2x2};
+
+/// Angle
+pub type Angle = f64;
+
+/// Cylindrical coordinates in ℝ^3
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct CylCoords {
+    pub r : f64,
+    pub angle : Angle,
+    pub z : f64
+}
+
+impl CylCoords {
+    #[inline]
+    pub fn to_cartesian(&self) -> Loc<f64, 3> {
+        let &CylCoords{r, angle, z} = self;
+        [r * angle.cos(), r * angle.sin(), z].into()
+    }
+
+    #[inline]
+    #[allow(dead_code)]
+    pub fn from_cartesian(coords : Loc<f64, 3>) -> Self {
+        let [x, y, z] = coords.into();
+        let r = (x*x + y*y).sqrt();
+        let angle = y.atan2(x);
+        CylCoords {r, angle, z}
+    }
+}
+
+/// Coordinates on a cap
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct CapPoint {
+    pub r : f64,
+    pub angle : Angle
+}
+
+#[inline]
+fn rotate(φ : f64, Loc([x, y]) : Loc<f64, 2>) -> Loc<f64, 2> {
+    let sin_φ = φ.sin();
+    let cos_φ = φ.cos();
+    [cos_φ * x - sin_φ * y, sin_φ * x + cos_φ * y].into()
+}
+
+impl CapPoint {
+    #[inline]
+    /// Convert to cylindrical coordinates given z coordinate
+    fn cyl_coords(&self, z : f64) -> CylCoords {
+        let CapPoint { r, angle } = *self;
+        CylCoords { r, angle, z }
+    }
+
+    #[inline]
+    /// Convert to cylindrical coordinates given z coordinate
+    fn cartesian_coords(&self) -> Loc<f64, 2> {
+        let CapPoint { r, angle } = *self;
+        [r * angle.cos(), r * angle.sin()].into()
+    }
+
+    #[inline]
+    #[allow(dead_code)]
+    /// Convert to cylindrical coordinates given z coordinate
+    fn from_cartesian(coords : Loc<f64, 2>) -> Self {
+        let [x, y] = coords.into();
+        let r = (x*x + y*y).sqrt();
+        let angle = y.atan2(x);
+        CapPoint { r, angle }
+    }
+
+    #[inline]
+    /// Calculate the vector between two points on the cap
+    fn log(&self, other : &CapPoint) -> Loc<f64, 2> {
+        other.cartesian_coords() - self.cartesian_coords()
+    }
+
+    #[inline]
+    /// Calculate partial exponential map until boundary.
+    /// Returns the final point within the cap as well as a remaining tangent on
+    /// the side of the cylinder, if `t` wasn't fully used.
+    fn partial_exp(&self, r : f64, t : &Tangent) -> (CapPoint, Option<Tangent>) {
+        let q = self.cartesian_coords() + t;
+        let n = q.norm2();
+        if n <= r {
+            (Self::from_cartesian(q), None)
+        } else {
+            let p = q * r / n;
+            (Self::from_cartesian(p), Some(q - p))
+        }
+    }
+
+    #[inline]
+    /// Convert tangent from side tangent to cap tangent
+    fn tangent_from_side(&self, top : bool, t : Tangent) -> Tangent {
+        if top {
+            // The angle is such that down would be rotated to self.angle, counterclockwise
+            // The new tangent is R[down +  t] - R[down] = Rt.
+            rotate(self.angle+f64::PI/2.0, t)
+        } else {
+            // The angle is such that up would be rotated to self.angle, clockwise
+            rotate(self.angle-f64::PI/2.0, t)
+        }
+    }
+
+    #[inline]
+    /// Convert tangent from cap tangent to tangent tangent
+    fn tangent_to_side(&self, top : bool, t : Tangent) -> Tangent {
+        if top {
+            // The angle is such that self.angle would be rotated to down, clockwise
+            rotate(-self.angle-f64::PI/2.0, t)
+        } else {
+            // The angle is such that self.angle would be rotated to up, counterclockwise
+            rotate(f64::PI/2.0-self.angle, t)
+        }
+    }
+}
+
+/// Coordinates on a side
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct SidePoint {
+    pub z : f64,
+    pub angle : Angle
+}
+
+#[inline]
+fn anglediff(mut φ1 : f64, mut φ2 : f64) -> f64 {
+    let π = f64::PI;
+    φ1 = normalise_angle(φ1);
+    φ2 = normalise_angle(φ2);
+    let α = φ2 - φ1;
+    if α >= 0.0 {
+        if α <= π {
+            α
+        } else {
+            α - 2.0 * π
+        }
+    } else {
+        if α >= -π {
+            α
+        } else {
+            2.0 * π + α
+        }
+    }
+}
+
+#[inline]
+pub fn normalise_angle(φ : f64) -> f64 {
+    let π = f64::PI;
+    φ.rem_euclid(2.0 * π)
+}
+
+impl SidePoint {
+    #[inline]
+    /// Convert to cylindrical coordinates given radius
+    fn cyl_coords(&self, r : f64) -> CylCoords {
+        let SidePoint { z, angle } = *self;
+        CylCoords { r, angle, z }
+    }
+
+    #[inline]
+    /// Calculate tangent vector between two points on the side, given radius
+    fn log(&self, r : f64, other : &SidePoint) -> Loc<f64, 2> {
+        let SidePoint{ z : z1, angle : angle1 } = *self;
+        let SidePoint{ z : z2, angle : angle2 } = *other;
+        let φ = anglediff(angle1, angle2);
+        // TODO: is this correct?
+        [r*φ, z2-z1].into()
+    }
+
+    #[inline]
+    /// Calculate partial exponential map under boundary
+    /// Returns a point on the next face, as well as a remaining tangent on
+    /// the side of the cylinder, if `t` wasn't fully used.
+    fn partial_exp(&self, r : f64, (a, b) : (f64, f64), t : &Tangent)
+        -> (SidePoint, Option<Tangent>)
+    {
+        assert!(a <= self.z && self.z <= b);
+        let Loc([_, h]) = *t;
+        let s = if h > 0.0 {
+            ((b - self.z)/h).min(1.0)
+        } else if h < 0.0 {
+            ((a - self.z)/h).min(1.0)
+        } else {
+            1.0
+        };
+        let d = t * s;
+        let p = self.unflatten(r, d);
+        if s < 1.0 {
+            (p, Some(t - d))
+        } else {
+            (p, None)
+        }
+    }
+
+    #[inline]
+    /// Unflattens another point in the local coordinate system of self
+    fn unflatten(&self, r : f64, Loc([v, h]) : Loc<f64, 2>) -> Self {
+        SidePoint{ z : self.z + h, angle : normalise_angle(self.angle + v / r) }
+    }
+
+}
+
+/// Point on a [`Cylinder`]
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub enum Point {
+    Top(CapPoint),
+    Bottom(CapPoint),
+    Side(SidePoint),
+}
+
+/// Face on a [`Cylinder`]
+#[derive(Copy, Clone, Debug, Eq, PartialEq, Serialize_repr, Deserialize_repr)]
+#[repr(u8)]
+pub enum Face {
+    Top,
+    Bottom,
+    Side,
+}
+
+#[derive(Clone, Debug, PartialEq)]
+pub struct OnCylinder<'a> {
+    cylinder : &'a Cylinder,
+    point : Point,
+}
+
+
+/// Cylinder configuration
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct CylinderConfig {
+    pub newton_iters : usize,
+}
+
+impl Default for CylinderConfig {
+    fn default() -> Self {
+        CylinderConfig { newton_iters : 10 }
+    }
+}
+
+/// A cylinder
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Cylinder {
+    /// Radius of the cylinder
+    pub radius : f64,
+    /// Height of the cylinder
+    pub height : f64,
+    /// Configuration for numerical methods
+    pub config : CylinderConfig
+}
+
+
+impl std::fmt::Display for Face {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        use Face::*;
+        let s = match *self {
+            Top => "Top",
+            Bottom => "Bottom",
+            Side => "Side",
+        };
+        write!(f, "{}", s)
+    }
+}
+
+impl Point {
+    fn face(&self) -> Face {
+        match *self {
+            Point::Top(..) => Face::Top,
+            Point::Bottom(_) => Face::Bottom,
+            Point::Side(_) => Face::Side,
+        }
+    }
+}
+
+impl<'a> FacedManifoldPoint for OnCylinder<'a> {
+    type Face = Face;
+    /// Returns the face of this point.
+    fn face(&self) -> Face {
+        self.point.face()
+    }
+}
+
+// Tangent vector
+type Tangent = Loc<f64, 2>;
+
+#[inline]
+fn best_tangent<I>(tangents : I) -> (Tangent, f64)
+where I : IntoIterator<Item = Tangent> {
+    tangents.into_iter()
+            .map(|t| (t, t.norm(L2)))
+            .reduce(|a@(_, la), b@(_, lb)| if la < lb { a } else { b })
+            .unwrap()
+}
+
+/// Swap the elements of a two-tuple
+#[inline]
+fn swap<A>((a, b) : (A, A)) -> (A, A) {
+    (b, a)
+}
+
+#[inline]
+fn indistinguishable(a : f64, b : f64) -> bool {
+    a > b - f64::EPSILON && a < b + f64::EPSILON
+}
+
+impl Cylinder {
+    /// Return the cylindrical coordinates of `point` on this cylinder
+    fn cyl_coords(&self, point : &Point) -> CylCoords {
+        match *point {
+            Point::Top(cap)    => { cap.cyl_coords(self.top_z()) },
+            Point::Bottom(cap) => { cap.cyl_coords(self.bottom_z()) },
+            Point::Side(side)  => { side.cyl_coords(self.radius) },
+        }
+    }
+
+    #[inline]
+    pub fn top_z(&self) -> f64 {
+        self.height / 2.0
+    }
+
+    #[inline]
+    pub fn bottom_z(&self) -> f64 {
+        -self.height / 2.0
+    }
+
+    /// Find angle where a geodesic from `side` to `(cap, z)` crosses the cap edge.
+    ///
+    /// Uses `newton_sym1x1`.
+    fn side_cap_crossing(
+        &self,
+        side : &SidePoint,
+        cap : &CapPoint, z : f64, // `z` is the z coordinate of cap
+    ) -> Angle {
+        let &SidePoint { z : z_1, angle : φ_1 } = side;
+        let &CapPoint { r : r_2, angle : φ_2 } = cap;
+        assert_ne!(z, z_1);
+        let d_1 = z - z_1;
+        let d2_1 = d_1 * d_1;
+        let r = self.radius;
+        let r2 = r * r;
+        let r2_2 = r_2 * r_2;
+        let rr_2 = r * r_2;
+
+        let g = |α_1 : f64| {
+            let ψ = φ_2 - φ_1 - α_1;
+            let ψ_cos = ψ.cos();
+            let ψ_sin = ψ.sin();
+            let ψ_sin2 = ψ_sin * ψ_sin;
+            let ψ_cos2 = ψ_cos * ψ_cos;
+            let α2_1 = α_1 * α_1;
+            let c = d2_1 + r2 * α2_1;
+            let e = r2 + r2_2 - 2.0 * rr_2 * ψ_cos;
+            let g = r2 * α_1 / c.sqrt() - rr_2 * ψ_sin / e.sqrt();
+            let h = r2 * d2_1 / c.powf(3.0/2.0)
+                    + r * r_2 * ((r2 + r2_2) * ψ_cos - rr_2 * (ψ_sin2 - 2.0 * ψ_cos2))
+                        / e.powf(3.0/2.0);
+            (h, g)
+        };
+
+        let α_1 = newton_sym1x1(g, 0.0, self.config.newton_iters);
+        normalise_angle(φ_1 + α_1)
+
+    }
+
+    /// Find angles where the geodesic passing through a cap at height `z` from `side1` to `side2`
+    /// crosses the cap edge. **Panics if `side2.angle < side1.angle`.**
+    ///
+    /// Uses `newton_sym2x2`.
+    fn side_cap_side_crossing(
+        &self,
+        side1 : &SidePoint,
+        z : f64,
+        side2 : &SidePoint
+    ) -> (Angle, Angle) {
+        let &SidePoint { z : z_1, angle : φ_1 } = side1;
+        let &SidePoint { z : z_2, angle : φ_2 } = side2;
+        assert!(φ_2 >= φ_1);
+        assert_ne!(z_1, z);
+        assert_ne!(z_2, z);
+        let r = self.radius;
+        let r2 = r * r;
+        let d_1 = z - z_1;
+        let d_2 = z - z_2;
+        let d2_1 = d_1 * d_1;
+        let d2_2 = d_2 * d_2;
+        let g = |α_1 : f64, α_2 : f64| {
+            let α2_1 = α_1 * α_1;
+            let α2_2 = α_2 * α_2;
+            let ψ = (α_1 + α_2 + φ_1 - φ_2) / 2.0;
+            let ψ_sin = ψ.sin();
+            let ψ_cos = ψ.cos();
+            let c_1 = d2_1 + r2 * α2_1;
+            let c_2 = d2_2 + r2 * α2_2;
+            let g_1 = r2 * α_1 / c_1.sqrt() - r * ψ_cos;
+            let g_2 = r2 * α_2 / c_2.sqrt() - r * ψ_cos;
+            let h_12 = (r / 2.0) * ψ_sin;
+            let h_11 = r2 * d2_1 / c_1.powf(3.0 / 2.0) + h_12;
+            let h_22 = r2 * d2_2 / c_2.powf(3.0 / 2.0) + h_12;
+            ([h_11, h_12, h_22], [g_1, g_2])
+        };
+
+        let [α_1, α_2] = newton_sym2x2(g, [0.0, 0.0], self.config.newton_iters);
+        (normalise_angle(φ_1 + α_1), normalise_angle(φ_2 + α_2))
+
+    }
+
+    /// Find angles where the geodesic passing through the side from `cap1` at height `z_1`
+    /// to `cap2` at height `z_2` crosses the cap edges.
+    /// **Panics if `cap2.angle < cap1.angle`.**
+    ///
+    /// Uses `newton_sym2x2`.
+    fn cap_side_cap_crossing(
+        &self,
+        cap1 : &CapPoint, z_1 : f64,
+        cap2 : &CapPoint, z_2 : f64,
+        init_by_cap2 : bool
+    ) -> (Angle, Angle) {
+        let r = self.radius;
+        let &CapPoint { r : r_1, angle : φ_1 } = cap1;
+        let &CapPoint { r : r_2, angle : φ_2 } = cap2;
+        assert!(φ_2 >= φ_1);
+        assert_ne!(r_1, r);
+        assert_ne!(r_2, r);
+        assert_ne!(z_2, z_1);
+        if r_1 == 0.0 && r_2 == 0.0 {
+            // Singular case: both points are in the middle of the caps.
+            return (φ_1, φ_1)
+        }
+        let r2 = r * r;
+        let d = (z_2 - z_1).abs();
+        let d2 = d * d;
+        let r2_1 = r_1 * r_1;
+        let r2_2 = r_2 * r_2;
+        let rr_1 = r * r_1;
+        let rr_2 = r * r_2;
+        let f = |α_1 : f64, α_2 : f64| {
+            let cos_α1 = α_1.cos();
+            let sin_α1 = α_1.sin();
+            let cos_α2 = α_2.cos();
+            let sin_α2 = α_2.sin();
+            //let cos2_α1 = cos_α1 * cos_α1;
+            let sin2_α1 = sin_α1 * sin_α1;
+            //let cos2_α2 = cos_α2 * cos_α2;
+            let sin2_α2 = sin_α2 * sin_α2;
+            let ψ = φ_2 - φ_1 - α_1 - α_2;
+            let ψ2 = ψ * ψ;
+            let ψ2r2 = ψ2 * r2;
+            //let r4 = r2 * r2;
+            let b = d2 + ψ2r2;
+            let c = r2 * ψ / b.sqrt();
+            let e_1 = r2 + r2_1 - 2.0 * rr_1 * cos_α1;
+            let e_2 = r2 + r2_2 - 2.0 * rr_2 * cos_α2;
+            let g_1 = rr_1 * sin_α1 / e_1.sqrt() - c;
+            let g_2 = rr_2 * sin_α2 / e_2.sqrt() - c;
+            let h_12 = r2 * (1.0  - ψ2r2 / b) / b.sqrt();
+            // let h_11 = rr_1 * ( (r2 + r2_1) * cos_α1 - rr_1 * ( sin2_α1 + 2.0 * cos2_α1) )
+            //                 / e_1.powf(3.0/2.0) + h_12;
+            // let h_22 = rr_2 * ( (r2 + r2_2) * cos_α2 - rr_2 * ( sin2_α2 + 2.0 * cos2_α2) )
+            //                 / e_2.powf(3.0/2.0) + h_12;
+            // let h_11 = rr_1 * cos_α1 / e_1.sqrt() - rr_1*rr_1 * sin2_α1 / e_1.powf(3.0/2.0) + h_12;
+            // let h_22 = rr_2 * cos_α2 / e_2.sqrt() - rr_2*rr_2 * sin2_α2 / e_2.powf(3.0/2.0) + h_12;
+            let h_11 = rr_1 * (cos_α1 - rr_1 * sin2_α1 / e_1) / e_1.sqrt() + h_12;
+            let h_22 = rr_2 * (cos_α2 - rr_2 * sin2_α2 / e_2) / e_2.sqrt() + h_12;
+            ([h_11, h_12, h_22], [g_1, g_2])
+        };
+
+        let α_init = if init_by_cap2 {
+            [φ_2 - φ_1, 0.0]
+        } else {
+            [0.0, φ_2 - φ_1]
+        };
+        let [α_1, α_2] = newton_sym2x2(f, α_init, self.config.newton_iters);
+        (normalise_angle(φ_1 + α_1), normalise_angle(φ_2 - α_2))
+    }
+
+    fn cap_side_log(
+        &self,
+        cap : &CapPoint, (z, top) : (f64, bool),
+        side : &SidePoint
+    ) -> Tangent {
+        let r = self.radius;
+        if indistinguishable(side.z, z) {
+            // Degenerate case
+            let capedge = CapPoint{ angle : side.angle, r };
+            cap.log(&capedge)
+        } else if indistinguishable(r, cap.r)
+                  && anglediff(side.angle, cap.angle).abs() < f64::PI/2.0 {
+            // Degenerate case
+            let sideedge = SidePoint{ angle : cap.angle, z};
+            cap.tangent_from_side(top, sideedge.log(r, side))
+        } else {
+            let φ = self.side_cap_crossing(side, cap, z);
+            let capedge = CapPoint{ angle : φ, r };
+            let sideedge = SidePoint{ angle : φ, z };
+            let t1 = cap.log(&capedge);
+            let t2 = sideedge.log(r, side);
+            // Either option should give the same result, but the first one avoids division.
+            t1 + capedge.tangent_from_side(top, t2)
+            // let n = t1.norm(L2);
+            // (t1/n)*(n + t2.norm(L2))
+        }
+    }
+
+    fn side_cap_log(
+        &self,
+        side : &SidePoint,
+        cap : &CapPoint, (z, top) : (f64, bool),
+    ) -> Tangent {
+        let r = self.radius;
+        if indistinguishable(side.z, z) {
+            // Degenerate case
+            let capedge = CapPoint{ angle : side.angle, r };
+            capedge.tangent_to_side(top, capedge.log(cap))
+        } else if indistinguishable(r, cap.r)
+                  && anglediff(side.angle, cap.angle).abs() < f64::PI/2.0 {
+            // Degenerate case
+            side.log(r, &SidePoint{ z, angle : cap.angle })
+        } else {
+            let φ = self.side_cap_crossing(side, cap, z);
+            let capedge = CapPoint{ angle : φ, r };
+            let sideedge = SidePoint{ angle : φ, z };
+            let t1 = side.log(r, &sideedge);
+            let t2 = capedge.log(cap);
+            // Either option should give the same result, but the first one avoids division.
+            t1 + capedge.tangent_to_side(top, t2)
+            // let n = t1.norm(L2);
+            // (t1/n)*(n + t2.norm(L2))
+        }
+    }
+
+    fn side_cap_side_log(
+        &self,
+        side1 : &SidePoint,
+        (z, top) : (f64, bool),
+        side2 : &SidePoint
+    ) -> Tangent {
+        let r = self.radius;
+        if indistinguishable(side1.z, z) {
+            // Degenerate case
+            let capedge1 = CapPoint{ angle : side1.angle, r };
+            capedge1.tangent_to_side(top, self.cap_side_log(&capedge1, (z, top), side2))
+        } else if indistinguishable(side2.z, z) {
+            // Degenerate case
+            let capedge2 = CapPoint{ angle : side2.angle, r };
+            self.side_cap_log(side1, &capedge2, (z, top))
+        } else {
+            let (φ1, φ2) = if side2.angle >= side1.angle {
+                self.side_cap_side_crossing(side1, z, side2)
+            } else {
+                swap(self.side_cap_side_crossing(side2, z, side1))
+            };
+            let capedge1 = CapPoint{ angle : φ1, r };
+            let sideedge1 = SidePoint{ angle : φ1, z };
+            let capedge2 = CapPoint{ angle : φ2, r };
+            let sideedge2 = SidePoint{ angle : φ2, z };
+            let t1 = side1.log(r, &sideedge1);
+            let t2 = capedge1.log(&capedge2);
+            let t3 = sideedge2.log(r, &side2);
+            // Any option should give the same result, but the first one avoids division.
+            // t1 + capedge1.tangent_to_side(top, t2 + capedge2.tangent_from_side(top, t3))
+            //
+            // let n = t2.norm(L2);
+            // t1 + capedge1.tangent_to_side(top, (t2/n)*(n + t3.norm(L2)))
+            //
+            let n = t1.norm(L2);
+            (t1/n)*(n + t2.norm(L2) + t3.norm(L2))
+            //
+            // let n = t1.norm(L2);
+            // let t23 = t2 + capedge2.tangent_from_side(top, t3);
+            // (t1/n)*(n + t23.norm(L2))
+        }
+    }
+
+    fn cap_side_cap_log(
+        &self,
+        cap1 : &CapPoint, (z1, top1) : (f64, bool),
+        cap2 : &CapPoint, (z2, top2) : (f64, bool),
+        init_by_cap2 : bool,
+    ) -> Tangent {
+        let r = self.radius;
+        if indistinguishable(cap1.r, r) {
+            // Degenerate case
+            let sideedge1 = SidePoint{ angle : cap1.angle, z : z1 };
+            cap1.tangent_from_side(top1, self.side_cap_log(&sideedge1, cap2, (z2, top2)))
+        } else if indistinguishable(cap2.r, r) {
+            // Degenerate case
+            let sideedge2 = SidePoint{ angle : cap2.angle, z : z2 };
+            self.cap_side_log(cap1, (z1, top1), &sideedge2)
+        } else {
+            let (φ1, φ2) = if cap2.angle >= cap1.angle {
+                self.cap_side_cap_crossing(cap1, z1, cap2, z2, init_by_cap2)
+            } else {
+                swap(self.cap_side_cap_crossing(cap2, z2, cap1, z1, !init_by_cap2))
+            };
+            let sideedge1 = SidePoint{ angle : φ1, z : z1 };
+            let capedge1 = CapPoint{ angle : φ1, r };
+            let sideedge2 = SidePoint{ angle : φ2, z : z2};
+            let capedge2 = CapPoint{ angle : φ2, r };
+            let t1 = cap1.log(&capedge1);
+            let t2 = sideedge1.log(r, &sideedge2);
+            let t3 = capedge2.log(cap2);
+            // Either option should give the same result, but the first one avoids division.
+            t1 + capedge1.tangent_from_side(top1, t2 + capedge2.tangent_to_side(top2, t3))
+            //let n = t1.norm(L2);
+            //(t1/n)*(n + t2.norm(L2) + t3.norm(L2))
+        }
+    }
+
+    /// Calculates both the logarithmic map and distance to another point
+    fn log_dist(&self, source : &Point, destination : &Point) -> (Tangent, f64) {
+        use Point::*;
+        match (source, destination) {
+            (Top(cap1), Top(cap2)) => {
+                best_tangent([cap1.log(cap2)])
+            },
+            (Bottom(cap1), Bottom(cap2)) => {
+                best_tangent([cap1.log(cap2)])
+            },
+            (Bottom(cap), Side(side)) => {
+                best_tangent([self.cap_side_log(cap, (self.bottom_z(), false), side)])
+            },
+            (Top(cap), Side(side)) => {
+                best_tangent([self.cap_side_log(cap, (self.top_z(), true), side)])
+            },
+            (Side(side), Bottom(cap)) => {
+                best_tangent([self.side_cap_log(side, cap, (self.bottom_z(), false))])
+            },
+            (Side(side), Top(cap)) => {
+                best_tangent([self.side_cap_log(side, cap, (self.top_z(), true))])
+            },
+            (Side(side1), Side(side2)) => {
+                best_tangent([
+                    side1.log(self.radius, side2),
+                    self.side_cap_side_log(side1, (self.top_z(), true), side2),
+                    self.side_cap_side_log(side1, (self.bottom_z(), false), side2),
+                ])
+            },
+            (Top(cap1), Bottom(cap2)) => {
+                best_tangent([
+                    // We try a few possible initialisations for Newton
+                    self.cap_side_cap_log(
+                        cap1, (self.top_z(), true),
+                        cap2, (self.bottom_z(), false),
+                        false
+                    ),
+                    self.cap_side_cap_log(
+                        cap1, (self.top_z(), true),
+                        cap2, (self.bottom_z(), false),
+                        true
+                    ),
+                ])
+            },
+            (Bottom(cap1), Top(cap2)) => {
+                best_tangent([
+                    // We try a few possible initialisations for Newton
+                    self.cap_side_cap_log(
+                        cap1, (self.bottom_z(), false),
+                        cap2, (self.top_z(), true),
+                        false
+                    ),
+                    self.cap_side_cap_log(
+                        cap1, (self.bottom_z(), false),
+                        cap2, (self.top_z(), true),
+                        true
+                    ),
+                ])
+            },
+        }
+    }
+
+    #[allow(unreachable_code)]
+    #[allow(unused_variables)]
+    fn partial_exp(&self, point : Point, t  : Tangent) -> (Point, Option<Tangent>) {
+        match point {
+            Point::Top(cap) => {
+                let (cap_new, t_new_basis) = cap.partial_exp(self.radius, &t);
+                match t_new_basis {
+                    None => (Point::Top(cap_new), None),
+                    Some(t_new) => {
+                        let side_new = SidePoint{ angle : cap_new.angle, z : self.top_z() };
+                        (Point::Side(side_new), Some(cap_new.tangent_to_side(true, t_new)))
+                    }
+                }
+            },
+            Point::Bottom(cap) => {
+                let (cap_new, t_new_basis) = cap.partial_exp(self.radius, &t);
+                match t_new_basis {
+                    None => (Point::Bottom(cap_new), None),
+                    Some(t_new) => {
+                        let side_new = SidePoint{ angle : cap_new.angle, z : self.bottom_z() };
+                        (Point::Side(side_new), Some(cap_new.tangent_to_side(false, t_new)))
+                    }
+                }
+            },
+            Point::Side(side) => {
+                let lims = (self.bottom_z(), self.top_z());
+                let (side_new, t_new_basis) = side.partial_exp(self.radius, lims, &t);
+                match t_new_basis {
+                    None => (Point::Side(side_new), None),
+                    Some(t_new) => {
+                        if side_new.z >= self.top_z() - f64::EPSILON {
+                            let cap_new = CapPoint { angle : side_new.angle, r : self.radius };
+                            (Point::Top(cap_new), Some(cap_new.tangent_from_side(true, t_new)))
+                        } else {
+                            let cap_new = CapPoint { angle : side_new.angle, r : self.radius };
+                            (Point::Bottom(cap_new), Some(cap_new.tangent_from_side(false, t_new)))
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    fn exp(&self, point : &Point, tangent : &Tangent) -> Point {
+        let mut p = *point;
+        let mut t = *tangent;
+        loop {
+            (p, t) = match self.partial_exp(p, t) {
+                (p, None) => break p,
+                (p, Some(t)) => (p, t),
+            };
+        }
+    }
+
+    /// Check that `point` has valid coordinates, and normalise angles
+    pub fn normalise(&self, point : Point) -> Option<Point> {
+        match point {
+            Point::Side(side) => {
+                let a = self.bottom_z();
+                let b = self.top_z();
+                (a <= side.z && side.z <= b).then(|| {
+                    Point::Side(SidePoint{ angle : normalise_angle(side.angle), .. side })
+                })
+            },
+            Point::Bottom(cap) => {
+                (cap.r <= self.radius).then(|| {
+                    Point::Bottom(CapPoint{ angle : normalise_angle(cap.angle), .. cap })
+                })
+            },
+           Point::Top(cap) => {
+                (cap.r <= self.radius).then(|| {
+                    Point::Top(CapPoint{ angle : normalise_angle(cap.angle), .. cap })
+                })
+            },
+        }
+    }
+
+    /// Convert `p` into a a point associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn point_on(&self, point : Point) -> OnCylinder<'_> {
+        match self.normalise(point) {
+            None => panic!("{point:?} not on cylinder {self:?}"),
+            Some(point) => OnCylinder { cylinder : self, point }
+        }
+    }
+
+    /// Convert `p` into a a point on side associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn point_on_side(&self, side : SidePoint) -> OnCylinder<'_> {
+        self.point_on(Point::Side(side))
+    }
+
+    /// Convert `p` into a a point on top associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn point_on_top(&self, cap : CapPoint) -> OnCylinder<'_> {
+        self.point_on(Point::Top(cap))
+    }
+
+    /// Convert `p` into a a point on bottom associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn point_on_bottom(&self, cap : CapPoint) -> OnCylinder<'_> {
+        self.point_on(Point::Bottom(cap))
+    }
+
+    /// Convert `p` into a a point on side associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn on_side(&self, angle : Angle, z : f64) -> OnCylinder<'_> {
+        self.point_on_side(SidePoint{ angle, z })
+    }
+
+    /// Convert `p` into a a point on top associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn on_top(&self, angle : Angle, r : f64) -> OnCylinder<'_> {
+        self.point_on_top(CapPoint{ angle, r })
+    }
+
+    /// Convert `p` into a a point on bottom associated with the cylinder.
+    ///
+    /// May panic if the coordinates are invalid.
+    pub fn on_bottom(&self, angle : Angle, r : f64) -> OnCylinder<'_> {
+        self.point_on_bottom(CapPoint{ angle, r })
+    }
+}
+
+impl<'a> OnCylinder<'a> {
+    /// Return the cylindrical coordinates of this point
+    pub fn cyl_coords(&self) -> CylCoords {
+        self.cylinder.cyl_coords(&self.point)
+    }
+}
+
+impl<'a> EmbeddedManifoldPoint for OnCylinder<'a> {
+    type EmbeddedCoords = Loc<f64, 3>;
+
+    /// Get embedded 3D coordinates
+    fn embedded_coords(&self) -> Loc<f64, 3> {
+        self.cyl_coords().to_cartesian()
+    }
+}
+
+impl<'a> ManifoldPoint for OnCylinder<'a> {
+    type Tangent = Tangent;
+
+    fn exp(self, tangent : &Self::Tangent) -> Self {
+        let cylinder = self.cylinder;
+        let point = cylinder.exp(&self.point, tangent);
+        OnCylinder { cylinder, point }
+    }
+
+    fn log(&self, other : &Self) -> Self::Tangent {
+        assert!(self.cylinder == other.cylinder);
+        self.cylinder.log_dist(&self.point, &other.point).0
+    }
+
+    fn dist_to(&self, other : &Self) -> f64 {
+        assert!(self.cylinder == other.cylinder);
+        self.cylinder.log_dist(&self.point, &other.point).1
+    }
+
+    fn tangent_origin(&self) -> Self::Tangent {
+        Loc([0.0, 0.0])
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    static CYL : Cylinder = Cylinder {
+        radius : 1.0,
+        height : 1.0,
+        config : CylinderConfig { newton_iters : 20 },
+    };
+
+    fn check_distance(distance : f64, expected : f64) {
+        let tol = 1e-10;
+        assert!(
+            (distance-expected).abs() < tol,
+            "distance = {distance}, expected = {expected}"
+        );
+    }
+
+    // fn check_distance_less(distance : f64, expected : f64) {
+    //     let tol = 1e-10;
+    //     assert!(
+    //         distance < expected + tol,
+    //         "distance = {distance}, expected = {expected}"
+    //     );
+    // }
+
+    #[test]
+    fn intra_cap_log_dist() {
+        let π = f64::PI;
+        let p1 = CYL.on_top(0.0, 0.5);
+        let p2 = CYL.on_top(π, 0.5);
+        let p3 = CYL.on_top(π/2.0, 0.5);
+
+        check_distance(p1.dist_to(&p2), 1.0);
+        check_distance(p2.dist_to(&p3), 0.5_f64.sqrt());
+        check_distance(p3.dist_to(&p1), 0.5_f64.sqrt());
+    }
+
+    #[test]
+    fn intra_side_log_dist() {
+        let π = f64::PI;
+        let p1 = CYL.on_side(0.0, 0.0);
+        let p2 = CYL.on_side(0.0, 0.4);
+        let p3 = CYL.on_side(π/2.0, 0.0);
+
+        check_distance(p1.dist_to(&p2), 0.4);
+        check_distance(p1.dist_to(&p3), π/2.0*CYL.radius);
+    }
+
+    #[test]
+    fn intra_side_over_cap_log_dist() {
+        let π = f64::PI;
+        let off = 0.05;
+        let z = CYL.top_z() - off;
+        let p1 = CYL.on_side(0.0, z);
+        let p2 = CYL.on_side(π, z);
+
+        check_distance(p1.dist_to(&p2), 2.0 * (CYL.radius + off));
+    }
+
+    #[test]
+    fn top_bottom_log_dist() {
+        let π = f64::PI;
+        let p1 = CYL.on_top(0.0, 0.0);
+        let p2 = CYL.on_bottom(0.0, 0.0);
+
+        check_distance(p1.dist_to(&p2), 2.0 * CYL.radius + CYL.height);
+
+        let p1 = CYL.on_top(0.0, CYL.radius / 2.0);
+        let p2 = CYL.on_bottom(0.0, CYL.radius / 2.0);
+        let p3 = CYL.on_bottom(π, CYL.radius / 2.0);
+
+        check_distance(p1.dist_to(&p2), CYL.radius + CYL.height);
+        check_distance(p1.dist_to(&p3), 2.0 * CYL.radius + CYL.height);
+    }
+
+    #[test]
+    fn top_side_log_dist() {
+        let p1 = CYL.on_top(0.0, 0.0);
+        let p2 = CYL.on_side(0.0, 0.0);
+
+        check_distance(p1.dist_to(&p2), CYL.radius + CYL.height / 2.0);
+    }
+}

--- a/src/main.rs	Thu Nov 07 16:52:05 2024 -0500
+++ b/src/main.rs	Wed Dec 04 23:19:46 2024 -0500
@@ -11,9 +11,11 @@
 mod manifold;
 mod fb;
 mod cube;
+mod cylinder;
 mod dist;
 mod zero;
 mod scaled;
+mod newton;
 
 use serde::{Serialize, Deserialize};
 use alg_tools::logger::Logger;
@@ -29,55 +31,61 @@
 
 use dist::{DistTo, DistToSquaredDiv2};
 use fb::{forward_backward, IterInfo, Desc, Prox};
-use manifold::EmbeddedManifoldPoint;
-use cube::*;
-use Face::*;
+use manifold::{EmbeddedManifoldPoint, ManifoldPoint, FacedManifoldPoint};
+use cube::OnCube;
+use cylinder::{CylCoords, Cylinder, CylinderConfig, OnCylinder, normalise_angle};
 #[allow(unused_imports)]
 use zero::ZeroFn;
 use scaled::Scaled;
 
+/// Location for saving results
+static PREFIX : &str = "res";
+
 /// Program entry point
 fn main() {
-    simple_cube_test().unwrap()
+    simple_cube_test(format!("{PREFIX}/cube").as_str()).unwrap();
+    simple_cylinder_test(format!("{PREFIX}/cylinder").as_str()).unwrap();
 }
 
 /// Helper structure for saving a point on a cube into a CSV file
 #[derive(Serialize,Deserialize,Debug)]
-struct CSVPoint {
+struct CSVPoint<Face> {
     face : Face,
     x : f64,
     y : f64,
     z : f64
 }
 
-impl From<&OnCube> for CSVPoint {
-    fn from(point : &OnCube) -> Self {
+impl<M> From<&M> for CSVPoint<M::Face>
+where M : EmbeddedManifoldPoint<EmbeddedCoords = Loc<f64, 3>> + FacedManifoldPoint {
+    fn from(point : &M) -> Self {
         let Loc([x,y,z]) = point.embedded_coords();
         let face = point.face();
         CSVPoint { face, x,  y, z }
     }
 }
 
-/// Helper structure for saving the log into a CSV file
+/// Helper structure for saving a point on a cylinder into a CSV file
 #[derive(Serialize,Deserialize,Debug)]
-struct CSVLog {
-    iter : usize,
-    value : f64,
-    // serde is junk
-    //#[serde(flatten)]
-    //point : CSVPoint
+struct CSVCylPoint<Face> {
     face : Face,
-    x : f64,
-    y : f64,
+    angle : f64,
+    r : f64,
     z : f64
 }
 
-
-/// Location for saving results
-static PREFIX : &str = "res";
+impl<'a> From<&'a crate::cylinder::OnCylinder<'a>> for CSVCylPoint<crate::cylinder::Face> {
+    fn from(point : &'a OnCylinder<'a>) -> Self {
+        let CylCoords {r, angle, z} = point.cyl_coords();
+        let face = point.face();
+        CSVCylPoint { face, r, angle, z }
+    }
+}
 
 /// A simple test on the cube
-fn simple_cube_test() -> DynError {
+fn simple_cube_test(dir : &str) -> DynError {
+    use crate::cube::Face;
+    use crate::cube::Face::*;
     
     let points = [
         OnCube::new(F1, Loc([0.5, 0.7])),
@@ -90,37 +98,34 @@
 
     let origin = OnCube::new(F4, Loc([0.5, 0.5]));
 
-    write_points(format!("{PREFIX}/data"), points.iter())?;
-    write_points(format!("{PREFIX}/origin"), std::iter::once(&origin))?;
+    std::fs::create_dir_all(dir)?;
+    write_csv(points.iter().map(CSVPoint::from), format!("{dir}/data.csv"))?;
+    write_csv(std::iter::once(&origin).map(CSVPoint::from), format!("{dir}/origin.csv"))?;
 
     let f = Sum::new(points.into_iter().map(DistToSquaredDiv2));
     //let g = ZeroFn::new();
     let g = Scaled::new(0.5, DistTo(origin));
     let τ = 0.1;
     
-    std::fs::create_dir_all(PREFIX)?;
     for face in Face::all() {
-        write_face_csv(format!("{PREFIX}/{face}"), face, 32, |x| f.apply(x) + g.apply(x))?;
+        write_cube_face_csv(format!("{dir}/{face}"), face, 32, |x| f.apply(x) + g.apply(x))?;
     }
-    write_face_imgs(128, |x| f.apply(x) + g.apply(x))?;
-
-    run_and_save("x1", &f, &g, OnCube::new(F3, Loc([0.1, 0.7])), τ)?;
-    run_and_save("x2", &f, &g, OnCube::new(F2, Loc([0.3, 0.1])), τ)?;
-    run_and_save("x3", &f, &g, OnCube::new(F6, Loc([0.6, 0.2])), τ)
-}
+    write_cube_face_imgs(dir, 128, |x| f.apply(x) + g.apply(x))?;
 
-/// Runs [forward_backward] and saves the results.
-pub fn run_and_save<F, G>(
-    name : &str,
-    f : &F,
-    g : &G,
-    x : OnCube,
-    τ : f64,
-) -> DynError
-where F : Desc<OnCube> +  Mapping<OnCube, Codomain = f64>,
-      G : Prox<OnCube> +  Mapping<OnCube, Codomain = f64> {
-    
-    let mut logger = Logger::new();
+    /// Helper structure for saving the log into a CSV file
+    #[derive(Serialize,Deserialize,Debug)]
+    struct CSVLog<Face> {
+        iter : usize,
+        value : f64,
+        // serde is junk
+        //#[serde(flatten)]
+        //point : CSVPoint
+        face : Face,
+        x : f64,
+        y : f64,
+        z : f64
+    }
+
     let logmap = |iter, IterInfo { value, point } : IterInfo<OnCube>| {
         let CSVPoint {x , y, z, face} = CSVPoint::from(&point);
         CSVLog {
@@ -128,6 +133,93 @@
             x, y, z, face
          }
     };
+
+    run_and_save(dir, "x1", &f, &g, OnCube::new(F3, Loc([0.1, 0.7])), τ, logmap)?;
+    run_and_save(dir, "x2", &f, &g, OnCube::new(F2, Loc([0.3, 0.1])), τ, logmap)?;
+    run_and_save(dir, "x3", &f, &g, OnCube::new(F6, Loc([0.6, 0.2])), τ, logmap)
+}
+
+
+/// A simple test on the cube
+fn simple_cylinder_test(dir : &str) -> DynError {
+    let π = f64::PI;
+
+    let cyl = Cylinder {
+        radius : 1.0,
+        height : 1.0,
+        config : CylinderConfig { newton_iters : 100 },
+        //Default::default(),
+    };
+    
+    let points = [
+        cyl.on_side(π/2.0, 0.4),
+        cyl.on_side(π*3.0/2.0, -0.1),
+        cyl.on_side(π, 0.2),
+        cyl.on_side(-π/5.0, -0.2),
+        cyl.on_side(π*2.0/3.0, -0.45),
+        cyl.on_top(π*3.0/4.0, 0.7),
+        cyl.on_bottom(π*5.0/4.0, 0.3),
+    ];
+
+    let origin = cyl.on_side(0.0, 0.0);
+
+    std::fs::create_dir_all(dir)?;
+    write_csv(points.iter().map(CSVCylPoint::from), format!("{dir}/data.csv"))?;
+    write_csv(std::iter::once(&origin).map(CSVCylPoint::from), format!("{dir}/origin.csv"))?;
+
+    let f = Sum::new(points.into_iter().map(DistToSquaredDiv2));
+    let g = Scaled::new(4.0, DistTo(origin));
+    let τ = 0.1;
+    
+    write_cylinder_faces_csv(dir, &cyl, 32, 32, 32, |x| f.apply(x) + g.apply(x))?;
+
+    /// Helper structure for saving the log into a CSV file
+    #[derive(Serialize,Deserialize,Debug)]
+    struct CSVLog<Face> {
+        iter : usize,
+        value : f64,
+        // serde is junk
+        //#[serde(flatten)]
+        //point : CSVCylPoint
+        face : Face,
+        angle : f64,
+        r : f64,
+        z : f64,
+    }
+
+    let logmap = |iter, IterInfo { value, point } : IterInfo<OnCylinder<'_>>| {
+        let CSVCylPoint {r, angle, z, face} = CSVCylPoint::from(&point);
+        CSVLog {
+            iter, value, //point : CSVPoint::from(&point)
+            r, angle, z, face
+         }
+    };
+
+    run_and_save(dir, "x1", &f, &g, cyl.on_top(0.0, 0.0), τ, logmap)?;
+    run_and_save(dir, "x2", &f, &g, cyl.on_side(-π*2.0/5.0, 0.0), τ, logmap)?;
+    run_and_save(dir, "x3", &f, &g,  cyl.on_bottom(-π*0.5, 0.8), τ, logmap)
+
+}
+
+
+/// Runs [forward_backward] and saves the results.
+pub fn run_and_save<F, G, M, I>(
+    dir : &str,
+    name : &str,
+    f : &F,
+    g : &G,
+    x : M,
+    τ : f64,
+    logmap : impl Fn(usize, IterInfo<M>) -> I
+) -> DynError
+where M : ManifoldPoint
+          + EmbeddedManifoldPoint<EmbeddedCoords = Loc<f64, 3>>
+          + FacedManifoldPoint,
+      F : Desc<M> +  Mapping<M, Codomain = f64>,
+      G : Prox<M> +  Mapping<M, Codomain = f64>,
+      I : Serialize {
+    
+    let mut logger = Logger::new();
     let iter = AlgIteratorOptions{
         max_iter : 20,
         verbose_iter : Verbose::Every(1),
@@ -138,14 +230,20 @@
     let x̂ = forward_backward(f, g, x, τ, iter);
     println!("result = {}\n{:?}", x̂.embedded_coords(), &x̂);
 
-    logger.write_csv(format!("{PREFIX}/{name}_log.csv"))?;
+    logger.write_csv(format!("{dir}/{name}_log.csv"))?;
 
     Ok(())
 }
 
 /// Writes the values of `f` on `face` of a [`OnCube`] into a PNG file
 /// with resolution `n × n`.
-fn write_face_imgs(n : usize, mut f : impl FnMut(&OnCube) -> f64) -> DynError {
+fn write_cube_face_imgs(
+    dir : &str,
+    n : usize,
+    mut f : impl FnMut(&OnCube) -> f64
+) -> DynError {
+    use crate::cube::Face;
+
     let grid = LinSpace {
         start : Loc([0.0, 0.0]),
         end : Loc([1.0, 1.0]),
@@ -175,7 +273,7 @@
                 *p = Rgb(rgb.map(|v| (v*(u8::RANGE_MAX as f64)) as u8))
             });
 
-        img.save_with_format(format!("{PREFIX}/{face}.png"), ImageFormat::Png)?;
+        img.save_with_format(format!("{dir}/{face}.png"), ImageFormat::Png)?;
     }
 
     Ok(())
@@ -183,7 +281,12 @@
 
 /// Writes the values of `f` on `face` of a [`OnCube`] into a CSV file
 /// with resolution `n × n`.
-fn write_face_csv(filename : String, face : Face, n : usize, mut f : impl FnMut(&OnCube) -> f64) -> DynError {
+fn write_cube_face_csv(
+    filename : String,
+    face : crate::cube::Face,
+    n : usize,
+    mut f : impl FnMut(&OnCube) -> f64
+) -> DynError {
 
     #[derive(Serialize)]
     struct CSVFace { u : f64, v : f64, value : f64 }
@@ -202,10 +305,70 @@
     Ok(())
 }
 
-/// Writes a list of points on a [`OnCube`] into a CSV file
-fn write_points<'a, I : Iterator<Item=&'a OnCube>>(filename : String, list : I) -> DynError {
+/// Writes the values of `f` on the faces of a `cylinder`.
+fn write_cylinder_faces_csv<'a>(
+    dir : &str,
+    cyl : &'a Cylinder,
+    n_height : usize,
+    n_radius : usize,
+    n_angle : usize,
+    mut f : impl FnMut(&OnCylinder<'a>) -> f64
+) -> DynError {
+    let π = f64::PI;
+
+    // Side front is [a, b] for the TikZ configuration.
+    let a = -π*5.0/6.0;
+    let b = π/6.0;
+
+    #[derive(Serialize)]
+    struct CSVFace { angle : f64, /*cos_angle : f64, sin_angle : f64,*/ v : f64, value : f64 }
+
+    let mkf = |angle, v, value| CSVFace {
+        angle,
+        //cos_angle : angle.cos(),
+        //sin_angle : angle.sin(),
+        v,
+        value
+    };
+
+    let side_half_grid = LinSpace {
+        start : Loc([0.0, cyl.bottom_z()]),
+        end : Loc([π, cyl.top_z()]),
+        count : [n_angle / 2, n_height]
+    };
 
-    write_csv(list.map(CSVPoint::from), format!("{filename}.csv"))?;
-    Ok(())
+    let side_grid = LinSpace {
+        start : Loc([0.0, cyl.bottom_z()]),
+        end : Loc([2.0*π, cyl.top_z()]),
+        count : [n_angle, n_height]
+    };
+    let cap_grid = LinSpace {
+        start : Loc([0.0, 0.0]),
+        end : Loc([2.0*π, cyl.radius]),
+        count : [n_angle, n_radius]
+    };
+
+    let side_front = side_grid.into_iter()
+        .map(|Loc([angle, v])| mkf(angle, v, f(&cyl.on_side(angle, v))));
+    write_csv(side_front, format!("{dir}/side.csv"))?;
+
+    let side_front = side_half_grid.into_iter()
+        .map(|Loc([angle, v])| (normalise_angle(angle + a), v))
+        .map(|(angle, v)| mkf(angle, v, f(&cyl.on_side(angle, v))));
+    write_csv(side_front, format!("{dir}/side_front.csv"))?;
+
+    let side_back = side_half_grid.into_iter()
+        .map(|Loc([angle, v])| (normalise_angle(angle + b), v))
+        .map(|(angle, v)| mkf(angle + π, v, f(&cyl.on_side(angle + π, v))));
+    write_csv(side_back, format!("{dir}/side_back.csv"))?;
+
+    let top = cap_grid.into_iter()
+        .map(|Loc([angle, v])| mkf(angle, v, f(&cyl.on_top(angle, v))));
+    write_csv(top, format!("{dir}/top.csv"))?;
+
+    let bottom = cap_grid.into_iter()
+        .map(|Loc([angle, v])| mkf(angle, v, f(&cyl.on_bottom(angle, v))));
+    write_csv(bottom, format!("{dir}/bottom.csv"))
+    
 }
 

--- a/src/manifold.rs	Thu Nov 07 16:52:05 2024 -0500
+++ b/src/manifold.rs	Wed Dec 04 23:19:46 2024 -0500
@@ -2,12 +2,13 @@
 Abstract traits for manifolds.
 */
 
+use serde::Serialize;
 use alg_tools::euclidean::Euclidean;
 
 /// A point on a manifold
 pub trait ManifoldPoint : Clone + PartialEq {
     // Type of tangent factors
-    type Tangent : Euclidean<f64, Output=Self::Tangent> + std::fmt::Debug;
+    type Tangent : Euclidean<f64, Output=Self::Tangent> + std::fmt::Debug + Serialize;
 
     /// Exponential map
     fn exp(self, tangent : &Self::Tangent) -> Self;
@@ -24,8 +25,16 @@
 
 /// Point on a manifold that possesses displayable embedded coordinates.
 pub trait EmbeddedManifoldPoint : ManifoldPoint + std::fmt::Debug {
-    type EmbeddedCoords : std::fmt::Display;
+    type EmbeddedCoords : std::fmt::Display + Serialize;
 
     /// Convert a point on a manifold into embedded coordinates
     fn embedded_coords(&self) -> Self::EmbeddedCoords;
 }
+
+/// Point on a manifold that possesses faces
+pub trait FacedManifoldPoint : ManifoldPoint + std::fmt::Debug {
+    type Face : std::fmt::Display + Serialize;
+
+    /// Convert a point on a manifold into embedded coordinates
+    fn face(&self) -> Self::Face;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/newton.rs	Wed Dec 04 23:19:46 2024 -0500
@@ -0,0 +1,34 @@
+/*!
+Newton method in 2D.
+*/
+
+use alg_tools::types::*;
+
+#[inline]
+pub fn newton_sym1x1<F : Float>(
+    g : impl Fn(F) -> (F, F),
+    mut x : F,
+    iters : usize
+) -> F {
+    for _i in 0..iters {
+        let (a, b) = g(x);
+        x -= b / a
+    }
+    x
+}
+
+#[inline]
+pub fn newton_sym2x2<F : Float>(
+    g : impl Fn(F, F) -> ([F; 3], [F; 2]),
+    [mut x1, mut x2] : [F; 2],
+    iters : usize
+) -> [F; 2] {
+    for _i in 0..iters {
+        let ([a11, a12, a22], [b1, b2]) = g(x1, x2);
+        let q = a11 * a22 - a12 * a12;
+        x1 -= (a22 * b1 - a12 * b2) / q;
+        x2 -= (a11 * b2 - a12 * b1) / q;
+    }
+    [x1, x2]
+}
+

--- a/visualisation/cube.tex	Thu Nov 07 16:52:05 2024 -0500
+++ b/visualisation/cube.tex	Wed Dec 04 23:19:46 2024 -0500
@@ -1,93 +1,39 @@
 \documentclass[tikz]{standalone}
-\usepackage{pgfplots}
-\usepackage{tikz-3dplot}
-\usepackage[svgnames]{xcolor}
-\usepgfplotslibrary{colorbrewer}
 
-\def\datapath{../res/}
+\input{plotsetup}
+
+\def\datapath{../res/cube}
 
 \begin{document}
 
-\tdplotsetmaincoords{50}{25}
-
-\def\onlyfront#1{
-    \ifnum#1=2\else\ifnum#1=4\else\ifnum#1=6\else
-    \def\pgfmathresult{}
-    \fi\fi\fi}
-
-\def\onlyback#1{
-    \ifnum#1=1\else\ifnum#1=3\else\ifnum#1=5\else
-    \def\pgfmathresult{}
-    \fi\fi\fi}
+\begin{tikzpicture}
+    \pgfplotsset{
+        onlyfront/.code = \pgfplotsset{x filter/.expression={%
+            or(\thisrow{face} == 2, or(\thisrow{face} == 4, \thisrow{face} == 6)) ? x : nan%
+        },},
+        onlyback/.code = \pgfplotsset{x filter/.expression={%
+            or(\thisrow{face} == 1, or(\thisrow{face} == 3, \thisrow{face} == 5)) ? x : nan%
+        },},
+    }
+    
+    \begin{axis}[illustr3d]
+        % Hidden edges
+        \addplot3[edge] coordinates {(0, 0, 0) (0, 1, 0) (0, 1, 1) };
+        \addplot3[edge] coordinates {(0, 1, 0) (1, 1, 0) };
 
-\pgfplotsset{
-    % forget plot = no legend
-    cube/.style = {thick, gray!50!black,forget plot},
-    onlyfront/.code = \pgfplotsset{x filter/.code={\onlyfront{\thisrow{face}}}},
-    onlyback/.code = \pgfplotsset{x filter/.code={\onlyback{\thisrow{face}}}},
-    data/.style = {mark=*, blue, only marks},
-    backdata/.style = {mark=*, blue!30!white, only marks, onlyback, forget plot},
-    iter1/.style = {mark=x, red},
-    iter2/.style = {mark=star, BlueViolet},
-    iter3/.style = {mark=asterisk, purple},
-    backiter1/.style = {mark=x, red!30!white, onlyback, forget plot},
-    backiter2/.style = {mark=star, BlueViolet!30!white, onlyback, forget plot},
-    backiter3/.style = {mark=asterisk, purple!30!white, onlyback, forget plot},
-    origin/.style = {mark=o, green, only marks, thick},
-    surfstyle/.style = {very nearly opaque, forget plot},
-    legend style = {
-        inner sep = 0pt,
-        outer xsep = 5pt,
-        outer ysep = 0pt,
-        legend cell align = left,
-        align = left,
-        draw = none,
-        fill = none,
-        font = \small
-    },
-    illustr3d/.style = {
-        axis equal,
-        width=5cm,
-        height=5cm,
-        scale only axis,
-        enlargelimits=false,
-        colormap access=map,
-        colormap/Paired,
-        colorbar,
-        point meta rel=axis wide,
-        shader = interp,
-        xlabel = {$x$},
-        ylabel = {$y$},
-        zlabel = {$z$},
-        ticks = none,
-        axis line style = {draw=none},
-        %axis x line = none,
-        %axis y line = none,
-        %axis z line = none,
-        legend columns = 3,
-        legend style = {
-            at = {(0.5, 1.1)},
-            anchor = north,
-            column sep = 1ex,
-        },
-        mark size=1.5pt,
-    }
-}
+        % Data on hidden faces
+        \addplot3[backdata] table[x=x,y=y,z=z] {\datapath/data.csv};
 
-\begin{tikzpicture}
-    \begin{axis}[illustr3d]
-        \addplot3[cube] coordinates {(0, 0, 0) (0, 1, 0) (0, 1, 1) };
-        \addplot3[cube] coordinates {(0, 1, 0) (1, 1, 0) };
-        \addplot3[backdata] table[x=x,y=y,z=z] {\datapath/data.csv};
+        % Iterates on hidden faces
         \addplot3[backiter1] table[x=x,y=y,z=z] {\datapath/x1_log.csv};
         \addplot3[backiter2] table[x=x,y=y,z=z] {\datapath/x2_log.csv};
         \addplot3[backiter3] table[x=x,y=y,z=z] {\datapath/x3_log.csv};
 
+        % Surface
         \addplot3[
             surf,
-            mesh/ordering=x varies,
-            mesh/cols=32,
-            mesh/rows=32,
+            mesh/ordering=x varies, % the second input coordinate stays fixed while first varies
+            mesh/cols=32, % number of first input coordinate points untils second changes
             surfstyle,
         ] table [
             x = u,
@@ -97,9 +43,8 @@
          ] {\datapath/F4.csv};
         \addplot3[
             surf,
-            mesh/ordering=y varies,
-            mesh/cols=32,
-            mesh/rows=32,
+            mesh/ordering=x varies, % the second input coordinate stays fixed while first varies
+            mesh/cols=32, % number of first input coordinate points untils second changes
             surfstyle,
         ] table [
             x expr = 1.0,
@@ -109,9 +54,8 @@
          ] {\datapath/F2.csv};
         \addplot3[
             surf,
-            mesh/ordering=x varies,
-            mesh/cols=32,
-            mesh/rows=32,
+            mesh/ordering=x varies, % the second input coordinate stays fixed while first varies
+            mesh/cols=32, % number of first input coordinate points untils second changes
             surfstyle,
         ] table [
             x = u,
@@ -121,16 +65,18 @@
         ] {\datapath/F6.csv};
     
         % Edges
-        \addplot3[cube] coordinates {(0, 0, 0) (1, 0, 0) (1, 0, 1) (0, 0, 1) (0, 0, 0)};
-        \addplot3[cube] coordinates {(0, 0, 1) (0, 1, 1) (1, 1, 1) (1, 0, 1)};
-        \addplot3[cube] coordinates {(1, 0, 0) (1, 1, 0) (1, 1, 1)};
+        \addplot3[edge] coordinates {(0, 0, 0) (1, 0, 0) (1, 0, 1) (0, 0, 1) (0, 0, 0)};
+        \addplot3[edge] coordinates {(0, 0, 1) (0, 1, 1) (1, 1, 1) (1, 0, 1)};
+        \addplot3[edge] coordinates {(1, 0, 0) (1, 1, 0) (1, 1, 1)};
 
+        % Data
         \addplot3[data,onlyfront] table[x=x,y=y,z=z] {\datapath/data.csv};
         \addlegendentry{Data}
 
         \addplot3[origin,onlyfront] table[x=x,y=y,z=z] {\datapath/origin.csv};
         \addlegendentry{Origin}
 
+        % Iterates
         \addplot3[iter1,onlyfront] table[x=x,y=y,z=z] {\datapath/x1_log.csv};
         \addlegendentry{Iterates 1}
         \addplot3[iter2,onlyfront] table[x=x,y=y,z=z] {\datapath/x2_log.csv};
@@ -140,4 +86,4 @@
     \end{axis}
 \end{tikzpicture}
 
-\end{document}
\ No newline at end of file
+\end{document}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/visualisation/cylinder.tex	Wed Dec 04 23:19:46 2024 -0500
@@ -0,0 +1,122 @@
+\documentclass[tikz]{standalone}
+
+\input{plotsetup}
+
+\def\datapath{../res/cylinder}
+
+\begin{document}
+
+\begin{tikzpicture}
+    \pgfplotsset{
+        onlyfront/.code = \pgfplotsset{x filter/.expression={%
+            \thisrow{face} == 0 || (%
+                and(\thisrow{face}==2, or(2*pi+\a <= \thisrow{angle}, \thisrow{angle} <= \b))
+            ) ? x : nan%
+        },},
+        onlyback/.code = \pgfplotsset{x filter/.expression={%
+            \thisrow{face} == 1 || (%
+                and(\thisrow{face}==2, or(2*pi+\a > \thisrow{angle}, \thisrow{angle} > \b))
+            ) ? x : nan%
+        },},
+    }
+
+    \pgfmathsetmacro{\a}{-pi*5/6}
+    \pgfmathsetmacro{\b}{pi/6}
+
+    \begin{axis}[illustr3d]
+        % Hidden edges
+        \addplot3[edge, domain=\b:(2*pi+\a), samples = 90, samples y = 0] ({cos(deg(x))}, {sin(deg(x))}, {-0.5});
+
+        % Data on hidden faces
+        \addplot3[backdata] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/data.csv};
+
+        % Iterates on hidden faces
+        \addplot3[backiter1,onlyback] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/x1_log.csv};
+        \addplot3[backiter2,onlyback] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/x2_log.csv};
+        \addplot3[backiter3,onlyback] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/x3_log.csv};
+
+        % Surface
+        \addplot3[
+            surf,
+            mesh/ordering=x varies, % the second input coordinate stays fixed while first varies
+            mesh/cols=16, % number of first input coordinate points untils second changes
+            surfstyle,
+            %x filter/.expression = {or(2*pi+\a <= \thisrow{angle}, \thisrow{angle} <= \b)},
+        ] table [
+            x expr = cos(deg(\thisrow{angle})),
+            y expr = sin(deg(\thisrow{angle})),
+            z = v,
+            point meta = \thisrow{value},
+         ] {\datapath/side_front.csv};
+        \addplot3[
+            surf,
+            mesh/ordering=x varies, % the second input coordinate stays fixed while first varies
+            mesh/cols=32, % number of first input coordinate points untils second changes
+            surfstyle,
+        ] table [
+            x expr = \thisrow{v} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{v} * sin(deg(\thisrow{angle})),
+            z expr = 0.5,
+            point meta = \thisrow{value},
+         ] {\datapath/top.csv};
+     
+        % Edges
+        \addplot3[edge, domain={0:2*pi}, samples = 60, samples y = 0] ({cos(deg(x))}, {sin(deg(x))}, {0.5});
+        \addplot3[edge, domain={\a:\b}, samples = 60, samples y = 0] ({cos(deg(x))}, {sin(deg(x))}, {-0.5});
+        \addplot3[edge] coordinates { (cos(deg(\a)), sin(deg(\a)), -0.5) (cos(deg(\a)), sin(deg(\a)), 0.5) };
+        \addplot3[edge] coordinates { (cos(deg(\b)), sin(deg(\b)), -0.5) (cos(deg(\b)), sin(deg(\b)), 0.5) };
+
+        % Data
+        \addplot3[data, onlyfront] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/data.csv};
+        \addlegendentry{Data}
+
+        \addplot3[origin, onlyfront] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/origin.csv};
+        \addlegendentry{Origin}
+
+        % Iterates
+        \addplot3[iter1,onlyfront] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/x1_log.csv};
+        \addlegendentry{Iterates 1}
+        \addplot3[iter2,onlyfront] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/x2_log.csv};
+        \addlegendentry{Iterates 2}
+        \addplot3[iter3,onlyfront] table[
+            x expr = \thisrow{r} * cos(deg(\thisrow{angle})),
+            y expr = \thisrow{r} * sin(deg(\thisrow{angle})),
+            z = z,
+        ] {\datapath/x3_log.csv};
+        \addlegendentry{Iterates 3}
+    \end{axis}
+\end{tikzpicture}
+
+\end{document}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/visualisation/plotsetup.tex	Wed Dec 04 23:19:46 2024 -0500
@@ -0,0 +1,56 @@
+\usepackage{pgfplots}
+\usepackage[svgnames]{xcolor}
+\usepgfplotslibrary{colorbrewer}
+
+\pgfplotsset{
+    compat=1.18,
+    % forget plot = no legend
+    edge/.style = {thick, gray!50!black,forget plot},
+    data/.style = {mark=*, blue, only marks},
+    backdata/.style = {mark=*, blue!30!white, only marks, onlyback, forget plot},
+    iter1/.style = {mark=x, red},
+    iter2/.style = {mark=star, BlueViolet},
+    iter3/.style = {mark=asterisk, purple},
+    backiter1/.style = {mark=x, red!30!white, onlyback, forget plot},
+    backiter2/.style = {mark=star, BlueViolet!30!white, onlyback, forget plot},
+    backiter3/.style = {mark=asterisk, purple!30!white, onlyback, forget plot},
+    origin/.style = {mark=o, green, only marks, thick},
+    surfstyle/.style = {very nearly opaque, forget plot},
+    legend style = {
+        inner sep = 0pt,
+        outer xsep = 5pt,
+        outer ysep = 0pt,
+        legend cell align = left,
+        align = left,
+        draw = none,
+        fill = none,
+        font = \small
+    },
+    illustr3d/.style = {
+        axis equal,
+        width=5cm,
+        height=5cm,
+        scale only axis,
+        enlargelimits=false,
+        colormap access=map,
+        colormap/Blues,
+        colorbar,
+        point meta rel=axis wide,
+        shader = interp,
+        xlabel = {$x$},
+        ylabel = {$y$},
+        zlabel = {$z$},
+        ticks = none,
+        axis line style = {draw=none},
+        %axis x line = none,
+        %axis y line = none,
+        %axis z line = none,
+        legend columns = 3,
+        legend style = {
+            at = {(0.5, 1.1)},
+            anchor = north,
+            column sep = 1ex,
+        },
+        mark size=1.5pt,
+    }
+}
\ No newline at end of file