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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

//! Fractional Brownian motion and fractional Gaussian noise.

use probability::distribution::{Gaussian, Sample};
use probability::source::Source;

use {Process, Stationary};
use gaussian::circulant_embedding;

macro_rules! hurst(
    ($value:expr) => ({
        let value = $value;
        debug_assert!(value > 0.0 && value < 1.0);
        value
    });
);

macro_rules! step(
    ($value:expr) => ({
        let value = $value;
        debug_assert!(value > 0.0);
        value
    });
);

/// A fractional Brownian motion.
pub struct Motion {
    hurst: f64,
}

/// A fractional Gaussian noise.
pub struct Noise {
    hurst: f64,
    step: f64,
}

impl Motion {
    /// Create a fractional Brownian motion.
    #[inline]
    pub fn new(hurst: f64) -> Motion {
        Motion { hurst: hurst!(hurst) }
    }

    /// Generate a sample path.
    pub fn sample<S>(&self, points: usize, step: f64, source: &mut S) -> Vec<f64>
        where S: Source
    {
        match points {
            0 => vec![],
            1 => vec![0.0],
            _ => {
                let mut data = vec![0.0];
                data.extend(Noise::new(self.hurst, step).sample(points - 1, source));
                for i in 2..points {
                    data[i] += data[i - 1];
                }
                data
            },
        }
    }
}

impl Noise {
    /// Create a fractional Gaussian noise.
    #[inline]
    pub fn new(hurst: f64, step: f64) -> Noise {
        Noise { hurst: hurst!(hurst), step: step!(step) }
    }

    /// Generate a sample path.
    pub fn sample<S>(&self, points: usize, source: &mut S) -> Vec<f64>
        where S: Source
    {
        match points {
            0 => vec![],
            1 => vec![Gaussian::new(0.0, Stationary::var(self).sqrt()).sample(source)],
            _ => {
                let n = points - 1;
                let gaussian = Gaussian::new(0.0, 1.0);
                let data = circulant_embedding(self, n, || gaussian.sample(source));
                data.iter().take(points).map(|point| point.re).collect()
            },
        }
    }
}

impl Process for Motion {
    type Index = f64;
    type State = f64;

    fn cov(&self, t: f64, s: f64) -> f64 {
        debug_assert!(t >= 0.0 && s >= 0.0);
        let power = 2.0 * self.hurst;
        0.5 * (t.powf(power) + s.powf(power) - (t - s).abs().powf(power))
    }
}

impl Process for Noise {
    type Index = usize;
    type State = f64;

    #[inline]
    fn cov(&self, t: usize, s: usize) -> f64 {
        Stationary::cov(self, if t < s { s - t } else { t - s })
    }
}

impl Stationary for Noise {
    type Distance = usize;

    fn cov(&self, delta: usize) -> f64 {
        let delta = delta as f64;
        let power = 2.0 * self.hurst;
        let var = self.step.powf(power);
        0.5 * var * ((delta + 1.0).powf(power) - 2.0 * delta.powf(power) +
                     (delta - 1.0).abs().powf(power))
    }
}

#[cfg(test)]
mod tests {
    use super::{Motion, Noise};

    #[test]
    fn motion_var() {
        use Process;
        let process = Motion::new(0.25);
        let variances = (0..3).map(|i| process.var(i as f64)).collect::<Vec<_>>();
        assert_eq!(&variances, &[0.0, 1.0, 2f64.sqrt()]);
    }

    #[test]
    fn noise_var() {
        let process = Noise::new(0.42, 1.0);
        {
            use Process;
            let variances = (0..3).map(|i| process.var(i)).collect::<Vec<_>>();
            assert_eq!(&variances, &[1.0, 1.0, 1.0]);
        }
        {
            use Stationary;
            assert_eq!(process.var(), 1.0);
        }
        let process = Noise::new(0.25, 0.01);
        {
            use Process;
            let variances = (0..3).map(|i| process.var(i)).collect::<Vec<_>>();
            assert_eq!(&variances, &[0.1, 0.1, 0.1]);
        }
        {
            use Stationary;
            assert_eq!(process.var(), 0.1);
        }
    }
}