diff --git a/Cargo.toml b/Cargo.toml index fcee07c..89305c3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,7 +11,7 @@ license = "GPL-3.0+" edition = "2021" [dependencies] -sys = { path = "gsl-sys", package = "GSL-sys", version = "3.0.0" } +sys = { path = "gsl-sys", package = "GSL-sys", version = "3.0.0", features = ["v2_1"] } paste = "1.0" [features] diff --git a/examples/Cargo.toml b/examples/Cargo.toml index ada63a1..7604638 100644 --- a/examples/Cargo.toml +++ b/examples/Cargo.toml @@ -1,126 +1,130 @@ -[package] -name = "gsl-examples" -version = "0.4.6" -authors = ["Guillaume Gomez "] - -[dependencies.GSL] -path = "../" - -[dependencies] -libc = "~0.2" - -[[bin]] -name = "blas" -path = "./blas.rs" - -[[bin]] -name = "bspline" -path = "./bspline.rs" - -[[bin]] -name = "cblas" -path = "./cblas.rs" - -[[bin]] -name = "cdf" -path = "./cdf.rs" - -[[bin]] -name = "chebyshev_approximation" -path = "./chebyshev_approximation.rs" - -[[bin]] -name = "combination" -path = "./combination.rs" - -[[bin]] -name = "diff" -path = "./diff.rs" - -[[bin]] -name = "eigen" -path = "./eigen.rs" - -[[bin]] -name = "eigen_nonsymm" -path = "./eigen_nonsymm.rs" - -[[bin]] -name = "fft" -path = "./fft.rs" - -[[bin]] -name = "fftmr" -path = "./fftmr.rs" - -[[bin]] -name = "filt_edge" -path = "./filt_edge.rs" -required-features = ["GSL/v2_5"] - -[[bin]] -name = "fitreg" -path = "./fitreg.rs" -required-features = ["GSL/v2_1"] - -[[bin]] -name = "fitreg2" -path = "./fitreg2.rs" -required-features = ["GSL/v2_1"] - -[[bin]] -name = "fitting" -path = "./fitting.rs" -required-features = ["GSL/v2_5"] - -[[bin]] -name = "fitting3" -path = "./fitting3.rs" - -[[bin]] -name = "gaussfilt" -path = "./gaussfilt.rs" -required-features = ["GSL/v2_5"] - -[[bin]] -name = "gaussfilt2" -path = "./gaussfilt2.rs" -required-features = ["GSL/v2_5"] - -[[bin]] -name = "histogram2d" -path = "./histogram2d.rs" - -[[bin]] -name = "impulse" -path = "./impulse.rs" -required-features = ["GSL/v2_5"] - -[[bin]] -name = "integration" -path = "./integration.rs" - -[[bin]] -name = "integration2" -path = "./integration2.rs" - -[[bin]] -name = "intro" -path = "./intro.rs" - -[[bin]] -name = "largefit" -path = "./largefit.rs" -required-features = ["GSL/v2_2"] - -[[bin]] -name = "rng" -path = "./rng.rs" - -[[bin]] -name = "statistics" -path = "./statistics.rs" - -[[bin]] -name = "vectors_and_matrices" -path = "./vectors_and_matrices.rs" +[package] +name = "gsl-examples" +version = "0.4.6" +authors = ["Guillaume Gomez "] + +[dependencies.GSL] +path = "../" + +[dependencies] +libc = "~0.2" + +[[bin]] +name = "blas" +path = "./blas.rs" + +[[bin]] +name = "bspline" +path = "./bspline.rs" + +[[bin]] +name = "cblas" +path = "./cblas.rs" + +[[bin]] +name = "cdf" +path = "./cdf.rs" + +[[bin]] +name = "chebyshev_approximation" +path = "./chebyshev_approximation.rs" + +[[bin]] +name = "combination" +path = "./combination.rs" + +[[bin]] +name = "diff" +path = "./diff.rs" + +[[bin]] +name = "eigen" +path = "./eigen.rs" + +[[bin]] +name = "eigen_nonsymm" +path = "./eigen_nonsymm.rs" + +[[bin]] +name = "fft" +path = "./fft.rs" + +[[bin]] +name = "fftmr" +path = "./fftmr.rs" + +[[bin]] +name = "filt_edge" +path = "./filt_edge.rs" +required-features = ["GSL/v2_5"] + +[[bin]] +name = "fitreg" +path = "./fitreg.rs" +required-features = ["GSL/v2_1"] + +[[bin]] +name = "fitreg2" +path = "./fitreg2.rs" +required-features = ["GSL/v2_1"] + +[[bin]] +name = "fitting" +path = "./fitting.rs" +required-features = ["GSL/v2_5"] + +[[bin]] +name = "fitting3" +path = "./fitting3.rs" + +[[bin]] +name = "gaussfilt" +path = "./gaussfilt.rs" +required-features = ["GSL/v2_5"] + +[[bin]] +name = "gaussfilt2" +path = "./gaussfilt2.rs" +required-features = ["GSL/v2_5"] + +[[bin]] +name = "histogram2d" +path = "./histogram2d.rs" + +[[bin]] +name = "impulse" +path = "./impulse.rs" +required-features = ["GSL/v2_5"] + +[[bin]] +name = "integration" +path = "./integration.rs" + +[[bin]] +name = "integration2" +path = "./integration2.rs" + +[[bin]] +name = "intro" +path = "./intro.rs" + +[[bin]] +name = "largefit" +path = "./largefit.rs" +required-features = ["GSL/v2_2"] + +[[bin]] +name = "multifit_nlinear" +path = "./multifit_nlinear.rs" + +[[bin]] +name = "rng" +path = "./rng.rs" + +[[bin]] +name = "statistics" +path = "./statistics.rs" + +[[bin]] +name = "vectors_and_matrices" +path = "./vectors_and_matrices.rs" diff --git a/examples/multifit_nlinear.rs b/examples/multifit_nlinear.rs new file mode 100644 index 0000000..c5baa2c --- /dev/null +++ b/examples/multifit_nlinear.rs @@ -0,0 +1,84 @@ +extern crate rgsl; + +use rgsl::gsl_multifit_nlinear_basic; +use rgsl::gsl_multifit_nlinear_basic_df; + +use rgsl::types::rng::Rng; +use rgsl::types::rng::RngType; + + +fn expb_f(params: Vec, t: f64, _args: Vec) -> f64 { + + let a = params.get(0).unwrap(); + let lambda = params.get(1).unwrap(); + let b = params.get(2).unwrap(); + + a * f64::exp(-lambda * t) + b +} + +fn expb_df_a(params: Vec, t: f64, _args: Vec) -> f64 { + let lambda = params.get(1).unwrap(); + f64::exp(-lambda * t) +} + +fn expb_df_lambda(params: Vec, t: f64, _args: Vec) -> f64 { + let a = params.get(0).unwrap(); + let lambda = params.get(1).unwrap(); + -t * a * f64::exp(-lambda * t) +} + +fn expb_df_b(_params: Vec, _t: f64, _args: Vec) -> f64 { + 1.0 +} + +fn main() { + + let params_out = vec![1.0, 1.0, 0.0]; + let mut ts = Vec::new(); + let mut ys = Vec::new(); + let args = vec![]; + let expb_dfs = vec![expb_df_a, expb_df_lambda, expb_df_b]; + + let mut rng = Rng::new(RngType::default()).unwrap(); + + for i in 0..100 { + + let rand_flt = rng.uniform(); + + let ti = (i as f64) * 3.0 / (100.0 - 1.0); + let yi = 1.0 + 5.0 * f64::exp(-1.5 * ti); + let si = 0.1 * yi; + let dy = si * rand_flt; + + ts.push(ti); + ys.push(yi + dy); + } + + let (params, parerr, status) = match unsafe { + gsl_multifit_nlinear_basic(expb_f, params_out.clone(), &ts, &ys, &args, 100) + } { + Ok(value) => value, + Err(_) => { + eprintln!("Error encountered during solve!"); + return; + } + }; + + println!("{:?}", params); + println!("{:?}", parerr); + println!("{}", status); + + let (params, parerr, status) = match unsafe { + gsl_multifit_nlinear_basic_df(expb_f, &expb_dfs, params_out.clone(), &ts, &ys, &args, 100) + } { + Ok(value) => value, + Err(_) => { + eprintln!("Error encountered during solve!"); + return; + } + }; + + println!("{:?}", params); + println!("{:?}", parerr); + println!("{}", status); +} diff --git a/gsl-sys/src/lib.rs b/gsl-sys/src/lib.rs index 5b88564..7763a0c 100644 --- a/gsl-sys/src/lib.rs +++ b/gsl-sys/src/lib.rs @@ -10,8 +10,814 @@ #![allow(clippy::redundant_static_lifetimes)] #![allow(clippy::upper_case_acronyms)] + pub extern crate libc; mod auto; - pub use auto::*; + +use std::ptr; + + +macro_rules! result_handler { + ($ret:ident, $value:expr) => {{ + if $ret == GSL_SUCCESS { + Ok($value) + } else { + Err($ret) + } + }}; +} + + +#[repr(C)] +#[derive(Debug, Clone)] +struct multifit_nlinear_data { + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: Vec, f64, Vec) -> f64>, + params_len: usize, + ts: Vec, + ys: Vec, + args: Vec +} + +#[repr(C)] +#[derive(Debug, Copy, Clone)] +pub struct gsl_multifit_nlinear_parameters { + trs: *const gsl_multifit_nlinear_trs, + scale: *const gsl_multifit_nlinear_scale, + solver: *const gsl_multifit_nlinear_solver, + fdtype: gsl_multifit_nlinear_fdtype, + factor_up: f64, + factor_down: f64, + avmax: f64, + h_df: f64, + h_fvv: f64 +} + +#[repr(C)] +#[derive(Debug, Copy, Clone)] +pub struct gsl_multifit_nlinear_fdf { + f: fn(*const gsl_vector, *mut multifit_nlinear_data, *mut gsl_vector) -> i32, + df: *const fn(*const gsl_vector, *mut multifit_nlinear_data, *mut gsl_matrix) -> i32, + fvv: *const fn(*const gsl_vector, *const gsl_vector, *mut multifit_nlinear_data, *mut gsl_vector) -> i32, + n: usize, + p: usize, + params: *mut multifit_nlinear_data, + nevalf: usize, + nevaldf: usize, + nevalfvv: usize +} + +extern "C" { + pub fn gsl_multifit_nlinear_default_parameters() -> self::gsl_multifit_nlinear_parameters; +} +extern "C" { + pub fn gsl_multifit_nlinear_alloc( + T: *const gsl_multifit_nlinear_type, + params: *const self::gsl_multifit_nlinear_parameters, + n: usize, + p: usize, + ) -> *mut gsl_multifit_nlinear_workspace; +} + +extern "C" { + pub fn gsl_multifit_nlinear_init( + x: *const gsl_vector, + fdf: *mut self::gsl_multifit_nlinear_fdf, + w: *mut gsl_multifit_nlinear_workspace, + ) -> ::std::os::raw::c_int; +} + +extern "C" { + pub fn gsl_multifit_nlinear_winit( + x: *const gsl_vector, + wts: *const gsl_vector, + fdf: *mut self::gsl_multifit_nlinear_fdf, + w: *mut gsl_multifit_nlinear_workspace, + ) -> ::std::os::raw::c_int; +} + +fn call_multifit_f(x: *const gsl_vector, data: *mut multifit_nlinear_data, f: *mut gsl_vector) -> i32 { + + let n = unsafe { (*data).ts.len() }; + + let ts: &Vec = unsafe { &(*data).ts }; + let ys: &Vec = unsafe { &(*data).ys }; + let args: &Vec = unsafe { &(*data).args }; + + let mut params: Vec = Vec::new(); + let params_len = unsafe { (*data).params_len }; + + for i in 0..params_len { + let param = unsafe { gsl_vector_get(x, i) }; + params.push(param); + } + + for i in 0..n { + let func_f = unsafe { (*data).func_f }; + let yi = func_f(params.clone(), ts[i], args.clone()); + unsafe { gsl_vector_set(f, i, yi - ys[i]); } + } + + GSL_SUCCESS +} + +fn call_multifit_df(x: *const gsl_vector, data: *mut multifit_nlinear_data, jacobian: *mut gsl_matrix) -> i32 { + + let n = unsafe { (*data).ts.len() }; + + let ts: &Vec = unsafe { &(*data).ts }; + let func_dfs: &Vec, f64, Vec) -> f64> = unsafe { &(*data).func_dfs }; + let args: &Vec = unsafe { &(*data).args }; + + let mut params: Vec = Vec::new(); + let params_len = unsafe { (*data).params_len }; + + for i in 0..params_len { + let param = unsafe { gsl_vector_get(x, i) }; + params.push(param); + } + + for i in 0..n { + for j in 0..params_len { + unsafe { gsl_matrix_set(jacobian, i, j, func_dfs[j](params.clone(), ts[i], args.clone())); } + } + } + + GSL_SUCCESS +} + +fn call_multifit_fvv(x: *const gsl_vector, v: *const gsl_vector, params: *mut multifit_nlinear_data, fvv: *mut gsl_vector) -> i32 { + GSL_SUCCESS +} + +pub unsafe fn gsl_multifit_nlinear_basic( + func_f: fn(Vec, f64, Vec) -> f64, + params_in: Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + gsl_multifit_nlinear_basic_df( + func_f, + &Vec::new(), + params_in, + ts, + ys, + args, + max_iters + ) +} + +pub unsafe fn gsl_multifit_nlinear_basic_df( + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: &Vec, f64, Vec) -> f64>, + params_in: Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + let mut params: Vec = params_in.clone(); + let mut parerr: Vec = Vec::with_capacity(params_in.len()); + let mut status: i32 = 0; + let mut err_code: i32 = GSL_FAILURE; + + parerr.resize(params_in.len(), 0.0); + + if ts.len() != ys.len() { + eprintln!("Time length does not match Ys length!"); + } else if params_in.len() > ts.len() { + eprintln!("Time length is shorter than parameter length!"); + } else { + run_gsl_multifit_nlinear_df( + func_f, + func_dfs, + params.as_mut_ptr(), + parerr.as_mut_ptr(), + params.len(), + &vec![], + ts, + ys, + args, + max_iters, + &mut status + ); + err_code = GSL_SUCCESS; + } + + result_handler!(err_code, (params, parerr, status)) +} + +pub unsafe fn gsl_multifit_nlinear_weighted( + func_f: fn(Vec, f64, Vec) -> f64, + params_in: Vec, + weights_in: &Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + gsl_multifit_nlinear_weighted_df( + func_f, + &Vec::new(), + params_in, + weights_in, + ts, + ys, + args, + max_iters + ) +} + +pub unsafe fn gsl_multifit_nlinear_weighted_df( + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: &Vec, f64, Vec) -> f64>, + params_in: Vec, + weights_in: &Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + let mut params: Vec = params_in.clone(); + let mut parerr: Vec = Vec::with_capacity(params_in.len()); + let mut status: i32 = 0; + let mut err_code: i32 = GSL_FAILURE; + + parerr.resize(params_in.len(), 0.0); + + if ts.len() != ys.len() { + eprintln!("Time length does not match Ys length!"); + } else if ts.len() != weights_in.len() { + eprintln!("Time length does not match weights length!"); + } else if params_in.len() > ts.len() { + eprintln!("Time length is shorter than parameter length!"); + } else { + run_gsl_multifit_nlinear_df( + func_f, + func_dfs, + params.as_mut_ptr(), + parerr.as_mut_ptr(), + params.len(), + weights_in, + ts, + ys, + args, + max_iters, + &mut status + ); + err_code = GSL_SUCCESS; + } + + result_handler!(err_code, (params, parerr, status)) +} +fn run_gsl_multifit_nlinear_df( + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: &Vec, f64, Vec) -> f64>, + params: *mut f64, + parerr: *mut f64, + params_len: usize, + weights_in: &Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64, + status: &mut i32 +) { + + let vars_len = ts.len(); + + let trust: *const gsl_multifit_nlinear_type = unsafe { gsl_multifit_nlinear_trust }; + let fdf_params: self::gsl_multifit_nlinear_parameters = unsafe { self::gsl_multifit_nlinear_default_parameters() }; + + let f: *mut gsl_vector; + let jacobian: *mut gsl_matrix; + let covar: *mut gsl_matrix = unsafe { gsl_matrix_alloc(params_len, params_len) }; + + let mut data_struct = multifit_nlinear_data { + func_f: func_f, + func_dfs: func_dfs.to_vec(), + params_len: params_len, + ts: ts.to_vec(), + ys: ys.to_vec(), + args: args.to_vec() + }; + let data_ptr: *mut multifit_nlinear_data = &mut data_struct as *mut multifit_nlinear_data; + let x = unsafe { gsl_vector_view_array(params, params_len) }; + let x_ptr: *const gsl_vector = &x.vector as *const gsl_vector; + let dof: f64 = (vars_len - params_len) as f64; + + let mut chisq: f64 = 0.0; + let mut chisq0: f64 = 0.0; + let mut info: i32 = 0; + + const xtol: f64 = 1e-8; + const gtol: f64 = 1e-8; + const ftol: f64 = 1e-8; + + let call_multifit_df_ptr = if func_dfs.len() == params_len { + call_multifit_df as *const fn(*const gsl_vector, *mut multifit_nlinear_data, *mut gsl_matrix) -> i32 + } else { + ptr::null() + }; + + let mut fdf = self::gsl_multifit_nlinear_fdf { + f: call_multifit_f, + df: call_multifit_df_ptr, + fvv: ptr::null(), + n: vars_len, + p: params_len, + params: data_ptr, + nevalf: 0, + nevaldf: 0, + nevalfvv: 0 + }; + let fdf_ptr: *mut self::gsl_multifit_nlinear_fdf = &mut fdf as *mut self::gsl_multifit_nlinear_fdf; + + /* allocate workspace with default parameters */ + let w: *mut gsl_multifit_nlinear_workspace = unsafe { + self::gsl_multifit_nlinear_alloc( + trust, + &fdf_params as *const self::gsl_multifit_nlinear_parameters, + vars_len, + params_len + ) + }; + + /* initialize solver with starting point */ + if weights_in.len() == ts.len() { + let weights = unsafe { gsl_vector_view_array(weights_in.clone().as_mut_ptr(), weights_in.len()) }; + let weights_ptr: *const gsl_vector = &weights.vector as *const gsl_vector; + unsafe { + self::gsl_multifit_nlinear_winit ( + x_ptr, + weights_ptr, + fdf_ptr, + w + ); + } + } else { + unsafe { + self::gsl_multifit_nlinear_init ( + x_ptr, + fdf_ptr, + w + ); + } + } + /* compute initial cost function */ + let f = unsafe { gsl_multifit_nlinear_residual(w) }; + unsafe { + gsl_blas_ddot(f, f, &mut chisq0); + } + + /* solve the system within a maximum of max_iters iterations */ + for i in 0..max_iters { + + let mut rcond: f64 = 0.0; + + unsafe { gsl_multifit_nlinear_iterate(w); } + unsafe { gsl_multifit_nlinear_test(xtol, gtol, ftol, &mut *status, w); } + + if *status != 0 { + + if *status == 1 { + println!("Reason for stopping: Small Step Size"); + } else if *status == 2 { + println!("Reason for stopping: Small Gradient"); + } + + break; + } + + unsafe { gsl_multifit_nlinear_rcond(&mut rcond, w); } + if rcond.is_nan() && i != 0 { + println!("Reason for stopping: Invalid Status"); + *status = -1; + + break; + } + } + + /* compute covariance of best fit parameters */ + let jacobian = unsafe { gsl_multifit_nlinear_jac(w) }; + unsafe { gsl_multifit_nlinear_covar(jacobian, 0.0, covar); } + + /* compute final cost */ + let residuals = unsafe { gsl_multifit_nlinear_residual(w) }; + unsafe { gsl_blas_ddot(residuals, residuals, &mut chisq); } + let chisq_dof = (chisq / dof).sqrt(); + + let x_gsl_vec: *const gsl_vector = unsafe { gsl_multifit_nlinear_position(w) }; + for i in 0..params_len { + unsafe { params.wrapping_add(i).write(gsl_vector_get(x_gsl_vec, i)); } + unsafe { parerr.wrapping_add(i).write(chisq_dof * gsl_matrix_get(covar, i, i).sqrt()); } + } + + unsafe { gsl_multifit_nlinear_free(w); } + unsafe { gsl_matrix_free(covar); } +} + +#[repr(C)] +#[derive(Debug, Clone)] +struct multilarge_nlinear_data { + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: Vec, f64, Vec) -> f64>, + params_len: usize, + ts: Vec, + ys: Vec, + args: Vec +} + +#[repr(C)] +#[derive(Debug, Copy, Clone)] +pub struct gsl_multilarge_nlinear_parameters { + trs: *const gsl_multilarge_nlinear_trs, + scale: *const gsl_multilarge_nlinear_scale, + solver: *const gsl_multilarge_nlinear_solver, + fdtype: gsl_multilarge_nlinear_fdtype, + factor_up: f64, + factor_down: f64, + avmax: f64, + h_df: f64, + h_fvv: f64 +} + +#[repr(C)] +#[derive(Debug, Copy, Clone)] +pub struct gsl_multilarge_nlinear_fdf { + f: fn(*const gsl_vector, *mut multilarge_nlinear_data, *mut gsl_vector) -> i32, + df: *const fn(*const gsl_vector, *mut multilarge_nlinear_data, *mut gsl_matrix) -> i32, + fvv: *const fn(*const gsl_vector, *const gsl_vector, *mut multilarge_nlinear_data, *mut gsl_vector) -> i32, + n: usize, + p: usize, + params: *mut multilarge_nlinear_data, + nevalf: usize, + nevaldfu: usize, + nevaldf2: usize, + nevalfvv: usize +} + +extern "C" { + pub fn gsl_multilarge_nlinear_default_parameters() -> self::gsl_multilarge_nlinear_parameters; +} +extern "C" { + pub fn gsl_multilarge_nlinear_alloc( + T: *const gsl_multilarge_nlinear_type, + params: *const self::gsl_multilarge_nlinear_parameters, + n: usize, + p: usize, + ) -> *mut gsl_multilarge_nlinear_workspace; +} + +extern "C" { + pub fn gsl_multilarge_nlinear_init( + x: *const gsl_vector, + fdf: *mut self::gsl_multilarge_nlinear_fdf, + w: *mut gsl_multilarge_nlinear_workspace, + ) -> ::std::os::raw::c_int; +} + +extern "C" { + pub fn gsl_multilarge_nlinear_winit( + x: *const gsl_vector, + wts: *const gsl_vector, + fdf: *mut self::gsl_multilarge_nlinear_fdf, + w: *mut gsl_multilarge_nlinear_workspace, + ) -> ::std::os::raw::c_int; +} + +fn call_multilarge_f(x: *const gsl_vector, data: *mut multilarge_nlinear_data, f: *mut gsl_vector) -> i32 { + + let n = unsafe { (*data).ts.len() }; + + let ts: &Vec = unsafe { &(*data).ts }; + let ys: &Vec = unsafe { &(*data).ys }; + let args: &Vec = unsafe { &(*data).args }; + + let mut params: Vec = Vec::new(); + let params_len = unsafe { (*data).params_len }; + + for i in 0..params_len { + let param = unsafe { gsl_vector_get(x, i) }; + params.push(param); + } + + for i in 0..n { + let func_f = unsafe { (*data).func_f }; + let yi = func_f(params.clone(), ts[i], args.clone()); + unsafe { gsl_vector_set(f, i, yi - ys[i]); } + } + + GSL_SUCCESS +} + +fn call_multilarge_df(x: *const gsl_vector, data: *mut multilarge_nlinear_data, jacobian: *mut gsl_matrix) -> i32 { + + let n = unsafe { (*data).ts.len() }; + + let ts: &Vec = unsafe { &(*data).ts }; + let func_dfs: &Vec, f64, Vec) -> f64> = unsafe { &(*data).func_dfs }; + let args: &Vec = unsafe { &(*data).args }; + + let mut params: Vec = Vec::new(); + let params_len = unsafe { (*data).params_len }; + + for i in 0..params_len { + let param = unsafe { gsl_vector_get(x, i) }; + params.push(param); + } + + for i in 0..n { + for j in 0..params_len { + unsafe { gsl_matrix_set(jacobian, i, j, func_dfs[j](params.clone(), ts[i], args.clone())); } + } + } + + GSL_SUCCESS +} + +fn call_multilarge_fvv(x: *const gsl_vector, v: *const gsl_vector, params: *mut multilarge_nlinear_data, fvv: *mut gsl_vector) -> i32 { + GSL_SUCCESS +} + +pub unsafe fn gsl_multilarge_nlinear_basic( + func_f: fn(Vec, f64, Vec) -> f64, + params_in: Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + gsl_multilarge_nlinear_basic_df( + func_f, + &Vec::new(), + params_in, + ts, + ys, + args, + max_iters + ) +} + +pub unsafe fn gsl_multilarge_nlinear_basic_df( + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: &Vec, f64, Vec) -> f64>, + params_in: Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + let mut params: Vec = params_in.clone(); + let mut parerr: Vec = Vec::with_capacity(params_in.len()); + let mut status: i32 = 0; + let mut err_code: i32 = GSL_FAILURE; + + parerr.resize(params_in.len(), 0.0); + + if ts.len() != ys.len() { + eprintln!("Time length does not match Ys length!"); + } else if params_in.len() > ts.len() { + eprintln!("Time length is shorter than parameter length!"); + } else { + run_gsl_multilarge_nlinear_df( + func_f, + func_dfs, + params.as_mut_ptr(), + parerr.as_mut_ptr(), + params.len(), + &vec![], + ts, + ys, + args, + max_iters, + &mut status + ); + err_code = GSL_SUCCESS; + } + + result_handler!(err_code, (params, parerr, status)) +} + +pub unsafe fn gsl_multilarge_nlinear_weighted( + func_f: fn(Vec, f64, Vec) -> f64, + params_in: Vec, + weights_in: &Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + gsl_multilarge_nlinear_weighted_df( + func_f, + &Vec::new(), + params_in, + weights_in, + ts, + ys, + args, + max_iters + ) +} + +pub unsafe fn gsl_multilarge_nlinear_weighted_df( + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: &Vec, f64, Vec) -> f64>, + params_in: Vec, + weights_in: &Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64 +) -> Result<(Vec, Vec, i32), i32> { + + let mut params: Vec = params_in.clone(); + let mut parerr: Vec = Vec::with_capacity(params_in.len()); + let mut status: i32 = 0; + let mut err_code: i32 = GSL_FAILURE; + + parerr.resize(params_in.len(), 0.0); + + if ts.len() != ys.len() { + eprintln!("Time length does not match Ys length!"); + } else if weights_in.len() != ts.len() { + eprintln!("Time length does not match weights length!"); + } else if params_in.len() > ts.len() { + eprintln!("Time length is shorter than parameter length!"); + } else { + run_gsl_multilarge_nlinear_df( + func_f, + func_dfs, + params.as_mut_ptr(), + parerr.as_mut_ptr(), + params.len(), + weights_in, + ts, + ys, + args, + max_iters, + &mut status + ); + err_code = GSL_SUCCESS; + } + + result_handler!(err_code, (params, parerr, status)) +} + +fn run_gsl_multilarge_nlinear_df( + func_f: fn(Vec, f64, Vec) -> f64, + func_dfs: &Vec, f64, Vec) -> f64>, + params: *mut f64, + parerr: *mut f64, + params_len: usize, + weights_in: &Vec, + ts: &Vec, + ys: &Vec, + args: &Vec, + max_iters: u64, + status: &mut i32 +) { + + let vars_len = ts.len(); + + let trust: *const gsl_multilarge_nlinear_type = unsafe { gsl_multilarge_nlinear_trust }; + let fdf_params: self::gsl_multilarge_nlinear_parameters = unsafe { self::gsl_multilarge_nlinear_default_parameters() }; + + let f: *mut gsl_vector; + let jacobian: *mut gsl_matrix; + let covar: *mut gsl_matrix = unsafe { gsl_matrix_alloc(params_len, params_len) }; + + let mut data_struct = multilarge_nlinear_data { + func_f: func_f, + func_dfs: func_dfs.to_vec(), + params_len: params_len, + ts: ts.to_vec(), + ys: ys.to_vec(), + args: args.to_vec() + }; + let data_ptr: *mut multilarge_nlinear_data = &mut data_struct as *mut multilarge_nlinear_data; + let x = unsafe { gsl_vector_view_array(params, params_len) }; + let x_ptr: *const gsl_vector = &x.vector as *const gsl_vector; + let dof: f64 = (vars_len - params_len) as f64; + + let mut chisq: f64 = 0.0; + let mut chisq0: f64 = 0.0; + let mut info: i32 = 0; + + const xtol: f64 = 1e-8; + const gtol: f64 = 1e-8; + const ftol: f64 = 1e-8; + + let call_multilarge_df_ptr = if func_dfs.len() == params_len { + call_multilarge_df as *const fn(*const gsl_vector, *mut multilarge_nlinear_data, *mut gsl_matrix) -> i32 + } else { + ptr::null() + }; + + let mut fdf = self::gsl_multilarge_nlinear_fdf { + f: call_multilarge_f, + df: call_multilarge_df_ptr, + fvv: ptr::null(), + n: vars_len, + p: params_len, + params: data_ptr, + nevalf: 0, + nevaldfu: 0, + nevaldf2: 0, + nevalfvv: 0 + }; + let fdf_ptr: *mut self::gsl_multilarge_nlinear_fdf = &mut fdf as *mut self::gsl_multilarge_nlinear_fdf; + + /* allocate workspace with default parameters */ + let w: *mut gsl_multilarge_nlinear_workspace = unsafe { + self::gsl_multilarge_nlinear_alloc( + trust, + &fdf_params as *const self::gsl_multilarge_nlinear_parameters, + vars_len, + params_len + ) + }; + + /* initialize solver with starting point */ + if weights_in.len() == ts.len() { + let weights = unsafe { gsl_vector_view_array(weights_in.clone().as_mut_ptr(), weights_in.len()) }; + let weights_ptr: *const gsl_vector = &weights.vector as *const gsl_vector; + unsafe { + self::gsl_multilarge_nlinear_winit ( + x_ptr, + weights_ptr, + fdf_ptr, + w + ); + } + } else { + unsafe { + self::gsl_multilarge_nlinear_init ( + x_ptr, + fdf_ptr, + w + ); + } + } + /* compute initial cost function */ + let f = unsafe { gsl_multilarge_nlinear_residual(w) }; + unsafe { + gsl_blas_ddot(f, f, &mut chisq0); + } + + /* solve the system within a maximum of max_iters iterations */ + for i in 0..max_iters { + + let mut rcond: f64 = 0.0; + + unsafe { gsl_multilarge_nlinear_iterate(w); } + unsafe { gsl_multilarge_nlinear_test(xtol, gtol, ftol, &mut *status, w); } + + if *status != 0 { + + if *status == 1 { + println!("Reason for stopping: Small Step Size"); + } else if *status == 2 { + println!("Reason for stopping: Small Gradient"); + } + + break; + } + + unsafe { gsl_multilarge_nlinear_rcond(&mut rcond, w); } + if rcond.is_nan() && i != 0 { + println!("Reason for stopping: Invalid Status"); + *status = -1; + + break; + } + } + + /* compute covariance of best fit parameters */ + unsafe { gsl_multilarge_nlinear_covar(covar, w); } + + /* compute final cost */ + let residuals = unsafe { gsl_multilarge_nlinear_residual(w) }; + unsafe { gsl_blas_ddot(residuals, residuals, &mut chisq); } + let chisq_dof = (chisq / dof).sqrt(); + + let x_gsl_vec: *const gsl_vector = unsafe { gsl_multilarge_nlinear_position(w) }; + for i in 0..params_len { + unsafe { params.wrapping_add(i).write(gsl_vector_get(x_gsl_vec, i)); } + unsafe { parerr.wrapping_add(i).write(chisq_dof * gsl_matrix_get(covar, i, i).sqrt()); } + } + + unsafe { gsl_multilarge_nlinear_free(w); } + unsafe { gsl_matrix_free(covar); } +} diff --git a/src/lib.rs b/src/lib.rs index 66d687e..1ae838f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,127 +1,136 @@ -// -// A rust binding for the GSL library by Guillaume Gomez (guillaume1.gomez@gmail.com) -// - -#![crate_name = "rgsl"] -#![cfg_attr(feature = "dox", feature(doc_cfg))] -#![allow(non_camel_case_types)] -#![allow(non_snake_case)] -#![allow(clippy::too_many_arguments)] -#![allow(clippy::excessive_precision)] -#![allow(clippy::many_single_char_names)] -#![allow(clippy::should_implement_trait)] -#![allow(clippy::type_complexity)] -#![doc = include_str!("../README.md")] - -pub use self::types::*; - -pub use self::elementary::Elementary; -pub use self::pow::Pow; -pub use self::trigonometric::Trigonometric; -pub use self::types::rng; -pub use self::utilities::IOStream; -pub use self::view::View; - -// enums part -pub use self::enums::*; - -mod enums; -mod macros; -mod utilities; -mod view; - -#[doc(hidden)] -pub mod ffi; - -pub mod randist; -pub mod types; - -pub mod airy; -pub mod bessel; -pub mod blas; -pub mod cblas; -pub mod clausen; -pub mod coulomb; -pub mod coupling_coefficients; -pub mod dawson; -pub mod debye; -pub mod dilogarithm; -pub mod eigen; -pub mod elementary; -pub mod elementary_operations; -pub mod elliptic; -pub mod error; -pub mod exponential; -pub mod exponential_integrals; -pub mod fermi_dirac; -pub mod fft; -#[cfg(feature = "v2_5")] -#[cfg_attr(feature = "dox", doc(cfg(feature = "v2_5")))] -pub mod filter; -pub mod fit; -pub mod gamma_beta; -pub mod gegenbauer; -pub mod hypergeometric; -pub mod integration; -pub mod interpolation; -pub mod jacobian_elliptic; -pub mod laguerre; -pub mod lambert_w; -pub mod legendre; -pub mod linear_algebra; -pub mod logarithm; -pub mod minimizer; -pub mod multifit; -#[cfg(feature = "v2_1")] -#[cfg_attr(feature = "dox", doc(cfg(feature = "v2_1")))] -pub mod multilarge; -pub mod multilinear; -pub mod multimin; -pub mod multiroot; -pub mod numerical_differentiation; -pub mod physical_constant; -pub mod polynomials; -pub mod pow; -pub mod power; -pub mod psi; -pub mod roots; -pub mod sort; -pub mod statistics; -pub mod stats; -pub mod synchrotron; -pub mod transport; -pub mod trigonometric; -pub mod util; -pub mod wavelet_transforms; -pub mod zeta; - -/// The maximum x such that gamma(x) is not considered an overflow. -pub static SF_GAMMA_XMAX: f64 = 171.0; -/// The maximum n such that gsl_sf_fact(n) does not give an overflow. -pub static SF_FACT_NMAX: f64 = 170.0; -/// The maximum n such that gsl_sf_doublefact(n) does not give an overflow. -pub static SF_DOUBLEFACT_NMAX: f64 = 297.0; - -pub static SF_MATHIEU_COEFF: u32 = 100; - -pub static DBL_EPSILON: f64 = 2.220_446_049_250_313_1e-16; -pub static SQRT_DBL_EPSILON: f64 = 1.490_116_119_384_765_6e-08; -pub static ROOT3_DBL_EPSILON: f64 = 6.055_454_452_393_342_9e-06; -pub static ROOT4_DBL_EPSILON: f64 = 1.220_703_125_000_000_0e-04; -pub static ROOT5_DBL_EPSILON: f64 = 7.400_959_797_414_050_5e-04; -pub static ROOT6_DBL_EPSILON: f64 = 2.460_783_300_575_925_1e-03; - -pub static DBL_MIN: f64 = 2.225_073_858_507_201_4e-308; -pub static SQRT_DBL_MIN: f64 = 1.491_668_146_240_041_3e-154; -pub static ROOT3_DBL_MIN: f64 = 2.812_644_285_236_299_6e-103; -pub static ROOT4_DBL_MIN: f64 = 1.221_338_669_755_462_0e-77; -pub static ROOT5_DBL_MIN: f64 = 2.947_602_296_969_176_3e-62; -pub static ROOT6_DBL_MIN: f64 = 5.303_436_890_579_821_8e-52; - -pub static DBL_MAX: f64 = std::f64::MAX; //1.7976931348623156e+308; -pub static SQRT_DBL_MAX: f64 = 1.340_780_792_994_259_6e+154; -pub static ROOT3_DBL_MAX: f64 = 5.643_803_094_122_289_7e+102; -pub static ROOT4_DBL_MAX: f64 = 1.157_920_892_373_162_0e+77; -pub static ROOT5_DBL_MAX: f64 = 4.476_546_622_757_270_7e+61; -pub static ROOT6_DBL_MAX: f64 = 2.375_668_978_229_561_2e+51; -pub static LOG_DBL_MAX: f64 = 7.097_827_128_933_839_7e+02; +// +// A rust binding for the GSL library by Guillaume Gomez (guillaume1.gomez@gmail.com) +// + +#![crate_name = "rgsl"] +#![cfg_attr(feature = "dox", feature(doc_cfg))] +#![allow(non_camel_case_types)] +#![allow(non_snake_case)] +#![allow(clippy::too_many_arguments)] +#![allow(clippy::excessive_precision)] +#![allow(clippy::many_single_char_names)] +#![allow(clippy::should_implement_trait)] +#![allow(clippy::type_complexity)] +#![doc = include_str!("../README.md")] + +pub use self::types::*; + +pub use self::elementary::Elementary; +pub use self::pow::Pow; +pub use self::trigonometric::Trigonometric; +pub use self::types::rng; +pub use self::utilities::IOStream; +pub use self::view::View; + +pub use sys::gsl_multifit_nlinear_basic; +pub use sys::gsl_multifit_nlinear_basic_df; +pub use sys::gsl_multifit_nlinear_weighted; +pub use sys::gsl_multifit_nlinear_weighted_df; +pub use sys::gsl_multilarge_nlinear_basic; +pub use sys::gsl_multilarge_nlinear_basic_df; +pub use sys::gsl_multilarge_nlinear_weighted; +pub use sys::gsl_multilarge_nlinear_weighted_df; + +// enums part +pub use self::enums::*; + +mod enums; +mod macros; +mod utilities; +mod view; + +#[doc(hidden)] +pub mod ffi; + +pub mod randist; +pub mod types; + +pub mod airy; +pub mod bessel; +pub mod blas; +pub mod cblas; +pub mod clausen; +pub mod coulomb; +pub mod coupling_coefficients; +pub mod dawson; +pub mod debye; +pub mod dilogarithm; +pub mod eigen; +pub mod elementary; +pub mod elementary_operations; +pub mod elliptic; +pub mod error; +pub mod exponential; +pub mod exponential_integrals; +pub mod fermi_dirac; +pub mod fft; +#[cfg(feature = "v2_5")] +#[cfg_attr(feature = "dox", doc(cfg(feature = "v2_5")))] +pub mod filter; +pub mod fit; +pub mod gamma_beta; +pub mod gegenbauer; +pub mod hypergeometric; +pub mod integration; +pub mod interpolation; +pub mod jacobian_elliptic; +pub mod laguerre; +pub mod lambert_w; +pub mod legendre; +pub mod linear_algebra; +pub mod logarithm; +pub mod minimizer; +pub mod multifit; +#[cfg(feature = "v2_1")] +#[cfg_attr(feature = "dox", doc(cfg(feature = "v2_1")))] +pub mod multilarge; +pub mod multilinear; +pub mod multimin; +pub mod multiroot; +pub mod numerical_differentiation; +pub mod physical_constant; +pub mod polynomials; +pub mod pow; +pub mod power; +pub mod psi; +pub mod roots; +pub mod sort; +pub mod statistics; +pub mod stats; +pub mod synchrotron; +pub mod transport; +pub mod trigonometric; +pub mod util; +pub mod wavelet_transforms; +pub mod zeta; + +/// The maximum x such that gamma(x) is not considered an overflow. +pub static SF_GAMMA_XMAX: f64 = 171.0; +/// The maximum n such that gsl_sf_fact(n) does not give an overflow. +pub static SF_FACT_NMAX: f64 = 170.0; +/// The maximum n such that gsl_sf_doublefact(n) does not give an overflow. +pub static SF_DOUBLEFACT_NMAX: f64 = 297.0; + +pub static SF_MATHIEU_COEFF: u32 = 100; + +pub static DBL_EPSILON: f64 = 2.220_446_049_250_313_1e-16; +pub static SQRT_DBL_EPSILON: f64 = 1.490_116_119_384_765_6e-08; +pub static ROOT3_DBL_EPSILON: f64 = 6.055_454_452_393_342_9e-06; +pub static ROOT4_DBL_EPSILON: f64 = 1.220_703_125_000_000_0e-04; +pub static ROOT5_DBL_EPSILON: f64 = 7.400_959_797_414_050_5e-04; +pub static ROOT6_DBL_EPSILON: f64 = 2.460_783_300_575_925_1e-03; + +pub static DBL_MIN: f64 = 2.225_073_858_507_201_4e-308; +pub static SQRT_DBL_MIN: f64 = 1.491_668_146_240_041_3e-154; +pub static ROOT3_DBL_MIN: f64 = 2.812_644_285_236_299_6e-103; +pub static ROOT4_DBL_MIN: f64 = 1.221_338_669_755_462_0e-77; +pub static ROOT5_DBL_MIN: f64 = 2.947_602_296_969_176_3e-62; +pub static ROOT6_DBL_MIN: f64 = 5.303_436_890_579_821_8e-52; + +pub static DBL_MAX: f64 = std::f64::MAX; //1.7976931348623156e+308; +pub static SQRT_DBL_MAX: f64 = 1.340_780_792_994_259_6e+154; +pub static ROOT3_DBL_MAX: f64 = 5.643_803_094_122_289_7e+102; +pub static ROOT4_DBL_MAX: f64 = 1.157_920_892_373_162_0e+77; +pub static ROOT5_DBL_MAX: f64 = 4.476_546_622_757_270_7e+61; +pub static ROOT6_DBL_MAX: f64 = 2.375_668_978_229_561_2e+51; +pub static LOG_DBL_MAX: f64 = 7.097_827_128_933_839_7e+02;