# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1708450457 18000
# Node ID 8aaa22fcd30268c7f031c0ec0e4420231d95b845
# Parent  2bee40c5e8349952943e39ea1089b645ca7df863# Parent  d14c877e14b70dd4c1ddc9436e4095e84964ef20
merge default

diff -r d14c877e14b7 -r 8aaa22fcd302 Cargo.lock
--- a/Cargo.lock	Tue Feb 20 12:33:16 2024 -0500
+++ b/Cargo.lock	Tue Feb 20 12:34:17 2024 -0500
@@ -4,7 +4,7 @@
 
 [[package]]
 name = "alg_tools"
-version = "0.1.0"
+version = "0.2.0-dev"
 dependencies = [
  "colored",
  "cpu-time",
diff -r d14c877e14b7 -r 8aaa22fcd302 Cargo.toml
--- a/Cargo.toml	Tue Feb 20 12:33:16 2024 -0500
+++ b/Cargo.toml	Tue Feb 20 12:34:17 2024 -0500
@@ -1,6 +1,6 @@
 [package]
 name = "alg_tools"
-version = "0.1.0"
+version = "0.2.0-dev"
 edition = "2021"
 rust-version = "1.67"
 authors = ["Tuomo Valkonen <tuomov@iki.fi>"]
diff -r d14c877e14b7 -r 8aaa22fcd302 src/bisection_tree/btfn.rs
--- a/src/bisection_tree/btfn.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/bisection_tree/btfn.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -4,7 +4,7 @@
 use std::marker::PhantomData;
 use std::sync::Arc;
 use crate::types::Float;
-use crate::mapping::{Apply, Mapping};
+use crate::mapping::{Apply, Mapping, Differentiable};
 //use crate::linops::{Apply, Linear};
 use crate::sets::Set;
 use crate::sets::Cube;
@@ -386,10 +386,8 @@
 
 make_btfn_unaryop!(Neg, neg);
 
-
-
 //
-// Mapping
+// Apply, Mapping, Differentiate
 //
 
 impl<'a, F : Float, G, BT, V, const N : usize> Apply<&'a Loc<F, N>>
@@ -422,6 +420,42 @@
     }
 }
 
+impl<'a, F : Float, G, BT, V, const N : usize> Differentiable<&'a Loc<F, N>>
+for BTFN<F, G, BT, N>
+where BT : BTImpl<F, N>,
+      G : SupportGenerator<F, N, Id=BT::Data>,
+      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Differentiable<&'a Loc<F, N>, Output = V>,
+      V : Sum {
+
+    type Output = V;
+
+    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.bt.iter_at(x)
+            .map(|&d| self.generator.support_for(d).differential(x))
+            .sum()
+    }
+}
+
+impl<F : Float, G, BT, V, const N : usize> Differentiable<Loc<F, N>>
+for BTFN<F, G, BT, N>
+where BT : BTImpl<F, N>,
+      G : SupportGenerator<F, N, Id=BT::Data>,
+      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Differentiable<Loc<F, N>, Output = V>,
+      V : Sum {
+
+    type Output = V;
+
+    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+        self.bt.iter_at(&x)
+            .map(|&d| self.generator.support_for(d).differential(x))
+            .sum()
+    }
+}
+
+//
+// GlobalAnalysis
+//
+
 impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg>
 for BTFN<F, G, BT, N>
 where BT : BTImpl<F, N>,
diff -r d14c877e14b7 -r 8aaa22fcd302 src/bisection_tree/either.rs
--- a/src/bisection_tree/either.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/bisection_tree/either.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -3,8 +3,8 @@
 use std::sync::Arc;
 
 use crate::types::*;
-use crate::mapping::Apply;
-use crate::iter::{Mappable,MapF,MapZ};
+use crate::mapping::{Apply, Differentiable};
+use crate::iter::{Mappable, MapF, MapZ};
 use crate::sets::Cube;
 use crate::loc::Loc;
 
@@ -190,6 +190,19 @@
     }
 }
 
+impl<F, S1, S2, X> Differentiable<X> for EitherSupport<S1, S2>
+where S1 : Differentiable<X, Output=F>,
+      S2 : Differentiable<X, Output=F> {
+    type Output = F;
+    #[inline]
+    fn differential(&self, x : X) -> F {
+        match self {
+            EitherSupport::Left(ref a) => a.differential(x),
+            EitherSupport::Right(ref b) => b.differential(x),
+        }
+    }
+}
+
 macro_rules! make_either_scalarop_rhs {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl<F : Float, G1, G2>
diff -r d14c877e14b7 -r 8aaa22fcd302 src/bisection_tree/support.rs
--- a/src/bisection_tree/support.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/bisection_tree/support.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -6,7 +6,7 @@
 use std::ops::{MulAssign,DivAssign,Neg};
 use crate::types::{Float, Num};
 use crate::maputil::map2;
-use crate::mapping::Apply;
+use crate::mapping::{Apply, Differentiable};
 use crate::sets::Cube;
 use crate::loc::Loc;
 use super::aggregator::Bounds;
@@ -145,6 +145,24 @@
     }
 }
 
+impl<'a, T, V, F : Float, const N : usize> Differentiable<&'a Loc<F, N>> for Shift<T,F,N>
+where T : Differentiable<Loc<F, N>, Output=V> {
+    type Output = V;
+    #[inline]
+    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.base_fn.differential(x - &self.shift)
+    }
+}
+
+impl<'a, T, V, F : Float, const N : usize> Differentiable<Loc<F, N>> for Shift<T,F,N>
+where T : Differentiable<Loc<F, N>, Output=V> {
+    type Output = V;
+    #[inline]
+    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+        self.base_fn.differential(x - &self.shift)
+    }
+}
+
 impl<'a, T, F : Float, const N : usize> Support<F,N> for Shift<T,F,N>
 where T : Support<F, N> {
     #[inline]
@@ -232,6 +250,28 @@
     }
 }
 
+impl<'a, T, V, F : Float, C, const N : usize> Differentiable<&'a Loc<F, N>> for Weighted<T, C>
+where T : for<'b> Differentiable<&'b Loc<F, N>, Output=V>,
+      V : std::ops::Mul<F, Output=V>,
+      C : Constant<Type=F> {
+    type Output = V;
+    #[inline]
+    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.base_fn.differential(x) * self.weight.value()
+    }
+}
+
+impl<'a, T, V, F : Float, C, const N : usize> Differentiable<Loc<F, N>> for Weighted<T, C>
+where T : Differentiable<Loc<F, N>, Output=V>,
+      V : std::ops::Mul<F, Output=V>,
+      C : Constant<Type=F> {
+    type Output = V;
+    #[inline]
+    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+        self.base_fn.differential(x) * self.weight.value()
+    }
+}
+
 impl<'a, T, F : Float, C, const N : usize> Support<F,N> for Weighted<T, C>
 where T : Support<F, N>,
       C : Constant<Type=F> {
diff -r d14c877e14b7 -r 8aaa22fcd302 src/iterate.rs
--- a/src/iterate.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/iterate.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -126,6 +126,9 @@
 
     /// Returns the current iteration count.
     fn iteration(&self) -> usize;
+
+    /// Indicates whether the iterator is quiet
+    fn is_quiet(&self) -> bool;
 }
 
 /// Result of a step of an [`AlgIterator`]
@@ -467,6 +470,8 @@
     verbose : bool,
     /// Whether results should be calculated.
     calc : bool,
+    /// Indicates whether the iteration is quiet
+    quiet : bool,
 }
 
 /// [`AlgIteratorFactory`] for [`BasicAlgIterator`]
@@ -585,10 +590,11 @@
     fn state(&self) -> BasicState {
         let iter = self.iter;
         let verbose = self.options.verbose_iter.is_verbose(iter);
-        BasicState{
+        BasicState {
             iter : iter,
             verbose : verbose,
             calc : verbose,
+            quiet : self.options.quiet
         }
     }
 }
@@ -605,7 +611,12 @@
 
     #[inline]
     fn iteration(&self) -> usize {
-        return self.iter;
+        self.iter
+    }
+
+    #[inline]
+    fn is_quiet(&self) -> bool {
+        self.quiet
     }
 }
 
diff -r d14c877e14b7 -r 8aaa22fcd302 src/loc.rs
--- a/src/loc.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/loc.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -135,6 +135,14 @@
     }
 }
 
+impl<F> Loc<F, 1> {
+    #[inline]
+    pub fn flatten1d(self) -> F {
+        let Loc([v]) = self;
+        v
+    }
+}
+
 impl<F, const N : usize> From<Loc<F, N>> for [F; N] {
     #[inline]
     fn from(other : Loc<F, N>) -> [F; N] {
@@ -237,6 +245,20 @@
 make_binop!(Add, add, AddAssign, add_assign);
 make_binop!(Sub, sub, SubAssign, sub_assign);
 
+impl<F : Float, const N : usize> std::iter::Sum for Loc<F, N> {
+    fn sum<I: Iterator<Item = Loc<F, N>>>(mut iter: I) -> Self {
+        match iter.next() {
+            None => Self::ORIGIN,
+            Some(mut v) => {
+                for w in iter {
+                    v += w
+                }
+                v
+            }
+        }
+    }
+}
+
 macro_rules! make_scalarop_rhs {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl<F : Num, const N : usize> $trait<F> for Loc<F, N> {
diff -r d14c877e14b7 -r 8aaa22fcd302 src/mapping.rs
--- a/src/mapping.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/mapping.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -42,39 +42,59 @@
     type Codomain = Codomain;
 }
 
-
-/// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to `F` a [`Float`].
-pub trait RealMapping<F : Float, const N : usize> : Mapping<Loc<F, N>, Codomain = F> {}
+/// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<F, N>`] to `F`.
+pub trait RealMapping<F : Float, const N : usize>
+: Mapping<Loc<F, N>, Codomain = F> {}
 
 impl<F : Float, T, const N : usize> RealMapping<F, N> for T
 where T : Mapping<Loc<F, N>, Codomain = F> {}
 
+/// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from
+/// [`Loc<F, N>`] to `F`.
+pub trait DifferentiableRealMapping<F : Float, const N : usize>
+: DifferentiableMapping<Loc<F, N>, Codomain = F, Differential=Loc<F, N>> {}
+
+impl<F : Float, T, const N : usize> DifferentiableRealMapping<F, N> for T
+where T : DifferentiableMapping<Loc<F, N>, Codomain = F, Differential=Loc<F, N>> {}
+
+
+/// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to [`Loc<F, M>`].
+pub trait RealVectorField<F : Float, const N : usize, const M : usize>
+: Mapping<Loc<F, N>, Codomain = Loc<F, M>> {}
+
+impl<F : Float, T, const N : usize, const M : usize> RealVectorField<F, N, M> for T
+where T : Mapping<Loc<F, N>, Codomain = Loc<F, M>> {}
+
 
 /// Trait for calculation the differential of `Self` as a mathematical function on `X`.
-pub trait Differentiate<X> {
+pub trait Differentiable<X> : Sized {
     type Output;
 
     /// Compute the differential of `self` at `x`.
     fn differential(&self, x : X) -> Self::Output;
 }
 
-
 /// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials
 /// `Differential`.
 ///
 /// This is automatically implemented when the relevant [`Differentiate`] are implemented.
 pub trait DifferentiableMapping<Domain>
 : Mapping<Domain>
-  + Differentiate<Domain, Output=Self::Differential>
-  + for<'a> Differentiate<&'a Domain, Output=Self::Differential>{
+  + Differentiable<Domain, Output=Self::Differential>
+  + for<'a> Differentiable<&'a Domain, Output=Self::Differential> {
     type Differential;
+
+    /// Form the differential mapping of `self`.
+    fn diff(self) -> Differential<Domain, Self> {
+        Differential{ g : self, _space : PhantomData }
+    }
 }
 
 
 impl<Domain, Differential, T> DifferentiableMapping<Domain> for T
 where T : Mapping<Domain>
-          + Differentiate<Domain, Output=Differential>
-          + for<'a> Differentiate<&'a Domain, Output=Differential> {
+          + Differentiable<Domain, Output=Differential>
+          + for<'a> Differentiable<&'a Domain, Output=Differential> {
     type Differential = Differential;
 }
 
@@ -103,7 +123,7 @@
     }
 }
 
-impl<Domain, M> Differentiate<Domain> for Sum<Domain, M>
+impl<Domain, M> Differentiable<Domain> for Sum<Domain, M>
 where M : DifferentiableMapping<Domain>,
       M :: Codomain : std::iter::Sum,
       M :: Differential : std::iter::Sum,
@@ -115,3 +135,101 @@
         self.components.iter().map(|c| c.differential(x)).sum()
     }
 }
+
+/// Container for the differential [`Mapping`] of a [`Differentiable`] mapping.
+pub struct Differential<X, G : DifferentiableMapping<X>> {
+    g : G,
+    _space : PhantomData<X>
+}
+
+impl<X, G : DifferentiableMapping<X>> Apply<X> for Differential<X, G> {
+    type Output = G::Differential;
+
+    #[inline]
+    fn apply(&self, x : X) -> G::Differential {
+        self.g.differential(x)
+    }
+}
+
+impl<'a, X, G : DifferentiableMapping<X>> Apply<&'a X> for Differential<X, G> {
+    type Output = G::Differential;
+
+    #[inline]
+    fn apply(&self, x : &'a X) -> G::Differential {
+        self.g.differential(x)
+    }
+}
+
+
+/// Container for flattening [`Loc`]`<F, 1>` codomain of a [`Mapping`] to `F`.
+pub struct FlattenedCodomain<X, F, G : Mapping<X, Codomain=Loc<F, 1>>> {
+    g : G,
+    _phantoms : PhantomData<(X, F)>
+}
+
+impl<X, F, G : Mapping<X, Codomain=Loc<F, 1>>> Apply<X> for FlattenedCodomain<X, F, G> {
+    type Output = F;
+
+    #[inline]
+    fn apply(&self, x : X) -> Self::Output {
+        self.g.apply(x).flatten1d()
+    }
+}
+
+/// An auto-trait for constructing a [`FlattenCodomain`] structure for
+/// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`.
+pub trait FlattenCodomain<X, F> : Mapping<X, Codomain=Loc<F, 1>> + Sized {
+    /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
+    fn flatten_codomain(self) -> FlattenedCodomain<X, F, Self> {
+        FlattenedCodomain{ g : self, _phantoms : PhantomData }
+    }
+}
+
+impl<X, F, G : Sized + Mapping<X, Codomain=Loc<F, 1>>> FlattenCodomain<X, F> for G {}
+
+
+/// Container for dimensional slicing [`Loc`]`<F, N>` codomain of a [`Mapping`] to `F`.
+pub struct SlicedCodomain<X, F, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> {
+    g : G,
+    slice : usize,
+    _phantoms : PhantomData<(X, F)>
+}
+
+impl<X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> Apply<X>
+for SlicedCodomain<X, F, G, N> {
+    type Output = F;
+
+    #[inline]
+    fn apply(&self, x : X) -> Self::Output {
+        let tmp : [F; N] = self.g.apply(x).into();
+        // Safety: `slice_codomain` below checks the range.
+        unsafe { *tmp.get_unchecked(self.slice) }
+    }
+}
+
+impl<'a, X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> Apply<&'a X>
+for SlicedCodomain<X, F, G, N> {
+    type Output = F;
+
+    #[inline]
+    fn apply(&self, x : &'a X) -> Self::Output {
+        let tmp : [F; N] = self.g.apply(x).into();
+        // Safety: `slice_codomain` below checks the range.
+        unsafe { *tmp.get_unchecked(self.slice) }
+    }
+}
+
+/// An auto-trait for constructing a [`FlattenCodomain`] structure for
+/// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`.
+pub trait SliceCodomain<X, F : Copy, const N : usize> : Mapping<X, Codomain=Loc<F, N>> + Sized {
+    /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
+    fn slice_codomain(self, slice : usize) -> SlicedCodomain<X, F, Self, N> {
+        assert!(slice < N);
+        SlicedCodomain{ g : self, slice, _phantoms : PhantomData }
+    }
+}
+
+impl<X, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>>, const N : usize>
+SliceCodomain<X, F, N>
+for G {}
+
diff -r d14c877e14b7 -r 8aaa22fcd302 src/norms.rs
--- a/src/norms.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/norms.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -12,7 +12,7 @@
 
 /// An exponent for norms.
 ///
-// Just a collection of desirabl attributes for a marker type
+// Just a collection of desirable attributes for a marker type
 pub trait NormExponent : Copy + Send + Sync + 'static {}
 
 
diff -r d14c877e14b7 -r 8aaa22fcd302 src/types.rs
--- a/src/types.rs	Tue Feb 20 12:33:16 2024 -0500
+++ b/src/types.rs	Tue Feb 20 12:34:17 2024 -0500
@@ -84,6 +84,7 @@
     const SQRT_2 : Self;
     const INFINITY : Self;
     const NEG_INFINITY : Self;
+    const NAN : Self;
     const FRAC_2_SQRT_PI : Self;
 }
 
@@ -133,6 +134,7 @@
     const SQRT_2 : Self = std::f64::consts::SQRT_2;
     const INFINITY : Self = std::f64::INFINITY;
     const NEG_INFINITY : Self = std::f64::NEG_INFINITY;
+    const NAN : Self = std::f64::NAN;
     const FRAC_2_SQRT_PI : Self = std::f64::consts::FRAC_2_SQRT_PI;
 }
 
@@ -150,6 +152,7 @@
     const SQRT_2 : Self = std::f32::consts::SQRT_2;
     const INFINITY : Self = std::f32::INFINITY;
     const NEG_INFINITY : Self = std::f32::NEG_INFINITY;
+    const NAN : Self = std::f32::NAN;
     const FRAC_2_SQRT_PI : Self = std::f32::consts::FRAC_2_SQRT_PI;
 }