Skip to main content

basin/problems/
sphere.rs

1//! N-dimensional Sphere function.
2//!
3//! `f(x) = Σᵢ xᵢ²`
4//!
5//! Smooth, convex, separable, unimodal. Global minimum at `x = (0, …, 0)`
6//! with `f = 0`. The trivial canary problem — every solver should solve it
7//! cleanly; failure indicates the implementation is broken.
8
9use core::marker::PhantomData;
10
11use super::spec::{Dimensionality, HasSpec, ProblemSpec, Properties, Reference};
12use crate::{CostFunction, Gradient};
13
14/// Evaluates the Sphere function at `x`.
15pub fn sphere(x: &[f64]) -> f64 {
16    x.iter().map(|v| v * v).sum()
17}
18
19/// Writes the Sphere gradient at `x` into `out`. Lengths must match.
20pub fn sphere_gradient(x: &[f64], out: &mut [f64]) {
21    debug_assert_eq!(x.len(), out.len());
22    for (g, &v) in out.iter_mut().zip(x.iter()) {
23        *g = 2.0 * v;
24    }
25}
26
27/// Pre-wrapped Sphere problem. Generic over the parameter backend `P`;
28/// the default `P = Vec<f64>` lets you write `Sphere::default()` for the
29/// common case. Backend impls (`nalgebra::DVector<f64>`, `ndarray::Array1<f64>`,
30/// `faer::Col<f64>`) are gated behind their respective features.
31pub struct Sphere<P = Vec<f64>>(PhantomData<fn() -> P>);
32
33impl<P> Sphere<P> {
34    /// Build a freshly typed problem instance. Pair with one of the
35    /// backend-specific impl blocks (Vec, nalgebra, ndarray, faer).
36    pub const fn new() -> Self {
37        Self(PhantomData)
38    }
39}
40
41impl<P> Default for Sphere<P> {
42    fn default() -> Self {
43        Self::new()
44    }
45}
46
47/// Catalogue entry for this problem.
48pub static SPHERE_SPEC: ProblemSpec = ProblemSpec {
49    name: "Sphere",
50    dim: Dimensionality::NDimensional { min: 1 },
51    properties: Properties {
52        smooth: true,
53        differentiable: true,
54        convex: true,
55        unimodal: true,
56        separable: true,
57        scalable: true,
58    },
59    references: &[Reference {
60        citation: "De Jong (1975)",
61        title: "An Analysis of the Behavior of a Class of Genetic Adaptive Systems",
62        source: "PhD thesis, University of Michigan",
63        doi: None,
64        url: Some("https://hdl.handle.net/2027.42/4507"),
65    }],
66    description: "Sum of squares: f(x) = Σ xᵢ². Convex, separable, unimodal. \
67                  Global minimum at x = (0, …, 0), value 0. The canonical \
68                  trivial canary — every solver should solve it cleanly.",
69};
70
71impl<P> HasSpec for Sphere<P> {
72    const SPEC: &'static ProblemSpec = &SPHERE_SPEC;
73}
74
75impl CostFunction for Sphere<Vec<f64>> {
76    type Param = Vec<f64>;
77    type Output = f64;
78    type Error = std::convert::Infallible;
79    fn cost(&self, x: &Vec<f64>) -> Result<f64, std::convert::Infallible> {
80        Ok(sphere(x))
81    }
82}
83
84impl Gradient for Sphere<Vec<f64>> {
85    type Gradient = Vec<f64>;
86    fn gradient(&self, x: &Vec<f64>) -> Result<Vec<f64>, std::convert::Infallible> {
87        let mut out = vec![0.0; x.len()];
88        sphere_gradient(x, &mut out);
89        Ok(out)
90    }
91}
92
93#[cfg(feature = "nalgebra")]
94mod nalgebra_impl {
95    use super::{Sphere, sphere, sphere_gradient};
96    use crate::{CostFunction, Gradient};
97    use nalgebra::DVector;
98
99    impl CostFunction for Sphere<DVector<f64>> {
100        type Param = DVector<f64>;
101        type Output = f64;
102        type Error = std::convert::Infallible;
103        fn cost(&self, x: &DVector<f64>) -> Result<f64, std::convert::Infallible> {
104            Ok(sphere(x.as_slice()))
105        }
106    }
107
108    impl Gradient for Sphere<DVector<f64>> {
109        type Gradient = DVector<f64>;
110        fn gradient(&self, x: &DVector<f64>) -> Result<DVector<f64>, std::convert::Infallible> {
111            let mut out = DVector::zeros(x.len());
112            sphere_gradient(x.as_slice(), out.as_mut_slice());
113            Ok(out)
114        }
115    }
116}
117
118#[cfg(feature = "ndarray")]
119mod ndarray_impl {
120    use super::{Sphere, sphere, sphere_gradient};
121    use crate::{CostFunction, Gradient};
122    use ndarray::Array1;
123
124    impl CostFunction for Sphere<Array1<f64>> {
125        type Param = Array1<f64>;
126        type Output = f64;
127        type Error = std::convert::Infallible;
128        fn cost(&self, x: &Array1<f64>) -> Result<f64, std::convert::Infallible> {
129            Ok(sphere(x.as_slice().expect("Array1 is contiguous")))
130        }
131    }
132
133    impl Gradient for Sphere<Array1<f64>> {
134        type Gradient = Array1<f64>;
135        fn gradient(&self, x: &Array1<f64>) -> Result<Array1<f64>, std::convert::Infallible> {
136            let mut out = Array1::zeros(x.len());
137            sphere_gradient(
138                x.as_slice().expect("Array1 is contiguous"),
139                out.as_slice_mut().expect("Array1 is contiguous"),
140            );
141            Ok(out)
142        }
143    }
144}
145
146#[cfg(feature = "faer")]
147mod faer_impl {
148    use super::Sphere;
149    use crate::{CostFunction, Gradient};
150    use faer::Col;
151
152    impl CostFunction for Sphere<Col<f64>> {
153        type Param = Col<f64>;
154        type Output = f64;
155        type Error = std::convert::Infallible;
156        fn cost(&self, x: &Col<f64>) -> Result<f64, std::convert::Infallible> {
157            let n = x.nrows();
158            let mut s = 0.0;
159            for i in 0..n {
160                s += x[i] * x[i];
161            }
162            Ok(s)
163        }
164    }
165
166    impl Gradient for Sphere<Col<f64>> {
167        type Gradient = Col<f64>;
168        fn gradient(&self, x: &Col<f64>) -> Result<Col<f64>, std::convert::Infallible> {
169            let n = x.nrows();
170            Ok(Col::<f64>::from_fn(n, |i| 2.0 * x[i]))
171        }
172    }
173}
174
175#[cfg(test)]
176mod tests {
177    use super::*;
178
179    #[test]
180    fn sphere_minimum_is_zero_at_origin() {
181        assert_eq!(sphere(&[0.0]), 0.0);
182        assert_eq!(sphere(&[0.0, 0.0, 0.0, 0.0]), 0.0);
183    }
184
185    #[test]
186    fn sphere_known_value() {
187        assert_eq!(sphere(&[1.0, 2.0, 3.0]), 14.0);
188    }
189
190    #[test]
191    fn sphere_gradient_zero_at_origin() {
192        let mut g = vec![0.0; 5];
193        sphere_gradient(&[0.0; 5], &mut g);
194        for v in g {
195            assert_eq!(v, 0.0);
196        }
197    }
198
199    #[test]
200    fn sphere_gradient_matches_finite_difference() {
201        let x = [-1.2, 1.0, 0.7, 0.4];
202        let mut g = vec![0.0; x.len()];
203        sphere_gradient(&x, &mut g);
204        let h = 1e-6;
205        for i in 0..x.len() {
206            let mut xp = x;
207            let mut xm = x;
208            xp[i] += h;
209            xm[i] -= h;
210            let fd = (sphere(&xp) - sphere(&xm)) / (2.0 * h);
211            assert!((g[i] - fd).abs() < 1e-6, "i={i}, g={}, fd={fd}", g[i]);
212        }
213    }
214
215    #[test]
216    fn spec_is_wired_up_via_has_spec_trait() {
217        let spec = <Sphere<Vec<f64>> as HasSpec>::SPEC;
218        assert_eq!(spec.name, "Sphere");
219        assert!(spec.properties.convex);
220        assert!(spec.properties.separable);
221        assert!(spec.properties.unimodal);
222        assert!(matches!(spec.dim, Dimensionality::NDimensional { min: 1 }));
223        assert!(!spec.references.is_empty());
224    }
225}