feature heap

Author: kmolan

Cargo.toml

@@ -19,4 +19,8 @@ num-complex = "0.4.6"
 [dev-dependencies]
 rand = "0.8"
 
+[features]
+default = []
+heap = []
+
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

README.md

@@ -43,6 +43,7 @@ Rust scientific computing for single and multi-variable calculus
 - [13. Line and Flux integrals](#13-line-and-flux-integrals)
 - [14. Curl and Divergence](#14-curl-and-divergence)
 - [15. Error Handling](#15-error-handling)
+- [16. Experimental](#16-experimental)
 
 ## 1. Single total derivatives
 ```rust
@@ -308,6 +309,40 @@ Wherever possible, "safe" versions of functions are provided that fill in the de
 However, that is not always possible either because no default argument can be assumed, or for functions that deliberately give users the freedom to tweak the parameters.
 In such cases, a `Result<T, ErrorCode>` object is returned instead, where all possible `ErrorCode`s can be viewed at [error_codes](./src/utils/error_codes.rs).
 
+##  16. Experimental
+Enable feature "heap" to access `std::Vec` based methods in certain modules. Currently this is only supported for the _Jacobian_ module via `get_on_heap()` and `get_on_heap_custom()` methods. The output is a dynamically allocated `Vec<Vec<T>>`. This is to support large datasets that might otherwise get a stack overflow with static arrays. Future plans might include adding such support for the approximation module.
+```rust
+//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 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);
+    }
+}
+```
 
 ## Contributions Guide
 See [CONTRIBUTIONS.md](./CONTRIBUTIONS.md)

src/lib.rs

@@ -3,6 +3,9 @@
 
 #![no_std]
 
+#[cfg(feature = "heap")]
+extern crate std;
+
 pub mod utils;
 pub mod numerical_integration;
 pub mod numerical_derivative;

src/numerical_derivative/jacobian.rs

@@ -3,6 +3,10 @@ use crate::numerical_derivative::mode as mode;
 use crate::utils::error_codes::ErrorCode;
 use num_complex::ComplexFloat;
 
+#[cfg(feature = "heap")]
+use std::{boxed::Box, vec::Vec};
+
+
 /// Returns the jacobian matrix for a given vector of functions
 /// Can handle multivariable functions of any order or complexity
 /// 
@@ -12,8 +16,6 @@ use num_complex::ComplexFloat;
 /// 
 /// where 'N' is the total number of variables, and 'M' is the total number of functions
 /// 
-/// consult the helper utility [`multicalc::utils::helper::transpose`] to transpose the matrix shape if required
-/// 
 /// NOTE: Returns a Result<T, ErrorCode>
 /// Possible ErrorCode are:
 /// VectorOfFunctionsCannotBeEmpty -> if function_matrix argument is an empty array
@@ -88,5 +90,101 @@ pub fn get_custom<T: ComplexFloat, const NUM_FUNCS: usize, const NUM_VARS: usize
         }
     }
 
+    return Ok(result);
+}
+
+
+/// This is a variant of [`get()`] that uses heap allocations.
+/// Useful when handling large datasets to avoid a stack overflow
+/// 
+/// Can handle multivariable functions of any order or complexity
+/// 
+/// The 2-D matrix returned has the structure [[df1/dvar1, df1/dvar2, ... , df1/dvarN],
+///                                            [         ...              ], 
+///                                            [dfM/dvar1, dfM/dvar2, ... , dfM/dvarN]]
+/// 
+/// where 'N' is the total number of variables, and 'M' is the total number of functions
+/// 
+/// NOTE: Returns a Result<T, ErrorCode>
+/// Possible ErrorCode are:
+/// VectorOfFunctionsCannotBeEmpty -> if function_matrix argument is an empty array
+/// 
+/// assume our function vector is (x*y*z ,  x^2 + y^2). First define both the functions
+/// ```
+/// use multicalc::numerical_derivative::jacobian;
+///    let my_func1 = | args: &[f64; 3] | -> f64 
+///    { 
+///        return args[0]*args[1]*args[2]; //x*y*z
+///    };
+/// 
+///    let my_func2 = | args: &[f64; 3] | -> f64 
+///    { 
+///        return args[0].powf(2.0) + args[1].powf(2.0); //x^2 + y^2
+///    };
+/// 
+/// //define the function vector
+/// let function_matrix: Vec<Box<dyn Fn(&[f64; 3]) -> f64>> = vec![Box::new(my_func1), Box::new(my_func2)];
+/// let points = [1.0, 2.0, 3.0]; //the point around which we want the jacobian matrix
+/// 
+/// let result = jacobian::get_on_heap(&function_matrix, &points).unwrap();
+/// ```
+/// 
+/// the above example can also be extended to complex numbers:
+///```
+/// use multicalc::numerical_derivative::jacobian;
+///    let my_func1 = | args: &[num_complex::Complex64; 3] | -> num_complex::Complex64 
+///    { 
+///        return args[0]*args[1]*args[2]; //x*y*z
+///    };
+/// 
+///    let my_func2 = | args: &[num_complex::Complex64; 3] | -> num_complex::Complex64 
+///    { 
+///        return args[0].powf(2.0) + args[1].powf(2.0); //x^2 + y^2
+///    };
+/// 
+/// //define the function vector
+/// let function_matrix: Vec<Box<dyn Fn(&[num_complex::Complex64; 3]) -> num_complex::Complex64>> = vec![Box::new(my_func1), Box::new(my_func2)];
+/// 
+/// //the point around which we want the jacobian matrix
+/// let points = [num_complex::c64(1.0, 3.0), num_complex::c64(2.0, 3.5), num_complex::c64(3.0, 0.0)];
+/// 
+/// let result = jacobian::get_on_heap(&function_matrix, &points).unwrap();
+///``` 
+/// 
+#[cfg(feature = "heap")]
+pub fn get_on_heap<T: ComplexFloat, const NUM_VARS: usize>(function_matrix: &Vec<Box<dyn Fn(&[T; NUM_VARS]) -> T>>, vector_of_points: &[T; NUM_VARS]) -> Result<Vec<Vec<T>>, ErrorCode>
+{    
+    return get_on_heap_custom(function_matrix, vector_of_points, mode::DEFAULT_STEP_SIZE, mode::DiffMode::CentralFixedStep);
+}
+
+
+///same as [get_on_heap()] but with the option to change the differentiation mode used, reserved for more advanced users
+/// NOTE: Returns a Result<T, ErrorCode>
+/// Possible ErrorCode are:
+/// VectorOfFunctionsCannotBeEmpty -> if function_matrix argument is an empty array
+/// NumberOfStepsCannotBeZero -> if the derivative step size is zero
+#[cfg(feature = "heap")]
+pub fn get_on_heap_custom<T: ComplexFloat, const NUM_VARS: usize>(function_matrix: &Vec<Box<dyn Fn(&[T; NUM_VARS]) -> T>>, vector_of_points: &[T; NUM_VARS], step_size: f64, mode: mode::DiffMode) -> Result<Vec<Vec<T>>, ErrorCode>
+{    
+    if function_matrix.is_empty()
+    {
+        return Err(ErrorCode::VectorOfFunctionsCannotBeEmpty);
+    }
+
+    let num_funcs = function_matrix.len();
+
+    let mut result: Vec<Vec<T>> = Vec::new();
+
+    for row_index in 0..num_funcs
+    {
+        let mut cur_row: Vec<T> = Vec::new();
+        for col_index in 0..NUM_VARS
+        {
+            cur_row.push(single_derivative::get_partial_custom(&function_matrix[row_index], col_index, vector_of_points, step_size, mode)?);
+        }
+
+        result.push(cur_row);
+    }
+    
     return Ok(result);
 }

src/numerical_derivative/test.rs

@@ -4,9 +4,12 @@ use crate::numerical_derivative::single_derivative as single_derivative;
 use crate::numerical_derivative::double_derivative as double_derivative;
 use crate::numerical_derivative::triple_derivative as triple_derivative;
 use crate::numerical_derivative::jacobian as jacobian;
-
 use crate::numerical_derivative::hessian as hessian;
 
+#[cfg(feature = "heap")]
+use std::{boxed::Box, vec::Vec};
+
+
 #[test]
 fn test_single_derivative_forward_difference() 
 {
@@ -723,6 +726,82 @@ fn test_jacobian_1_complex()
     }
 }
 
+#[test]
+#[cfg(feature = "heap")]
+fn test_jacobian_2() 
+{
+    //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 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);
+        }
+    }
+}
+
+#[test]
+#[cfg(feature = "heap")]
+fn test_jacobian_2_complex() 
+{
+    //function is x*y*z
+    let func1 = | args: &[num_complex::Complex64; 3] | -> num_complex::Complex64 
+    { 
+        return args[0]*args[1]*args[2];
+    };
+
+    //function is x^2 + y^2
+    let func2 = | args: &[num_complex::Complex64; 3] | -> num_complex::Complex64 
+    { 
+        return args[0].powf(2.0) + args[1].powf(2.0);
+    };
+
+    let function_matrix: Vec<Box<dyn Fn(&[num_complex::Complex64; 3]) -> num_complex::Complex64>> = std::vec![Box::new(func1), Box::new(func2)];
+
+    //the point around which we want the jacobian matrix
+    let points = [num_complex::c64(1.0, 3.0), num_complex::c64(2.0, 3.5), num_complex::c64(3.0, 0.0)];
+
+    let result = 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 = [[points[1]*points[2], points[0]*points[2], points[0]*points[1]], 
+                                                   [2.0*points[0], 2.0*points[1], 0.0*points[2]]];
+
+
+    for i in 0..function_matrix.len()
+    {
+        for j in 0..points.len()
+        {
+            assert!(num_complex::ComplexFloat::abs(result[i][j].re - expected_result[i][j].re) < 0.000001);
+            assert!(num_complex::ComplexFloat::abs(result[i][j].im - expected_result[i][j].im) < 0.000001);
+        }
+    }
+}
+
 #[test]
 fn test_jacobian_1_error() 
 {