From 5bdde3ef70b75feeca6ed4f63d7193a4d3951934 Mon Sep 17 00:00:00 2001 From: Stefan Kroboth Date: Sat, 2 Mar 2024 12:10:06 +0100 Subject: [PATCH] Be even more generic --- crates/finitediff/src/lib.rs | 28 ++-- crates/finitediff/src/ndarray_m/diff.rs | 47 ++++--- crates/finitediff/src/ndarray_m/jacobian.rs | 139 +++++++++++++------- crates/finitediff/src/vec/diff.rs | 8 +- crates/finitediff/src/vec/hessian.rs | 4 - crates/finitediff/src/vec/jacobian.rs | 3 +- 6 files changed, 136 insertions(+), 93 deletions(-) diff --git a/crates/finitediff/src/lib.rs b/crates/finitediff/src/lib.rs index 32248b71c..7d78f526e 100644 --- a/crates/finitediff/src/lib.rs +++ b/crates/finitediff/src/lib.rs @@ -651,19 +651,19 @@ where type OperatorOutput = ndarray::Array1; fn forward_diff(&self, f: &dyn Fn(&Self) -> f64) -> Self { - forward_diff_ndarray_f64(self, f) + forward_diff_ndarray(self, f) } fn central_diff(&self, f: &dyn Fn(&ndarray::Array1) -> f64) -> Self { - central_diff_ndarray_f64(self, f) + central_diff_ndarray(self, f) } fn forward_jacobian(&self, fs: &dyn Fn(&Self) -> Self::OperatorOutput) -> Self::Jacobian { - forward_jacobian_ndarray_f64(self, fs) + forward_jacobian_ndarray(self, fs) } fn central_jacobian(&self, fs: &dyn Fn(&Self) -> Self::OperatorOutput) -> Self::Jacobian { - central_jacobian_ndarray_f64(self, fs) + central_jacobian_ndarray(self, fs) } fn forward_jacobian_vec_prod( @@ -671,7 +671,7 @@ where fs: &dyn Fn(&Self) -> Self::OperatorOutput, p: &Self, ) -> Self { - forward_jacobian_vec_prod_ndarray_f64(self, fs, p) + forward_jacobian_vec_prod_ndarray(self, fs, p) } fn central_jacobian_vec_prod( @@ -679,7 +679,7 @@ where fs: &dyn Fn(&Self) -> Self::OperatorOutput, p: &Self, ) -> Self { - central_jacobian_vec_prod_ndarray_f64(self, fs, p) + central_jacobian_vec_prod_ndarray(self, fs, p) } fn forward_jacobian_pert( @@ -687,7 +687,7 @@ where fs: &dyn Fn(&Self) -> Self::OperatorOutput, pert: &PerturbationVectors, ) -> Self::Jacobian { - forward_jacobian_pert_ndarray_f64(self, fs, pert) + forward_jacobian_pert_ndarray(self, fs, pert) } fn central_jacobian_pert( @@ -695,7 +695,7 @@ where fs: &dyn Fn(&Self) -> Self::OperatorOutput, pert: &PerturbationVectors, ) -> Self::Jacobian { - central_jacobian_pert_ndarray_f64(self, fs, pert) + central_jacobian_pert_ndarray(self, fs, pert) } fn forward_hessian(&self, g: &dyn Fn(&Self) -> Self::OperatorOutput) -> Self::Jacobian { @@ -825,7 +825,7 @@ mod tests_vec { #[test] fn test_forward_diff_vec_f64_trait() { let grad = x1().forward_diff(&f1); - let res = vec![1.0f64, 2.0]; + let res = [1.0f64, 2.0]; for i in 0..2 { assert!((res[i] - grad[i]).abs() < COMP_ACC) @@ -833,7 +833,7 @@ mod tests_vec { let p = vec![1.0f64, 2.0f64]; let grad = p.forward_diff(&f1); - let res = vec![1.0f64, 4.0]; + let res = [1.0f64, 4.0]; for i in 0..2 { assert!((res[i] - grad[i]).abs() < COMP_ACC) @@ -843,7 +843,7 @@ mod tests_vec { #[test] fn test_central_diff_vec_f64_trait() { let grad = x1().central_diff(&f1); - let res = vec![1.0f64, 2.0]; + let res = [1.0f64, 2.0]; for i in 0..2 { assert!((res[i] - grad[i]).abs() < COMP_ACC) @@ -851,7 +851,7 @@ mod tests_vec { let p = vec![1.0f64, 2.0f64]; let grad = p.central_diff(&f1); - let res = vec![1.0f64, 4.0]; + let res = [1.0f64, 4.0]; for i in 0..2 { assert!((res[i] - grad[i]).abs() < COMP_ACC) @@ -958,7 +958,7 @@ mod tests_vec { #[test] fn test_forward_hessian_vec_prod_vec_f64_trait() { let hessian = x3().forward_hessian_vec_prod(&g, &p2()); - let res = vec![0.0, 6.0, 10.0, 18.0]; + let res = [0.0, 6.0, 10.0, 18.0]; // println!("hessian:\n{:#?}", hessian); // println!("diff:\n{:#?}", diff); for i in 0..4 { @@ -969,7 +969,7 @@ mod tests_vec { #[test] fn test_central_hessian_vec_prod_vec_f64_trait() { let hessian = x3().central_hessian_vec_prod(&g, &p2()); - let res = vec![0.0, 6.0, 10.0, 18.0]; + let res = [0.0, 6.0, 10.0, 18.0]; // println!("hessian:\n{:#?}", hessian); // println!("diff:\n{:#?}", diff); for i in 0..4 { diff --git a/crates/finitediff/src/ndarray_m/diff.rs b/crates/finitediff/src/ndarray_m/diff.rs index 0f4412731..50b4909c1 100644 --- a/crates/finitediff/src/ndarray_m/diff.rs +++ b/crates/finitediff/src/ndarray_m/diff.rs @@ -5,33 +5,44 @@ // http://opensource.org/licenses/MIT>, at your option. This file may not be // copied, modified, or distributed except according to those terms. +use num::{Float, FromPrimitive}; + use crate::utils::*; -use crate::EPS_F64; -pub fn forward_diff_ndarray_f64( - x: &ndarray::Array1, - f: &dyn Fn(&ndarray::Array1) -> f64, -) -> ndarray::Array1 { +pub fn forward_diff_ndarray( + x: &ndarray::Array1, + f: &dyn Fn(&ndarray::Array1) -> F, +) -> ndarray::Array1 +where + F: Float, +{ + let eps_sqrt = F::epsilon().sqrt(); + let fx = (f)(x); let mut xt = x.clone(); (0..x.len()) .map(|i| { - let fx1 = mod_and_calc(&mut xt, f, i, EPS_F64.sqrt()); - (fx1 - fx) / (EPS_F64.sqrt()) + let fx1 = mod_and_calc(&mut xt, f, i, eps_sqrt); + (fx1 - fx) / eps_sqrt }) .collect() } -pub fn central_diff_ndarray_f64( - x: &ndarray::Array1, - f: &dyn Fn(&ndarray::Array1) -> f64, -) -> ndarray::Array1 { +pub fn central_diff_ndarray( + x: &ndarray::Array1, + f: &dyn Fn(&ndarray::Array1) -> F, +) -> ndarray::Array1 +where + F: Float + FromPrimitive, +{ + let eps_sqrt = F::epsilon().sqrt(); + let mut xt = x.clone(); (0..x.len()) .map(|i| { - let fx1 = mod_and_calc(&mut xt, f, i, EPS_F64.sqrt()); - let fx2 = mod_and_calc(&mut xt, f, i, -EPS_F64.sqrt()); - (fx1 - fx2) / (2.0 * EPS_F64.sqrt()) + let fx1 = mod_and_calc(&mut xt, f, i, eps_sqrt); + let fx2 = mod_and_calc(&mut xt, f, i, -eps_sqrt); + (fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_sqrt) }) .collect() } @@ -51,7 +62,7 @@ mod tests { fn test_forward_diff_ndarray_f64() { let p = ndarray::Array1::from(vec![1.0f64, 1.0f64]); - let grad = forward_diff_ndarray_f64(&p, &f); + let grad = forward_diff_ndarray(&p, &f); let res = vec![1.0f64, 2.0]; (0..2) @@ -59,7 +70,7 @@ mod tests { .count(); let p = ndarray::Array1::from(vec![1.0f64, 2.0f64]); - let grad = forward_diff_ndarray_f64(&p, &f); + let grad = forward_diff_ndarray(&p, &f); let res = vec![1.0f64, 4.0]; (0..2) @@ -70,7 +81,7 @@ mod tests { fn test_central_diff_ndarray_f64() { let p = ndarray::Array1::from(vec![1.0f64, 1.0f64]); - let grad = central_diff_ndarray_f64(&p, &f); + let grad = central_diff_ndarray(&p, &f); let res = vec![1.0f64, 2.0]; (0..2) @@ -78,7 +89,7 @@ mod tests { .count(); let p = ndarray::Array1::from(vec![1.0f64, 2.0f64]); - let grad = central_diff_ndarray_f64(&p, &f); + let grad = central_diff_ndarray(&p, &f); let res = vec![1.0f64, 4.0]; (0..2) diff --git a/crates/finitediff/src/ndarray_m/jacobian.rs b/crates/finitediff/src/ndarray_m/jacobian.rs index 966627ffa..db2026c41 100644 --- a/crates/finitediff/src/ndarray_m/jacobian.rs +++ b/crates/finitediff/src/ndarray_m/jacobian.rs @@ -5,38 +5,54 @@ // http://opensource.org/licenses/MIT>, at your option. This file may not be // copied, modified, or distributed except according to those terms. -use ndarray::Array2; +use std::ops::AddAssign; -use crate::pert::*; -use crate::utils::*; -use crate::EPS_F64; +use ndarray::Array2; +use ndarray::ScalarOperand; +use num::Float; +use num::FromPrimitive; + +use crate::pert::PerturbationVectors; +use crate::utils::mod_and_calc; + +pub fn forward_jacobian_ndarray( + x: &ndarray::Array1, + fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, +) -> ndarray::Array2 +where + F: Float, +{ + let eps_sqrt = F::epsilon().sqrt(); -pub fn forward_jacobian_ndarray_f64( - x: &ndarray::Array1, - fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, -) -> ndarray::Array2 { let fx = (fs)(x); let mut xt = x.clone(); let rn = fx.len(); let n = x.len(); let mut out = Array2::zeros((n, rn)); for i in 0..n { - let fx1 = mod_and_calc(&mut xt, fs, i, EPS_F64.sqrt()); + let fx1 = mod_and_calc(&mut xt, fs, i, eps_sqrt); for j in 0..rn { - out[(i, j)] = (fx1[j] - fx[j]) / EPS_F64.sqrt(); + out[(i, j)] = (fx1[j] - fx[j]) / eps_sqrt; } } out } -pub fn central_jacobian_ndarray_f64( - x: &ndarray::Array1, - fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, -) -> ndarray::Array2 { +pub fn central_jacobian_ndarray( + x: &ndarray::Array1, + fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, +) -> ndarray::Array2 +where + F: Float + FromPrimitive, +{ + let eps_sqrt = F::epsilon().sqrt(); + let mut xt = x.clone(); // TODO: get rid of this! fx is only needed to calculate rn in order to be able to allocate the // array for the jacobian. + // A question, many years after this comment was written: Shouldn't that be the same length as + // x anyways? let fx = (fs)(x); let rn = fx.len(); @@ -44,49 +60,63 @@ pub fn central_jacobian_ndarray_f64( let mut out = Array2::zeros((n, rn)); for i in 0..n { - let fx1 = mod_and_calc(&mut xt, fs, i, EPS_F64.sqrt()); - let fx2 = mod_and_calc(&mut xt, fs, i, -EPS_F64.sqrt()); + let fx1 = mod_and_calc(&mut xt, fs, i, eps_sqrt); + let fx2 = mod_and_calc(&mut xt, fs, i, -eps_sqrt); for j in 0..rn { - out[(i, j)] = (fx1[j] - fx2[j]) / (2.0 * EPS_F64.sqrt()); + out[(i, j)] = (fx1[j] - fx2[j]) / (F::from_f64(2.0).unwrap() * eps_sqrt); } } out } -pub fn forward_jacobian_vec_prod_ndarray_f64( - x: &ndarray::Array1, - fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, - p: &ndarray::Array1, -) -> ndarray::Array1 { +pub fn forward_jacobian_vec_prod_ndarray( + x: &ndarray::Array1, + fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, + p: &ndarray::Array1, +) -> ndarray::Array1 +where + F: Float + ScalarOperand, +{ + let eps_sqrt = F::epsilon().sqrt(); let fx = (fs)(x); - let x1 = x + &p.mapv(|pi| EPS_F64.sqrt() * pi); + // TODO: Why not use mod_and_calc? + let x1 = x + &p.mapv(|pi| eps_sqrt * pi); let fx1 = (fs)(&x1); - (fx1 - fx) / EPS_F64.sqrt() + (fx1 - fx) / eps_sqrt } -pub fn central_jacobian_vec_prod_ndarray_f64( - x: &ndarray::Array1, - fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, - p: &ndarray::Array1, -) -> ndarray::Array1 { - let x1 = x + &p.mapv(|pi| EPS_F64.sqrt() * pi); - let x2 = x + &p.mapv(|pi| -EPS_F64.sqrt() * pi); +pub fn central_jacobian_vec_prod_ndarray( + x: &ndarray::Array1, + fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, + p: &ndarray::Array1, +) -> ndarray::Array1 +where + F: Float + FromPrimitive + ScalarOperand, +{ + let eps_sqrt = F::epsilon().sqrt(); + let x1 = x + &p.mapv(|pi| eps_sqrt * pi); + let x2 = x + &p.mapv(|pi| -eps_sqrt * pi); let fx1 = (fs)(&x1); let fx2 = (fs)(&x2); - (fx1 - fx2) / (2.0 * EPS_F64.sqrt()) + (fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_sqrt) } -pub fn forward_jacobian_pert_ndarray_f64( - x: &ndarray::Array1, - fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, +pub fn forward_jacobian_pert_ndarray( + x: &ndarray::Array1, + fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, pert: &PerturbationVectors, -) -> ndarray::Array2 { +) -> ndarray::Array2 +where + F: Float + AddAssign, +{ + let eps_sqrt = F::epsilon().sqrt(); + let fx = (fs)(x); let mut xt = x.clone(); let mut out = ndarray::Array2::zeros((fx.len(), x.len())); for pert_item in pert.iter() { for j in pert_item.x_idx.iter() { - xt[*j] += EPS_F64.sqrt(); + xt[*j] += eps_sqrt; } let fx1 = (fs)(&xt); @@ -97,29 +127,34 @@ pub fn forward_jacobian_pert_ndarray_f64( for (k, x_idx) in pert_item.x_idx.iter().enumerate() { for j in pert_item.r_idx[k].iter() { - out[(*x_idx, *j)] = (fx1[*j] - fx[*j]) / EPS_F64.sqrt(); + out[(*x_idx, *j)] = (fx1[*j] - fx[*j]) / eps_sqrt; } } } out } -pub fn central_jacobian_pert_ndarray_f64( - x: &ndarray::Array1, - fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, +pub fn central_jacobian_pert_ndarray( + x: &ndarray::Array1, + fs: &dyn Fn(&ndarray::Array1) -> ndarray::Array1, pert: &PerturbationVectors, -) -> ndarray::Array2 { +) -> ndarray::Array2 +where + F: Float + FromPrimitive + AddAssign, +{ + let eps_sqrt = F::epsilon().sqrt(); + let mut out = ndarray::Array2::zeros((1, 1)); let mut xt = x.clone(); for (i, pert_item) in pert.iter().enumerate() { for j in pert_item.x_idx.iter() { - xt[*j] += EPS_F64.sqrt(); + xt[*j] += eps_sqrt; } let fx1 = (fs)(&xt); for j in pert_item.x_idx.iter() { - xt[*j] = x[*j] - EPS_F64.sqrt(); + xt[*j] = x[*j] - eps_sqrt; } let fx2 = (fs)(&xt); @@ -134,7 +169,7 @@ pub fn central_jacobian_pert_ndarray_f64( for (k, x_idx) in pert_item.x_idx.iter().enumerate() { for j in pert_item.r_idx[k].iter() { - out[(*x_idx, *j)] = (fx1[*j] - fx2[*j]) / (2.0 * EPS_F64.sqrt()); + out[(*x_idx, *j)] = (fx1[*j] - fx2[*j]) / (F::from_f64(2.0).unwrap() * eps_sqrt); } } } @@ -143,6 +178,8 @@ pub fn central_jacobian_pert_ndarray_f64( #[cfg(test)] mod tests { + use crate::PerturbationVector; + use super::*; use ndarray; use ndarray::{array, Array1}; @@ -199,7 +236,7 @@ mod tests { #[test] fn test_forward_jacobian_ndarray_f64() { - let jacobian = forward_jacobian_ndarray_f64(&x(), &f); + let jacobian = forward_jacobian_ndarray(&x(), &f); let res = res1(); // println!("{:?}", jacobian); for i in 0..6 { @@ -211,7 +248,7 @@ mod tests { #[test] fn test_central_jacobian_ndarray_f64() { - let jacobian = central_jacobian_ndarray_f64(&x(), &f); + let jacobian = central_jacobian_ndarray(&x(), &f); let res = res1(); // println!("{:?}", jacobian); for i in 0..6 { @@ -223,7 +260,7 @@ mod tests { #[test] fn test_forward_jacobian_vec_prod_ndarray_f64() { - let jacobian = forward_jacobian_vec_prod_ndarray_f64(&x(), &f, &p()); + let jacobian = forward_jacobian_vec_prod_ndarray(&x(), &f, &p()); let res = res2(); // println!("{:?}", jacobian); // the accuracy for this is pretty bad!! @@ -234,7 +271,7 @@ mod tests { #[test] fn test_central_jacobian_vec_prod_ndarray_f64() { - let jacobian = central_jacobian_vec_prod_ndarray_f64(&x(), &f, &p()); + let jacobian = central_jacobian_vec_prod_ndarray(&x(), &f, &p()); let res = res2(); // println!("{:?}", jacobian); for i in 0..6 { @@ -244,7 +281,7 @@ mod tests { #[test] fn test_forward_jacobian_pert_ndarray_f64() { - let jacobian = forward_jacobian_pert_ndarray_f64(&x(), &f, &pert()); + let jacobian = forward_jacobian_pert_ndarray(&x(), &f, &pert()); let res = res1(); // println!("jacobian:\n{:?}", jacobian); // println!("res:\n{:?}", res); @@ -257,7 +294,7 @@ mod tests { #[test] fn test_central_jacobian_pert_ndarray_f64() { - let jacobian = central_jacobian_pert_ndarray_f64(&x(), &f, &pert()); + let jacobian = central_jacobian_pert_ndarray(&x(), &f, &pert()); let res = res1(); // println!("jacobian:\n{:?}", jacobian); // println!("res:\n{:?}", res); diff --git a/crates/finitediff/src/vec/diff.rs b/crates/finitediff/src/vec/diff.rs index a6c2ac1ab..f32c381c0 100644 --- a/crates/finitediff/src/vec/diff.rs +++ b/crates/finitediff/src/vec/diff.rs @@ -54,7 +54,7 @@ mod tests { fn test_forward_diff_vec_f64() { let p = vec![1.0f64, 1.0f64]; let grad = forward_diff_vec(&p, &f); - let res = vec![1.0f64, 2.0]; + let res = [1.0f64, 2.0]; (0..2) .map(|i| assert!((res[i] - grad[i]).abs() < COMP_ACC)) @@ -62,7 +62,7 @@ mod tests { let p = vec![1.0f64, 2.0f64]; let grad = forward_diff_vec(&p, &f); - let res = vec![1.0f64, 4.0]; + let res = [1.0f64, 4.0]; (0..2) .map(|i| assert!((res[i] - grad[i]).abs() < COMP_ACC)) @@ -73,7 +73,7 @@ mod tests { fn test_central_diff_vec_f64() { let p = vec![1.0f64, 1.0f64]; let grad = central_diff_vec(&p, &f); - let res = vec![1.0f64, 2.0]; + let res = [1.0f64, 2.0]; (0..2) .map(|i| assert!((res[i] - grad[i]).abs() < COMP_ACC)) @@ -81,7 +81,7 @@ mod tests { let p = vec![1.0f64, 2.0f64]; let grad = central_diff_vec(&p, &f); - let res = vec![1.0f64, 4.0]; + let res = [1.0f64, 4.0]; (0..2) .map(|i| assert!((res[i] - grad[i]).abs() < COMP_ACC)) diff --git a/crates/finitediff/src/vec/hessian.rs b/crates/finitediff/src/vec/hessian.rs index cb4fcfd00..ae00aa459 100644 --- a/crates/finitediff/src/vec/hessian.rs +++ b/crates/finitediff/src/vec/hessian.rs @@ -10,10 +10,6 @@ use std::ops::AddAssign; use num::{Float, FromPrimitive}; use crate::utils::{mod_and_calc, restore_symmetry_vec, KV}; -use crate::EPS_F64; - -/// I wish this wasn't necessary! -const EPS_F64_NOGRAD: f64 = EPS_F64 * 2.0; pub fn forward_hessian_vec(x: &Vec, grad: &dyn Fn(&Vec) -> Vec) -> Vec> where diff --git a/crates/finitediff/src/vec/jacobian.rs b/crates/finitediff/src/vec/jacobian.rs index 0331561b4..1c0957436 100644 --- a/crates/finitediff/src/vec/jacobian.rs +++ b/crates/finitediff/src/vec/jacobian.rs @@ -7,8 +7,7 @@ use std::ops::AddAssign; -use num::Float; -use num::FromPrimitive; +use num::{Float, FromPrimitive}; use crate::pert::PerturbationVectors; use crate::utils::mod_and_calc;