use argmin::prelude::*;
use argmin::solver::gaussnewton::GaussNewton;
use nalgebra::{DMatrix, DVector};
type Rate = f64;
type S = f64;
type Measurement = (S, Rate);
struct Problem {
data: Vec<Measurement>,
}
impl ArgminOp for Problem {
type Param = DVector<f64>;
type Output = DVector<f64>;
type Hessian = ();
type Jacobian = DMatrix<f64>;
type Float = f64;
fn apply(&self, p: &Self::Param) -> Result<Self::Output, Error> {
Ok(DVector::from_vec(
self.data
.iter()
.map(|(s, rate)| rate - (p[0] * s) / (p[1] + s))
.collect(),
))
}
fn jacobian(&self, p: &Self::Param) -> Result<Self::Jacobian, Error> {
Ok(DMatrix::from_fn(7, 2, |si, i| {
if i == 0 {
-self.data[si].0 / (p[1] + self.data[si].0)
} else {
p[0] * self.data[si].0 / (p[1] + self.data[si].0).powi(2)
}
}))
}
}
fn run() -> Result<(), Error> {
let cost = Problem {
data: vec![
(0.038, 0.050),
(0.194, 0.127),
(0.425, 0.094),
(0.626, 0.2122),
(1.253, 0.2729),
(2.5, 0.2665),
(3.74, 0.3317),
],
};
let init_param: DVector<f64> = DVector::from_vec(vec![0.9, 0.2]);
let solver: GaussNewton<f64> = GaussNewton::new();
let res = Executor::new(cost, solver, init_param)
.add_observer(ArgminSlogLogger::term(), ObserverMode::Always)
.max_iters(10)
.run()?;
std::thread::sleep(std::time::Duration::from_secs(1));
println!("{}", res);
Ok(())
}
fn main() {
if let Err(ref e) = run() {
println!("{}", e);
std::process::exit(1);
}
}