1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

use ndarray::prelude::*;

/// `structs` implementing this trait represent regions in n-dimensional space.
pub trait Set {
    /// Returns whether or not a given vector `p` lies in the instance.
    fn contains(&self, p: &Array1<f64>) -> bool;

    /// Returns a bounding box for the set.
    /// This function is used for point process simulation by rejection, but can also be used to implement a Monte Carlo estimation of the set"s measure.
    fn bounding_box(&self) -> Array2<f64>;
}

/// General n-dimensional hyperrectangle.
pub struct Rectangle(Array1<f64>, Array1<f64>);

impl Rectangle {
    /// `close` is the point with smaller coordinates,
    /// `far` is the further point delimiting the rectangle,
    /// with the largest coordinates.
    pub fn new(close: Array1<f64>, far: Array1<f64>) -> Rectangle {
        Rectangle(close, far)
    }
}


impl Set for Rectangle {
    fn contains(&self, p: &Array1<f64>) -> bool {
        assert_eq!(p.len(), self.0.shape()[0]);
        
        // check if p is further away than the closer point
        let further = self.0.iter().zip(p.iter())
            .fold(true, |acc: bool, (v,w)| {
                acc & (w > v)
            });
        
        // check if p is closer than the far point
        let closer = self.1.iter().zip(p.iter())
            .fold(true, |acc: bool, (v,w)| {
                acc & (w < v)
            });

        // if both conditions are true then we're in the rectangle
        further & closer
    }

    fn bounding_box(&self) -> Array<f64, Ix2> {
        let mut result = unsafe {
            Array::uninitialized((self.0.shape()[0], 2))
        };
        
        result.slice_mut(s![0,..]).assign(&self.0);
        result.slice_mut(s![1,..]).assign(&self.1);
        result
    }
}

/// The n-dimensional ball.
pub struct Ball {
    center: Array1<f64>,
    radius: f64
}

impl Ball {
    pub fn new(center: Array1<f64>, radius: f64) -> Ball {
        assert!(radius > 0.0);

        Ball {
            center, radius
        }
    }
}

impl Set for Ball {
    fn contains(&self, p: &Array1<f64>) -> bool {
        let diff = &self.center - p;
        let distance = diff.dot(&diff).sqrt();
        distance <= self.radius
    }

    fn bounding_box(&self) -> Array<f64, Ix2> {
        // Dimension of current space
        let n = self.center.shape()[0];

        let mut res = unsafe {
            Array::uninitialized((n, 2))
        };

        for i in 0..n {
            res[[0,i]] = self.center[i] - self.radius;
            res[[1,i]] = self.center[i] + self.radius;
        }

        res
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn rectangle_test() {
        let close = array![0.0, 1.0];
        let far = array![1.0, 4.0];
        
        let rect = Rectangle::new(close, far);

        let p = array![0.5, 1.5];
        assert!(rect.contains(&p));

        let p = array![-1.0,2.0];
        assert!(!rect.contains(&p));

    }
}