geop_core_sketch/
dual.rs

1//! Forward-mode automatic differentiation for constraint residuals.
2//!
3//! Every residual is written once, generically over [`Scalar`] (the same
4//! trait the rest of the kernel runs on), and evaluated either with a plain
5//! [`Scalar`] value (the value only) or with [`Dual`] (the value plus its
6//! exact partial derivatives with respect to the constraint's own
7//! variables). A constraint touches only a handful of variables, so a dual
8//! number carries a small fixed-size gradient rather than one entry per
9//! sketch variable.
10//!
11//! Building `Dual` over an interval [`Scalar`] rather than a bare `f64`
12//! means a residual's value is a rigorous enclosure of the true real number
13//! it stands for, not a float that silently absorbs whatever rounding `sin`,
14//! `sqrt` or `PI` introduced along the way — the sketch solver's own
15//! uncertainty about an irrational quantity shows up as interval width
16//! instead of vanishing.
17
18use geop_core_math::{
19    geop_error::GeopResult,
20    scalars::{Field, Ring, Scalar},
21};
22
23/// The most variables a single constraint may depend on. The largest
24/// constraints (tangency or equality between two arcs) touch two arcs of
25/// five variables each.
26pub const MAX_LOCAL_VARS: usize = 10;
27
28/// A value together with its gradient with respect to up to
29/// [`MAX_LOCAL_VARS`] seeded variables.
30#[derive(Clone, Copy, Debug)]
31pub struct Dual<S: Scalar> {
32    pub v: S,
33    pub d: [S; MAX_LOCAL_VARS],
34}
35
36impl<S: Scalar> core::fmt::Display for Dual<S> {
37    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
38        core::fmt::Display::fmt(&self.v, f)
39    }
40}
41
42impl<S: Scalar> Default for Dual<S> {
43    fn default() -> Self {
44        Dual::cst(S::default())
45    }
46}
47
48impl<S: Scalar> Dual<S> {
49    /// A constant, with no dependence on any variable.
50    pub fn cst(v: S) -> Self {
51        Dual {
52            v,
53            d: [S::ZERO; MAX_LOCAL_VARS],
54        }
55    }
56
57    /// The `slot`-th independent variable, with value `v`.
58    pub fn var(v: S, slot: usize) -> Self {
59        let mut d = [S::ZERO; MAX_LOCAL_VARS];
60        d[slot] = S::ONE;
61        Dual { v, d }
62    }
63
64    /// `f(self)` for a scalar function with value `v` and derivative `df =
65    /// f'(self.v)`: `d <- d * df` by the chain rule.
66    fn chain(self, v: S, df: S) -> Self {
67        let mut d = self.d;
68        for x in &mut d {
69            *x = x.mul(df);
70        }
71        Dual { v, d }
72    }
73}
74
75impl<S: Scalar> Ring for Dual<S> {
76    fn add(self, o: Self) -> Self {
77        let mut d = self.d;
78        for (x, y) in d.iter_mut().zip(o.d) {
79            *x = x.add(y);
80        }
81        Dual {
82            v: self.v.add(o.v),
83            d,
84        }
85    }
86
87    fn sub(self, o: Self) -> Self {
88        let mut d = self.d;
89        for (x, y) in d.iter_mut().zip(o.d) {
90            *x = x.sub(y);
91        }
92        Dual {
93            v: self.v.sub(o.v),
94            d,
95        }
96    }
97
98    fn mul(self, o: Self) -> Self {
99        let mut d = [S::ZERO; MAX_LOCAL_VARS];
100        for (i, x) in d.iter_mut().enumerate() {
101            *x = self.d[i].mul(o.v).add(o.d[i].mul(self.v));
102        }
103        Dual {
104            v: self.v.mul(o.v),
105            d,
106        }
107    }
108
109    fn neg(self) -> Self {
110        self.chain(self.v.neg(), S::ZERO.sub(S::ONE))
111    }
112}
113
114impl<S: Scalar> Field for Dual<S> {
115    fn div(self, o: Self) -> GeopResult<Self> {
116        let v = self.v.div(o.v)?;
117        let mut d = [S::ZERO; MAX_LOCAL_VARS];
118        for (i, x) in d.iter_mut().enumerate() {
119            // Quotient rule: (self' - v * o') / o.v.
120            *x = self.d[i].sub(v.mul(o.d[i])).div(o.v)?;
121        }
122        Ok(Dual { v, d })
123    }
124}
125
126// `could_be_equal`/`is_infinite`/etc. below act on the value alone, ignoring
127// the gradient — a residual's dual number is compared and reported on the
128// same terms a plain [`Scalar`] would be, and these interval-lattice
129// operations (`intersect`, `midpoint`, `width`, ...) exist on [`Dual`] only
130// because [`Scalar`] requires them, not because a constraint ever calls
131// them: the sketch solver only ever adds, multiplies, divides, and takes
132// `sqrt`/`sin`/`cos`/`abs` of a residual.
133impl<S: Scalar> Scalar for Dual<S> {
134    const ZERO: Self = Dual {
135        v: S::ZERO,
136        d: [S::ZERO; MAX_LOCAL_VARS],
137    };
138    const ONE: Self = Dual {
139        v: S::ONE,
140        d: [S::ZERO; MAX_LOCAL_VARS],
141    };
142    const TWO: Self = Dual {
143        v: S::TWO,
144        d: [S::ZERO; MAX_LOCAL_VARS],
145    };
146    const E: Self = Dual {
147        v: S::E,
148        d: [S::ZERO; MAX_LOCAL_VARS],
149    };
150    const PI: Self = Dual {
151        v: S::PI,
152        d: [S::ZERO; MAX_LOCAL_VARS],
153    };
154    const INFINITY: Self = Dual {
155        v: S::INFINITY,
156        d: [S::ZERO; MAX_LOCAL_VARS],
157    };
158    const ENTIRE: Self = Dual {
159        v: S::ENTIRE,
160        d: [S::ZERO; MAX_LOCAL_VARS],
161    };
162
163    fn from_i64(v: i64) -> Self {
164        Dual::cst(S::from_i64(v))
165    }
166
167    fn from_f64(v: f64) -> Self {
168        Dual::cst(S::from_f64(v))
169    }
170
171    fn from_ratio(num: i64, den: i64) -> GeopResult<Self> {
172        Ok(Dual::cst(S::from_ratio(num, den)?))
173    }
174
175    fn to_f64(self) -> f64 {
176        self.v.to_f64()
177    }
178
179    fn abs(self) -> Self {
180        // Not differentiable at a sign change; 0 is the subgradient (the
181        // straddling case falls through to the `+1` branch, matching
182        // `Scalar::abs`'s own choice to keep the value non-negative there).
183        if self.v.definitely_less(S::ZERO) {
184            self.chain(self.v.abs(), S::ZERO.sub(S::ONE))
185        } else {
186            self.chain(self.v.abs(), S::ONE)
187        }
188    }
189
190    fn sqrt(self) -> GeopResult<Self> {
191        let s = self.v.sqrt()?;
192        // `sqrt` is not differentiable at 0. Residuals only take roots of
193        // squared lengths, whose gradient vanishes there too, so 0 is the
194        // right subgradient: it leaves the other terms in charge.
195        let df = if s.definitely_greater(S::ZERO) {
196            S::ONE.div(S::TWO.mul(s))?
197        } else {
198            S::ZERO
199        };
200        Ok(self.chain(s, df))
201    }
202
203    fn sin(self) -> Self {
204        self.chain(self.v.sin(), self.v.cos())
205    }
206
207    fn cos(self) -> Self {
208        let df = self.v.sin().neg();
209        self.chain(self.v.cos(), df)
210    }
211
212    fn could_be_equal(self, other: Self) -> bool {
213        self.v.could_be_equal(other.v)
214    }
215
216    fn definitely_not_equal(self, other: Self) -> bool {
217        self.v.definitely_not_equal(other.v)
218    }
219
220    fn could_be_greater(self, other: Self) -> bool {
221        self.v.could_be_greater(other.v)
222    }
223
224    fn definitely_greater(self, other: Self) -> bool {
225        self.v.definitely_greater(other.v)
226    }
227
228    fn could_be_less(self, other: Self) -> bool {
229        self.v.could_be_less(other.v)
230    }
231
232    fn definitely_less(self, other: Self) -> bool {
233        self.v.definitely_less(other.v)
234    }
235
236    fn is_infinite(self) -> bool {
237        self.v.is_infinite()
238    }
239
240    fn is_finite(self) -> bool {
241        self.v.is_finite()
242    }
243
244    fn midpoint(self) -> Self {
245        Dual {
246            v: self.v.midpoint(),
247            d: self.d,
248        }
249    }
250
251    fn is_sharp(self) -> bool {
252        self.v.is_sharp()
253    }
254
255    fn width(self) -> Self {
256        Dual::cst(self.v.width())
257    }
258
259    fn lower(self) -> Self {
260        Dual {
261            v: self.v.lower(),
262            d: self.d,
263        }
264    }
265
266    fn upper(self) -> Self {
267        Dual {
268            v: self.v.upper(),
269            d: self.d,
270        }
271    }
272
273    fn intersect(self, other: Self) -> Self {
274        Dual {
275            v: self.v.intersect(other.v),
276            d: self.d,
277        }
278    }
279
280    fn union(self, other: Self) -> Self {
281        Dual {
282            v: self.v.union(other.v),
283            d: self.d,
284        }
285    }
286
287    fn is_subset_of(self, other: Self) -> bool {
288        self.v.is_subset_of(other.v)
289    }
290}
291
292#[cfg(test)]
293mod tests {
294    use super::*;
295    use geop_core_math::scalars::scal_in_f64::ScalInF64;
296
297    /// Every operation's derivative against a central difference.
298    #[test]
299    fn dual_matches_finite_differences() {
300        fn f<S: Scalar>(x: S, y: S) -> S {
301            let numerator = x.mul(y).add(x.sin().mul(y.cos()));
302            let denom = x.mul(x).add(S::ONE).sqrt().unwrap();
303            numerator.div(denom).unwrap().sub(y.sub(x).abs())
304        }
305        let (x, y) = (0.7, -1.3);
306        let d = f(
307            Dual::<ScalInF64>::var(ScalInF64::from_f64(x), 0),
308            Dual::<ScalInF64>::var(ScalInF64::from_f64(y), 1),
309        );
310        let h = 1e-6;
311        let fv = |x: f64, y: f64| f(ScalInF64::from_f64(x), ScalInF64::from_f64(y)).to_f64();
312        let dx = (fv(x + h, y) - fv(x - h, y)) / (2.0 * h);
313        let dy = (fv(x, y + h) - fv(x, y - h)) / (2.0 * h);
314        assert!((d.v.to_f64() - fv(x, y)).abs() < 1e-15);
315        assert!((d.d[0].to_f64() - dx).abs() < 1e-8, "{} vs {dx}", d.d[0]);
316        assert!((d.d[1].to_f64() - dy).abs() < 1e-8, "{} vs {dy}", d.d[1]);
317    }
318}
geop_core_sketch/dual.rs

geop_core_sketch/
dual.rs
