Function jacobian 

pub fn jacobian<F, G>(
    f: &G,
    x: &Array1<F>,
    epsilon: F,
) -> LinalgResult<Array2<F>>
where
    F: Float + NumAssign + Sum + ScalarOperand,
    G: Fn(&Array1<F>) -> Array1<F>,

Calculate the Jacobian matrix for a function that maps from R^n to R^m

Computes a numerical approximation of the Jacobian matrix for a function that takes an n-dimensional input and produces an m-dimensional output.

Arguments

• f - A function that maps from R^n to R^m
• x - The point at which to evaluate the Jacobian
• epsilon - The step size for the finite difference approximation

Returns

• The Jacobian matrix of shape (m, n)

Examples

use scirs2_core::ndarray::{array, Array1};
use scirs2_linalg::gradient::jacobian;

// Define a simple function R^2 -> R^3
// f(x,y) = [x^2 + y, 2*x + y^2, x*y]
let f = |v: &Array1<f64>| -> Array1<f64> {
    let x = v[0];
    let y = v[1];
    array![x*x + y, 2.0*x + y*y, x*y]
};

let x = array![2.0, 3.0];  // Point at which to evaluate the Jacobian
let epsilon = 1e-5;

let jac = jacobian(&f, &x, epsilon).unwrap();

// Analytical Jacobian at (2,3) is:
// [2x, 1]     [4, 1]
// [2, 2y]  =  [2, 6]
// [y, x]      [3, 2]

assert!((jac[[0, 0]] - 4.0).abs() < 1e-4);
assert!((jac[[0, 1]] - 1.0).abs() < 1e-4);
assert!((jac[[1, 0]] - 2.0).abs() < 1e-4);
assert!((jac[[1, 1]] - 6.0).abs() < 1e-4);
assert!((jac[[2, 0]] - 3.0).abs() < 1e-4);
assert!((jac[[2, 1]] - 2.0).abs() < 1e-4);