Skip to main content

optimization_engine/constraints/
ball2.rs

1use super::Constraint;
2use crate::FunctionCallResult;
3use num::Float;
4use std::iter::Sum;
5
6#[derive(Copy, Clone)]
7/// A Euclidean ball, that is, a set given by $B_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert \leq r\\}$
8/// or a Euclidean ball centered at a point $x_c$, that is, $B_2^{x_c, r} = \\{x \in \mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert \leq r\\}$
9pub struct Ball2<'a, T = f64> {
10    center: Option<&'a [T]>,
11    radius: T,
12}
13
14impl<'a, T: Float> Ball2<'a, T> {
15    /// Construct a new Euclidean ball with given center and radius
16    /// If no `center` is given, then it is assumed to be in the origin
17    ///
18    /// # Example
19    ///
20    /// ```
21    /// use optimization_engine::constraints::{Ball2, Constraint};
22    ///
23    /// let ball = Ball2::new(None, 1.0);
24    /// let mut x = [2.0, 0.0];
25    /// ball.project(&mut x).unwrap();
26    /// ```
27    pub fn new(center: Option<&'a [T]>, radius: T) -> Self {
28        assert!(radius > T::zero());
29
30        Ball2 { center, radius }
31    }
32}
33
34impl<'a, T> Constraint<T> for Ball2<'a, T>
35where
36    T: Float + Sum<T>,
37{
38    fn project(&self, x: &mut [T]) -> FunctionCallResult {
39        if let Some(center) = &self.center {
40            assert_eq!(
41                x.len(),
42                center.len(),
43                "x and xc have incompatible dimensions"
44            );
45            let mut norm_difference = T::zero();
46            x.iter().zip(center.iter()).for_each(|(a, b)| {
47                let diff_ = *a - *b;
48                norm_difference = norm_difference + diff_ * diff_
49            });
50
51            norm_difference = norm_difference.sqrt();
52
53            if norm_difference > self.radius {
54                x.iter_mut().zip(center.iter()).for_each(|(x, c)| {
55                    *x = *c + self.radius * (*x - *c) / norm_difference;
56                });
57            }
58        } else {
59            let norm_x = crate::matrix_operations::norm2(x);
60            if norm_x > self.radius {
61                let norm_over_radius = norm_x / self.radius;
62                x.iter_mut().for_each(|x_| *x_ = *x_ / norm_over_radius);
63            }
64        }
65        Ok(())
66    }
67
68    fn is_convex(&self) -> bool {
69        true
70    }
71}