Skip to content

Commit

Permalink
Merge pull request #144 from light-curve/rand-0.8
Browse files Browse the repository at this point in the history
Update rand to 0.8
  • Loading branch information
hombit authored Oct 10, 2023
2 parents 5b535fa + a1c348b commit 1590b18
Show file tree
Hide file tree
Showing 8 changed files with 88 additions and 56 deletions.
4 changes: 2 additions & 2 deletions Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ hyperdual = "1.1"
light-curve-common = "0.1.0"
plotters = { version = "0.3.5", default-features = false, features = ["errorbar", "line_series", "ttf"] }
plotters-bitmap = "0.3.3"
rand = "0.7"
rand_distr = "0.2"
rand = "0.8"
rand_distr = "0.4"
rayon = "1.5"
realfft= "3.1"
rustfft = "6.1"
Expand Down
2 changes: 1 addition & 1 deletion benches/fit.rs
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ pub fn bench_fit_straight_line(c: &mut Criterion) {
let y: Vec<_> = x.iter().map(|&x| x + thread_rng().gen::<f64>()).collect();
let w: Vec<_> = x
.iter()
.map(|_| thread_rng().gen_range(10.0, 100.0))
.map(|_| thread_rng().gen_range(10.0..100.0))
.collect();
let ts = TimeSeries::new(&x, &y, &w);

Expand Down
15 changes: 8 additions & 7 deletions src/features/bazin_fit.rs
Original file line number Diff line number Diff line change
Expand Up @@ -445,20 +445,21 @@ mod tests {
})
.collect();
let w: Vec<_> = model.iter().copied().map(f64::recip).collect();
println!("{:?}\n{:?}\n{:?}\n{:?}", t, model, m, w);
println!("t = {:?}\nmodel = {:?}\nm = {:?}\nw = {:?}", t, model, m, w);
let mut ts = TimeSeries::new(&t, &m, &w);

// curve_fit(lambda t, a, b, t0, rise, fall: b + a * np.exp(-(t-t0)/fall) / (1 + np.exp(-(t-t0) / rise)), xdata=t, ydata=m, sigma=np.array(w)**-0.5, p0=[1e4, 1e3, 30, 10, 30])
let desired = [
9.89658673e+03,
1.11312724e+03,
3.06401284e+01,
9.75027284e+00,
2.86714363e+01,
1.01014622e+04,
9.76594899e+02,
3.02450768e+01,
1.00640456e+01,
2.99357915e+01,
];

let values = eval.eval(&mut ts).unwrap();
assert_relative_eq!(&values[..5], &desired[..], max_relative = 0.01);
assert_relative_eq!(&values[..5], &param_true[..], max_relative = 0.03);
assert_relative_eq!(&values[..5], &desired[..], max_relative = 0.02);
}

#[cfg(any(feature = "ceres-source", feature = "ceres-system"))]
Expand Down
7 changes: 4 additions & 3 deletions src/features/linexp_fit.rs
Original file line number Diff line number Diff line change
Expand Up @@ -423,14 +423,15 @@ mod tests {
})
.collect();
let w: Vec<_> = model.iter().copied().map(f64::recip).collect();
println!("{:?}\n{:?}\n{:?}\n{:?}", t, model, m, w);
println!("t = {:?}\nmodel = {:?}\nm = {:?}\nw = {:?}", t, model, m, w);
let mut ts = TimeSeries::new(&t, &m, &w);

// curve_fit(lambda t, a, t0, fall, b : b + a * ((t - t0) / fall) * np.exp(-(t - t0) / fall), xdata=t, ydata=m, sigma=np.array(w)**-0.5, p0=[800, 10, 15, 30])
let desired = [982.42262317, -15.08873767, 20.42325543, 13.23259373];
let desired = [986.62990444, -15.1956711, 20.05763093, 15.54839175];

let values = eval.eval(&mut ts).unwrap();
assert_relative_eq!(&values[..4], &desired[..], max_relative = 0.05);
assert_relative_eq!(&values[..NPARAMS], &param_true[..], max_relative = 0.07);
assert_relative_eq!(&values[..NPARAMS], &desired[..], max_relative = 0.04);
}

#[cfg(any(feature = "ceres-source", feature = "ceres-system"))]
Expand Down
41 changes: 21 additions & 20 deletions src/features/villar_fit.rs
Original file line number Diff line number Diff line change
Expand Up @@ -644,13 +644,13 @@ mod tests {

// curve_fit(lambda t, a, c, t0, tau_rise, tau_fall, nu, gamma: c + a * (1 - nu * np.where(t - t0 < gamma, (t - t0) / gamma, 1)) / (1 + np.exp(-(t-t0) / tau_rise)) * np.where(t > t0 + gamma, np.exp(-(t-t0-gamma) / tau_fall), 1.0), xdata=t, ydata=m, sigma=np.array(w)**-0.5, p0=[1e4, 1e3, 20, 5, 30, 0.3, 30])
let desired = [
1.00899754e+04,
1.00689083e+03,
2.02385802e+01,
5.04283765e+00,
3.00679883e+01,
3.00400783e-01,
2.94149214e+01,
9.94199441e+03,
1.00094110e+03,
1.98928845e+01,
4.96049541e+00,
2.96149274e+01,
2.97679866e-01,
3.05544920e+01,
];

let values = eval.eval(&mut ts).unwrap();
Expand All @@ -660,25 +660,26 @@ mod tests {
t, model, m, w, values[NPARAMS]
);

assert_relative_eq!(&values[..NPARAMS], &desired[..], max_relative = 0.01);
assert_relative_eq!(&values[..NPARAMS], &param_true[..], max_relative = 0.1);
assert_relative_eq!(&values[..NPARAMS], &desired[..], max_relative = 0.05);
}

// It doesn't converge to the right place
// #[cfg(any(feature = "ceres-source", feature = "ceres-system"))]
// #[test]
// fn villar_fit_noisy_ceres() {
// villar_fit_noisy(VillarFit::new(
// CeresCurveFit::new().into(),
// LnPrior::none(),
// VillarInitsBounds::Default,
// ));
// }
#[cfg(any(feature = "ceres-source", feature = "ceres-system"))]
#[test]
fn villar_fit_noisy_ceres() {
let ceres = CeresCurveFit::default();
villar_fit_noisy(VillarFit::new(
ceres.into(),
LnPrior::none(),
VillarInitsBounds::Default,
));
}

#[cfg(any(feature = "ceres-source", feature = "ceres-system"))]
#[test]
fn villar_fit_noizy_mcmc_plus_ceres() {
let ceres = CeresCurveFit::default();
let mcmc = McmcCurveFit::new(512, Some(ceres.into()));
let mcmc = McmcCurveFit::new(1024, Some(ceres.into()));
villar_fit_noisy(VillarFit::new(
mcmc.into(),
LnPrior::none(),
Expand Down Expand Up @@ -721,7 +722,7 @@ mod tests {
LnPrior1D::log_normal(f64::ln(30.0), 1.0),
]);
let lmsder = LmsderCurveFit::new(1);
let mcmc = McmcCurveFit::new(1024, Some(lmsder.into()));
let mcmc = McmcCurveFit::new(2048, Some(lmsder.into()));
villar_fit_noisy(VillarFit::new(
mcmc.into(),
prior,
Expand Down
20 changes: 15 additions & 5 deletions src/nl_fit/ceres.rs
Original file line number Diff line number Diff line change
Expand Up @@ -168,9 +168,6 @@ mod tests {
const N: usize = 300;
const NOISE: f64 = 0.5;

// curve_fit(lambda x, a, b, c: b * np.exp(-a * x) * x**2 + c, xdata=t, ydata=y, p0=[1, 1, 1], xtol=1e-6)
let desired = [0.7450598836400693, 1.981911479079224, 0.5094446163866907];

let param_true = [0.75, 2.0, 0.5];
let param_init = [1.0, 1.0, 1.0];

Expand All @@ -182,6 +179,12 @@ mod tests {
nonlinear_func(x, &param_true) + NOISE * eps
});
let inv_err: Array1<_> = vec![1.0 / NOISE; N].into();
println!(
"t = {:?}\ny = {:?}\ninv_err = {:?}",
t.as_slice().unwrap(),
y.as_slice().unwrap(),
inv_err.as_slice().unwrap()
);
let ts = Rc::new(Data { t, m: y, inv_err });

let fitter = CeresCurveFit::new(14, None);
Expand All @@ -194,8 +197,15 @@ mod tests {
nonlinear_func_dump_ln_prior,
);

// curve_fit(lambda x, a, b, c: b * np.exp(-a * x) * x**2 + c, xdata=t, ydata=y, sigma=1/np.array(inv_err), p0=[1, 1, 1], xtol=1e-6)
let desired = [0.76007721, 2.0225076, 0.49238112];

// Not as good as for LMSDER
assert_abs_diff_eq!(&result.x[..], &param_true[..], epsilon = 0.04);
assert_abs_diff_eq!(&result.x[..], &desired[..], epsilon = 0.02);
assert_abs_diff_eq!(
&result.x[..],
&param_true[..],
epsilon = NOISE / (N as f64).sqrt()
);
assert_abs_diff_eq!(&result.x[..], &desired[..], epsilon = 0.04);
}
}
53 changes: 36 additions & 17 deletions src/nl_fit/lmsder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,8 @@ mod tests {
use super::*;
use crate::straight_line_fit::StraightLineFitterResult;

use light_curve_common::{all_close, linspace};
use approx::assert_relative_eq;
use light_curve_common::linspace;
use rand::prelude::*;
use rand_distr::StandardNormal;
use std::ops::Mul;
Expand Down Expand Up @@ -422,7 +423,11 @@ mod tests {
assert_eq!(result.status, Value::Success);
let param = result.x();

all_close(param.as_slice().unwrap(), &[-247.0, 39.0], 1e-9);
assert_relative_eq!(
param.as_slice().unwrap(),
[-247.0, 39.0].as_slice(),
max_relative = 1e-9
);
}

#[test]
Expand Down Expand Up @@ -466,7 +471,11 @@ mod tests {
assert_eq!(result.status, Value::Success);
let param = result.x();

all_close(param.as_slice().unwrap(), &[0.0, 1.0, 1.0], 1e-8);
assert_relative_eq!(
param.as_slice().unwrap(),
[0.0, 1.0, 1.0].as_slice(),
epsilon = 1e-8
);
}

#[inline]
Expand All @@ -485,9 +494,6 @@ mod tests {
const NOISE: f64 = 0.5;
const RTOL: f64 = 1e-6;

// curve_fit(lambda x, a, b, c: b * np.exp(-a * x) * x**2 + c, xdata=t, ydata=y, p0=[1, 1, 1], xtol=1e-6)
let desired = [0.7450598836400693, 1.981911479079224, 0.5094446163866907];

let param_true = [0.75, 2.0, 0.5];
let param_init = [1.0, 1.0, 1.0];
let param_init = VectorF64::from_slice(&param_init).unwrap();
Expand All @@ -502,6 +508,7 @@ mod tests {
nonlinear_func(&param_true, x) + NOISE * eps
})
.collect();
println!("t = {:?}\ny = {:?}", t, y);
let data = Rc::new(Data { t, y, err: None });

let function = {
Expand Down Expand Up @@ -546,12 +553,19 @@ mod tests {
assert_eq!(result.status, Value::Success);
let param = result.x();

all_close(
// curve_fit(lambda x, a, b, c: b * np.exp(-a * x) * x**2 + c, xdata=t, ydata=y, p0=[1, 1, 1], xtol=1e-6)
let desired = [0.76007721, 2.0225076, 0.49238112];

assert_relative_eq!(
param.as_slice().unwrap(),
&param_true,
NOISE / (N as f64).sqrt(),
param_true.as_slice(),
max_relative = NOISE / (N as f64).sqrt(),
);
assert_relative_eq!(
param.as_slice().unwrap(),
desired.as_slice(),
max_relative = RTOL
);
all_close(param.as_slice().unwrap(), &desired, RTOL);
}

#[test]
Expand All @@ -560,9 +574,6 @@ mod tests {
const NOISE: f64 = 0.5;
const RTOL: f64 = 1e-6;

// curve_fit(lambda x, a, b, c: b * np.exp(-a * x) * x**2 + c, xdata=t, ydata=y, p0=[1, 1, 1], xtol=1e-6)
let desired = [0.7450598836400693, 1.981911479079224, 0.5094446163866907];

let param_true = [0.75, 2.0, 0.5];
let param_init = [1.0, 1.0, 1.0];
let param_init = VectorF64::from_slice(&param_init).unwrap();
Expand All @@ -577,6 +588,7 @@ mod tests {
nonlinear_func(&param_true, x) + NOISE * eps
})
.collect();
println!("t = {:?}\ny = {:?}", t, y);
let data = Rc::new(Data { t, y, err: None });
let data_real = data.clone();
let data_dual = data;
Expand Down Expand Up @@ -608,11 +620,18 @@ mod tests {
assert_eq!(result.status, Value::Success);
let param = result.x();

all_close(
// curve_fit(lambda x, a, b, c: b * np.exp(-a * x) * x**2 + c, xdata=t, ydata=y, p0=[1, 1, 1], xtol=1e-6)
let desired = [0.76007721, 2.0225076, 0.49238112];

assert_relative_eq!(
param.as_slice().unwrap(),
param_true.as_slice(),
max_relative = NOISE / ((N - param_true.len()) as f64).sqrt(),
);
assert_relative_eq!(
param.as_slice().unwrap(),
&param_true,
NOISE / (N as f64).sqrt(),
desired.as_slice(),
max_relative = RTOL
);
all_close(param.as_slice().unwrap(), &desired, RTOL);
}
}
2 changes: 1 addition & 1 deletion src/periodogram/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,7 @@ mod tests {

#[test]
fn direct_vs_fft_unevenly_sin_cos() {
const OMEGA1: f64 = 0.472;
const OMEGA1: f64 = 0.222;
const OMEGA2: f64 = 1.222;
const AMPLITUDE2: f64 = 2.0;
const NOISE_AMPLITUDE: f64 = 1.0;
Expand Down

0 comments on commit 1590b18

Please sign in to comment.