Mercurial > repos > non-riemannian-opt / changeset

--- a/Cargo.toml	Fri Dec 06 14:27:14 2024 -0500
+++ b/Cargo.toml	Fri Dec 06 14:57:11 2024 -0500
@@ -24,4 +24,5 @@
 image = "~0.24.3"
 serde_repr = "0.1"
 
-
+[build-dependencies]
+regex = "~1.7.0"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/build.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -0,0 +1,38 @@
+use std::env;
+use regex::{Regex, Captures};
+
+fn proc<A : AsRef<str>>(re : &str, str : A) -> String {
+    let need_to_escape = Regex::new(r"([_*\\])").unwrap();
+    Regex::new(re).unwrap().replace_all(str.as_ref(), |caps : &Captures| {
+        format!("{}{}{}",
+            caps.get(1).unwrap().as_str(),
+            need_to_escape.replace_all(caps.get(2).unwrap().as_str(), "\\$1"),
+            caps.get(3).unwrap().as_str()
+        )
+    }).to_string()
+}
+
+fn main() {
+    let out_dir = env::var("OUT_DIR").unwrap();
+    
+    // Since rust is stuck in 80's 7-bit gringo ASCII world, so that rustdoc does not support
+    // markdown KaTeX mathematics, we have to process the README to include horrible horrible
+    // horrible escapes for the math, and then use an vomit-inducingly ugly javasccript
+    // workaround to process the math on the fly.
+
+    println!("cargo:rerun-if-changed=README.md");
+
+    let readme = std::fs::read_to_string("README.md")
+        .expect("Error reading README");
+
+    // Escape _, *, and \ in equations.
+    let readme_uglified = proc(r"(?m)([^$]\$)([^$]+)(\$[^$])",
+                               proc(r"([^$]\$\$)([^$]+)(\$\$[^$])", readme));
+    // Remove the instructions for building the documentation
+    let readme_cut = Regex::new("## Internals(.*|\n)*")
+        .unwrap()
+        .replace_all(&readme_uglified, "");
+
+    std::fs::write(out_dir + "/README_uglified.md", readme_cut.as_bytes())
+        .expect("Error saving uglified README");
+}

--- a/src/cube.rs	Fri Dec 06 14:27:14 2024 -0500
+++ b/src/cube.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -14,26 +14,40 @@
 pub enum Face {F1 = 1, F2 = 2, F3 = 3, F4 = 4, F5 = 5, F6 = 6}
 use Face::*;
 
+/// General point in 2D
 pub type Point = Loc<f64, 2>;
 
+/// Types for faces adjacent to a given face.
 pub type AdjacentFaces = [Face; 4];
 
+/// Types of paths on a cube
 #[derive(Clone, Debug, Serialize)]
 pub enum Path {
+    /// Direct path from an unindicated source face to a `destination` face.
     Direct { destination : Face },
+    /// Indirect path from an unindicated source face to a `destination` face,
+    /// via an `intermediate` face.
     Indirect { destination : Face, intermediate : Face },
 }
 
 /// An iterator over paths on a cube, from a source face to a destination face.
 #[derive(Clone, Debug)]
 pub enum PathIter {
+    /// Direct path to a destination.
     Same {
+        /// Deistination face
         destination : Face,
+        /// Indicator whether the only possible [`Path::Direct`] has already been returned.
         exhausted : bool
     },
+    /// Path via several possible intermedite faces.
+    /// This is used to generate several [`Path::Indirect`].
     Indirect {
+        /// Destination face
         destination : Face,
+        /// Possible intermediate faces
         intermediate : AdjacentFaces,
+        /// Intermediate face index counter.
         current : usize
     }
 }
@@ -197,7 +211,7 @@
     }
 
     /// Indicates whether an unfolded point `p` is on this face, i.e.,
-    /// has coordinates in [0,1]².
+    /// has coordinates in $\[0,1\]^2$.
     pub fn is_in_face(&self, p: &Point) -> bool {
         p.iter().all(|t| 0.0 <= *t && *t <= 1.0)
     }
@@ -264,6 +278,7 @@
     }
 }
 
+/// Point on a the surface of the unit cube $\[0,1\]^3$.
 #[derive(Clone, Debug, PartialEq, Serialize)]
 pub struct OnCube {
     face : Face,
@@ -348,6 +363,7 @@
 mod tests {
     use super::*;
 
+    /// Tests that the distancse between the centers of all the faces are correctly calculated.
     #[test]
     fn center_distance() {
         let center = Loc([0.5, 0.5]);
@@ -367,6 +383,8 @@
         }
     }
 
+    /// Tests that the distances between points on the boundaries of distinct faces are
+    /// correctly calculated.
     #[test]
     fn boundary_distance() {
         let left = Loc([0.0, 0.5]);
@@ -399,6 +417,7 @@
     }
 
 
+    /// Tests that the conversions between the coordinate systems of each face is working correctly.
     #[test]
     fn convert_adjacent() {
         let point = Loc([0.4, 0.6]);
@@ -440,6 +459,7 @@
     //     }
     // }
 
+    /// Tests that the logarithmic map is working correctly between adjacent faces.
     #[test]
     fn log_adjacent() {
         let p1 = OnCube{ face : F1, point : Loc([0.5, 0.5])};
@@ -448,6 +468,7 @@
         assert_eq!(p1.log(&p2).norm(L2), 1.0);
     }
 
+    /// Tests that the logarithmic map is working correctly between opposing faces.
     #[test]
     fn log_opposing_equal() {
         let p1 = OnCube{ face : F1, point : Loc([0.5, 0.5])};
@@ -456,6 +477,8 @@
         assert_eq!(p1.log(&p2).norm(L2), 2.0);
     }
 
+    /// Tests that the logarithmic map is working correctly between opposing faces when there
+    /// is a unique shortest geodesic.
     #[test]
     fn log_opposing_unique_shortest() {
         let p1 = OnCube{ face : F1, point : Loc([0.3, 0.25])};

--- a/src/cylinder.rs	Fri Dec 06 14:27:14 2024 -0500
+++ b/src/cylinder.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -23,12 +23,14 @@
 }
 
 impl CylCoords {
+    /// Convert to cartesian coordinates.
     #[inline]
     pub fn to_cartesian(&self) -> Loc<f64, 3> {
         let &CylCoords{r, angle, z} = self;
         [r * angle.cos(), r * angle.sin(), z].into()
     }
 
+    /// Form cylindrical coordinates from cartesian coordinates.
     #[inline]
     #[allow(dead_code)]
     pub fn from_cartesian(coords : Loc<f64, 3>) -> Self {
@@ -131,6 +133,7 @@
     pub angle : Angle
 }
 
+/// Calculates the difference between two angles, normalised to [–π, π].
 #[inline]
 fn anglediff(mut φ1 : f64, mut φ2 : f64) -> f64 {
     let π = f64::PI;
@@ -152,6 +155,7 @@
     }
 }
 
+/// Normalises an angle to [0, 2π].
 #[inline]
 pub fn normalise_angle(φ : f64) -> f64 {
     let π = f64::PI;
@@ -209,7 +213,7 @@
 
 }
 
-/// Point on a [`Cylinder`]
+/// Point on some [`Cylinder`]
 #[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
 pub enum Point {
     Top(CapPoint),
@@ -226,16 +230,20 @@
     Side,
 }
 
+/// Point on a specific [`Cylinder`]
 #[derive(Clone, Debug, PartialEq)]
 pub struct OnCylinder<'a> {
+    /// Face and coordinates of the point
+    point : Point,
+    /// The specific cylinder the `point` lies on.
     cylinder : &'a Cylinder,
-    point : Point,
 }
 
 
 /// Cylinder configuration
 #[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
 pub struct CylinderConfig {
+    /// Number of Newton iterates two take to solve geodesics.
     pub newton_iters : usize,
 }
 
@@ -270,6 +278,7 @@
 }
 
 impl Point {
+    /// Returns the face of the point
     fn face(&self) -> Face {
         match *self {
             Point::Top(..) => Face::Top,
@@ -290,6 +299,8 @@
 // Tangent vector
 type Tangent = Loc<f64, 2>;
 
+/// Goes through list of potential tangents (corresponding to several geodesics taking
+/// different paths), and retuns the shortest tangent vector and the corresponding length.
 #[inline]
 fn best_tangent<I>(tangents : I) -> (Tangent, f64)
 where I : IntoIterator<Item = Tangent> {
@@ -305,6 +316,7 @@
     (b, a)
 }
 
+/// Indicates whether the `a` and `b` are a distance of less than [`f64::EPSILON`].
 #[inline]
 fn indistinguishable(a : f64, b : f64) -> bool {
     a > b - f64::EPSILON && a < b + f64::EPSILON
@@ -320,11 +332,13 @@
         }
     }
 
+    /// Returns the z coordinate of the top (cap) of the cylinder.
     #[inline]
     pub fn top_z(&self) -> f64 {
         self.height / 2.0
     }
 
+    /// Returns the z coordinate of the bottom (cap) of the cylinder.
     #[inline]
     pub fn bottom_z(&self) -> f64 {
         -self.height / 2.0
@@ -332,7 +346,7 @@
 
     /// Find angle where a geodesic from `side` to `(cap, z)` crosses the cap edge.
     ///
-    /// Uses `newton_sym1x1`.
+    /// Uses [`newton_sym1x1`].
     fn side_cap_crossing(
         &self,
         side : &SidePoint,
@@ -372,7 +386,7 @@
     /// 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`.
+    /// Uses [`newton_sym2x2`].
     fn side_cap_side_crossing(
         &self,
         side1 : &SidePoint,
@@ -415,7 +429,7 @@
     /// to `cap2` at height `z_2` crosses the cap edges.
     /// **Panics if `cap2.angle < cap1.angle`.**
     ///
-    /// Uses `newton_sym2x2`.
+    /// Uses [`newton_sym2x2`].
     fn cap_side_cap_crossing(
         &self,
         cap1 : &CapPoint, z_1 : f64,
@@ -480,6 +494,9 @@
         (normalise_angle(φ_1 + α_1), normalise_angle(φ_2 - α_2))
     }
 
+    /// Calculates the log between points on a `cap` and a `side` of the cylinder, in this direction.
+    /// The `z` coordinate of the cap and an indication whether it is a `top` or bottom cap must
+    /// also be given.
     fn cap_side_log(
         &self,
         cap : &CapPoint, (z, top) : (f64, bool),
@@ -514,6 +531,9 @@
         }
     }
 
+    /// Calculates the log between points on a `side` and a `cap` of the cylinder, in this direction.
+    /// The `z` coordinate of the cap and an indication whether it is a `top` or bottom cap must
+    /// also be given.
     fn side_cap_log(
         &self,
         side : &SidePoint,
@@ -547,6 +567,9 @@
         }
     }
 
+    /// Calculates the log between points on two sides of the cylinder, assuming the corresonding
+    /// geodesic crosses a cap at height `z`. The `top` parameter indicates whether this is a
+    /// top cap.
     fn side_cap_side_log(
         &self,
         side1 : &SidePoint,
@@ -596,6 +619,9 @@
         }
     }
 
+    /// Calculates the log between points on two caps of the cylinder, assuming the corresonding
+    /// geodesic crosses the side. The heights of the caps and indications whether they are
+    /// top caps, must also be given.
     fn cap_side_cap_log(
         &self,
         cap1 : &CapPoint, (z1, top1) : (f64, bool),
@@ -644,7 +670,7 @@
         }
     }
 
-    /// Calculates both the logarithmic map and distance to another point
+    /// Calculates both the logarithmic map and distance between two points.
     fn log_dist(&self, source : &Point, destination : &Point) -> (Tangent, f64) {
         use Point::*;
         match (source, destination) {
@@ -706,6 +732,9 @@
         }
     }
 
+    /// Calculates the exponential map of `point` in the direction `t` until the next edge.
+    /// If `t` was fully exhausted before reaching an edge, the second return value is [`None`],
+    /// otherwise it is a [`Some`] of the remaining tangent, translated at the returned point.
     #[allow(unreachable_code)]
     #[allow(unused_variables)]
     fn partial_exp(&self, point : Point, t  : Tangent) -> (Point, Option<Tangent>) {
@@ -751,6 +780,7 @@
         }
     }
 
+    /// Calculates the exponential map from `point` in the direction of `tangent`.
     fn exp(&self, point : &Point, tangent : &Tangent) -> Point {
         let mut p = *point;
         let mut t = *tangent;
@@ -931,8 +961,9 @@
         config : CylinderConfig { newton_iters : 20 },
     };
 
+    /// Tests for correct distance calculation beween two points on the same cap.
     #[test]
-    fn intra_cap_log_dist() {
+    fn intra_cap_dist() {
         let π = f64::PI;
         let p1 = CYL.on_top(0.0, 0.5);
         let p2 = CYL.on_top(π, 0.5);
@@ -943,8 +974,9 @@
         check_distance!(p3.dist_to(&p1), 0.5_f64.sqrt());
     }
 
+    /// Tests for correct distance calculation beween two points on the side.
     #[test]
-    fn intra_side_log_dist() {
+    fn intra_side_dist() {
         let π = f64::PI;
         let p1 = CYL.on_side(0.0, 0.0);
         let p2 = CYL.on_side(0.0, 0.4);
@@ -954,8 +986,10 @@
         check_distance!(p1.dist_to(&p3), π/2.0*CYL.radius);
     }
 
+    /// Tests for correct distance calculation beween two points on the side, when the corresponding
+    /// minimal geodesic has to cross a cap.
     #[test]
-    fn intra_side_over_cap_log_dist() {
+    fn intra_side_over_cap_dist() {
         let π = f64::PI;
         let off = 0.05;
         let z = CYL.top_z() - off;
@@ -965,8 +999,9 @@
         check_distance!(p1.dist_to(&p2), 2.0 * (CYL.radius + off));
     }
 
+    /// Tests for correct distance calculation between points on top and bottom caps.
     #[test]
-    fn top_bottom_log_dist() {
+    fn top_bottom_dist() {
         let π = f64::PI;
         let p1 = CYL.on_top(0.0, 0.0);
         let p2 = CYL.on_bottom(0.0, 0.0);
@@ -981,14 +1016,16 @@
         check_distance!(p1.dist_to(&p3), 2.0 * CYL.radius + CYL.height);
     }
 
+    /// Test for correct distance calculation between points on the side and a cap.
     #[test]
-    fn top_side_log_dist() {
+    fn top_side_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);
     }
 
+    /// Tests for correct tangent calculation from a cap to the side.
     #[test]
     fn cap_side_partial_exp() {
         let π = f64::PI;
@@ -1023,6 +1060,7 @@
         }
     }
 
+    /// Tests for correct tangent calculation from the side to a cap.
     #[test]
     fn side_top_exp() {
         let π = f64::PI;

--- a/src/fb.rs	Fri Dec 06 14:27:14 2024 -0500
+++ b/src/fb.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -11,11 +11,13 @@
 
 /// Trait for function objects that implement gradients
 pub trait Grad<M : ManifoldPoint> {
+    /// Calculates the gradient of `self` at `x`.
     fn grad(&self, x : &M) -> M::Tangent;
 }
 
 /// Trait for function objects that implement gradient steps
 pub trait Desc<M : ManifoldPoint> {
+    /// Calculates the gradient steps of `self` at `x` for the step length `τ`.
     fn desc(&self, τ : f64, x : M) -> M;
 }
 
@@ -39,6 +41,7 @@
 
 /// Trait for function objects that implement proximal steps
 pub trait Prox<M : ManifoldPoint> {
+    /// Calculates the proximap map of `self` at `x` for the step length `τ`.
     fn prox(&self, τ : f64, x : M) -> M;
 }
 

--- a/src/main.rs	Fri Dec 06 14:27:14 2024 -0500
+++ b/src/main.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -1,6 +1,5 @@
-/*!
-Optimisation on non-Riemannian manifolds.
-*/
+// Optimisation on non-Riemannian manifolds.
+#![doc = include_str!(concat!(env!("OUT_DIR"), "/README_uglified.md"))]
 
 // We use unicode. We would like to use much more of it than Rust allows.
 // Live with it. Embrace it.
@@ -140,7 +139,7 @@
 }
 
 
-/// A simple test on the cube
+/// A simple test on a [`Cylinder`].
 fn simple_cylinder_test(dir : &str) -> DynError {
     let π = f64::PI;
 
@@ -202,7 +201,7 @@
 }
 
 
-/// Runs [forward_backward] and saves the results.
+/// Runs [`forward_backward`] and saves the results.
 pub fn run_and_save<F, G, M, I>(
     dir : &str,
     name : &str,

--- a/src/newton.rs	Fri Dec 06 14:27:14 2024 -0500
+++ b/src/newton.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -4,6 +4,10 @@
 
 use alg_tools::types::*;
 
+/// Approximately solves $f(x)=b$ using Newton's method in 1D.
+///
+/// The function `g` should return $(f'(x), f(x))$.
+/// The initial iterate will be `x`, and exactly `iters` iterations are taken.
 #[inline]
 pub fn newton_sym1x1<F : Float>(
     g : impl Fn(F) -> (F, F),
@@ -17,6 +21,11 @@
     x
 }
 
+/// Approximately solves $f(x)=b$ using Newton's method in "D.
+///
+/// The function `g` should return $(∇f(x), f(x))$.
+/// The Hessian $A=∇f(x)$ should be symmetric and given in the form $[A_{11}, A_{12}, A_{22}]$.
+/// The initial iterate will be `[x1, x2]`, and exactly `iters` iterations are taken.
 #[inline]
 pub fn newton_sym2x2<F : Float>(
     g : impl Fn(F, F) -> ([F; 3], [F; 2]),

--- a/src/scaled.rs	Fri Dec 06 14:27:14 2024 -0500
+++ b/src/scaled.rs	Fri Dec 06 14:57:11 2024 -0500
@@ -8,11 +8,14 @@
 
 /// Structure for a function of type `G`, scaled by a scalar.
 pub struct Scaled<G> {
+    /// Scaling factor
     α : f64,
+    /// The base function
     g : G,
 }
 
 impl<G> Scaled<G> {
+    /// Creates a new scaled function.
     #[allow(dead_code)]
     pub fn new(α : f64, g : G) -> Self {
         Scaled{ α, g }