From a1c348b984dfaf454692c34879c91884a305b33d Mon Sep 17 00:00:00 2001 From: Konstantin Malanchev Date: Tue, 10 Oct 2023 08:19:31 -0400 Subject: [PATCH] Update rand to 0.8 --- Cargo.toml | 4 +-- benches/fit.rs | 2 +- src/features/bazin_fit.rs | 15 ++++++----- src/features/linexp_fit.rs | 7 ++--- src/features/villar_fit.rs | 41 +++++++++++++++-------------- src/nl_fit/ceres.rs | 20 ++++++++++---- src/nl_fit/lmsder.rs | 53 ++++++++++++++++++++++++++------------ src/periodogram/mod.rs | 2 +- 8 files changed, 88 insertions(+), 56 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 169f0182..5fd4dca6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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" diff --git a/benches/fit.rs b/benches/fit.rs index 374209a7..df36a693 100644 --- a/benches/fit.rs +++ b/benches/fit.rs @@ -17,7 +17,7 @@ pub fn bench_fit_straight_line(c: &mut Criterion) { let y: Vec<_> = x.iter().map(|&x| x + thread_rng().gen::()).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); diff --git a/src/features/bazin_fit.rs b/src/features/bazin_fit.rs index dc0a9caf..a4c77df6 100644 --- a/src/features/bazin_fit.rs +++ b/src/features/bazin_fit.rs @@ -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], ¶m_true[..], max_relative = 0.03); + assert_relative_eq!(&values[..5], &desired[..], max_relative = 0.02); } #[cfg(any(feature = "ceres-source", feature = "ceres-system"))] diff --git a/src/features/linexp_fit.rs b/src/features/linexp_fit.rs index 17a988bf..eea9198c 100644 --- a/src/features/linexp_fit.rs +++ b/src/features/linexp_fit.rs @@ -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], ¶m_true[..], max_relative = 0.07); + assert_relative_eq!(&values[..NPARAMS], &desired[..], max_relative = 0.04); } #[cfg(any(feature = "ceres-source", feature = "ceres-system"))] diff --git a/src/features/villar_fit.rs b/src/features/villar_fit.rs index 1c1286fa..c3d030e9 100644 --- a/src/features/villar_fit.rs +++ b/src/features/villar_fit.rs @@ -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(); @@ -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], ¶m_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(), @@ -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, diff --git a/src/nl_fit/ceres.rs b/src/nl_fit/ceres.rs index 3acb6929..2b70c0e3 100644 --- a/src/nl_fit/ceres.rs +++ b/src/nl_fit/ceres.rs @@ -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]; @@ -182,6 +179,12 @@ mod tests { nonlinear_func(x, ¶m_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); @@ -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[..], ¶m_true[..], epsilon = 0.04); - assert_abs_diff_eq!(&result.x[..], &desired[..], epsilon = 0.02); + assert_abs_diff_eq!( + &result.x[..], + ¶m_true[..], + epsilon = NOISE / (N as f64).sqrt() + ); + assert_abs_diff_eq!(&result.x[..], &desired[..], epsilon = 0.04); } } diff --git a/src/nl_fit/lmsder.rs b/src/nl_fit/lmsder.rs index accab0bc..af6426a6 100644 --- a/src/nl_fit/lmsder.rs +++ b/src/nl_fit/lmsder.rs @@ -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; @@ -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] @@ -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] @@ -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(¶m_init).unwrap(); @@ -502,6 +508,7 @@ mod tests { nonlinear_func(¶m_true, x) + NOISE * eps }) .collect(); + println!("t = {:?}\ny = {:?}", t, y); let data = Rc::new(Data { t, y, err: None }); let function = { @@ -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(), - ¶m_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] @@ -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(¶m_init).unwrap(); @@ -577,6 +588,7 @@ mod tests { nonlinear_func(¶m_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; @@ -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(), - ¶m_true, - NOISE / (N as f64).sqrt(), + desired.as_slice(), + max_relative = RTOL ); - all_close(param.as_slice().unwrap(), &desired, RTOL); } } diff --git a/src/periodogram/mod.rs b/src/periodogram/mod.rs index af7a29df..7ce3c5cf 100644 --- a/src/periodogram/mod.rs +++ b/src/periodogram/mod.rs @@ -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;