Skip to content

Commit

Permalink
Merge pull request #98 from light-curve/ceres
Browse files Browse the repository at this point in the history
Ceres as a new gradient optimizer for *Fit features
hombit authored Feb 17, 2023

Partially verified

This commit is signed with the committer’s verified signature.
gsmet’s contribution has been verified via GPG key.
We cannot verify signatures from co-authors, and some of the co-authors attributed to this commit require their commits to be signed.
2 parents 327f9ba + 22e3b1b commit e4028c7
Showing 14 changed files with 383 additions and 25 deletions.
17 changes: 10 additions & 7 deletions .github/workflows/test.yml
Original file line number Diff line number Diff line change
@@ -18,6 +18,8 @@ jobs:
- run: cargo build --all-targets --features gsl,fftw-source --no-default-features
- run: cargo build --all-targets --features gsl,fftw-system --no-default-features
- run: cargo build --all-targets --features gsl,fftw-mkl --no-default-features
- run: cargo build --all-targets --features gsl,fftw-source,ceres-source --no-default-features
- run: cargo build --all-targets --features fftw-source,ceres-source --no-default-features

msrv-build:
runs-on: ubuntu-latest
@@ -37,7 +39,7 @@ jobs:
with:
toolchain: ${{ steps.get_msrv.outputs.msrv }}
# Build "normal" target only, see https://github.com/light-curve/light-curve-feature/issues/74
- run: cargo +${{ steps.get_msrv.outputs.msrv }} build --no-default-features --features gsl,fftw-source
- run: cargo +${{ steps.get_msrv.outputs.msrv }} build --no-default-features --features ceres-source,fftw-source,gsl

test:
runs-on: ubuntu-latest
@@ -50,7 +52,7 @@ jobs:
run: |
sudo apt-get update
sudo apt-get install -y libfftw3-dev libgsl-dev
- run: cargo test --no-default-features --features=fftw-source,gsl
- run: cargo test --no-default-features --features=ceres-source,fftw-source,gsl

test-release:
runs-on: ubuntu-latest
@@ -63,7 +65,7 @@ jobs:
run: |
sudo apt-get update
sudo apt-get install -y libfftw3-dev libgsl-dev
- run: cargo test --profile=release-with-debug --no-default-features --features=fftw-source,gsl
- run: cargo test --profile=release-with-debug --no-default-features --features=ceres-source,fftw-source,gsl

examples:
runs-on: ubuntu-latest
@@ -76,7 +78,7 @@ jobs:
run: |
sudo apt-get update
sudo apt-get install -y libfftw3-dev libgsl-dev
- run: cargo run --example plot_snia_curve_fits --no-default-features --features=fftw-source,gsl -- -n=1
- run: cargo run --example plot_snia_curve_fits --no-default-features --features=ceres-source,fftw-source,gsl -- -n=1

fmt:
runs-on: ubuntu-latest
@@ -98,7 +100,7 @@ jobs:
run: |
sudo apt-get update
sudo apt-get install -y libfftw3-dev libgsl-dev
- run: cargo clippy --all-targets --no-default-features --features=fftw-source,gsl -- -D warnings
- run: cargo clippy --all-targets --no-default-features --features=ceres-source,fftw-source,gsl -- -D warnings

fmt-test-util:
runs-on: ubuntu-latest
@@ -134,6 +136,7 @@ jobs:
with:
submodules: true
- name: build
# ceres-source doesn't build on ARM yet
run: cargo build --all-targets --no-default-features --features=fftw-source,gsl

macos:
@@ -145,7 +148,7 @@ jobs:
submodules: true
- name: Install gsl
run: brew install gsl
- run: cargo build --all-targets --no-default-features --features=fftw-source,gsl
- run: cargo build --all-targets --no-default-features --features=ceres-source,fftw-source,gsl

windows:
runs-on: windows-latest
@@ -154,4 +157,4 @@ jobs:
- uses: actions/checkout@v3
with:
submodules: true
- run: cargo build --all-targets --no-default-features --features=fftw-source
- run: cargo build --all-targets --no-default-features --features=ceres-source,fftw-source
2 changes: 1 addition & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

### Added

--
- `CeresCurveFit` and optional `ceres-source` Cargo feature to use Ceres Solver for curve fitting

### Changed

4 changes: 3 additions & 1 deletion Cargo.toml
Original file line number Diff line number Diff line change
@@ -16,6 +16,7 @@ bench = false

[features]
default = ["fftw-source"]
ceres-source = ["ceres-solver/source"]
fftw-system = ["fftw/system"]
fftw-source = ["fftw/source"]
fftw-mkl = ["fftw/intel-mkl"]
@@ -33,6 +34,7 @@ debug = true
# it is a dependency of intel-mkl-tool, we pin it to temporary solve
# https://github.com/rust-math/intel-mkl-src/issues/68
anyhow = "<1.0.49"
ceres-solver = { version = "0.2.0", optional = true }
conv = "^0.3.3"
emcee = "^0.3.0"
emcee_rand = { version = "^0.3.15", package = "rand" }
@@ -83,4 +85,4 @@ rustdoc-args = [
"katex-header.html",
]
no-default-features = true
features = ["fftw-system", "gsl"]
features = ["ceres-source", "fftw-system", "gsl"]
42 changes: 39 additions & 3 deletions examples/plot_snia_curve_fits.rs
Original file line number Diff line number Diff line change
@@ -3,8 +3,8 @@ use light_curve_feature::{
features::VillarLnPrior, prelude::*, BazinFit, Feature, FeatureEvaluator, McmcCurveFit,
TimeSeries, VillarFit,
};
#[cfg(feature = "gsl")]
use light_curve_feature::{LmsderCurveFit, LnPrior};
#[cfg(all(feature = "ceres-source", feature = "gsl"))]
use light_curve_feature::{CeresCurveFit, LmsderCurveFit, LnPrior};
use light_curve_feature_test_util::iter_sn1a_flux_ts;
use ndarray::{Array1, ArrayView1};
use plotters::prelude::*;
@@ -20,7 +20,7 @@ fn main() {
#[allow(clippy::vec_init_then_push)]
let features = {
let mut features: Vec<(&str, Feature<_>)> = vec![];
#[cfg(feature = "gsl")]
#[cfg(all(feature = "gsl", feature = "ceres-source"))]
{
features.push((
"BazinFit LMSDER",
@@ -40,6 +40,24 @@ fn main() {
)
.into(),
));
features.push((
"BazinFit Ceres",
BazinFit::new(
CeresCurveFit::default().into(),
LnPrior::none(),
BazinFit::default_inits_bounds(),
)
.into(),
));
features.push((
"BazinFit MCMC+Ceres",
BazinFit::new(
McmcCurveFit::new(1024, Some(CeresCurveFit::default().into())).into(),
LnPrior::none(),
BazinFit::default_inits_bounds(),
)
.into(),
));
features.push((
"VillarFit LMSDER",
VillarFit::new(
@@ -58,6 +76,24 @@ fn main() {
)
.into(),
));
features.push((
"VillarFit Ceres",
VillarFit::new(
CeresCurveFit::default().into(),
LnPrior::none(),
VillarFit::default_inits_bounds(),
)
.into(),
));
features.push((
"VillarFit MCMC+Ceres",
VillarFit::new(
McmcCurveFit::new(1024, Some(CeresCurveFit::default().into())).into(),
LnPrior::none(),
VillarFit::default_inits_bounds(),
)
.into(),
));
}
features.push((
"VillarFit MCMC+prior",
62 changes: 60 additions & 2 deletions src/features/bazin_fit.rs
Original file line number Diff line number Diff line change
@@ -42,9 +42,10 @@ impl BazinFit {
/// New [BazinFit] instance
///
/// `algorithm` specifies which optimization method is used, it is an instance of the
/// [CurveFitAlgorithm], currently supported algorithms are [MCMC](McmcCurveFit) and
/// [CurveFitAlgorithm], currently supported algorithms are [MCMC](McmcCurveFit),
/// [LMSDER](crate::nl_fit::LmsderCurveFit) (a Levenberg–Marquard algorithm modification,
/// requires `gsl` Cargo feature).
/// requires `gsl` Cargo feature), and [Ceres](crate::nl_fit::CeresCurveFit) (trust-region
/// algorithm, requires `ceres` Cargo feature).
///
/// `ln_prior` is an instance of [BazinLnPrior] and specifies the natural logarithm of the prior
/// to use. Some curve-fit algorithms doesn't support this and ignores the prior
@@ -405,6 +406,8 @@ mod tests {
use super::*;
use crate::nl_fit::LnPrior1D;
use crate::tests::*;
#[cfg(feature = "ceres-source")]
use crate::CeresCurveFit;
#[cfg(feature = "gsl")]
use crate::LmsderCurveFit;
use crate::TimeSeries;
@@ -456,6 +459,61 @@ mod tests {
assert_relative_eq!(&values[..5], &desired[..], max_relative = 0.01);
}

#[cfg(feature = "ceres-source")]
#[test]
fn bazin_fit_noisy_ceres() {
bazin_fit_noisy(BazinFit::new(
CeresCurveFit::default().into(),
LnPrior::none(),
BazinInitsBounds::Default,
));
}

#[cfg(feature = "ceres-source")]
#[test]
fn bazin_fit_noizy_mcmc_plus_ceres() {
let ceres = CeresCurveFit::default();
let mcmc = McmcCurveFit::new(512, Some(ceres.into()));
bazin_fit_noisy(BazinFit::new(
mcmc.into(),
LnPrior::none(),
BazinInitsBounds::Default,
));
}

// Currently fails, we need better support of bounds from ceres-solver Rust crate
// #[cfg(feature = "ceres-source")]
// #[test]
// fn bazin_fit_noizy_ceres_with_bounds() {
// let lmsder = CeresCurveFit::new();
// let mcmc = McmcCurveFit::new(512, Some(lmsder.into()));
// bazin_fit_noisy(BazinFit::new(
// CeresCurveFit::new().into(),
// LnPrior::none(),
// BazinInitsBounds::option_arrays(
// [None; 5],
// [None; 5],
// [None, None, Some(50.0), Some(50.0), Some(50.0)],
// ),
// ));
// }

#[cfg(feature = "ceres-source")]
#[test]
fn bazin_fit_noizy_mcmc_plus_ceres_and_bounds() {
let ceres = CeresCurveFit::default();
let mcmc = McmcCurveFit::new(10, Some(ceres.into()));
bazin_fit_noisy(BazinFit::new(
mcmc.into(),
LnPrior::none(),
BazinInitsBounds::option_arrays(
[None; 5],
[None; 5],
[None, None, Some(50.0), Some(50.0), Some(50.0)],
),
));
}

#[cfg(feature = "gsl")]
#[test]
fn bazin_fit_noisy_lmsder() {
28 changes: 27 additions & 1 deletion src/features/villar_fit.rs
Original file line number Diff line number Diff line change
@@ -593,6 +593,9 @@ impl From<LnPrior<NPARAMS>> for VillarLnPrior {
#[allow(clippy::excessive_precision)]
mod tests {
use super::*;

#[cfg(feature = "ceres-source")]
use crate::nl_fit::CeresCurveFit;
#[cfg(feature = "gsl")]
use crate::nl_fit::LnPrior1D;
use crate::tests::*;
@@ -613,7 +616,7 @@ mod tests {
[0.0; 11],
);

#[cfg(feature = "gsl")]
#[cfg(any(feature = "gsl", feature = "ceres-source"))]
fn villar_fit_noisy(eval: VillarFit) {
const N: usize = 50;

@@ -655,6 +658,29 @@ mod tests {
assert_relative_eq!(&values[..NPARAMS], &desired[..], max_relative = 0.01);
}

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

#[cfg(feature = "ceres-source")]
#[test]
fn villar_fit_noizy_mcmc_plus_ceres() {
let ceres = CeresCurveFit::default();
let mcmc = McmcCurveFit::new(512, Some(ceres.into()));
villar_fit_noisy(VillarFit::new(
mcmc.into(),
LnPrior::none(),
VillarInitsBounds::Default,
));
}

#[cfg(feature = "gsl")]
#[test]
fn villar_fit_noisy_lmsder() {
2 changes: 2 additions & 0 deletions src/lib.rs
Original file line number Diff line number Diff line change
@@ -29,6 +29,8 @@ mod lnerfc;

mod nl_fit;
pub use nl_fit::evaluator::FitFeatureEvaluatorGettersTrait;
#[cfg(feature = "ceres-source")]
pub use nl_fit::CeresCurveFit;
#[cfg(feature = "gsl")]
pub use nl_fit::LmsderCurveFit;
pub use nl_fit::{prior, LnPrior, LnPrior1D};
15 changes: 15 additions & 0 deletions src/nl_fit/bounds.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
pub(super) fn within_bounds<T, const NPARAMS: usize>(
x: &[T; NPARAMS],
lower: &[T; NPARAMS],
upper: &[T; NPARAMS],
) -> bool
where
T: PartialOrd,
{
for i in 0..NPARAMS {
if x[i] < lower[i] || x[i] > upper[i] {
return false;
}
}
true
}
201 changes: 201 additions & 0 deletions src/nl_fit/ceres.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,201 @@
use crate::nl_fit::constants::PARAMETER_TOLERANCE;
use crate::nl_fit::curve_fit::{CurveFitResult, CurveFitTrait};
use crate::nl_fit::data::Data;

use ceres_solver::{CurveFitProblem1D, CurveFunctionType, LossFunction, SolverOptions};
use ndarray::Zip;
use schemars::JsonSchema;
use serde::{Deserialize, Serialize};
use std::rc::Rc;

/// Ceres-Solver non-linear least-squares wrapper
///
/// Requires `ceres` Cargo feature
///
/// Non-linear squares-based light-curve fitters. It requires the function Jacobean. It supports
/// boundaries, but not priors.
#[derive(Clone, Debug, Serialize, Deserialize, JsonSchema)]
#[serde(rename = "Ceres")]
pub struct CeresCurveFit {
niterations: u16,
loss_factor: Option<f64>,
}

impl CeresCurveFit {
/// Create a new [CeresCurveFit].
///
/// # Arguments
/// - `niterations`: number of iterations
/// - `loss_factor`: if specified, use Huber loss function with the given factor to transform
/// the squared norm of the residuals. This is useful to reduce the influence of outliers.
pub fn new(niterations: u16, loss_factor: Option<f64>) -> Self {
if let Some(loss_factor) = loss_factor {
assert!(loss_factor > 0.0, "loss_factor must be positive");
}
Self {
niterations,
loss_factor,
}
}

#[inline]
pub fn default_niterations() -> u16 {
10
}

#[inline]
pub fn default_loss_factor() -> Option<f64> {
None
}
}

impl Default for CeresCurveFit {
fn default() -> Self {
Self::new(Self::default_niterations(), Self::default_loss_factor())
}
}

impl CurveFitTrait for CeresCurveFit {
fn curve_fit<F, DF, LP, const NPARAMS: usize>(
&self,
ts: Rc<Data<f64>>,
x0: &[f64; NPARAMS],
bounds: (&[f64; NPARAMS], &[f64; NPARAMS]),
model: F,
derivatives: DF,
_ln_prior: LP,
) -> CurveFitResult<f64, NPARAMS>
where
F: 'static + Clone + Fn(f64, &[f64; NPARAMS]) -> f64,
DF: 'static + Clone + Fn(f64, &[f64; NPARAMS], &mut [f64; NPARAMS]),
LP: Clone + Fn(&[f64; NPARAMS]) -> f64,
{
let func: CurveFunctionType = {
let model = model.clone();
Box::new(move |t, parameters, y, jacobians| {
let parameters = parameters.try_into().unwrap();
*y = model(t, parameters);
if !y.is_finite() {
*y = f64::MAX.sqrt();
return false;
}
if let Some(jacobians) = jacobians {
let jacobians: &mut [_; NPARAMS] = jacobians.try_into().unwrap();
let der = {
let mut der = [0.0; NPARAMS];
derivatives(t, parameters, &mut der);
der
};
for (input, output) in der.into_iter().zip(jacobians.iter_mut()) {
if let Some(output) = output {
if !input.is_finite() {
return false;
}
*output = input;
}
}
}
true
})
};

let lower_bounds: Vec<_> = bounds.0.iter().map(|&v| Some(v)).collect();
let upper_bounds: Vec<_> = bounds.1.iter().map(|&v| Some(v)).collect();

let options = SolverOptions::builder()
.parameter_tolerance(PARAMETER_TOLERANCE)
.max_num_iterations(self.niterations as i32)
.build()
.unwrap();

let mut problem_builder = CurveFitProblem1D::builder()
.x(ts.t.as_slice().unwrap())
.y(ts.m.as_slice().unwrap())
.inverse_error(ts.inv_err.as_slice().unwrap())
.func(func)
.parameters(x0)
.lower_bounds(&lower_bounds)
.upper_bounds(&upper_bounds);
if let Some(loss_factor) = self.loss_factor {
problem_builder = problem_builder.loss(LossFunction::cauchy(loss_factor));
};
let solution = problem_builder.build().unwrap().solve(&options);
let x = solution.parameters.try_into().unwrap();
let success = solution.summary.is_solution_usable();

let reduced_chi2 = Zip::from(&ts.t)
.and(&ts.m)
.and(&ts.inv_err)
.fold(0.0, |acc, &t, &m, &inv_err| {
acc + ((model(t, &x) - m) * inv_err).powi(2)
})
/ (ts.t.len() - NPARAMS) as f64;
CurveFitResult {
x,
reduced_chi2,
success,
}
}
}

#[cfg(test)]
#[allow(clippy::unreadable_literal)]
#[allow(clippy::excessive_precision)]
mod tests {
use super::*;

use approx::assert_abs_diff_eq;
use ndarray::Array1;
use rand::prelude::*;
use rand_distr::StandardNormal;

fn nonlinear_func(t: f64, param: &[f64; 3]) -> f64 {
param[1] * f64::exp(-param[0] * t) * t.powi(2) + param[2]
}

fn nonlinear_func_derivatives(t: f64, param: &[f64; 3], derivatives: &mut [f64; 3]) {
derivatives[0] = -param[1] * f64::exp(-param[2] * t) * t.powi(3);
derivatives[1] = f64::exp(-param[0] * t) * t.powi(2);
derivatives[2] = 1.0;
}

fn nonlinear_func_dump_ln_prior(_param: &[f64; 3]) -> f64 {
0.0
}

#[test]
fn nonlinear() {
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];

let mut rng = StdRng::seed_from_u64(0);

let t = Array1::linspace(0.0, 10.0, N);
let y = t.mapv(|x| {
let eps: f64 = rng.sample(StandardNormal);
nonlinear_func(x, &param_true) + NOISE * eps
});
let inv_err: Array1<_> = vec![1.0 / NOISE; N].into();
let ts = Rc::new(Data { t, m: y, inv_err });

let fitter = CeresCurveFit::new(14, None);
let result = fitter.curve_fit(
ts,
&param_init,
(&[f64::NEG_INFINITY; 3], &[f64::INFINITY; 3]),
nonlinear_func,
nonlinear_func_derivatives,
nonlinear_func_dump_ln_prior,
);

// 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);
}
}
1 change: 1 addition & 0 deletions src/nl_fit/constants.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
pub(super) const PARAMETER_TOLERANCE: f64 = 1e-4;
4 changes: 4 additions & 0 deletions src/nl_fit/curve_fit.rs
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
#[cfg(feature = "ceres-source")]
use crate::nl_fit::ceres::CeresCurveFit;
use crate::nl_fit::data::Data;
#[cfg(feature = "gsl")]
use crate::nl_fit::lmsder::LmsderCurveFit;
@@ -39,6 +41,8 @@ pub trait CurveFitTrait: Clone + Debug + Serialize + DeserializeOwned {
#[derive(Clone, Debug, Serialize, Deserialize, JsonSchema)]
#[non_exhaustive]
pub enum CurveFitAlgorithm {
#[cfg(feature = "ceres-source")]
Ceres(CeresCurveFit),
#[cfg(feature = "gsl")]
Lmsder(LmsderCurveFit),
Mcmc(McmcCurveFit),
5 changes: 3 additions & 2 deletions src/nl_fit/lmsder.rs
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#[cfg(test)]
use crate::float_trait::Float;
use crate::nl_fit::constants::PARAMETER_TOLERANCE;
use crate::nl_fit::curve_fit::{CurveFitResult, CurveFitTrait};
use crate::nl_fit::data::Data;

@@ -120,8 +121,8 @@ impl NlsProblem {
fn new(fit_function: MultiFitFunctionFdf) -> Self {
Self {
max_iter: Self::default_max_iter(),
atol: 0.0,
rtol: 1e-4,
atol: PARAMETER_TOLERANCE,
rtol: PARAMETER_TOLERANCE,
fit_function,
}
}
15 changes: 7 additions & 8 deletions src/nl_fit/mcmc.rs
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
use crate::float_trait::Float;
use crate::nl_fit::bounds::within_bounds;
use crate::nl_fit::curve_fit::{CurveFitAlgorithm, CurveFitResult, CurveFitTrait};
use crate::nl_fit::data::Data;

@@ -236,14 +237,12 @@ where
}

fn lnprior(&self, params: &Guess) -> f32 {
for (&p, (&lower, &upper)) in params
.values
.iter()
.zip(self.lower.iter().zip(self.upper.iter()))
{
if p < lower || p > upper {
return f32::NEG_INFINITY;
}
if !within_bounds(
(&params.values[..]).try_into().unwrap(),
self.lower,
self.upper,
) {
return f32::NEG_INFINITY;
}
(self.ln_prior)(params)
}
10 changes: 10 additions & 0 deletions src/nl_fit/mod.rs
Original file line number Diff line number Diff line change
@@ -1,3 +1,13 @@
pub(self) mod bounds;

#[cfg(feature = "ceres-source")]
pub mod ceres;
#[cfg(feature = "ceres-source")]
pub use ceres::CeresCurveFit;

#[cfg(any(feature = "gsl", feature = "ceres-source"))]
pub(self) mod constants;

pub mod curve_fit;
pub use curve_fit::{CurveFitAlgorithm, CurveFitResult, CurveFitTrait};

0 comments on commit e4028c7

Please sign in to comment.