// assume we want to differentiate f(x) = x^3
let func = | arg: f64 | -> f64 
{ 
    return arg*arg*arg;
};

let point = 2.0; //the point at which we want to differentiate

use multicalc::numerical_derivative::derivator::*;
use multicalc::numerical_derivative::finite_difference::*;
 
let derivator = SingleVariableSolver::default();
//alternatively, you can also create the derivator with custom parameters using SingleVariableSolver::from_parameters()

let val = derivator.get(1, &my_func, point).unwrap(); //single derivative
assert!(f64::abs(val - 12.0) < 1e-7);
let val = derivator.get(2, &my_func, point).unwrap(); //double derivative
assert!(f64::abs(val - 12.0) < 1e-5);
let val = derivator.get(3, &my_func, point).unwrap(); //triple derivative
assert!(f64::abs(val - 6.0) < 1e-3);

//for single and double derivatives, you can also just use the convenience wrappers
let val = derivator.get_single(&my_func, point).unwrap(); //single derivative
assert!(f64::abs(val - 12.0) < 1e-7);
let val = derivator.get_double(&my_func, point).unwrap(); //double derivative
assert!(f64::abs(val - 12.0) < 1e-5);

//function is f(x,y,z) = y*sin(x) + x*cos(y) + x*y*e^z
let func = | args: &[f64; 3] | -> f64 
{ 
    return args[1]*args[0].sin() + args[0]*args[1].cos() + args[0]*args[1]*args[2].exp();
};

let point = [1.0, 2.0, 3.0]; //the point of interest

use multicalc::numerical_derivative::derivator::*;
use multicalc::numerical_derivative::finite_difference::*;

let derivator = MultiVariableSolver::default();
//alternatively, you can also create the derivator with custom parameters using MultiVariableSolver::from_parameters()

let idx: [usize; 2] = [0, 1]; //mixed partial double derivate d(df/dx)/dy
let val = derivator.get(2, &my_func, &idx, &point).unwrap();
let expected_value = f64::cos(1.0) - f64::sin(2.0) + f64::exp(3.0);
assert!(f64::abs(val - expected_value) < 0.001);
 
let idx: [usize; 2] = [1, 1]; //partial double derivative d(df/dy)/dy 
let val = derivator.get(2, &my_func, &idx, &point).unwrap();
let expected_value = -1.0*f64::cos(2.0);
assert!(f64::abs(val - expected_value) < 0.0001);

//integrate 2*x . the function would be:
let my_func = | arg: f64 | -> f64 
{ 
    return 2.0*arg;
};

use multicalc::numerical_integration::integrator::*;
use multicalc::numerical_integration::iterative_integration;

let integrator = iterative_integration::SingleVariableSolver::default();  
//alternatively, create a custom integrator using SingleVariableSolver::from_parameters()

//instead of using iterative algorithms, you can also use gaussian quadratures instead:
//let integrator = gaussian_integration::SingleVariableSolver::default();  

let integration_limit = [[0.0, 2.0]; 1]; //desired integration limit 
let val = integrator.get(1, &my_func, &integration_limit).unwrap(); //single integration
assert!(f64::abs(val - 4.0) < 1e-6);

let integration_limit = [[0.0, 2.0], [-1.0, 1.0]]; //desired integration limits
let val = integrator.get(2, &my_func, &integration_limit).unwrap(); //double integration
assert!(f64::abs(val - 8.0) < 1e-6);
///
let integration_limit = [[0.0, 2.0], [0.0, 2.0], [0.0, 2.0]]; //desired integration limits
let val = integrator.get(3, &my_func, &integration_limit).unwrap(); //triple integration
assert!(f64::abs(val - 16.0) < 1e-6);


//you can also just use convenience wrappers for single and double integration
let integration_limit = [0.0, 2.0];
let val = integrator.get_single(&my_func, &integration_limit).unwrap(); //single integration
assert!(f64::abs(val - 4.0) < 1e-6);

let integration_limit = [[0.0, 2.0], [-1.0, 1.0]];
let val = integrator.get_double(&my_func, &integration_limit).unwrap(); //double integration
assert!(f64::abs(val - 8.0) < 1e-6);

//equation is 2.0*x + y*z
let func = | args: &[f64; 3] | -> f64 
{ 
    return 2.0*args[0] + args[1]*args[2];
};

//for multivariable integration to work, we must know the final value of all variables
//for variables that are constant, they will have their constant values
//for variables that are being integrated, their values will be the final upper limit of integration
//in this example, we are keeping 'z' constant at 3.0, and integrating x and y to upper limits 1.0 and 2.0 respectively
let point = [1.0, 2.0, 3.0];

let iterator = iterative_integration::MultiVariableSolver::default();
//you can also create a custom integrator using MultiVariableSolver::from_parameters()

//instead of using iterative algorithms, you can also use gaussian quadratures instead:
//let integrator = gaussian_integration::MultiVariableSolver::default();  

//integration for x, known to be x*x + x*y*z, expect a value of ~7.00
let integration_limit = [0.0, 1.0];
let val = iterator.get_single_partial(&func, 0, &integration_limit, &point).unwrap();
//alternative syntax: iterator.get(1, [0], &func, &[integration_limit; 1], &point);
assert!(f64::abs(val - 7.0) < 1e-7);


//integration for y, known to be 2.0*x*y + y*y*z/2.0, expect a value of ~10.00 
let integration_limit = [0.0, 2.0];
let val = iterator.get_single_partial(&func, 1, &integration_limit, &point).unwrap();
//alternative syntax: iterator.get(1, [1], &func, &[integration_limit; 1], &point);
assert!(f64::abs(val - 10.0) < 1e-7);


//double integration for first x then y, expect a value of ~8.0
let integration_limit = [[0.0, 1.0], [0.0, 2.0]];
let val = integrator.get_double_partial(&func, [0, 1], &integration_limit, &point).unwrap();
//alternative syntax: iterator.get(2, [0, 1], &func, &integration_limit, &point);
assert!(f64::abs(val - 8.0) < 1e-7);

//triple integration for x, expect a value of ~7.0
let integration_limit = [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]];
let val = integrator.get(3, [0, 0, 0], &func, &integration_limit, &point).unwrap();
assert!(f64::abs(val - 7.0) < 1e-7);

//function is x*y*z
let func1 = | args: &[f64; 3] | -> f64 
{ 
    return args[0]*args[1]*args[2];
};

//function is x^2 + y^2
let func2 = | args: &[f64; 3] | -> f64 
{ 
    return args[0].powf(2.0) + args[1].powf(2.0);
};

let function_matrix: [&dyn Fn(&[f64; 3]) -> f64; 2] = [&func1, &func2];

let points = [1.0, 2.0, 3.0]; //the point around which we want the jacobian matrix

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let jacobian = Jacobian::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = jacobian.get(&function_matrix, &points).unwrap();

let expected_result = [[6.0, 3.0, 2.0], [2.0, 4.0, 0.0]];
for i in 0..function_vector.len()
{
    for j in 0..points.len()
    {
        //numerical error less than 1e-6
        assert!(f64::abs(result[i][j] - expected_result[i][j]) < 1e-6);
    }
}

//function is y*sin(x) + 2*x*e^y
let func = | args: &[f64; 2] | -> f64 
{ 
    return args[1]*args[0].sin() + 2.0*args[0]*args[1].exp();
};

let points = [1.0, 2.0]; //the point around which we want the hessian matrix

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let hessian = Hessian::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = hessian.get(&func, &points).unwrap();

assert!(result.len() == points.len()); //number of rows
assert!(result[0].len() == points.len()); //number of columns

let expected_result = [[-2.0*f64::sin(1.0), f64::cos(1.0) + 2.0*f64::exp(2.0)], 
                                        [f64::cos(1.0) + 2.0*f64::exp(2.0), 2.0*f64::exp(2.0)]];

for i in 0..points.len()
{
    for j in 0..points.len()
    {
        assert!(f64::abs(result[i][j] - expected_result[i][j]) < 1e-5);
    }
}

//function is x + y^2 + z^3, which we want to linearize
let function_to_approximate = | args: &[f64; 3] | -> f64 
{ 
    return args[0] + args[1].powf(2.0) + args[2].powf(3.0);
};

let point = [1.0, 2.0, 3.0]; //the point we want to linearize around

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let approximator = LinearApproximator::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = approximator.get(&function_to_approximate, &point).unwrap();
assert!(f64::abs(function_to_approximate(&point) - result.get_prediction_value(&point)) < 1e-9);

//function is e^(x/2) + sin(y) + 2.0*z
let function_to_approximate = | args: &[f64; 3] | -> f64 
{ 
    return f64::exp(args[0]/2.0) + f64::sin(args[1]) + 2.0*args[2];
};

let point = [0.0, 3.14/2.0, 10.0]; //the point we want to approximate around

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let approximator = QuadraticApproximator::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = approximator.get(&function_to_approximate, &point).unwrap();

assert!(f64::abs(function_to_approximate(&point) - result.get_prediction_value(&point)) < 1e-9);

//vector field is (y, -x). On a 2D plane this would like a tornado rotating counter-clockwise
//curve is a unit circle, defined by (Cos(t), Sin(t))
//limit t goes from 0->2*pi

let vector_field_matrix: [&dyn Fn(&f64, &f64) -> f64; 2] = [&(|_:&f64, y:&f64|-> f64 { *y }), &(|x:&f64, _:&f64|-> f64 { -x })];

let transformation_matrix: [&dyn Fn(&f64) -> f64; 2] = [&(|t:&f64|->f64 { t.cos() }), &(|t:&f64|->f64 { t.sin() })];

let integration_limit = [0.0, 6.28];

//line integral of a unit circle curve on our vector field from 0 to 2*pi, expect an answer of -2.0*pi
let val = line_integral::get_2d(&vector_field_matrix, &transformation_matrix, &integration_limit).unwrap();
assert!(f64::abs(val + 6.28) < 0.01);

//flux integral of a unit circle curve on our vector field from 0 to 2*pi, expect an answer of 0.0
let val = flux_integral::get_2d(&vector_field_matrix, &transformation_matrix, &integration_limit).unwrap();
assert!(f64::abs(val - 0.0) < 0.01);

//vector field is (2*x*y, 3*cos(y))
//x-component
let vf_x = | args: &[f64; 2] | -> f64 
{ 
    return 2.0*args[0]*args[1];
};

//y-component
let vf_y = | args: &[f64; 2] | -> f64 
{ 
    return 3.0*args[1].cos()
};

let vector_field_matrix: [&dyn Fn(&[f64; 2]) -> f64; 2] = [&vf_x, &vf_y];

let point = [1.0, 3.14]; //the point of interest

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;
let derivator = MultiVariableSolver::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

//curl is known to be -2*x, expect and answer of -2.0
let val = curl::get_2d(derivator, &vector_field_matrix, &point).unwrap();
assert!(f64::abs(val + 2.0) < 0.000001); //numerical error less than 1e-6

//divergence is known to be 2*y - 3*sin(y), expect and answer of 6.27
let val = divergence::get_2d(derivator, &vector_field_matrix, &point).unwrap();
assert!(f64::abs(val - 6.27) < 0.01);

//function is x*y*z
let func1 = | args: &[f64; 3] | -> f64 
{ 
    return args[0]*args[1]*args[2];
};

//function is x^2 + y^2
let func2 = | args: &[f64; 3] | -> f64 
{ 
    return args[0].powf(2.0) + args[1].powf(2.0);
};

let function_matrix: Vec<Box<dyn Fn(&[f64; 3]) -> f64>> = std::vec![Box::new(func1), Box::new(func2)];

let points = [1.0, 2.0, 3.0]; //the point around which we want the jacobian matrix

let jacobian = Jacobian::<MultiVariableSolver>::default();

let result: Vec<Vec<f64>> = jacobian.get_on_heap(&function_matrix, &points).unwrap();

assert!(result.len() == function_matrix.len()); //number of rows
assert!(result[0].len() == points.len()); //number of columns

let expected_result = [[6.0, 3.0, 2.0], [2.0, 4.0, 0.0]];

for i in 0..function_matrix.len()
{
    for j in 0..points.len()
    {
        assert!(f64::abs(result[i][j] - expected_result[i][j]) < 0.000001);
    }
}