Play with rust and implement semi integration in rust

.gitignore

@@ -135,3 +135,18 @@ dmypy.json
 
 # pixi environments
 .pixi
+
+# Generated by Cargo
+# will have compiled files and executables
+debug/
+target/
+
+# Remove Cargo.lock from gitignore if creating an executable, leave it for libraries
+# More information here https://doc.rust-lang.org/cargo/guide/cargo-toml-vs-cargo-lock.html
+Cargo.lock
+
+# These are backup files generated by rustfmt
+**/*.rs.bk
+
+# MSVC Windows builds of rustc generate these, which store debugging information
+*.pdb

Cargo.toml

@@ -0,0 +1,15 @@
+[package]
+name = "ec-tools"
+version = "0.1.0"
+edition = "2024"
+
+# [package.metadata.maturin]
+# name = "ec_tools.semi_integration_rs"
+
+[dependencies]
+pyo3 = { version = "^0.24", features = ["extension-module"] }
+numpy = "^0.24"
+
+[lib]
+name = "ec_tools"
+crate-type = ["cdylib"]

ec_tools/helper.py

@@ -69,6 +69,31 @@ def determine_scan_rate(t, x):
     return np.abs(np.diff(x) / np.diff(t)).mean()
 
 
+def find_extrema_indeces(y, mode="all"):
+    """Return the indexes of limits of cyclic voltammetry.
+
+    Workaround for Windows platform:
+    list(array([...], dtype=int64))
+    TEST:
+    >>> E = np.array([0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.14, 0.13, 0.12, 0.11, 0.10, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14])
+    >>> list(find_extrema_indeces(E))
+    [5, 11]
+
+    >>> list(find_extrema_indeces(E, mode="pos"))
+    [5]
+    """
+    signs = np.diff(np.sign(np.diff(y)))
+
+    if mode == "all":
+        extrema = np.where((signs == 2) | (signs == -2))[0]
+    elif mode == "pos":
+        extrema = np.where(signs == -2)[0]
+    elif mode == "neg":
+        extrema = np.where(signs == 2)[0]
+
+    return extrema + 1
+
+
 def detect_step(t, x):
     """Returns the index of the step in given t and x arrays.
     Index is the where the changed value of t located.

ec_tools/semi_integration_rs/__init__.py

@@ -0,0 +1 @@
+from ec_tools.semi_integration_rs import *

ec_tools/semi_integration_rs/__init__.pyi

@@ -0,0 +1,4 @@
+def semi_integration_rs(param1: str, param2: int) -> int:
+    """
+    Function description
+    """

pyproject.toml

@@ -1,6 +1,6 @@
 [build-system]
-requires = ["setuptools"]
-build-backend = "setuptools.build_meta"
+requires = ["maturin>=1.5,<2.0"]
+build-backend = "maturin"
 
 [project]
 name = "ec-tools"
@@ -10,6 +10,8 @@ readme      = "README.md"
 requires-python = ">=3.10.0"
 dependencies    = [ "numpy>=2.0.3,<3", "transonic>=0.7.3,<0.9" ]
 version = "0.1.0"
+[tool.maturin]
+features = ["pyo3/extension-module"]
 
 [tool.setuptools]
 packages = ["ec_tools",]

src/lib.rs

@@ -0,0 +1,10 @@
+// use pyo3::prelude::*;
+
+mod semi_integration_rs;
+
+// #[pymodule]
+// fn _rs(_py: Python, m: &PyModule) -> PyResult<()> {
+//     semi_integration::init_linalg(_py, m)?;
+
+//     Ok(())
+// }

src/semi_integration_rs/mod.rs

@@ -0,0 +1,79 @@
+use pyo3::{pymodule, types::PyModule, PyResult, Python};
+use numpy::ndarray::{Array1, ArrayView1};
+use numpy::{IntoPyArray, PyReadonlyArray1, PyArray1};
+
+use std::f64::consts::PI;
+
+#[pymodule]
+fn semi_integration_rs(_py: Python<'_>, m: &PyModule) -> PyResult<()> {
+    fn prepare_kernel(q: f64, Δx: f64, N: u64, c1: u64, c2: u64) -> (Array1<f64>, Array1<f64>, Array1<f64>) {
+        let length = N as f64;
+        let τ0 = Δx * length.powf(0.5);
+        let a0 = (PI * q).sin() / (PI *q * τ0.powf(q));
+        // total number of filters
+        let n_filters = 2 * c1 * c2 + 1;
+        // dimension according to the number of filters
+        // filter weights
+        let mut w1 = Array1::<f64>::zeros(n_filters as usize);
+        let mut w2 = Array1::<f64>::zeros(n_filters as usize);
+        // auxiliary array
+        let s = Array1::<f64>::zeros(n_filters as usize);
+
+        for i in 0..2*c1 as usize *c2 as usize {
+            let j:f64 = i as f64 - c1 as f64 * c2 as f64;
+            let a_j = (a0 / c2 as f64) * f64::exp(j as f64 / c2 as f64);
+            let t_j = τ0 * f64::exp(-j as f64 / (q*c2 as f64) );
+            w1[i] = t_j / (Δx + t_j);
+            w2[i] = a_j * (1f64 - w1[i])};
+        (s, w1, w2)
+    }
+
+    fn semi_integration(y: ArrayView1<'_, f64>, q: Option<f64>, Δx: Option<f64>, c1: Option<u64>, c2: Option<u64>) -> Array1<f64> {
+        let N = y.len();
+        let mut R = Array1::<f64>::zeros(N);
+        let ( mut s, w1, w2) = prepare_kernel(
+            q.unwrap_or( -0.5),
+            Δx.unwrap_or(1.0),
+            N as u64,
+            c1.unwrap_or(2u64),
+            c2.unwrap_or(8u64)
+        );
+        for k in 1..N {
+            for i in 0..s.len() {
+                s[i] = s[i]*w1[i] + &y[k]*w2[i];
+                R[k] = R[k] + s[i];
+            }
+        }
+    R
+    } 
+// wrapper of `semi_integration_rs`
+#[pyfn(m)]
+#[pyo3(name = "semi_integration")]
+fn semi_integration_py<'py>(
+    py: Python<'py>,
+    y: PyReadonlyArray1<f64>,
+    q: Option<f64>,
+    Δx: Option<f64>,
+    c1: Option<u64>,
+    c2: Option<u64>,
+) -> &'py PyArray1<f64> {
+    semi_integration(y.as_array(), q, Δx, c1, c2).into_pyarray(py)
+    // Array1::<f64>::zeros(2).into_pyarray(py)
+}
+
+// wrapper of `prepare_kernel_rs`
+#[pyfn(m)]
+#[pyo3(name = "prepare_kernel")]
+fn prepare_kernel_py<'py>(
+    py: Python<'py>,
+    q: f64,
+    Δx: f64,
+    N: u64,
+    c1: u64,
+    c2: u64,
+) -> (&'py PyArray1<f64>, &'py PyArray1<f64>, &'py PyArray1<f64>) {
+    let (s , w1, w2) = prepare_kernel(q, Δx, N as u64, c1, c2);
+    (s.into_pyarray(py) , w1.into_pyarray(py), w2.into_pyarray(py))
+}
+    Ok(())
+}