compile --no-default-features --target wasm32-unknown-unknown

Cargo.toml

@@ -7,27 +7,51 @@ description = "Rust library for solving analytic and ode-defined pharmacometric
 license = "GPL-3.0"
 documentation = "https://lapkb.github.io/pharmsol/"
 
+
 [features]
-default = []
+
+
+default = [
+    "native"
+]
+
+# `native` now contains only crates that are known to break or behave poorly
+# on wasm (threads, native-optimized code, or platform-specific backends).
+native = [
+    "dep:rayon",
+    "ndarray/rayon",
+    "dep:argmin",
+    "dep:argmin-math",
+    "dep:getrandom",
+    "dep:statrs",
+]
+
+wasm = []
+
 exa = ["libloading"]
 
 [dependencies]
 cached = { version = "0.56.0" }
 csv = "1.3.0"
-diffsol = "=0.7.0"
+diffsol = { version = "=0.7.0",  features = ["faer", "nalgebra"] }
 libloading = { version = "0.8.6", optional = true, features = [] }
-nalgebra = "0.34.1"
-ndarray = { version = "0.16.1", features = ["rayon"] }
-rand = "0.9.0"
-rand_distr = "0.5.0"
-rayon = "1.10.0"
+nalgebra = { version = "0.34.1" }
+ndarray = { version = "0.16.1"}
+rand = { version = "0.9.0"}
+rand_distr = { version = "0.5.0"}
+rayon = { version = "1.10.0", optional = true }
 serde = { version = "1.0.201", features = ["derive"] }
 serde_json = "1.0.117"
-statrs = "0.18.0"
+statrs = { version = "0.18.0", optional = true }
 thiserror = "2.0.11"
-argmin = "0.11.0"
-argmin-math = "0.5.1"
+argmin = { version = "0.11.0", optional = true }
+argmin-math = { version = "0.5.1", optional = true }
 tracing = "0.1.41"
+getrandom = { version = "0.3", optional = true }
+
+[target.'cfg(target_arch = "wasm32")'.dependencies]
+
+getrandom = { version = "0.3", features = ["wasm_js"] }
 
 [dev-dependencies]
 criterion = { version = "0.7.0", features = ["html_reports"] }

src/lib.rs

@@ -10,7 +10,10 @@ pub use crate::data::builder::SubjectBuilderExt;
 pub use crate::data::Interpolation::*;
 pub use crate::data::*;
 pub use crate::equation::*;
+#[cfg(feature = "native")]
 pub use crate::optimize::effect::get_e2;
+
+#[cfg(feature = "native")]
 pub use crate::optimize::spp::SppOptimizer;
 pub use crate::simulator::equation::{self, ODE};
 pub use error::PharmsolError;

src/optimize/effect.rs

@@ -1,8 +1,10 @@
+#[cfg(feature = "native")]
 use argmin::{
     core::{CostFunction, Executor, TerminationReason, TerminationStatus},
     solver::neldermead::NelderMead,
 };
 
+#[cfg(feature = "native")]
 #[derive(Debug, Clone)]
 struct BestM0 {
     a: f64,
@@ -14,6 +16,7 @@ struct BestM0 {
 }
 
 /// We'll optimize over y = ln(xm), so Param = f64 (the log of xm)
+#[cfg(feature = "native")]
 impl CostFunction for BestM0 {
     type Param = f64; // this is ln(xm)
     type Output = f64;
@@ -59,12 +62,14 @@ impl CostFunction for BestM0 {
     }
 }
 
+#[cfg(feature = "native")]
 #[derive(thiserror::Error, Debug)]
 pub enum BestM0Error {
     #[error("Optimization error: {0}")]
     OptimizationError(String),
 }
 
+#[cfg(feature = "native")]
 impl BestM0 {
     /// start and step are in log-space (ln(x))
     fn get_best(&self, start_log: f64, step_log: f64) -> Result<(f64, f64, bool), BestM0Error> {
@@ -111,8 +116,7 @@ impl BestM0 {
     }
 }
 
-/// find_m0 left largely as-is, but consider returning Result<f64>
-/// Keep in mind it expects a,b,h1,h2 in valid ranges
+#[cfg(feature = "native")]
 fn find_m0(afinal: f64, b: f64, alpha: f64, h1: f64, h2: f64) -> f64 {
     let noint = 1000;
     let del_a = afinal / (noint as f64);
@@ -150,6 +154,7 @@ fn find_m0(afinal: f64, b: f64, alpha: f64, h1: f64, h2: f64) -> f64 {
     xm
 }
 
+#[cfg(feature = "native")]
 pub fn get_e2(a: f64, b: f64, w: f64, h1: f64, h2: f64, alpha_s: f64) -> f64 {
     // trivial cases
     if a.abs() < 1.0e-12 && b.abs() < 1.0e-12 {
@@ -237,3 +242,11 @@ pub fn get_e2(a: f64, b: f64, w: f64, h1: f64, h2: f64, alpha_s: f64) -> f64 {
         }
     }
 }
+
+// Non-native (wasm) stub implementation
+#[cfg(not(feature = "native"))]
+pub fn get_e2(_a: f64, _b: f64, _w: f64, _h1: f64, _h2: f64, _alpha_s: f64) -> f64 {
+    // Provide a conservative default for wasm builds. This is a best-effort
+    // stub — accurate behavior requires the native optimizer.
+    0.0
+}

src/optimize/spp.rs

@@ -1,19 +1,24 @@
+#[cfg(feature = "native")]
 use argmin::{
     core::{CostFunction, Error, Executor},
     solver::neldermead::NelderMead,
 };
 
+#[cfg(feature = "native")]
 use ndarray::{Array1, Axis};
 
+#[cfg(feature = "native")]
 use crate::{prelude::simulator::psi, Data, Equation, ErrorModels};
 
+#[cfg(feature = "native")]
 pub struct SppOptimizer<'a, E: Equation> {
     equation: &'a E,
     data: &'a Data,
     sig: &'a ErrorModels,
     pyl: &'a Array1<f64>,
 }
 
+#[cfg(feature = "native")]
 impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
     type Param = Vec<f64>;
     type Output = f64;
@@ -41,6 +46,7 @@ impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
     }
 }
 
+#[cfg(feature = "native")]
 impl<'a, E: Equation> SppOptimizer<'a, E> {
     pub fn new(
         equation: &'a E,
@@ -66,6 +72,7 @@ impl<'a, E: Equation> SppOptimizer<'a, E> {
     }
 }
 
+#[cfg(feature = "native")]
 fn create_initial_simplex(initial_point: &[f64]) -> Vec<Vec<f64>> {
     let num_dimensions = initial_point.len();
     let perturbation_percentage = 0.008;
@@ -91,3 +98,39 @@ fn create_initial_simplex(initial_point: &[f64]) -> Vec<Vec<f64>> {
 
     vertices
 }
+
+// WASM / non-native stubs
+#[cfg(not(feature = "native"))]
+use ndarray::Array1;
+
+#[cfg(not(feature = "native"))]
+use crate::{Data, Equation, ErrorModels};
+
+#[cfg(not(feature = "native"))]
+pub struct SppOptimizer<'a, E> {
+    _phantom: std::marker::PhantomData<&'a E>,
+}
+
+#[cfg(not(feature = "native"))]
+impl<'a, E: Equation> SppOptimizer<'a, E> {
+    pub fn new(
+        _equation: &'a E,
+        _data: &'a Data,
+        _sig: &'a ErrorModels,
+        _pyl: &'a Array1<f64>,
+    ) -> Self {
+        SppOptimizer {
+            _phantom: std::marker::PhantomData,
+        }
+    }
+
+    pub fn optimize_point(
+        self,
+        _spp: Array1<f64>,
+    ) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
+        Err(Box::new(std::io::Error::new(
+            std::io::ErrorKind::Other,
+            "SppOptimizer is disabled when feature `native` is not enabled",
+        )))
+    }
+}

src/simulator/equation/sde/mod.rs

@@ -4,6 +4,8 @@ use diffsol::{NalgebraContext, Vector};
 use nalgebra::DVector;
 use ndarray::{concatenate, Array2, Axis};
 use rand::{rng, Rng};
+
+#[cfg(feature = "native")]
 use rayon::prelude::*;
 
 use cached::proc_macro::cached;
@@ -142,9 +144,18 @@ impl State for Vec<DVector<f64>> {
     /// * `input` - Index of the input compartment
     /// * `amount` - Amount to add to the compartment
     fn add_bolus(&mut self, input: usize, amount: f64) {
-        self.par_iter_mut().for_each(|particle| {
-            particle[input] += amount;
-        });
+        #[cfg(feature = "native")]
+        {
+            self.par_iter_mut().for_each(|particle| {
+                particle[input] += amount;
+            });
+        }
+        #[cfg(not(feature = "native"))]
+        {
+            self.iter_mut().for_each(|particle| {
+                particle[input] += amount;
+            });
+        }
     }
 }
 
@@ -239,20 +250,40 @@ impl EquationPriv for SDE {
         ti: f64,
         tf: f64,
     ) -> Result<(), PharmsolError> {
-        state.par_iter_mut().for_each(|particle| {
-            *particle = simulate_sde_event(
-                &self.drift,
-                &self.diffusion,
-                particle.clone().into(),
-                support_point,
-                covariates,
-                infusions,
-                ti,
-                tf,
-            )
-            .inner()
-            .clone();
-        });
+        #[cfg(feature = "native")]
+        {
+            state.par_iter_mut().for_each(|particle| {
+                *particle = simulate_sde_event(
+                    &self.drift,
+                    &self.diffusion,
+                    particle.clone().into(),
+                    support_point,
+                    covariates,
+                    infusions,
+                    ti,
+                    tf,
+                )
+                .inner()
+                .clone();
+            });
+        }
+        #[cfg(not(feature = "native"))]
+        {
+            for particle in state.iter_mut() {
+                *particle = simulate_sde_event(
+                    &self.drift,
+                    &self.diffusion,
+                    particle.clone().into(),
+                    support_point,
+                    covariates,
+                    infusions,
+                    ti,
+                    tf,
+                )
+                .inner()
+                .clone();
+            }
+        }
         Ok(())
     }
     fn nparticles(&self) -> usize {
@@ -275,17 +306,34 @@ impl EquationPriv for SDE {
         output: &mut Self::P,
     ) -> Result<(), PharmsolError> {
         let mut pred = vec![Prediction::default(); self.nparticles];
-        pred.par_iter_mut().enumerate().for_each(|(i, p)| {
-            let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
-            (self.out)(
-                &x[i].clone().into(),
-                &V::from_vec(support_point.clone(), NalgebraContext),
-                observation.time(),
-                covariates,
-                &mut y,
-            );
-            *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
-        });
+        #[cfg(feature = "native")]
+        {
+            pred.par_iter_mut().enumerate().for_each(|(i, p)| {
+                let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
+                (self.out)(
+                    &x[i].clone().into(),
+                    &V::from_vec(support_point.clone(), NalgebraContext),
+                    observation.time(),
+                    covariates,
+                    &mut y,
+                );
+                *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
+            });
+        }
+        #[cfg(not(feature = "native"))]
+        {
+            for (i, p) in pred.iter_mut().enumerate() {
+                let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
+                (self.out)(
+                    &x[i].clone().into(),
+                    &V::from_vec(support_point.clone(), NalgebraContext),
+                    observation.time(),
+                    covariates,
+                    &mut y,
+                );
+                *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
+            }
+        }
         let out = Array2::from_shape_vec((self.nparticles, 1), pred.clone())?;
         *output = concatenate(Axis(1), &[output.view(), out.view()]).unwrap();
         //e = y[t] .- x[:,1]