Commit b50ab27a authored by Egor Larionov's avatar Egor Larionov

Added hrbf fitting to data. Added kernel tests.

parent 0b52aef2
......@@ -6,4 +6,7 @@ authors = ["Egor Larionov <[email protected]>"]
[dependencies]
num-traits = "0.1"
nalgebra = "*"
lazy_static = "*"
[dev-dependencies]
approx = "0.1"
rand = "0.3"
This diff is collapsed.
......@@ -27,13 +27,11 @@
//! It must be that g(x), h(x,3) -> 0 as x -> 0 for the HRBF derivatives to exist.
use num_traits::{Float};
use std::fmt::Debug;
use std::marker::PhantomData;
use Real;
/// Kernel trait declaring all of the necessary derivatives.
pub trait Kernel<T> : Copy + Clone + Debug + Default
where T: Real
pub trait Kernel<T> : Copy + Clone + Default
where T: Float
{
/// Main kernel function φ(x)
fn f(&self, x: T) -> T;
......@@ -70,18 +68,18 @@ pub trait Kernel<T> : Copy + Clone + Debug + Default
}
/// Local kernel trait defines the radial fall-off for appropriate kernels.
pub trait LocalKernel<T> where T: Real {
pub trait LocalKernel<T> where T: Float {
fn new(r: T) -> Self;
}
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Pow2<T>(PhantomData<T>);
/// Default constructor for `Pow2` kernel
impl<T: Real> Default for Pow2<T> { fn default() -> Self { Pow2(PhantomData) } }
impl<T: Float> Default for Pow2<T> { fn default() -> Self { Pow2(PhantomData) } }
impl<T: Real> Kernel<T> for Pow2<T> {
impl<T: Float> Kernel<T> for Pow2<T> {
fn f(&self, x: T) -> T { x*x }
fn df(&self, x: T) -> T { T::from(2.0).unwrap()*x }
fn ddf(&self, _: T) -> T { T::from(2.0).unwrap() }
......@@ -94,13 +92,13 @@ impl<T: Real> Kernel<T> for Pow2<T> {
}
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Pow3<T>(::std::marker::PhantomData<T>);
/// Default constructor for `Pow3` kernel
impl<T: Real> Default for Pow3<T> { fn default() -> Self { Pow3(PhantomData) } }
impl<T: Float> Default for Pow3<T> { fn default() -> Self { Pow3(PhantomData) } }
impl<T: Real> Kernel<T> for Pow3<T> {
impl<T: Float> Kernel<T> for Pow3<T> {
fn f(&self, x: T) -> T { x*x*x }
fn df(&self, x: T) -> T { T::from(3.0).unwrap()*x*x }
fn ddf(&self, x: T) -> T { T::from(6.0).unwrap()*x }
......@@ -119,13 +117,13 @@ impl<T: Real> Kernel<T> for Pow3<T> {
}
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Pow4<T>(::std::marker::PhantomData<T>);
/// Default constructor for `Pow4` kernel
impl<T: Real> Default for Pow4<T> { fn default() -> Self { Pow4(PhantomData) } }
impl<T: Float> Default for Pow4<T> { fn default() -> Self { Pow4(PhantomData) } }
impl<T: Real> Kernel<T> for Pow4<T> {
impl<T: Float> Kernel<T> for Pow4<T> {
fn f(&self, x: T) -> T { x*x*x*x }
fn df(&self, x: T) -> T { T::from(4.0).unwrap()*x*x*x }
fn ddf(&self, x: T) -> T { T::from(12.0).unwrap()*x*x }
......@@ -139,13 +137,13 @@ impl<T: Real> Kernel<T> for Pow4<T> {
/// x^5 kernel.
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Pow5<T>(::std::marker::PhantomData<T>);
/// Default constructor for `Pow5` kernel
impl<T: Real> Default for Pow5<T> { fn default() -> Self { Pow5(PhantomData) } }
impl<T: Float> Default for Pow5<T> { fn default() -> Self { Pow5(PhantomData) } }
impl<T: Real> Kernel<T> for Pow5<T> {
impl<T: Float> Kernel<T> for Pow5<T> {
fn f(&self, x: T) -> T { x*x*x*x*x }
fn df(&self, x: T) -> T { T::from(5.0).unwrap()*x*x*x*x }
fn ddf(&self, x: T) -> T { T::from(20.0).unwrap()*x*x*x }
......@@ -159,20 +157,20 @@ impl<T: Real> Kernel<T> for Pow5<T> {
/// Gaussian kernel.
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Gauss<T> {
r: T,
}
impl<T: Real> LocalKernel<T> for Gauss<T> {
impl<T: Float> LocalKernel<T> for Gauss<T> {
fn new(r: T) -> Self { Gauss { r } }
}
/// Default constructor for `Gauss` kernel
impl<T: Real> Default for Gauss<T> { fn default() -> Self { Gauss { r: T::one() } } }
impl<T: Float> Default for Gauss<T> { fn default() -> Self { Gauss { r: T::one() } } }
impl<T: Real> Kernel<T> for Gauss<T> {
fn f(&self, x: T) -> T { Float::exp(-x * x / self.r) }
impl<T: Float> Kernel<T> for Gauss<T> {
fn f(&self, x: T) -> T { T::exp(-x * x / self.r) }
fn df(&self, x: T) -> T { -(T::from(2.0).unwrap()/self.r) * x * self.f(x) }
fn ddf(&self, x: T) -> T {
let _2 = T::from(2.0).unwrap();
......@@ -209,19 +207,19 @@ impl<T: Real> Kernel<T> for Gauss<T> {
/// Quintic kernel. Generates a positive definite hrbf fitting matrix.
/// Third and fourth order hrbf derivatives don't exist at x = 0.
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Csrbf31<T> {
r: T,
}
impl<T: Real> LocalKernel<T> for Csrbf31<T> {
impl<T: Float> LocalKernel<T> for Csrbf31<T> {
fn new(r: T) -> Self { Csrbf31 { r } }
}
/// Default constructor for `Csrbf31` kernel
impl<T: Real> Default for Csrbf31<T> { fn default() -> Self { Csrbf31 { r: T::one() } } }
impl<T: Float> Default for Csrbf31<T> { fn default() -> Self { Csrbf31 { r: T::one() } } }
impl<T: Real> Kernel<T> for Csrbf31<T> {
impl<T: Float> Kernel<T> for Csrbf31<T> {
fn f(&self, x: T) -> T {
let _1 = T::one();
let x = x/self.r;
......@@ -322,19 +320,19 @@ impl<T: Real> Kernel<T> for Csrbf31<T> {
}
}
#[derive(Copy, Clone, Debug)]
#[derive(Copy, Clone)]
pub struct Csrbf42<T> {
r: T,
}
impl<T: Real> LocalKernel<T> for Csrbf42<T> {
impl<T: Float> LocalKernel<T> for Csrbf42<T> {
fn new(r: T) -> Self { Csrbf42 { r } }
}
/// Default constructor for `Csrbf42` kernel
impl<T: Real> Default for Csrbf42<T> { fn default() -> Self { Csrbf42 { r: T::one() } } }
impl<T: Float> Default for Csrbf42<T> { fn default() -> Self { Csrbf42 { r: T::one() } } }
impl<T: Real> Kernel<T> for Csrbf42<T> {
impl<T: Float> Kernel<T> for Csrbf42<T> {
fn f(&self, x: T) -> T {
let _1 = T::one();
let x = x/self.r;
......
This diff is collapsed.
This diff is collapsed.
extern crate num_traits;
extern crate rand;
mod autodiff;
use num_traits::Float;
use autodiff::{cst, diff, grad};
// NOTE: we don't need approximate equality here because we compare with the exact derivative
// expression, which means that the derivative and expected values should be identically computed.
#[test]
fn simple_test() {
assert_eq!(diff(|_| cst(1.0), 0.0), 0.0);
assert_eq!(diff(|x| x, 0.0), 1.0);
assert_eq!(diff(|x| x*x, 0.0), 0.0);
assert_eq!(diff(|x| x*x, 1.0), 2.0);
assert_eq!(diff(|x| Float::exp(-x*x/cst(2.0)), 0.0), 0.0);
}
#[test]
fn random_test() {
use self::rand::{SeedableRng, StdRng};
use self::rand::distributions::{IndependentSample, Range};
let seed: &[_] = &[1,2,3,4];
let mut rng: StdRng = SeedableRng::from_seed(seed);
let range = Range::new(-1.0, 1.0);
for _ in 0..99 {
let t = range.ind_sample(&mut rng);
assert_eq!(diff(|x| Float::exp(-x*x/cst(2.0)), t), -t*Float::exp(-t*t/2.0));
}
}
#[test]
fn grad_test() {
use self::rand::{SeedableRng, StdRng};
use self::rand::distributions::{IndependentSample, Range};
let seed: &[_] = &[1,2,3,4];
let mut rng: StdRng = SeedableRng::from_seed(seed);
let range = Range::new(-1.0, 1.0);
for _ in 0..99 {
let t = vec![range.ind_sample(&mut rng), range.ind_sample(&mut rng)];
let expected = vec![-0.5*t[1]*Float::exp(-t[0]*t[1]/2.0), -0.5*t[0]*Float::exp(-t[0]*t[1]/2.0)];
assert_eq!(grad(|x| Float::exp(-x[0]*x[1]/cst(2.0)), t), expected);
}
}
extern crate num_traits;
extern crate hrbf;
extern crate rand;
#[macro_use]
extern crate approx;
#[allow(dead_code)]
mod autodiff;
use hrbf::kernel::*;
use autodiff::{Num, cst};
const TEST_VALS: [f64;4] = [0.0, 1.0, 0.5, ::std::f64::consts::PI];
const TEST_RADIUS: f64 = 2.0;
fn test_kernel<F,K: Kernel<Num>>(ker: K, x0: f64, compare: F)
where F: Fn(f64, f64)
{
let x = Num { val: x0, eps: 1.0 };
let f = ker.f(x);
let df = ker.df(x);
let ddf = ker.ddf(x);
let dddf = ker.dddf(x);
let ddddf = ker.ddddf(x);
compare(f.eps, df.val);
compare(df.eps, ddf.val);
compare(ddf.eps, dddf.val);
compare(dddf.eps, ddddf.val);
if x0 != 0.0 {
let df_l = ker.df_l(x);
let g = ker.g(x);
let g_l = ker.g_l(x);
let h3 = ker.h(x, cst(3.0));
let h52 = ker.h(x, cst(5.0/2.0));
compare(x0*df_l.val, df.val);
compare(x0*x0*g.val, ddf.val*x0 - df.val);
compare(x0*g_l.val, g.val);
compare(x0*x0*x0*h3.val, x0*x0*dddf.val - 3.0*(x0*ddf.val - df.val));
compare(x0*x0*x0*h52.val, x0*x0*dddf.val - 0.5*5.0*(x0*ddf.val - df.val));
}
}
fn test_kernel_simple<K: Kernel<Num>>(kern: K) {
for &x in TEST_VALS.iter() { test_kernel(kern, x, easy_compare); }
}
fn test_kernel_random<K: Kernel<Num>>(kern: K) {
use self::rand::{SeedableRng, StdRng};
use self::rand::distributions::{IndependentSample, Range};
let seed: &[_] = &[1,2,3,4];
let mut rng: StdRng = SeedableRng::from_seed(seed);
let range = Range::new(-1.0, 1.0);
for _ in 0..999 {
let x = range.ind_sample(&mut rng);
test_kernel(kern, x, hard_compare);
}
}
fn easy_compare(a: f64, b: f64) {
assert_ulps_eq!(a,b, max_ulps=6);
}
fn hard_compare(a: f64, b: f64) {
assert_relative_eq!(a,b, max_relative=1e-13, epsilon=1e-14);
}
#[test]
fn pow2_test() {
let kern = Pow2::<Num>::default();
test_kernel_simple(kern);
test_kernel_random(kern);
}
#[test]
fn pow3_test() {
let kern = Pow3::<Num>::default();
test_kernel_simple(kern);
test_kernel_random(kern);
}
#[test]
fn pow4_test() {
let kern = Pow4::<Num>::default();
test_kernel_simple(kern);
test_kernel_random(kern);
}
#[test]
fn pow5_test() {
let kern = Pow5::<Num>::default();
test_kernel_simple(kern);
test_kernel_random(kern);
}
#[test]
fn gauss_test() {
let kern = Gauss::<Num>::new(cst(TEST_RADIUS));
test_kernel_simple(kern);
test_kernel_random(kern);
}
#[test]
fn csrbf31_test() {
let kern = Csrbf31::<Num>::new(cst(TEST_RADIUS));
test_kernel_simple(kern);
test_kernel_random(kern);
}
#[test]
fn csrbf42_test() {
let kern = Csrbf42::<Num>::new(cst(TEST_RADIUS));
test_kernel_simple(kern);
test_kernel_random(kern);
}
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