Commit 34a912a3 authored by Rui Vieira's avatar Rui Vieira 🍵
Browse files

Add monotonic cubic spline interpolation

parent b0117336
/target
**/*.rs.bk
.idea
* 0.0.3
* Add Monotonic Cubic Spline interpolation
* 0.0.2
* Add `pdf` and `logpdf` to the Normal distribution
\ No newline at end of file
This diff is collapsed.
[package]
name = "mentat"
version = "0.0.2"
version = "0.0.3"
authors = ["Rui Vieira <ruidevieira@googlemail.com>"]
exclude = [".idea/**/*"]
exclude = [".idea/**/*", "docs/**/*"]
description = "A Rust library for scientific computing."
homepage = "https://gitlab.com/ruivieira/mentat"
repository = "https://gitlab.com/ruivieira/mentat"
readme = "README.md"
keywords = ["science", "statistics"]
keywords = ["science", "statistics", "math"]
categories = ["science"]
license = "Apache-2.0"
[dependencies]
rand = "0.6.4"
\ No newline at end of file
rand = "0.6.4"
[dev-dependencies]
matplotrust = "0.1.5"
\ No newline at end of file
# mentat
A Rust library for scientific computing
\ No newline at end of file
A Rust library for scientific computing
# features
* Monotonic Cubic Spline interpolation
# examples
Full examples are available [here](docs/README.md).
\ No newline at end of file
# Examples
# Monotonic Cubic Spline interpolation
```rust
let x = vec![0.0, 2.0, 3.0, 10.0];
let y = vec![1.0, 4.0, 8.0, 10.5];
let mut spline = MonotonicCubicSpline::new(&x, &y);
let mut x_interp = Vec::new();
let mut y_interp = Vec::new();
for i in 0..100 {
let p = i as f64 / 10.0;
x_interp.push(p);
y_interp.push(spline.interpolate(p));
}
```
![plot](figures/monotonic_cubic_spline.png)
\ No newline at end of file
extern crate rand;
#[cfg(test)]
#[macro_use]
extern crate matplotrust;
use rand::prelude::*;
use rand::distributions::Normal as N;
use std::f64::consts::PI;
......@@ -42,6 +47,93 @@ impl Normal {
}
}
pub struct MonotonicCubicSpline {
m_x: Vec<f64>,
m_y: Vec<f64>,
m_m: Vec<f64>
}
impl MonotonicCubicSpline {
pub fn new(x : &Vec<f64>, y : &Vec<f64>) -> MonotonicCubicSpline {
assert!(x.len() == y.len() && x.len() >= 2 && y.len() >= 2, "Must have at least 2 control points.");
let n = x.len();
let mut secants = vec![0.0 ; n - 1];
let mut slopes = vec![0.0 ; n];
for i in 0..(n-1) {
let h = *x.get(i + 1).unwrap() - *x.get(i).unwrap();
assert!(h > 0.0, "Control points must be monotonically increasing.");
secants[i] = (*y.get(i + 1).unwrap() - *y.get(i).unwrap()) / h;
}
slopes[0] = secants[0];
for i in 1..(n-1) {
slopes[i] = (secants[i - 1] + secants[i]) * 0.5;
}
slopes[n - 1] = secants[n - 2];
for i in 0..(n-1) {
if secants[i] == 0.0 {
slopes[i] = 0.0;
slopes[i + 1] = 0.0;
} else {
let a = slopes[i] / secants[i];
let b = slopes[i + 1] / secants[i];
let h = a.hypot(b);
if h > 9.0 {
let t = 3.0 / h;
slopes[i] = t * a * secants[i];
slopes[i + 1] = t * b * secants[i];
}
}
}
let spline = MonotonicCubicSpline {
m_x: x.clone(),
m_y: y.clone(),
m_m: slopes
};
return spline;
}
pub fn hermite(point: f64, x : (f64, f64), y: (f64, f64), m: (f64, f64)) -> f64 {
let h: f64 = x.1 - x.0;
let t = (point - x.0) / h;
return (y.0 * (1.0 + 2.0 * t) + h * m.0 * t) * (1.0 - t) * (1.0 - t)
+ (y.1 * (3.0 - 2.0 * t) + h * m.1 * (t - 1.0)) * t * t;
}
pub fn interpolate(&mut self, point : f64) -> f64 {
let n = self.m_x.len();
if point <= *self.m_x.get(0).unwrap() {
return *self.m_y.get(0).unwrap();
}
if point >= *self.m_x.get(n - 1).unwrap() {
return *self.m_y.get(n - 1).unwrap();
}
let mut i = 0;
while point >= *self.m_x.get(i + 1).unwrap() {
i += 1;
if point == *self.m_x.get(i).unwrap() {
return *self.m_y.get(i).unwrap();
}
}
return MonotonicCubicSpline::hermite(point,
(*self.m_x.get(i).unwrap(), *self.m_x.get(i+1).unwrap()),
(*self.m_y.get(i).unwrap(), *self.m_y.get(i+1).unwrap()),
(*self.m_m.get(i).unwrap(), *self.m_m.get(i+1).unwrap()));
}
}
#[cfg(test)]
mod test_normal {
......@@ -71,4 +163,29 @@ mod test_normal {
assert_delta!(0.3989423, _pdf, 1e-5);
}
#[test]
fn interpolation() {
use matplotrust::*;
let x = vec![0.0, 2.0, 3.0, 10.0];
let y = vec![1.0, 4.0, 8.0, 10.5];
let mut spline = MonotonicCubicSpline::new(&x, &y);
let mut x_interp = Vec::new();
let mut y_interp = Vec::new();
for i in 0..100 {
let p = i as f64 / 10.0;
x_interp.push(p);
y_interp.push(spline.interpolate(p));
}
let mut figure = Figure::new();
let points = scatter_plot::<f64, f64>(x, y, None);
let interpolation = line_plot::<f64, f64>(x_interp, y_interp, None);
figure.add_plot(points);
figure.add_plot(interpolation);
figure.save("./docs/figures/monotonic_cubic_spline.png", None);
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment