diff --git a/CMakeLists.txt b/CMakeLists.txt index ec8400d8..d086d966 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,11 +143,11 @@ endif() # set(SRC_FILES ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolbng.c + ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolgraphics.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolboxes.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolcmd.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolcomparts.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolfilament.c - ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolgraphics.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolmolec.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolport.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smollattice.c @@ -155,8 +155,7 @@ set(SRC_FILES ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolrule.c ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolsim.cpp ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolsurface.c - ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolwall.c - ) + ${CMAKE_SOURCE_DIR}/source/Smoldyn/smolwall.c) set(MAIN_FILES ${CMAKE_SOURCE_DIR}/source/Smoldyn/smoldyn.cpp) diff --git a/rust/Cargo.toml b/rust/Cargo.toml new file mode 100644 index 00000000..29890292 --- /dev/null +++ b/rust/Cargo.toml @@ -0,0 +1,3 @@ +[workspace] +resolver = "2" +members = ["libsmoldyn", "smoldyn-rs"] diff --git a/rust/libsmoldyn/Cargo.toml b/rust/libsmoldyn/Cargo.toml new file mode 100644 index 00000000..dad29f56 --- /dev/null +++ b/rust/libsmoldyn/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "libsmoldyn" +version = "0.27.1" +edition = "2024" + +[dependencies] +cxx = { version = "1.0.202", features = ["c++17"] } +tracing = "0.1.44" + +[build-dependencies] +cmake = "0.1.58" +cxx-build = "1.0.202" diff --git a/rust/libsmoldyn/build.rs b/rust/libsmoldyn/build.rs new file mode 100644 index 00000000..9f842328 --- /dev/null +++ b/rust/libsmoldyn/build.rs @@ -0,0 +1,23 @@ +use cmake::Config; + +fn main() { + let dst = Config::new("../..") + .define("OPTION_PYTHON", "OFF") + .define("OPTION_STATIC", "ON") + .define("OPTION_USE_OPENGL", "OFF") + .define("OPTION_USE_LIBTIFF", "OFF") + .build(); + + cxx_build::bridge("src/lib.rs") + .include("../../source") + .include("include") + .include(format!("{}/build", dst.display())) + .compile("simulation"); + + println!("cargo:rerun-if-changed=include/smolrs.h"); + // cxx_build's rerun-if-changed replaces cargo's default, so list the C/C++ sources too. + println!("cargo:rerun-if-changed=../../source"); + println!("cargo:rerun-if-changed=../../CMakeLists.txt"); + println!("cargo:rustc-link-search=native={}/build", dst.display()); + println!("cargo:rustc-link-lib=static=smoldyn_static"); +} diff --git a/rust/libsmoldyn/include/smolrs.h b/rust/libsmoldyn/include/smolrs.h new file mode 100644 index 00000000..f0a264e2 --- /dev/null +++ b/rust/libsmoldyn/include/smolrs.h @@ -0,0 +1,45 @@ +// Read-only accessors for `simstruct` fields that libsmoldyn.h has no getter for. +// +// Rust sees `simstruct` as an opaque type, so reading fields goes through these +// inline functions; that keeps the struct layout in C++ where it is defined. +#pragma once + +#include +#include +#include + +#include "Smoldyn/smoldyn.h" + +inline int smolrs_dim(const simstruct &sim) { return sim.dim; } +inline double smolrs_time(const simstruct &sim) { return sim.time; } +inline double smolrs_time_start(const simstruct &sim) { return sim.tmin; } +inline double smolrs_time_stop(const simstruct &sim) { return sim.tmax; } +inline double smolrs_time_step(const simstruct &sim) { return sim.dt; } + +// Position of the low (highside=0) or high (highside=1) wall along dimension d. +// Returns NaN if the walls are not allocated or d/highside are out of range. +inline double smolrs_wall_pos(const simstruct &sim, int d, int highside) { + if (!sim.wlist || d < 0 || d >= sim.dim || highside < 0 || highside > 1) + return NAN; + return sim.wlist[2 * d + highside]->pos; +} + +// Length of the configuration file path ("" for simulations not loaded from a file). +inline size_t smolrs_filepath_len(const simstruct &sim) { + return sim.filepath ? strlen(sim.filepath) : 0; +} + +// Number of species slots, including the "empty" species at index 0. +inline int smolrs_nspecies(const simstruct &sim) { + return sim.mols ? sim.mols->nspecies : 0; +} + +// Name of species i, or nullptr if out of range. Owned by the simulation. +inline const char *smolrs_species_name(const simstruct &sim, int i) { + if (!sim.mols || i < 0 || i >= sim.mols->nspecies) + return nullptr; + return sim.mols->spname[i]; +} + +// Free an array that libsmoldyn allocated with calloc (smolGetOutputData). +inline void smolrs_free_doubles(double *array) { free(array); } diff --git a/rust/libsmoldyn/src/error.rs b/rust/libsmoldyn/src/error.rs new file mode 100644 index 00000000..4a9f0472 --- /dev/null +++ b/rust/libsmoldyn/src/error.rs @@ -0,0 +1,81 @@ +use std::ffi::{CStr, c_char}; +use std::fmt; + +use crate::ffi::{self, ErrorCode}; + +pub(crate) const STRCHARLONG: usize = 4096; + +pub type SmolResult = Result; + +#[derive(Debug, Clone)] +pub enum SmolError { + Smoldyn { + code: ErrorCode, + function: String, + message: String, + }, + InvalidArgument(String), +} + +impl SmolError { + pub(crate) fn last(fallback: ErrorCode) -> Self { + let mut function = [0 as c_char; STRCHARLONG]; + let mut message = [0 as c_char; STRCHARLONG]; + let code = unsafe { ffi::smolGetError(function.as_mut_ptr(), message.as_mut_ptr(), 1) }; + let code = if code == ErrorCode::Ok { + fallback + } else { + code + }; + let to_string = |buf: &[c_char]| { + unsafe { CStr::from_ptr(buf.as_ptr()) } + .to_string_lossy() + .into_owned() + }; + SmolError::Smoldyn { + code, + function: to_string(&function), + message: to_string(&message), + } + } +} + +impl fmt::Display for SmolError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + SmolError::Smoldyn { + code, + function, + message, + } => write!(f, "smoldyn error {code:?} in {function}: {message}"), + SmolError::InvalidArgument(msg) => write!(f, "invalid argument: {msg}"), + } + } +} + +impl std::error::Error for SmolError {} + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Status { + Ok, + Notify, + Warning, +} + +pub(crate) fn check(code: ErrorCode) -> SmolResult { + match code { + ErrorCode::Ok => Ok(Status::Ok), + ErrorCode::Notify => { + // clear it, otherwise libsmoldyn keeps returning ECnotify from later calls. + let note = SmolError::last(code); + tracing::debug!("{note}"); + Ok(Status::Notify) + } + ErrorCode::Warning => { + let err = SmolError::last(code); + tracing::warn!("{err}"); + Ok(Status::Warning) + } + _ => Err(SmolError::last(code)), + } +} diff --git a/rust/libsmoldyn/src/lib.rs b/rust/libsmoldyn/src/lib.rs new file mode 100644 index 00000000..066b7dec --- /dev/null +++ b/rust/libsmoldyn/src/lib.rs @@ -0,0 +1,679 @@ +mod error; +mod sim; + +pub use error::{SmolError, SmolResult}; +pub use ffi::{CmptLogic, DrawMode, MolecState, PanelFace, PanelShape, RevParam, SrfAction}; +pub use sim::{OutputData, Progress, RateKind, Side, Sim, SurfaceStyle}; + +#[allow(clippy::missing_safety_doc, clippy::too_many_arguments)] +#[cxx::bridge] +pub mod ffi { + #[repr(i32)] + #[derive(Debug)] + enum ErrorCode { + #[cxx_name = "ECok"] + Ok = 0, + #[cxx_name = "ECnotify"] + Notify = -1, + #[cxx_name = "ECwarning"] + Warning = -2, + #[cxx_name = "ECnonexist"] + Nonexist = -3, + #[cxx_name = "ECall"] + All = -4, + #[cxx_name = "ECmissing"] + Missing = -5, + #[cxx_name = "ECbounds"] + Bounds = -6, + #[cxx_name = "ECsyntax"] + Syntax = -7, + #[cxx_name = "ECerror"] + Error = -8, + #[cxx_name = "ECmemory"] + Memory = -9, + #[cxx_name = "ECbug"] + Bug = -10, + #[cxx_name = "ECsame"] + Same = -11, + #[cxx_name = "ECwildcard"] + Wildcard = -12, + } + + #[repr(i32)] + #[derive(Debug)] + enum MolecState { + #[cxx_name = "MSsoln"] + Soln = 0, + #[cxx_name = "MSfront"] + Front = 1, + #[cxx_name = "MSback"] + Back = 2, + #[cxx_name = "MSup"] + Up = 3, + #[cxx_name = "MSdown"] + Down = 4, + #[cxx_name = "MSbsoln"] + Bsoln = 5, + #[cxx_name = "MSall"] + All = 6, + #[cxx_name = "MSnone"] + None = 7, + #[cxx_name = "MSsome"] + Some = 8, + } + + #[repr(i32)] + #[derive(Debug)] + enum PanelFace { + #[cxx_name = "PFfront"] + Front = 0, + #[cxx_name = "PFback"] + Back = 1, + #[cxx_name = "PFnone"] + None = 2, + #[cxx_name = "PFboth"] + Both = 3, + } + + #[repr(i32)] + #[derive(Debug)] + enum PanelShape { + #[cxx_name = "PSrect"] + Rect = 0, + #[cxx_name = "PStri"] + Tri = 1, + #[cxx_name = "PSsph"] + Sph = 2, + #[cxx_name = "PScyl"] + Cyl = 3, + #[cxx_name = "PShemi"] + Hemi = 4, + #[cxx_name = "PSdisk"] + Disk = 5, + #[cxx_name = "PSall"] + All = 6, + #[cxx_name = "PSnone"] + None = 7, + } + + #[repr(i32)] + #[derive(Debug)] + enum SrfAction { + #[cxx_name = "SAreflect"] + Reflect = 0, + #[cxx_name = "SAtrans"] + Trans = 1, + #[cxx_name = "SAabsorb"] + Absorb = 2, + #[cxx_name = "SAjump"] + Jump = 3, + #[cxx_name = "SAport"] + Port = 4, + #[cxx_name = "SAmult"] + Mult = 5, + #[cxx_name = "SAno"] + No = 6, + #[cxx_name = "SAnone"] + None = 7, + #[cxx_name = "SAadsorb"] + Adsorb = 8, + #[cxx_name = "SArevdes"] + Revdes = 9, + #[cxx_name = "SAirrevdes"] + Irrevdes = 10, + #[cxx_name = "SAflip"] + Flip = 11, + } + + #[repr(i32)] + #[derive(Debug)] + enum DrawMode { + #[cxx_name = "DMno"] + No = 0, + #[cxx_name = "DMvert"] + Vert = 1, + #[cxx_name = "DMedge"] + Edge = 2, + #[cxx_name = "DMve"] + Ve = 3, + #[cxx_name = "DMface"] + Face = 4, + #[cxx_name = "DMvf"] + Vf = 5, + #[cxx_name = "DMef"] + Ef = 6, + #[cxx_name = "DMvef"] + Vef = 7, + #[cxx_name = "DMnone"] + None = 8, + } + + #[repr(i32)] + #[derive(Debug)] + enum CmptLogic { + #[cxx_name = "CLequal"] + Equal = 0, + #[cxx_name = "CLequalnot"] + Equalnot = 1, + #[cxx_name = "CLand"] + And = 2, + #[cxx_name = "CLor"] + Or = 3, + #[cxx_name = "CLxor"] + Xor = 4, + #[cxx_name = "CLandnot"] + Andnot = 5, + #[cxx_name = "CLornot"] + Ornot = 6, + #[cxx_name = "CLnone"] + None = 7, + } + + #[repr(i32)] + #[derive(Debug)] + enum RevParam { + #[cxx_name = "RPnone"] + None = 0, + #[cxx_name = "RPirrev"] + Irrev = 1, + #[cxx_name = "RPconfspread"] + Confspread = 2, + #[cxx_name = "RPbounce"] + Bounce = 3, + #[cxx_name = "RPpgem"] + Pgem = 4, + #[cxx_name = "RPpgemmax"] + Pgemmax = 5, + #[cxx_name = "RPpgemmaxw"] + Pgemmaxw = 6, + #[cxx_name = "RPratio"] + Ratio = 7, + #[cxx_name = "RPunbindrad"] + Unbindrad = 8, + #[cxx_name = "RPpgem2"] + Pgem2 = 9, + #[cxx_name = "RPpgemmax2"] + Pgemmax2 = 10, + #[cxx_name = "RPratio2"] + Ratio2 = 11, + #[cxx_name = "RPoffset"] + Offset = 12, + #[cxx_name = "RPfixed"] + Fixed = 13, + } + + unsafe extern "C++" { + include!("Smoldyn/libsmoldyn.h"); + + type simstruct; + type ErrorCode; + type MolecState; + type PanelFace; + type PanelShape; + type SrfAction; + type DrawMode; + type CmptLogic; + type RevParam; + + unsafe fn smolGetVersion() -> f64; + + unsafe fn smolGetError( + errorfunction: *mut c_char, + errorstring: *const c_char, + clearerror: i32, + ) -> ErrorCode; + unsafe fn smolClearError(); + unsafe fn smolSetDebugMode(debugmode: i32); + + unsafe fn smolNewSim(dim: i32, lowbounds: *mut f64, highbounds: *mut f64) + -> *mut simstruct; + unsafe fn smolUpdateSim(sim: *mut simstruct) -> ErrorCode; + unsafe fn smolRunTimeStep(sim: *mut simstruct) -> ErrorCode; + unsafe fn smolRunSim(sim: *mut simstruct) -> ErrorCode; + unsafe fn smolRunSimUntil(sim: *mut simstruct, breaktime: f64) -> ErrorCode; + unsafe fn smolFreeSim(sim: *mut simstruct) -> ErrorCode; + unsafe fn smolDisplaySim(sim: *mut simstruct) -> ErrorCode; + + unsafe fn smolPrepareSimFromFile( + filepath: *const c_char, + filename: *const c_char, + flags: *const c_char, + ) -> *mut simstruct; + unsafe fn smolLoadSimFromFile( + filepath: *const c_char, + filename: *const c_char, + simpointer: *mut *mut simstruct, + flags: *const c_char, + ) -> ErrorCode; + unsafe fn smolReadConfigString( + sim: *mut simstruct, + statement: *const c_char, + parameters: *mut c_char, + ) -> ErrorCode; + + unsafe fn smolSetSimFlags(sim: *mut simstruct, flags: *const c_char) -> ErrorCode; + unsafe fn smolSetSimTimes( + sim: *mut simstruct, + timestart: f64, + timestop: f64, + timestep: f64, + ) -> ErrorCode; + unsafe fn smolSetTimeStart(sim: *mut simstruct, timestart: f64) -> ErrorCode; + unsafe fn smolSetTimeStop(sim: *mut simstruct, timestop: f64) -> ErrorCode; + unsafe fn smolSetTimeNow(sim: *mut simstruct, timenow: f64) -> ErrorCode; + unsafe fn smolSetTimeStep(sim: *mut simstruct, timestep: f64) -> ErrorCode; + // `long int`, which is 64 bits on LP64 platforms. + unsafe fn smolSetRandomSeed(sim: *mut simstruct, seed: i64) -> ErrorCode; + unsafe fn smolSetPartitions( + sim: *mut simstruct, + method: *const c_char, + value: f64, + ) -> ErrorCode; + + unsafe fn smolSetGraphicsParams( + sim: *mut simstruct, + method: *const c_char, + timesteps: i32, + delay: i32, + ) -> ErrorCode; + unsafe fn smolSetTiffParams( + sim: *mut simstruct, + timesteps: i32, + tiffname: *const c_char, + lowcount: i32, + highcount: i32, + ) -> ErrorCode; + unsafe fn smolSetLightParams( + sim: *mut simstruct, + lightindex: i32, + ambient: *mut f64, + diffuse: *mut f64, + specular: *mut f64, + position: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetBackgroundStyle(sim: *mut simstruct, color: *mut f64) -> ErrorCode; + unsafe fn smolSetFrameStyle( + sim: *mut simstruct, + thickness: f64, + color: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetGridStyle( + sim: *mut simstruct, + thickness: f64, + color: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetTextStyle(sim: *mut simstruct, color: *mut f64) -> ErrorCode; + unsafe fn smolAddTextDisplay(sim: *mut simstruct, item: *mut c_char) -> ErrorCode; + + unsafe fn smolSetOutputPath(sim: *mut simstruct, path: *const c_char) -> ErrorCode; + unsafe fn smolAddOutputFile( + sim: *mut simstruct, + filename: *mut c_char, + suffix: i32, + append: i32, + ) -> ErrorCode; + unsafe fn smolAddOutputData(sim: *mut simstruct, dataname: *mut c_char) -> ErrorCode; + unsafe fn smolOpenOutputFiles(sim: *mut simstruct, overwrite: i32) -> ErrorCode; + unsafe fn smolAddCommand( + sim: *mut simstruct, + kind: c_char, + on: f64, + off: f64, + step: f64, + multiplier: f64, + commandstring: *const c_char, + ) -> ErrorCode; + unsafe fn smolAddCommandFromString(sim: *mut simstruct, string: *mut c_char) -> ErrorCode; + unsafe fn smolGetOutputData( + sim: *mut simstruct, + dataname: *mut c_char, + nrow: *mut i32, + ncol: *mut i32, + array: *mut *mut f64, + erase: i32, + ) -> ErrorCode; + unsafe fn smolRunCommand(sim: *mut simstruct, commandstring: *const c_char) -> ErrorCode; + + unsafe fn smolAddSpecies( + sim: *mut simstruct, + species: *const c_char, + mollist: *const c_char, + ) -> ErrorCode; + unsafe fn smolGetSpeciesIndex(sim: *mut simstruct, species: *const c_char) -> i32; + unsafe fn smolGetSpeciesName(sim: *mut simstruct, speciesindex: i32, species: *mut c_char); + unsafe fn smolSetSpeciesMobility( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + difc: f64, + drift: *mut f64, + difmatrix: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetMoleculeColor( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + color: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetMoleculeSize( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + size: f64, + ) -> ErrorCode; + unsafe fn smolAddMolList(sim: *mut simstruct, mollist: *const c_char) -> ErrorCode; + unsafe fn smolGetMolListIndex(sim: *mut simstruct, mollist: *const c_char) -> i32; + unsafe fn smolGetMolListName( + sim: *mut simstruct, + mollistindex: i32, + mollist: *mut c_char, + ) -> *mut c_char; + unsafe fn smolSetMolList( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + mollist: *const c_char, + ) -> ErrorCode; + unsafe fn smolSetMaxMolecules(sim: *mut simstruct, maxmolecules: i32) -> ErrorCode; + unsafe fn smolAddSolutionMolecules( + sim: *mut simstruct, + species: *const c_char, + number: i32, + lowposition: *mut f64, + highposition: *mut f64, + ) -> ErrorCode; + unsafe fn smolAddCompartmentMolecules( + sim: *mut simstruct, + species: *const c_char, + number: i32, + compartment: *const c_char, + ) -> ErrorCode; + unsafe fn smolAddSurfaceMolecules( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + number: i32, + surface: *const c_char, + panelshape: PanelShape, + panel: *const c_char, + position: *mut f64, + ) -> ErrorCode; + unsafe fn smolGetMoleculeCount( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + ) -> i32; + unsafe fn smolSetMoleculeStyle( + sim: *mut simstruct, + species: *const c_char, + state: MolecState, + size: f64, + color: *mut f64, + ) -> ErrorCode; + + unsafe fn smolSetBoundaryType( + sim: *mut simstruct, + dimension: i32, + highside: i32, + kind: c_char, + ) -> ErrorCode; + unsafe fn smolAddSurface(sim: *mut simstruct, surface: *const c_char) -> ErrorCode; + unsafe fn smolGetSurfaceIndex(sim: *mut simstruct, surface: *const c_char) -> i32; + unsafe fn smolGetSurfaceName( + sim: *mut simstruct, + surfaceindex: i32, + surface: *mut c_char, + ) -> *mut c_char; + unsafe fn smolSetReactionIntersurface( + sim: *mut simstruct, + reaction: *const c_char, + rulelist: *mut i32, + ) -> ErrorCode; + unsafe fn smolSetSurfaceAction( + sim: *mut simstruct, + surface: *const c_char, + face: PanelFace, + species: *const c_char, + state: MolecState, + action: SrfAction, + newspecies: *const c_char, + ) -> ErrorCode; + unsafe fn smolSetSurfaceRate( + sim: *mut simstruct, + surface: *const c_char, + species: *const c_char, + state: MolecState, + state1: MolecState, + state2: MolecState, + rate: f64, + newspecies: *const c_char, + isinternal: i32, + ) -> ErrorCode; + unsafe fn smolAddPanel( + sim: *mut simstruct, + surface: *const c_char, + panelshape: PanelShape, + panel: *const c_char, + axisstring: *const c_char, + params: *mut f64, + ) -> ErrorCode; + unsafe fn smolGetPanelIndex( + sim: *mut simstruct, + surface: *const c_char, + panelshapeptr: *mut PanelShape, + panel: *const c_char, + ) -> i32; + unsafe fn smolGetPanelName( + sim: *mut simstruct, + surface: *const c_char, + panelshape: PanelShape, + panelindex: i32, + panel: *mut c_char, + ) -> *mut c_char; + unsafe fn smolSetPanelJump( + sim: *mut simstruct, + surface: *const c_char, + panel1: *const c_char, + face1: PanelFace, + panel2: *const c_char, + face2: PanelFace, + isbidirectional: i32, + ) -> ErrorCode; + unsafe fn smolAddSurfaceUnboundedEmitter( + sim: *mut simstruct, + surface: *const c_char, + face: PanelFace, + species: *const c_char, + emitamount: f64, + emitposition: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetSurfaceSimParams( + sim: *mut simstruct, + parameter: *const c_char, + value: f64, + ) -> ErrorCode; + unsafe fn smolAddPanelNeighbor( + sim: *mut simstruct, + surface1: *const c_char, + panel1: *const c_char, + surface2: *const c_char, + panel2: *const c_char, + reciprocal: i32, + ) -> ErrorCode; + unsafe fn smolSetSurfaceStyle( + sim: *mut simstruct, + surface: *const c_char, + face: PanelFace, + mode: DrawMode, + thickness: f64, + color: *mut f64, + stipplefactor: i32, + stipplepattern: i32, + shininess: f64, + ) -> ErrorCode; + + unsafe fn smolAddCompartment(sim: *mut simstruct, compartment: *const c_char) -> ErrorCode; + unsafe fn smolGetCompartmentIndex(sim: *mut simstruct, compartment: *const c_char) -> i32; + unsafe fn smolGetCompartmentName( + sim: *mut simstruct, + compartmentindex: i32, + compartment: *mut c_char, + ) -> *mut c_char; + unsafe fn smolAddCompartmentSurface( + sim: *mut simstruct, + compartment: *const c_char, + surface: *const c_char, + ) -> ErrorCode; + unsafe fn smolAddCompartmentPoint( + sim: *mut simstruct, + compartment: *const c_char, + point: *mut f64, + ) -> ErrorCode; + unsafe fn smolAddCompartmentLogic( + sim: *mut simstruct, + compartment: *const c_char, + logic: CmptLogic, + compartment2: *const c_char, + ) -> ErrorCode; + + unsafe fn smolAddReaction( + sim: *mut simstruct, + reaction: *const c_char, + reactant1: *const c_char, + rstate1: MolecState, + reactant2: *const c_char, + rstate2: MolecState, + nproduct: i32, + productspecies: *mut *const c_char, + productstates: *mut MolecState, + rate: f64, + ) -> ErrorCode; + unsafe fn smolGetReactionIndex( + sim: *mut simstruct, + orderptr: *mut i32, + reaction: *const c_char, + ) -> i32; + unsafe fn smolGetReactionName( + sim: *mut simstruct, + order: i32, + reactionindex: i32, + reaction: *mut c_char, + ) -> *mut c_char; + unsafe fn smolSetReactionRate( + sim: *mut simstruct, + reaction: *const c_char, + rate: f64, + kind: i32, + ) -> ErrorCode; + unsafe fn smolGetReactionRate( + sim: *mut simstruct, + reaction: *const c_char, + rate: *mut f64, + ) -> ErrorCode; + unsafe fn smolSetReactionRegion( + sim: *mut simstruct, + reaction: *const c_char, + compartment: *const c_char, + surface: *const c_char, + ) -> ErrorCode; + unsafe fn smolSetReactionProducts( + sim: *mut simstruct, + reaction: *const c_char, + method: RevParam, + parameter: f64, + product: *const c_char, + position: *mut f64, + ) -> ErrorCode; + + unsafe fn smolAddPort( + sim: *mut simstruct, + port: *const c_char, + surface: *const c_char, + face: PanelFace, + ) -> ErrorCode; + unsafe fn smolGetPortIndex(sim: *mut simstruct, port: *const c_char) -> i32; + unsafe fn smolGetPortName( + sim: *mut simstruct, + portindex: i32, + port: *mut c_char, + ) -> *mut c_char; + unsafe fn smolAddPortMolecules( + sim: *mut simstruct, + port: *const c_char, + nmolec: i32, + species: *const c_char, + positions: *mut *mut f64, + ) -> ErrorCode; + unsafe fn smolGetPortMolecules( + sim: *mut simstruct, + port: *const c_char, + species: *const c_char, + state: MolecState, + remove: i32, + ) -> i32; + + unsafe fn smolAddLattice( + sim: *mut simstruct, + lattice: *const c_char, + min: *const f64, + max: *const f64, + dx: *const f64, + btype: *const c_char, + ) -> ErrorCode; + unsafe fn smolAddLatticePort( + sim: *mut simstruct, + lattice: *const c_char, + port: *const c_char, + ) -> ErrorCode; + unsafe fn smolAddLatticeSpecies( + sim: *mut simstruct, + lattice: *const c_char, + species: *const c_char, + ) -> ErrorCode; + unsafe fn smolGetLatticeIndex(sim: *mut simstruct, lattice: *const c_char) -> i32; + unsafe fn smolGetLatticeName( + sim: *mut simstruct, + latticeindex: i32, + lattice: *mut c_char, + ) -> *mut c_char; + unsafe fn smolAddLatticeMolecules( + sim: *mut simstruct, + lattice: *const c_char, + species: *const c_char, + number: i32, + lowposition: *mut f64, + highposition: *mut f64, + ) -> ErrorCode; + unsafe fn smolAddLatticeReaction( + sim: *mut simstruct, + lattice: *const c_char, + reaction: *const c_char, + mov: i32, + ) -> ErrorCode; + } + + unsafe extern "C++" { + include!("smolrs.h"); + + fn smolrs_dim(sim: &simstruct) -> i32; + fn smolrs_time(sim: &simstruct) -> f64; + fn smolrs_time_start(sim: &simstruct) -> f64; + fn smolrs_time_stop(sim: &simstruct) -> f64; + fn smolrs_time_step(sim: &simstruct) -> f64; + fn smolrs_wall_pos(sim: &simstruct, d: i32, highside: i32) -> f64; + fn smolrs_nspecies(sim: &simstruct) -> i32; + fn smolrs_species_name(sim: &simstruct, i: i32) -> *const c_char; + fn smolrs_filepath_len(sim: &simstruct) -> usize; + unsafe fn smolrs_free_doubles(array: *mut f64); + } +} + +pub fn version() -> String { + let version = unsafe { crate::ffi::smolGetVersion() }; + + format!("{version}") +} + +pub fn set_debug_mode(on: bool) { + unsafe { ffi::smolSetDebugMode(on as i32) }; +} diff --git a/rust/libsmoldyn/src/sim/compartments.rs b/rust/libsmoldyn/src/sim/compartments.rs new file mode 100644 index 00000000..d8a9f27f --- /dev/null +++ b/rust/libsmoldyn/src/sim/compartments.rs @@ -0,0 +1,54 @@ +use super::{Sim, cstring, index, name, ok, to_i32, vector}; +use crate::error::SmolResult; +use crate::ffi::{self, CmptLogic}; + +impl Sim { + pub fn add_compartment(&mut self, compartment: &str) -> SmolResult<()> { + let compartment = cstring(compartment)?; + ok(unsafe { ffi::smolAddCompartment(self.p(), compartment.as_ptr()) }) + } + + pub fn compartment_index(&self, compartment: &str) -> SmolResult { + let compartment = cstring(compartment)?; + index(unsafe { ffi::smolGetCompartmentIndex(self.p(), compartment.as_ptr()) }) + } + + pub fn compartment_name(&self, index: usize) -> SmolResult { + let index = to_i32(index, "compartment index")?; + name(|buf| unsafe { ffi::smolGetCompartmentName(self.p(), index, buf) }) + } + + pub fn add_compartment_surface(&mut self, compartment: &str, surface: &str) -> SmolResult<()> { + let compartment = cstring(compartment)?; + let surface = cstring(surface)?; + ok(unsafe { + ffi::smolAddCompartmentSurface(self.p(), compartment.as_ptr(), surface.as_ptr()) + }) + } + + pub fn add_compartment_point(&mut self, compartment: &str, point: &[f64]) -> SmolResult<()> { + let mut point = vector(point, self.dim(), "point")?; + let compartment = cstring(compartment)?; + ok(unsafe { + ffi::smolAddCompartmentPoint(self.p(), compartment.as_ptr(), point.as_mut_ptr()) + }) + } + + pub fn add_compartment_logic( + &mut self, + compartment: &str, + logic: CmptLogic, + compartment2: &str, + ) -> SmolResult<()> { + let compartment = cstring(compartment)?; + let compartment2 = cstring(compartment2)?; + ok(unsafe { + ffi::smolAddCompartmentLogic( + self.p(), + compartment.as_ptr(), + logic, + compartment2.as_ptr(), + ) + }) + } +} diff --git a/rust/libsmoldyn/src/sim/graphics.rs b/rust/libsmoldyn/src/sim/graphics.rs new file mode 100644 index 00000000..f55df5e8 --- /dev/null +++ b/rust/libsmoldyn/src/sim/graphics.rs @@ -0,0 +1,111 @@ +use std::ffi::c_char; + +use super::{Sim, cstring, invalid, mut_cstring, ok, opt_cstring, opt_ptr}; +use crate::error::SmolResult; +use crate::ffi; + +fn opt4(v: &mut Option<[f64; 4]>) -> *mut f64 { + v.as_mut().map_or(std::ptr::null_mut(), |v| v.as_mut_ptr()) +} + +fn opt_count(n: Option) -> i32 { + n.map_or(-1, |n| n.min(i32::MAX as u32) as i32) +} + +impl Sim { + pub fn set_graphics_params( + &mut self, + method: &str, + timesteps: Option, + delay: Option, + ) -> SmolResult<()> { + let method = cstring(method)?; + ok(unsafe { + ffi::smolSetGraphicsParams( + self.p(), + method.as_ptr(), + opt_count(timesteps).max(0), + opt_count(delay), + ) + }) + } + + pub fn set_tiff_params( + &mut self, + timesteps: Option, + name: Option<&str>, + low_count: Option, + high_count: Option, + ) -> SmolResult<()> { + // libsmoldyn strcpy's the model's file path into a STRCHAR (256) buffer + // before appending the name. + if name.is_some() && ffi::smolrs_filepath_len(self.raw()) >= 256 { + return Err(invalid("model file path too long for a TIFF name")); + } + let name = opt_cstring(name)?; + ok(unsafe { + ffi::smolSetTiffParams( + self.p(), + opt_count(timesteps).max(0), + opt_ptr(&name), + opt_count(low_count), + opt_count(high_count), + ) + }) + } + + pub fn set_light_params( + &mut self, + index: i32, + mut ambient: Option<[f64; 4]>, + mut diffuse: Option<[f64; 4]>, + mut specular: Option<[f64; 4]>, + mut position: Option<[f64; 4]>, + ) -> SmolResult<()> { + ok(unsafe { + ffi::smolSetLightParams( + self.p(), + index, + opt4(&mut ambient), + opt4(&mut diffuse), + opt4(&mut specular), + opt4(&mut position), + ) + }) + } + + pub fn set_background_style(&mut self, color: [f64; 4]) -> SmolResult<()> { + let mut color = color; + ok(unsafe { ffi::smolSetBackgroundStyle(self.p(), color.as_mut_ptr()) }) + } + + pub fn set_frame_style( + &mut self, + thickness: Option, + mut color: Option<[f64; 4]>, + ) -> SmolResult<()> { + ok( + unsafe { + ffi::smolSetFrameStyle(self.p(), thickness.unwrap_or(-1.0), opt4(&mut color)) + }, + ) + } + + pub fn set_grid_style( + &mut self, + thickness: Option, + mut color: Option<[f64; 4]>, + ) -> SmolResult<()> { + ok(unsafe { ffi::smolSetGridStyle(self.p(), thickness.unwrap_or(-1.0), opt4(&mut color)) }) + } + + pub fn set_text_style(&mut self, color: [f64; 4]) -> SmolResult<()> { + let mut color = color; + ok(unsafe { ffi::smolSetTextStyle(self.p(), color.as_mut_ptr()) }) + } + + pub fn add_text_display(&mut self, item: &str) -> SmolResult<()> { + let mut item = mut_cstring(item)?; + ok(unsafe { ffi::smolAddTextDisplay(self.p(), item.as_mut_ptr() as *mut c_char) }) + } +} diff --git a/rust/libsmoldyn/src/sim/lattices.rs b/rust/libsmoldyn/src/sim/lattices.rs new file mode 100644 index 00000000..142f8896 --- /dev/null +++ b/rust/libsmoldyn/src/sim/lattices.rs @@ -0,0 +1,100 @@ +use super::{Sim, cstring, index, invalid, name, ok, opt_mut_ptr, opt_vector, to_i32, vector}; +use crate::error::SmolResult; +use crate::ffi; + +impl Sim { + pub fn add_lattice( + &mut self, + lattice: &str, + min: &[f64], + max: &[f64], + dx: &[f64], + btype: &str, + ) -> SmolResult<()> { + let dim = self.dim(); + let min = vector(min, dim, "min")?; + let max = vector(max, dim, "max")?; + let dx = vector(dx, dim, "dx")?; + if btype.len() < dim { + return Err(invalid(format!("btype needs {dim} characters"))); + } + let lattice = cstring(lattice)?; + let btype = cstring(btype)?; + ok(unsafe { + ffi::smolAddLattice( + self.p(), + lattice.as_ptr(), + min.as_ptr(), + max.as_ptr(), + dx.as_ptr(), + btype.as_ptr(), + ) + }) + } + + pub fn add_lattice_port(&mut self, lattice: &str, port: &str) -> SmolResult<()> { + let lattice = cstring(lattice)?; + let port = cstring(port)?; + ok(unsafe { ffi::smolAddLatticePort(self.p(), lattice.as_ptr(), port.as_ptr()) }) + } + + pub fn add_lattice_species(&mut self, lattice: &str, species: &str) -> SmolResult<()> { + let lattice = cstring(lattice)?; + let species = cstring(species)?; + ok(unsafe { ffi::smolAddLatticeSpecies(self.p(), lattice.as_ptr(), species.as_ptr()) }) + } + + pub fn lattice_index(&self, lattice: &str) -> SmolResult { + let lattice = cstring(lattice)?; + index(unsafe { ffi::smolGetLatticeIndex(self.p(), lattice.as_ptr()) }) + } + + pub fn lattice_name(&self, index: usize) -> SmolResult { + let index = to_i32(index, "lattice index")?; + name(|buf| unsafe { ffi::smolGetLatticeName(self.p(), index, buf) }) + } + + pub fn add_lattice_molecules( + &mut self, + lattice: &str, + species: &str, + number: usize, + low: Option<&[f64]>, + high: Option<&[f64]>, + ) -> SmolResult<()> { + let dim = self.dim(); + let mut low = opt_vector(low, dim, "low position")?; + let mut high = opt_vector(high, dim, "high position")?; + let lattice = cstring(lattice)?; + let species = cstring(species)?; + let number = to_i32(number, "number")?; + ok(unsafe { + ffi::smolAddLatticeMolecules( + self.p(), + lattice.as_ptr(), + species.as_ptr(), + number, + opt_mut_ptr(&mut low), + opt_mut_ptr(&mut high), + ) + }) + } + + pub fn add_lattice_reaction( + &mut self, + lattice: &str, + reaction: &str, + move_products: bool, + ) -> SmolResult<()> { + let lattice = cstring(lattice)?; + let reaction = cstring(reaction)?; + ok(unsafe { + ffi::smolAddLatticeReaction( + self.p(), + lattice.as_ptr(), + reaction.as_ptr(), + move_products as i32, + ) + }) + } +} diff --git a/rust/libsmoldyn/src/sim/mod.rs b/rust/libsmoldyn/src/sim/mod.rs new file mode 100644 index 00000000..38057d49 --- /dev/null +++ b/rust/libsmoldyn/src/sim/mod.rs @@ -0,0 +1,349 @@ +mod compartments; +mod graphics; +mod lattices; +mod molecules; +mod output; +mod ports; +mod reactions; +mod surfaces; + +pub use output::OutputData; +pub use reactions::RateKind; +pub use surfaces::{Side, SurfaceStyle}; + +use std::ffi::{CStr, CString, c_char}; +use std::marker::PhantomData; +use std::path::Path; +use std::ptr::NonNull; +use std::sync::atomic::{AtomicBool, Ordering}; + +use crate::error::{STRCHARLONG, SmolError, SmolResult, Status, check}; +use crate::ffi::{self, ErrorCode}; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Progress { + Running, + Finished, +} + +/// ```compile_fail +/// fn is_clone() {} +/// is_clone::(); +/// ``` +/// ```compile_fail +/// fn is_sync() {} +/// is_sync::(); +/// ``` +pub struct Sim { + ptr: NonNull, + outputs_open: bool, + _not_sync: PhantomData<*mut ()>, +} + +// SAFETY: `Sim` exclusively owns its `simstruct`, and nothing else keeps a +// pointer to it, so moving it to another thread is fine. +unsafe impl Send for Sim {} + +impl Sim { + // TODO: flags should be an enum + pub fn from_file(path: impl AsRef, flags: &str) -> SmolResult { + let (fileroot, filename) = split_path(path.as_ref())?; + let flags = cstring(flags)?; + let ptr = unsafe { + ffi::smolPrepareSimFromFile(fileroot.as_ptr(), filename.as_ptr(), flags.as_ptr()) + }; + Self::from_raw(ptr) + } + + pub fn load_from_file(path: impl AsRef, flags: &str) -> SmolResult { + let (fileroot, filename) = split_path(path.as_ref())?; + let flags = cstring(flags)?; + let mut ptr = std::ptr::null_mut(); + let code = unsafe { + ffi::smolLoadSimFromFile( + fileroot.as_ptr(), + filename.as_ptr(), + &mut ptr, + flags.as_ptr(), + ) + }; + let sim = Self::from_raw(ptr); + check(code)?; + sim + } + + pub fn new(low: &[f64], high: &[f64]) -> SmolResult { + if low.len() != high.len() || !(1..=3).contains(&low.len()) { + return Err(SmolError::InvalidArgument(format!( + "low and high must have the same length (1 to 3), got {} and {}", + low.len(), + high.len() + ))); + } + let mut low = low.to_vec(); + let mut high = high.to_vec(); + let ptr = unsafe { ffi::smolNewSim(low.len() as i32, low.as_mut_ptr(), high.as_mut_ptr()) }; + Self::from_raw(ptr) + } + + fn from_raw(ptr: *mut ffi::simstruct) -> SmolResult { + NonNull::new(ptr) + .map(|ptr| Sim { + ptr, + outputs_open: false, + _not_sync: PhantomData, + }) + .ok_or_else(|| SmolError::last(ErrorCode::Error)) + } + + /// # Safety + /// Don't free the pointer or use it after `self` is dropped or during a `&mut self` call. + pub unsafe fn as_ptr(&self) -> *mut ffi::simstruct { + self.ptr.as_ptr() + } + + fn p(&self) -> *mut ffi::simstruct { + self.ptr.as_ptr() + } + + fn raw(&self) -> &ffi::simstruct { + // SAFETY: `ptr` is valid for the lifetime of `self`, and `&self` + // guarantees no `&mut self` method is changing it. + unsafe { self.ptr.as_ref() } + } + + pub fn dim(&self) -> usize { + ffi::smolrs_dim(self.raw()) as usize + } + + pub fn time(&self) -> f64 { + ffi::smolrs_time(self.raw()) + } + + pub fn time_start(&self) -> f64 { + ffi::smolrs_time_start(self.raw()) + } + + pub fn time_stop(&self) -> f64 { + ffi::smolrs_time_stop(self.raw()) + } + + pub fn time_step(&self) -> f64 { + ffi::smolrs_time_step(self.raw()) + } + + pub fn bounds(&self) -> (Vec, Vec) { + let sim = self.raw(); + (0..self.dim() as i32) + .map(|d| { + ( + ffi::smolrs_wall_pos(sim, d, 0), + ffi::smolrs_wall_pos(sim, d, 1), + ) + }) + .unzip() + } + + pub fn update(&mut self) -> SmolResult<()> { + check(unsafe { ffi::smolUpdateSim(self.ptr.as_ptr()) }).map(drop) + } + + pub fn step(&mut self) -> SmolResult { + match check(unsafe { ffi::smolRunTimeStep(self.ptr.as_ptr()) })? { + Status::Notify => Ok(Progress::Finished), + Status::Ok | Status::Warning => Ok(Progress::Running), + } + } + + pub fn run(&mut self) -> SmolResult<()> { + self.run_until(self.time_stop()) + } + + pub fn run_with( + &mut self, + stop: &AtomicBool, + mut on_step: impl FnMut(&Sim), + ) -> SmolResult { + if !self.outputs_open { + check(unsafe { ffi::smolOpenOutputFiles(self.ptr.as_ptr(), 1) })?; + self.outputs_open = true; + } + while !stop.load(Ordering::Relaxed) { + let progress = self.step()?; + on_step(self); + if progress == Progress::Finished { + return Ok(Progress::Finished); + } + } + Ok(Progress::Running) + } + + pub fn run_until(&mut self, breaktime: f64) -> SmolResult<()> { + check(unsafe { ffi::smolRunSimUntil(self.ptr.as_ptr(), breaktime) }).map(drop) + } + + pub fn display(&self) -> SmolResult<()> { + check(unsafe { ffi::smolDisplaySim(self.ptr.as_ptr()) }).map(drop) + } + + pub fn set_times(&mut self, start: f64, stop: f64, step: f64) -> SmolResult<()> { + check(unsafe { ffi::smolSetSimTimes(self.ptr.as_ptr(), start, stop, step) }).map(drop) + } + + pub fn set_time_start(&mut self, start: f64) -> SmolResult<()> { + ok(unsafe { ffi::smolSetTimeStart(self.p(), start) }) + } + + pub fn set_time_stop(&mut self, stop: f64) -> SmolResult<()> { + ok(unsafe { ffi::smolSetTimeStop(self.p(), stop) }) + } + + pub fn set_time_now(&mut self, now: f64) -> SmolResult<()> { + check(unsafe { ffi::smolSetTimeNow(self.ptr.as_ptr(), now) }).map(drop) + } + + pub fn set_time_step(&mut self, step: f64) -> SmolResult<()> { + check(unsafe { ffi::smolSetTimeStep(self.ptr.as_ptr(), step) }).map(drop) + } + + pub fn set_flags(&mut self, flags: &str) -> SmolResult<()> { + // smolSetSimFlags strncpy's into a STRCHAR (256) buffer without terminating it. + if flags.len() >= 256 { + return Err(SmolError::InvalidArgument( + "flags must be < 256 bytes".into(), + )); + } + let flags = cstring(flags)?; + ok(unsafe { ffi::smolSetSimFlags(self.p(), flags.as_ptr()) }) + } + + pub fn set_random_seed(&mut self, seed: i64) -> SmolResult<()> { + ok(unsafe { ffi::smolSetRandomSeed(self.p(), seed) }) + } + + pub fn set_partitions(&mut self, method: &str, value: f64) -> SmolResult<()> { + let method = cstring(method)?; + ok(unsafe { ffi::smolSetPartitions(self.p(), method.as_ptr(), value) }) + } + + pub fn read_config(&mut self, statement: &str, parameters: &str) -> SmolResult<()> { + let statement = cstring(statement)?; + // libsmoldyn may modify `parameters` in place, so hand it a private copy. + let mut parameters = cstring(parameters)?.into_bytes_with_nul(); + check(unsafe { + ffi::smolReadConfigString( + self.ptr.as_ptr(), + statement.as_ptr(), + parameters.as_mut_ptr() as *mut c_char, + ) + }) + .map(drop) + } +} + +impl Drop for Sim { + fn drop(&mut self) { + unsafe { ffi::smolFreeSim(self.ptr.as_ptr()) }; + } +} + +impl std::fmt::Debug for Sim { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("Sim").field("ptr", &self.ptr).finish() + } +} + +// TODO: simplify path handling inside libsmoldyn later. +fn split_path(path: &Path) -> SmolResult<(CString, CString)> { + let filename = path + .file_name() + .ok_or_else(|| SmolError::InvalidArgument(format!("no file name in {path:?}")))?; + // libsmoldyn concatenates filepath and filename, so filepath must end + // with a separator (or be empty for the current directory). + let mut fileroot = match path.parent() { + Some(parent) if !parent.as_os_str().is_empty() => path_to_bytes(parent), + _ => Vec::new(), + }; + if !fileroot.is_empty() { + fileroot.push(std::path::MAIN_SEPARATOR as u8); + } + Ok((cstring(fileroot)?, cstring(path_to_bytes(filename))?)) +} + +fn invalid(msg: impl Into) -> SmolError { + SmolError::InvalidArgument(msg.into()) +} + +fn cstring(s: impl Into>) -> SmolResult { + CString::new(s).map_err(|e| SmolError::InvalidArgument(e.to_string())) +} + +fn opt_cstring(s: Option<&str>) -> SmolResult> { + s.map(cstring).transpose() +} + +fn opt_ptr(s: &Option) -> *const c_char { + s.as_ref().map_or(std::ptr::null(), |s| s.as_ptr()) +} + +fn mut_cstring(s: &str) -> SmolResult> { + Ok(cstring(s)?.into_bytes_with_nul()) +} + +fn ok(code: ErrorCode) -> SmolResult<()> { + check(code).map(drop) +} + +fn index(ret: i32) -> SmolResult { + if ret < 0 { + Err(SmolError::last(ErrorCode { repr: ret })) + } else { + Ok(ret as usize) + } +} + +fn count(ret: i32) -> SmolResult { + index(ret) +} + +fn to_i32(n: usize, what: &str) -> SmolResult { + i32::try_from(n).map_err(|_| invalid(format!("{what} is too large: {n}"))) +} + +fn name(get: impl FnOnce(*mut c_char) -> *mut c_char) -> SmolResult { + let mut buf = vec![0 as c_char; STRCHARLONG]; + let ret = get(buf.as_mut_ptr()); + if ret.is_null() { + return Err(SmolError::last(ErrorCode::Nonexist)); + } + Ok(unsafe { CStr::from_ptr(buf.as_ptr()) } + .to_string_lossy() + .into_owned()) +} + +fn vector(v: &[f64], n: usize, what: &str) -> SmolResult> { + if v.len() != n { + return Err(invalid(format!("{what} needs {n} values, got {}", v.len()))); + } + Ok(v.to_vec()) +} + +fn opt_vector(v: Option<&[f64]>, n: usize, what: &str) -> SmolResult>> { + v.map(|v| vector(v, n, what)).transpose() +} + +fn opt_mut_ptr(v: &mut Option>) -> *mut f64 { + v.as_mut().map_or(std::ptr::null_mut(), |v| v.as_mut_ptr()) +} + +// Thanks +#[cfg(unix)] +fn path_to_bytes>(path: P) -> Vec { + use std::os::unix::ffi::OsStrExt; + path.as_ref().as_os_str().as_bytes().to_vec() +} + +#[cfg(not(unix))] +fn path_to_bytes>(path: P) -> Vec { + path.as_ref().to_string_lossy().to_string().into_bytes() +} diff --git a/rust/libsmoldyn/src/sim/molecules.rs b/rust/libsmoldyn/src/sim/molecules.rs new file mode 100644 index 00000000..ddba01cb --- /dev/null +++ b/rust/libsmoldyn/src/sim/molecules.rs @@ -0,0 +1,231 @@ +use std::ffi::{CStr, c_char}; + +use super::{ + Sim, count, cstring, index, name, ok, opt_cstring, opt_mut_ptr, opt_ptr, opt_vector, to_i32, + vector, +}; +use crate::error::{STRCHARLONG, SmolError, SmolResult}; +use crate::ffi::{self, ErrorCode, MolecState, PanelShape}; + +impl Sim { + pub fn add_species(&mut self, species: &str, mollist: Option<&str>) -> SmolResult<()> { + let species = cstring(species)?; + let mollist = opt_cstring(mollist)?; + ok(unsafe { ffi::smolAddSpecies(self.p(), species.as_ptr(), opt_ptr(&mollist)) }) + } + + pub fn species(&self) -> Vec { + let sim = self.raw(); + (1..ffi::smolrs_nspecies(sim)) + .filter_map(|i| { + let name = ffi::smolrs_species_name(sim, i); + // SAFETY: non-null names are NUL-terminated strings owned by the simulation. + (!name.is_null()).then(|| { + unsafe { CStr::from_ptr(name) } + .to_string_lossy() + .into_owned() + }) + }) + .collect() + } + + pub fn species_index(&self, species: &str) -> SmolResult { + let species = cstring(species)?; + index(unsafe { ffi::smolGetSpeciesIndex(self.p(), species.as_ptr()) }) + } + + pub fn species_name(&self, index: usize) -> SmolResult { + let index = to_i32(index, "species index")?; + let mut buf = vec![0 as c_char; STRCHARLONG]; + unsafe { + ffi::smolClearError(); + ffi::smolGetSpeciesName(self.p(), index, buf.as_mut_ptr()); + } + // smolGetSpeciesName returns nothing; an empty buffer means it failed. + if buf[0] == 0 { + return Err(SmolError::last(ErrorCode::Nonexist)); + } + Ok(unsafe { CStr::from_ptr(buf.as_ptr()) } + .to_string_lossy() + .into_owned()) + } + + pub fn set_species_mobility( + &mut self, + species: &str, + state: MolecState, + difc: Option, + drift: Option<&[f64]>, + difmatrix: Option<&[f64]>, + ) -> SmolResult<()> { + let dim = self.dim(); + let species = cstring(species)?; + let mut drift = opt_vector(drift, dim, "drift")?; + let mut difmatrix = opt_vector(difmatrix, dim * dim, "difmatrix")?; + ok(unsafe { + ffi::smolSetSpeciesMobility( + self.p(), + species.as_ptr(), + state, + difc.unwrap_or(-1.0), + opt_mut_ptr(&mut drift), + opt_mut_ptr(&mut difmatrix), + ) + }) + } + + pub fn set_molecule_color( + &mut self, + species: &str, + state: MolecState, + color: [f64; 3], + ) -> SmolResult<()> { + let species = cstring(species)?; + let mut color = color; + ok(unsafe { + ffi::smolSetMoleculeColor(self.p(), species.as_ptr(), state, color.as_mut_ptr()) + }) + } + + pub fn set_molecule_size( + &mut self, + species: &str, + state: MolecState, + size: f64, + ) -> SmolResult<()> { + let species = cstring(species)?; + ok(unsafe { ffi::smolSetMoleculeSize(self.p(), species.as_ptr(), state, size) }) + } + + pub fn set_molecule_style( + &mut self, + species: &str, + state: MolecState, + size: Option, + color: Option<[f64; 3]>, + ) -> SmolResult<()> { + let species = cstring(species)?; + let mut color = color; + let color_ptr = color + .as_mut() + .map_or(std::ptr::null_mut(), |c| c.as_mut_ptr()); + ok(unsafe { + ffi::smolSetMoleculeStyle( + self.p(), + species.as_ptr(), + state, + size.unwrap_or(-1.0), + color_ptr, + ) + }) + } + + pub fn add_mol_list(&mut self, mollist: &str) -> SmolResult<()> { + let mollist = cstring(mollist)?; + ok(unsafe { ffi::smolAddMolList(self.p(), mollist.as_ptr()) }) + } + + pub fn mol_list_index(&self, mollist: &str) -> SmolResult { + let mollist = cstring(mollist)?; + index(unsafe { ffi::smolGetMolListIndex(self.p(), mollist.as_ptr()) }) + } + + pub fn mol_list_name(&self, index: usize) -> SmolResult { + let index = to_i32(index, "molecule list index")?; + name(|buf| unsafe { ffi::smolGetMolListName(self.p(), index, buf) }) + } + + pub fn set_mol_list( + &mut self, + species: &str, + state: MolecState, + mollist: &str, + ) -> SmolResult<()> { + let species = cstring(species)?; + let mollist = cstring(mollist)?; + ok(unsafe { ffi::smolSetMolList(self.p(), species.as_ptr(), state, mollist.as_ptr()) }) + } + + pub fn set_max_molecules(&mut self, max: usize) -> SmolResult<()> { + let max = to_i32(max, "max molecules")?; + ok(unsafe { ffi::smolSetMaxMolecules(self.p(), max) }) + } + + pub fn add_solution_molecules( + &mut self, + species: &str, + number: usize, + low: Option<&[f64]>, + high: Option<&[f64]>, + ) -> SmolResult<()> { + let dim = self.dim(); + let species = cstring(species)?; + let number = to_i32(number, "number")?; + let mut low = opt_vector(low, dim, "low position")?; + let mut high = opt_vector(high, dim, "high position")?; + ok(unsafe { + ffi::smolAddSolutionMolecules( + self.p(), + species.as_ptr(), + number, + opt_mut_ptr(&mut low), + opt_mut_ptr(&mut high), + ) + }) + } + + pub fn add_compartment_molecules( + &mut self, + species: &str, + number: usize, + compartment: &str, + ) -> SmolResult<()> { + let species = cstring(species)?; + let number = to_i32(number, "number")?; + let compartment = cstring(compartment)?; + ok(unsafe { + ffi::smolAddCompartmentMolecules( + self.p(), + species.as_ptr(), + number, + compartment.as_ptr(), + ) + }) + } + + #[allow(clippy::too_many_arguments)] + pub fn add_surface_molecules( + &mut self, + species: &str, + state: MolecState, + number: usize, + surface: &str, + shape: PanelShape, + panel: &str, + position: Option<&[f64]>, + ) -> SmolResult<()> { + let dim = self.dim(); + let species = cstring(species)?; + let number = to_i32(number, "number")?; + let surface = cstring(surface)?; + let panel = cstring(panel)?; + let mut position = position.map(|p| vector(p, dim, "position")).transpose()?; + ok(unsafe { + ffi::smolAddSurfaceMolecules( + self.p(), + species.as_ptr(), + state, + number, + surface.as_ptr(), + shape, + panel.as_ptr(), + opt_mut_ptr(&mut position), + ) + }) + } + + pub fn molecule_count(&self, species: &str, state: MolecState) -> SmolResult { + let species = cstring(species)?; + count(unsafe { ffi::smolGetMoleculeCount(self.p(), species.as_ptr(), state) }) + } +} diff --git a/rust/libsmoldyn/src/sim/output.rs b/rust/libsmoldyn/src/sim/output.rs new file mode 100644 index 00000000..b0535e29 --- /dev/null +++ b/rust/libsmoldyn/src/sim/output.rs @@ -0,0 +1,120 @@ +use std::ffi::c_char; + +use super::{Sim, cstring, mut_cstring, ok}; +use crate::error::{SmolError, SmolResult}; +use crate::ffi; + +#[derive(Debug, Clone, Default, PartialEq)] +pub struct OutputData { + pub nrow: usize, + pub ncol: usize, + pub data: Vec, +} + +impl OutputData { + pub fn row(&self, i: usize) -> Option<&[f64]> { + (i < self.nrow).then(|| &self.data[i * self.ncol..(i + 1) * self.ncol]) + } + + pub fn rows(&self) -> impl Iterator { + self.data.chunks(self.ncol.max(1)).take(self.nrow) + } +} + +impl Sim { + pub fn set_output_path(&mut self, path: &str) -> SmolResult<()> { + let path = cstring(path)?; + ok(unsafe { ffi::smolSetOutputPath(self.p(), path.as_ptr()) }) + } + + pub fn add_output_file( + &mut self, + filename: &str, + suffix: Option, + append: bool, + ) -> SmolResult<()> { + let mut filename = mut_cstring(filename)?; + let suffix = suffix.map_or(-1, |s| s.min(i32::MAX as u32) as i32); + ok(unsafe { + ffi::smolAddOutputFile( + self.p(), + filename.as_mut_ptr() as *mut c_char, + suffix, + append as i32, + ) + }) + } + + pub fn add_output_data(&mut self, dataname: &str) -> SmolResult<()> { + let mut dataname = mut_cstring(dataname)?; + ok(unsafe { ffi::smolAddOutputData(self.p(), dataname.as_mut_ptr() as *mut c_char) }) + } + + pub fn open_output_files(&mut self, overwrite: bool) -> SmolResult<()> { + ok(unsafe { ffi::smolOpenOutputFiles(self.p(), overwrite as i32) })?; + self.outputs_open = true; + Ok(()) + } + + pub fn add_command(&mut self, command: &str) -> SmolResult<()> { + let mut command = mut_cstring(command)?; + ok(unsafe { ffi::smolAddCommandFromString(self.p(), command.as_mut_ptr() as *mut c_char) }) + } + + pub fn add_timed_command( + &mut self, + kind: char, + on: f64, + off: f64, + step: f64, + multiplier: f64, + command: &str, + ) -> SmolResult<()> { + if !kind.is_ascii() { + return Err(SmolError::InvalidArgument(format!( + "invalid command type {kind:?}" + ))); + } + let command = cstring(command)?; + ok(unsafe { + ffi::smolAddCommand( + self.p(), + kind as u8 as c_char, + on, + off, + step, + multiplier, + command.as_ptr(), + ) + }) + } + + pub fn run_command(&mut self, command: &str) -> SmolResult<()> { + let command = cstring(command)?; + ok(unsafe { ffi::smolRunCommand(self.p(), command.as_ptr()) }) + } + + pub fn output_data(&mut self, dataname: &str, erase: bool) -> SmolResult { + let mut dataname = mut_cstring(dataname)?; + let (mut nrow, mut ncol) = (0i32, 0i32); + let mut array: *mut f64 = std::ptr::null_mut(); + ok(unsafe { + ffi::smolGetOutputData( + self.p(), + dataname.as_mut_ptr() as *mut c_char, + &mut nrow, + &mut ncol, + &mut array, + erase as i32, + ) + })?; + if array.is_null() { + return Ok(OutputData::default()); + } + let (nrow, ncol) = (nrow.max(0) as usize, ncol.max(0) as usize); + // SAFETY: libsmoldyn calloc'd nrow*ncol doubles for us to free. + let data = unsafe { std::slice::from_raw_parts(array, nrow * ncol) }.to_vec(); + unsafe { ffi::smolrs_free_doubles(array) }; + Ok(OutputData { nrow, ncol, data }) + } +} diff --git a/rust/libsmoldyn/src/sim/ports.rs b/rust/libsmoldyn/src/sim/ports.rs new file mode 100644 index 00000000..85cb8621 --- /dev/null +++ b/rust/libsmoldyn/src/sim/ports.rs @@ -0,0 +1,77 @@ +use super::{Sim, count, cstring, index, invalid, name, ok, to_i32, vector}; +use crate::error::SmolResult; +use crate::ffi::{self, MolecState, PanelFace}; + +impl Sim { + pub fn add_port(&mut self, port: &str, surface: &str, face: PanelFace) -> SmolResult<()> { + let port = cstring(port)?; + let surface = cstring(surface)?; + ok(unsafe { ffi::smolAddPort(self.p(), port.as_ptr(), surface.as_ptr(), face) }) + } + + pub fn port_index(&self, port: &str) -> SmolResult { + let port = cstring(port)?; + index(unsafe { ffi::smolGetPortIndex(self.p(), port.as_ptr()) }) + } + + pub fn port_name(&self, index: usize) -> SmolResult { + let index = to_i32(index, "port index")?; + name(|buf| unsafe { ffi::smolGetPortName(self.p(), index, buf) }) + } + + pub fn add_port_molecules( + &mut self, + port: &str, + species: &str, + number: usize, + positions: Option<&[&[f64]]>, + ) -> SmolResult<()> { + let dim = self.dim(); + let port = cstring(port)?; + let species = cstring(species)?; + let nmolec = to_i32(number, "number")?; + let mut positions = match positions { + Some(p) if p.len() != number => { + return Err(invalid(format!( + "{number} molecules but {} positions", + p.len() + ))); + } + Some(p) => Some( + p.iter() + .map(|p| vector(p, dim, "position")) + .collect::>>()?, + ), + None => None, + }; + let mut rows: Option> = positions + .as_mut() + .map(|p| p.iter_mut().map(|row| row.as_mut_ptr()).collect()); + let rows_ptr = rows + .as_mut() + .map_or(std::ptr::null_mut(), |r| r.as_mut_ptr()); + ok(unsafe { + ffi::smolAddPortMolecules(self.p(), port.as_ptr(), nmolec, species.as_ptr(), rows_ptr) + }) + } + + pub fn port_molecules( + &mut self, + port: &str, + species: &str, + state: MolecState, + remove: bool, + ) -> SmolResult { + let port = cstring(port)?; + let species = cstring(species)?; + count(unsafe { + ffi::smolGetPortMolecules( + self.p(), + port.as_ptr(), + species.as_ptr(), + state, + remove as i32, + ) + }) + } +} diff --git a/rust/libsmoldyn/src/sim/reactions.rs b/rust/libsmoldyn/src/sim/reactions.rs new file mode 100644 index 00000000..e3d19f83 --- /dev/null +++ b/rust/libsmoldyn/src/sim/reactions.rs @@ -0,0 +1,143 @@ +use std::ffi::{CString, c_char}; + +use super::{ + Sim, cstring, index, invalid, name, ok, opt_cstring, opt_mut_ptr, opt_ptr, opt_vector, to_i32, +}; +use crate::error::SmolResult; +use crate::ffi::{self, MolecState, RevParam}; + +// MAXPRODUCT in smoldyn.h. +const MAX_PRODUCT: usize = 256; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum RateKind { + Rate, + Internal, + Probability, +} + +impl Sim { + pub fn add_reaction( + &mut self, + reaction: &str, + reactants: &[(&str, MolecState)], + products: &[(&str, MolecState)], + rate: Option, + ) -> SmolResult<()> { + if reactants.len() > 2 { + return Err(invalid("a reaction has at most 2 reactants")); + } + if products.len() > MAX_PRODUCT { + return Err(invalid(format!( + "a reaction has at most {MAX_PRODUCT} products" + ))); + } + let reaction = cstring(reaction)?; + let reactant = |i: usize| -> SmolResult<(Option, MolecState)> { + match reactants.get(i) { + Some(&(species, state)) => Ok((Some(cstring(species)?), state)), + None => Ok((None, MolecState::None)), + } + }; + let (r1, s1) = reactant(0)?; + let (r2, s2) = reactant(1)?; + let names = products + .iter() + .map(|&(species, _)| cstring(species)) + .collect::>>()?; + let mut name_ptrs: Vec<*const c_char> = names.iter().map(|n| n.as_ptr()).collect(); + let mut states: Vec = products.iter().map(|&(_, state)| state).collect(); + ok(unsafe { + ffi::smolAddReaction( + self.p(), + reaction.as_ptr(), + opt_ptr(&r1), + s1, + opt_ptr(&r2), + s2, + products.len() as i32, + name_ptrs.as_mut_ptr(), + states.as_mut_ptr(), + rate.unwrap_or(-1.0), + ) + }) + } + + pub fn reaction_index(&self, reaction: &str) -> SmolResult<(usize, usize)> { + let reaction = cstring(reaction)?; + let mut order = -1; + let i = + index(unsafe { ffi::smolGetReactionIndex(self.p(), &mut order, reaction.as_ptr()) })?; + Ok((order as usize, i)) + } + + pub fn reaction_name(&self, order: usize, index: usize) -> SmolResult { + let order = to_i32(order, "reaction order")?; + let index = to_i32(index, "reaction index")?; + name(|buf| unsafe { ffi::smolGetReactionName(self.p(), order, index, buf) }) + } + + pub fn set_reaction_rate( + &mut self, + reaction: &str, + rate: f64, + kind: RateKind, + ) -> SmolResult<()> { + let reaction = cstring(reaction)?; + let kind = match kind { + RateKind::Rate => 0, + RateKind::Internal => 1, + RateKind::Probability => 2, + }; + ok(unsafe { ffi::smolSetReactionRate(self.p(), reaction.as_ptr(), rate, kind) }) + } + + pub fn reaction_rate(&self, reaction: &str) -> SmolResult { + let reaction = cstring(reaction)?; + let mut rate = 0.0; + ok(unsafe { ffi::smolGetReactionRate(self.p(), reaction.as_ptr(), &mut rate) })?; + Ok(rate) + } + + pub fn set_reaction_region( + &mut self, + reaction: &str, + compartment: Option<&str>, + surface: Option<&str>, + ) -> SmolResult<()> { + let reaction = cstring(reaction)?; + let compartment = opt_cstring(compartment)?; + let surface = opt_cstring(surface)?; + ok(unsafe { + ffi::smolSetReactionRegion( + self.p(), + reaction.as_ptr(), + opt_ptr(&compartment), + opt_ptr(&surface), + ) + }) + } + + pub fn set_reaction_products( + &mut self, + reaction: &str, + method: RevParam, + parameter: f64, + product: Option<&str>, + position: Option<&[f64]>, + ) -> SmolResult<()> { + let mut position = opt_vector(position, self.dim(), "position")?; + let reaction = cstring(reaction)?; + let product = opt_cstring(product)?; + ok(unsafe { + ffi::smolSetReactionProducts( + self.p(), + reaction.as_ptr(), + method, + parameter, + opt_ptr(&product), + opt_mut_ptr(&mut position), + ) + }) + } +} diff --git a/rust/libsmoldyn/src/sim/surfaces.rs b/rust/libsmoldyn/src/sim/surfaces.rs new file mode 100644 index 00000000..b3f1277a --- /dev/null +++ b/rust/libsmoldyn/src/sim/surfaces.rs @@ -0,0 +1,301 @@ +use super::{Sim, cstring, index, invalid, name, ok, opt_cstring, opt_ptr, to_i32, vector}; +use crate::error::SmolResult; +use crate::ffi::{self, DrawMode, MolecState, PanelFace, PanelShape, SrfAction}; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Side { + Low, + High, + Both, +} + +#[derive(Debug, Clone, Copy)] +pub struct SurfaceStyle { + pub mode: DrawMode, + pub thickness: Option, + pub color: Option<[f64; 4]>, + pub stipple_factor: Option, + pub stipple_pattern: Option, + pub shininess: Option, +} + +impl Default for SurfaceStyle { + fn default() -> Self { + SurfaceStyle { + mode: DrawMode::None, + thickness: None, + color: None, + stipple_factor: None, + stipple_pattern: None, + shininess: None, + } + } +} + +// Enough for the parameters of any panel shape (at most 9, for 3D triangles, +// cylinders and hemispheres). smolAddPanel reads a shape-dependent number of +// values, so we always pass a buffer this large and missing values read as 0 +// instead of out of bounds. +const MAX_PANEL_PARAMS: usize = 16; + +// MAXPRODUCT in smoldyn.h: RxnSetIntersurfaceRules reads one rule per product. +const MAX_PRODUCT: usize = 256; + +impl Sim { + pub fn set_boundary_type(&mut self, dim: usize, side: Side, kind: char) -> SmolResult<()> { + if dim >= self.dim() { + return Err(invalid(format!("dimension {dim} out of range"))); + } + if !matches!(kind, 'r' | 'p' | 'a' | 't') { + return Err(invalid(format!("invalid boundary type {kind:?}"))); + } + let highside = match side { + Side::Low => 0, + Side::High => 1, + Side::Both => -1, + }; + ok(unsafe { ffi::smolSetBoundaryType(self.p(), dim as i32, highside, kind as u8 as _) }) + } + + pub fn add_surface(&mut self, surface: &str) -> SmolResult<()> { + let surface = cstring(surface)?; + ok(unsafe { ffi::smolAddSurface(self.p(), surface.as_ptr()) }) + } + + pub fn surface_index(&self, surface: &str) -> SmolResult { + let surface = cstring(surface)?; + index(unsafe { ffi::smolGetSurfaceIndex(self.p(), surface.as_ptr()) }) + } + + pub fn surface_name(&self, index: usize) -> SmolResult { + let index = to_i32(index, "surface index")?; + name(|buf| unsafe { ffi::smolGetSurfaceName(self.p(), index, buf) }) + } + + pub fn set_reaction_intersurface(&mut self, reaction: &str, rules: &[i32]) -> SmolResult<()> { + if rules.len() > MAX_PRODUCT { + return Err(invalid(format!("at most {MAX_PRODUCT} rules"))); + } + let reaction = cstring(reaction)?; + // libsmoldyn reads as many rules as the reaction has products. + let mut buf = [0i32; MAX_PRODUCT]; + if rules.is_empty() { + buf[0] = -1; + } else { + buf[..rules.len()].copy_from_slice(rules); + } + ok(unsafe { + ffi::smolSetReactionIntersurface(self.p(), reaction.as_ptr(), buf.as_mut_ptr()) + }) + } + + pub fn set_surface_action( + &mut self, + surface: &str, + face: PanelFace, + species: &str, + state: MolecState, + action: SrfAction, + new_species: Option<&str>, + ) -> SmolResult<()> { + let surface = cstring(surface)?; + let species = cstring(species)?; + let new_species = opt_cstring(new_species)?; + ok(unsafe { + ffi::smolSetSurfaceAction( + self.p(), + surface.as_ptr(), + face, + species.as_ptr(), + state, + action, + opt_ptr(&new_species), + ) + }) + } + + #[allow(clippy::too_many_arguments)] + pub fn set_surface_rate( + &mut self, + surface: &str, + species: &str, + state: MolecState, + state1: MolecState, + state2: MolecState, + rate: f64, + new_species: Option<&str>, + internal: bool, + ) -> SmolResult<()> { + let surface = cstring(surface)?; + let species = cstring(species)?; + let new_species = opt_cstring(new_species)?; + ok(unsafe { + ffi::smolSetSurfaceRate( + self.p(), + surface.as_ptr(), + species.as_ptr(), + state, + state1, + state2, + rate, + opt_ptr(&new_species), + internal as i32, + ) + }) + } + + pub fn add_panel( + &mut self, + surface: &str, + shape: PanelShape, + panel: Option<&str>, + axis: Option<&str>, + params: &[f64], + ) -> SmolResult<()> { + if params.len() > MAX_PANEL_PARAMS { + return Err(invalid(format!( + "at most {MAX_PANEL_PARAMS} panel parameters" + ))); + } + let surface = cstring(surface)?; + let panel = opt_cstring(panel)?; + let axis = opt_cstring(axis)?; + let mut buf = [0.0; MAX_PANEL_PARAMS]; + buf[..params.len()].copy_from_slice(params); + ok(unsafe { + ffi::smolAddPanel( + self.p(), + surface.as_ptr(), + shape, + opt_ptr(&panel), + opt_ptr(&axis), + buf.as_mut_ptr(), + ) + }) + } + + pub fn panel_index(&self, surface: &str, panel: &str) -> SmolResult<(PanelShape, usize)> { + let surface = cstring(surface)?; + let panel = cstring(panel)?; + let mut shape = PanelShape::None; + let i = index(unsafe { + ffi::smolGetPanelIndex(self.p(), surface.as_ptr(), &mut shape, panel.as_ptr()) + })?; + Ok((shape, i)) + } + + pub fn panel_name(&self, surface: &str, shape: PanelShape, index: usize) -> SmolResult { + let surface = cstring(surface)?; + let index = to_i32(index, "panel index")?; + name(|buf| unsafe { ffi::smolGetPanelName(self.p(), surface.as_ptr(), shape, index, buf) }) + } + + pub fn set_panel_jump( + &mut self, + surface: &str, + panel1: &str, + face1: PanelFace, + panel2: &str, + face2: PanelFace, + bidirectional: bool, + ) -> SmolResult<()> { + let surface = cstring(surface)?; + let panel1 = cstring(panel1)?; + let panel2 = cstring(panel2)?; + ok(unsafe { + ffi::smolSetPanelJump( + self.p(), + surface.as_ptr(), + panel1.as_ptr(), + face1, + panel2.as_ptr(), + face2, + bidirectional as i32, + ) + }) + } + + pub fn add_surface_unbounded_emitter( + &mut self, + surface: &str, + face: PanelFace, + species: &str, + amount: f64, + position: &[f64], + ) -> SmolResult<()> { + let mut position = vector(position, self.dim(), "position")?; + let surface = cstring(surface)?; + let species = cstring(species)?; + ok(unsafe { + ffi::smolAddSurfaceUnboundedEmitter( + self.p(), + surface.as_ptr(), + face, + species.as_ptr(), + amount, + position.as_mut_ptr(), + ) + }) + } + + pub fn set_surface_sim_params(&mut self, parameter: &str, value: f64) -> SmolResult<()> { + let parameter = cstring(parameter)?; + ok(unsafe { ffi::smolSetSurfaceSimParams(self.p(), parameter.as_ptr(), value) }) + } + + pub fn add_panel_neighbor( + &mut self, + surface1: &str, + panel1: &str, + surface2: &str, + panel2: &str, + reciprocal: bool, + ) -> SmolResult<()> { + let (surface1, panel1) = (cstring(surface1)?, cstring(panel1)?); + let (surface2, panel2) = (cstring(surface2)?, cstring(panel2)?); + ok(unsafe { + ffi::smolAddPanelNeighbor( + self.p(), + surface1.as_ptr(), + panel1.as_ptr(), + surface2.as_ptr(), + panel2.as_ptr(), + reciprocal as i32, + ) + }) + } + + pub fn set_surface_style( + &mut self, + surface: &str, + face: PanelFace, + style: &SurfaceStyle, + ) -> SmolResult<()> { + // libsmoldyn accepts "all" here but then indexes srflist[-5]. + if surface == "all" { + return Err(invalid( + "set_surface_style needs a specific surface, not \"all\"", + )); + } + let surface = cstring(surface)?; + let mut color = style.color; + let color_ptr = color + .as_mut() + .map_or(std::ptr::null_mut(), |c| c.as_mut_ptr()); + ok(unsafe { + ffi::smolSetSurfaceStyle( + self.p(), + surface.as_ptr(), + face, + style.mode, + style.thickness.unwrap_or(-1.0), + color_ptr, + style + .stipple_factor + .map_or(-1, |f| f.min(i32::MAX as u32) as i32), + style.stipple_pattern.map_or(-1, i32::from), + style.shininess.unwrap_or(-1.0), + ) + }) + } +} diff --git a/rust/libsmoldyn/tests/api.rs b/rust/libsmoldyn/tests/api.rs new file mode 100644 index 00000000..0acdecc1 --- /dev/null +++ b/rust/libsmoldyn/tests/api.rs @@ -0,0 +1,449 @@ +use std::fs; +use std::path::PathBuf; +use std::sync::{Mutex, MutexGuard}; + +use libsmoldyn::{ + CmptLogic, MolecState, PanelFace, PanelShape, RateKind, RevParam, Side, Sim, SmolError, + SrfAction, SurfaceStyle, +}; + +// libsmoldyn keeps global error state, so tests must not run concurrently. +static LOCK: Mutex<()> = Mutex::new(()); + +fn lock() -> MutexGuard<'static, ()> { + LOCK.lock().unwrap_or_else(|e| e.into_inner()) +} + +fn box_sim() -> Sim { + let mut sim = Sim::new(&[-1.0; 3], &[1.0; 3]).unwrap(); + sim.set_times(0.0, 0.01, 0.002).unwrap(); + sim.set_random_seed(1).unwrap(); + sim +} + +#[test] +fn settings() { + let _guard = lock(); + let mut sim = box_sim(); + sim.set_time_start(0.5).unwrap(); + sim.set_time_stop(2.0).unwrap(); + sim.set_time_step(0.25).unwrap(); + assert_eq!( + (sim.time_start(), sim.time_stop(), sim.time_step()), + (0.5, 2.0, 0.25) + ); + sim.set_flags("q").unwrap(); + assert!(matches!( + sim.set_flags(&"q".repeat(300)), + Err(SmolError::InvalidArgument(_)) + )); + sim.set_partitions("molperbox", 4.0).unwrap(); + assert!(sim.set_partitions("nonsense", 4.0).is_err()); + for side in [Side::Low, Side::High, Side::Both] { + sim.set_boundary_type(0, side, 'r').unwrap(); + } + assert!(sim.set_boundary_type(3, Side::Both, 'r').is_err()); + assert!(sim.set_boundary_type(0, Side::Both, 'x').is_err()); +} + +#[test] +fn molecules() { + let _guard = lock(); + let mut sim = box_sim(); + sim.add_species("A", None).unwrap(); + sim.add_mol_list("slow").unwrap(); + sim.add_species("B", Some("slow")).unwrap(); + + assert_eq!(sim.species_index("A").unwrap(), 1); + assert_eq!(sim.species_index("B").unwrap(), 2); + assert!(sim.species_index("C").is_err()); + assert_eq!(sim.species_name(2).unwrap(), "B"); + assert!(sim.species_name(10).is_err()); + + let slow = sim.mol_list_index("slow").unwrap(); + assert_eq!(sim.mol_list_name(slow).unwrap(), "slow"); + sim.set_mol_list("A", MolecState::Soln, "slow").unwrap(); + + sim.set_species_mobility( + "A", + MolecState::All, + Some(1.0), + Some(&[0.1, 0.0, 0.0]), + None, + ) + .unwrap(); + assert!(matches!( + sim.set_species_mobility("A", MolecState::All, None, Some(&[0.1]), None), + Err(SmolError::InvalidArgument(_)) + )); + sim.set_molecule_color("A", MolecState::All, [1.0, 0.0, 0.0]) + .unwrap(); + sim.set_molecule_size("A", MolecState::All, 2.0).unwrap(); + sim.set_molecule_style("B", MolecState::All, Some(3.0), None) + .unwrap(); + assert!( + sim.set_molecule_color("A", MolecState::All, [2.0, 0.0, 0.0]) + .is_err() + ); + + sim.set_max_molecules(1000).unwrap(); + sim.add_solution_molecules("A", 10, None, None).unwrap(); + sim.add_solution_molecules("B", 4, Some(&[0.0; 3]), Some(&[0.5; 3])) + .unwrap(); + assert!( + sim.add_solution_molecules("B", 1, Some(&[0.0]), None) + .is_err() + ); + sim.update().unwrap(); + + assert_eq!(sim.molecule_count("A", MolecState::Soln).unwrap(), 10); + assert_eq!(sim.molecule_count("B", MolecState::All).unwrap(), 4); + assert_eq!(sim.molecule_count("all", MolecState::All).unwrap(), 14); + sim.run().unwrap(); +} + +fn ball_sim() -> Sim { + let mut sim = box_sim(); + sim.add_species("A", None).unwrap(); + sim.add_species("B", None).unwrap(); + sim.add_surface("ball").unwrap(); + sim.add_panel( + "ball", + PanelShape::Sph, + Some("s1"), + None, + &[0.0, 0.0, 0.0, 0.5, 10.0, 10.0], + ) + .unwrap(); + sim.set_surface_action( + "ball", + PanelFace::Both, + "all", + MolecState::Soln, + SrfAction::Reflect, + None, + ) + .unwrap(); + sim.add_compartment("inside").unwrap(); + sim.add_compartment_surface("inside", "ball").unwrap(); + sim.add_compartment_point("inside", &[0.0, 0.0, 0.0]) + .unwrap(); + sim +} + +#[test] +fn surfaces_and_compartments() { + let _guard = lock(); + let mut sim = ball_sim(); + + assert_eq!(sim.surface_index("ball").unwrap(), 0); + assert_eq!(sim.surface_name(0).unwrap(), "ball"); + assert!(sim.surface_index("nope").is_err()); + assert_eq!(sim.panel_index("ball", "s1").unwrap(), (PanelShape::Sph, 0)); + assert_eq!(sim.panel_name("ball", PanelShape::Sph, 0).unwrap(), "s1"); + assert!(sim.panel_name("ball", PanelShape::Sph, 3).is_err()); + assert!(sim.add_compartment_point("inside", &[0.0]).is_err()); + + assert_eq!(sim.compartment_index("inside").unwrap(), 0); + assert_eq!(sim.compartment_name(0).unwrap(), "inside"); + sim.add_compartment("outside").unwrap(); + sim.add_compartment_logic("outside", CmptLogic::Equalnot, "inside") + .unwrap(); + + sim.set_surface_rate( + "ball", + "B", + MolecState::Soln, + MolecState::Soln, + MolecState::Front, + 0.1, + None, + false, + ) + .unwrap(); + sim.set_surface_sim_params("epsilon", 1e-6).unwrap(); + sim.set_surface_style( + "ball", + PanelFace::Both, + &SurfaceStyle { + thickness: Some(1.0), + color: Some([0.0, 0.0, 1.0, 1.0]), + ..Default::default() + }, + ) + .unwrap(); + assert!(matches!( + sim.set_surface_style("all", PanelFace::Both, &SurfaceStyle::default()), + Err(SmolError::InvalidArgument(_)) + )); + + sim.add_compartment_molecules("A", 20, "inside").unwrap(); + sim.add_surface_molecules( + "B", + MolecState::Front, + 5, + "ball", + PanelShape::All, + "all", + None, + ) + .unwrap(); + sim.update().unwrap(); + assert_eq!(sim.molecule_count("A", MolecState::Soln).unwrap(), 20); + assert_eq!(sim.molecule_count("B", MolecState::Front).unwrap(), 5); + + sim.set_species_mobility("A", MolecState::All, Some(1.0), None, None) + .unwrap(); + sim.run().unwrap(); + assert_eq!(sim.molecule_count("A", MolecState::Soln).unwrap(), 20); +} + +#[test] +fn panel_jump_and_neighbor() { + let _guard = lock(); + let mut sim = box_sim(); + sim.add_surface("walls").unwrap(); + sim.add_panel( + "walls", + PanelShape::Rect, + Some("r1"), + Some("+0"), + &[-0.5, -1.0, -1.0, 2.0, 2.0], + ) + .unwrap(); + sim.add_panel( + "walls", + PanelShape::Rect, + Some("r2"), + Some("-0"), + &[0.5, -1.0, -1.0, 2.0, 2.0], + ) + .unwrap(); + assert!( + sim.add_panel("walls", PanelShape::Rect, Some("r3"), None, &[0.0; 5]) + .is_err() + ); + assert_eq!( + sim.panel_index("walls", "r2").unwrap(), + (PanelShape::Rect, 1) + ); + sim.set_panel_jump("walls", "r1", PanelFace::Front, "r2", PanelFace::Back, true) + .unwrap(); + sim.add_panel_neighbor("walls", "r1", "walls", "r2", true) + .unwrap(); + sim.add_species("A", None).unwrap(); + sim.add_surface_unbounded_emitter("walls", PanelFace::Front, "A", 1.0, &[0.0; 3]) + .unwrap(); + sim.update().unwrap(); +} + +#[test] +fn reactions() { + let _guard = lock(); + let mut sim = box_sim(); + for s in ["A", "B", "C"] { + sim.add_species(s, None).unwrap(); + } + // a bimolecular reaction needs diffusing reactants to get a binding radius. + sim.set_species_mobility("all", MolecState::All, Some(1.0), None, None) + .unwrap(); + sim.add_reaction( + "bind", + &[("A", MolecState::Soln), ("B", MolecState::Soln)], + &[("C", MolecState::Soln)], + Some(10.0), + ) + .unwrap(); + sim.add_reaction("decay", &[("C", MolecState::Soln)], &[], None) + .unwrap(); + sim.add_reaction("make", &[], &[("A", MolecState::Soln)], Some(0.0)) + .unwrap(); + assert!( + sim.add_reaction("bad", &[("X", MolecState::Soln)], &[], None) + .is_err() + ); + assert!(matches!( + sim.add_reaction("three", &[("A", MolecState::Soln); 3], &[], None), + Err(SmolError::InvalidArgument(_)) + )); + + assert_eq!(sim.reaction_index("bind").unwrap(), (2, 0)); + assert_eq!(sim.reaction_index("decay").unwrap(), (1, 0)); + assert_eq!(sim.reaction_index("make").unwrap(), (0, 0)); + assert!(sim.reaction_index("nope").is_err()); + assert_eq!(sim.reaction_name(1, 0).unwrap(), "decay"); + assert_eq!(sim.reaction_rate("bind").unwrap(), 10.0); + + sim.set_reaction_rate("decay", 2.5, RateKind::Rate).unwrap(); + assert_eq!(sim.reaction_rate("decay").unwrap(), 2.5); + sim.set_reaction_products("bind", RevParam::Irrev, 0.0, None, None) + .unwrap(); + sim.set_reaction_intersurface("bind", &[1]).unwrap(); + sim.set_reaction_intersurface("bind", &[]).unwrap(); + + sim.add_compartment("everywhere").unwrap(); + sim.add_compartment_point("everywhere", &[0.0; 3]).unwrap(); + sim.set_reaction_region("decay", Some("everywhere"), None) + .unwrap(); + + sim.add_solution_molecules("A", 50, None, None).unwrap(); + sim.add_solution_molecules("B", 50, None, None).unwrap(); + sim.update().unwrap(); + sim.run().unwrap(); +} + +#[test] +fn output_data_and_commands() { + let _guard = lock(); + let mut sim = box_sim(); + sim.add_species("A", None).unwrap(); + sim.add_solution_molecules("A", 7, None, None).unwrap(); + sim.add_output_data("counts").unwrap(); + sim.add_command("E molcount counts").unwrap(); + sim.add_timed_command('@', 0.004, 0.0, 0.0, 0.0, "molcount counts") + .unwrap(); + assert!( + sim.add_timed_command('é', 0.0, 0.0, 0.0, 0.0, "stop") + .is_err() + ); + sim.update().unwrap(); + sim.run_command("molcount counts").unwrap(); + sim.run().unwrap(); + + let data = sim.output_data("counts", false).unwrap(); + assert_eq!(data.ncol, 2, "time + one species"); + assert!(data.nrow >= 7, "{data:?}"); + assert!(data.rows().all(|row| row[1] == 7.0), "{data:?}"); + assert_eq!(data.rows().count(), data.nrow); + assert!(data.row(data.nrow).is_none()); + + let erased = sim.output_data("counts", true).unwrap(); + assert_eq!(erased, data); + assert_eq!(sim.output_data("counts", false).unwrap().nrow, 0); + assert!(sim.output_data("nope", false).is_err()); +} + +#[test] +fn output_files() { + let _guard = lock(); + let dir = PathBuf::from(env!("CARGO_TARGET_TMPDIR")); + let _ = fs::remove_file(dir.join("api_out.txt")); + let mut sim = box_sim(); + sim.add_species("A", None).unwrap(); + sim.set_output_path(&format!("{}/", dir.display())).unwrap(); + sim.add_output_file("api_out.txt", None, false).unwrap(); + sim.add_command("E molcount api_out.txt").unwrap(); + sim.update().unwrap(); + sim.run().unwrap(); + drop(sim); + let lines = fs::read_to_string(dir.join("api_out.txt")) + .unwrap() + .lines() + .count(); + assert!(lines >= 5, "{lines} lines"); +} + +#[test] +fn ports() { + let _guard = lock(); + let mut sim = box_sim(); + sim.add_species("A", None).unwrap(); + sim.add_surface("edge").unwrap(); + sim.add_panel( + "edge", + PanelShape::Rect, + Some("r1"), + Some("+0"), + &[0.9, -1.0, -1.0, 2.0, 2.0], + ) + .unwrap(); + sim.set_surface_action( + "edge", + PanelFace::Front, + "A", + MolecState::Soln, + SrfAction::Port, + None, + ) + .unwrap(); + sim.add_port("out", "edge", PanelFace::Front).unwrap(); + assert_eq!(sim.port_index("out").unwrap(), 0); + assert_eq!(sim.port_name(0).unwrap(), "out"); + sim.update().unwrap(); + + sim.add_port_molecules("out", "A", 3, None).unwrap(); + let at = [0.95, 0.0, 0.0]; + sim.add_port_molecules("out", "A", 1, Some(&[&at])).unwrap(); + assert!(sim.add_port_molecules("out", "A", 2, Some(&[&at])).is_err()); + assert_eq!(sim.molecule_count("A", MolecState::All).unwrap(), 4); + assert_eq!( + sim.port_molecules("out", "A", MolecState::All, false) + .unwrap(), + 0 + ); +} + +#[test] +fn lattices() { + let _guard = lock(); + let mut sim = box_sim(); + sim.add_species("A", None).unwrap(); + // this used to crash: smolAddLattice read sim->latticess before it existed. + sim.add_lattice("lat", &[-1.0; 3], &[1.0; 3], &[0.5; 3], "rrr") + .unwrap(); + assert!( + sim.add_lattice("lat2", &[-1.0; 2], &[1.0; 3], &[0.5; 3], "rrr") + .is_err() + ); + assert_eq!(sim.lattice_index("lat").unwrap(), 0); + assert_eq!(sim.lattice_name(0).unwrap(), "lat"); + sim.add_lattice_species("lat", "A").unwrap(); + sim.add_lattice_molecules("lat", "A", 5, None, None) + .unwrap(); + sim.add_reaction("decay", &[("A", MolecState::Soln)], &[], Some(1.0)) + .unwrap(); + sim.add_lattice_reaction("lat", "decay", false).unwrap(); +} + +#[test] +fn graphics_settings() { + let _guard = lock(); + let mut sim = box_sim(); + sim.set_graphics_params("none", Some(1), Some(0)).unwrap(); + assert!(sim.set_graphics_params("nonsense", None, None).is_err()); + sim.set_background_style([0.0, 0.0, 0.0, 1.0]).unwrap(); + sim.set_frame_style(Some(2.0), Some([1.0, 1.0, 1.0, 1.0])) + .unwrap(); + sim.set_grid_style(Some(1.0), None).unwrap(); + sim.set_text_style([1.0, 1.0, 1.0, 1.0]).unwrap(); + sim.set_light_params(0, None, Some([1.0, 1.0, 1.0, 1.0]), None, None) + .unwrap(); + assert!( + sim.set_light_params(-1, None, Some([1.0; 4]), None, None) + .is_err() + ); + sim.add_text_display("time").unwrap(); + sim.set_tiff_params(Some(10), Some("snap"), Some(1), None) + .unwrap(); +} + +#[test] +fn load_from_file_then_modify() { + let _guard = lock(); + let model = PathBuf::from(env!("CARGO_TARGET_TMPDIR")).join("api_load.txt"); + fs::write( + &model, + "dim 2\nboundaries 0 0 1\nboundaries 1 0 1\ntime_start 0\ntime_stop 0.01\ntime_step 0.002\nspecies A\n", + ) + .unwrap(); + let mut sim = Sim::load_from_file(&model, "q").unwrap(); + assert_eq!(sim.dim(), 2); + sim.add_species("B", None).unwrap(); + sim.add_solution_molecules("B", 3, None, None).unwrap(); + sim.update().unwrap(); + assert_eq!(sim.species(), ["A", "B"]); + assert_eq!(sim.molecule_count("B", MolecState::All).unwrap(), 3); + sim.run().unwrap(); + + assert!(Sim::load_from_file("/nonexistent/model.txt", "q").is_err()); +} diff --git a/rust/libsmoldyn/tests/sim.rs b/rust/libsmoldyn/tests/sim.rs new file mode 100644 index 00000000..b791949e --- /dev/null +++ b/rust/libsmoldyn/tests/sim.rs @@ -0,0 +1,237 @@ +use std::fs; +use std::path::PathBuf; +use std::sync::Mutex; +use std::sync::atomic::{AtomicBool, Ordering}; + +use libsmoldyn::{MolecState, Progress, Sim, SmolError}; + +// libsmoldyn keeps global error state, so tests must not run concurrently. +static LOCK: Mutex<()> = Mutex::new(()); + +const TINY_MODEL: &str = r#" +dim 3 +boundaries 0 -1 1 +boundaries 1 -1 1 +boundaries 2 -1 1 +time_start 0 +time_stop 0.01 +time_step 0.002 + +species A +difc A 1 +"#; + +fn tiny_model(name: &str) -> PathBuf { + let dir = PathBuf::from(env!("CARGO_TARGET_TMPDIR")); + fs::create_dir_all(&dir).unwrap(); + let model = dir.join(name); + fs::write(&model, TINY_MODEL).unwrap(); + model +} + +fn lock() -> std::sync::MutexGuard<'static, ()> { + LOCK.lock().unwrap_or_else(|e| e.into_inner()) +} + +#[test] +fn from_file_and_run() { + let _guard = lock(); + let mut sim = Sim::from_file(tiny_model("sim_run.txt"), "q").unwrap(); + sim.run().unwrap(); +} + +#[test] +fn missing_file_is_smoldyn_error() { + let _guard = lock(); + let err = Sim::from_file("/nonexistent/model.txt", "q").unwrap_err(); + assert!(matches!(err, SmolError::Smoldyn { .. }), "{err}"); +} + +#[test] +fn nul_byte_is_invalid_argument() { + let _guard = lock(); + let err = Sim::from_file("bad\0name.txt", "").unwrap_err(); + assert!(matches!(err, SmolError::InvalidArgument(_)), "{err}"); +} + +#[test] +fn new_rejects_bad_dims() { + let _guard = lock(); + assert!(matches!( + Sim::new(&[0.0; 4], &[1.0; 4]), + Err(SmolError::InvalidArgument(_)) + )); + assert!(matches!( + Sim::new(&[0.0], &[1.0, 1.0]), + Err(SmolError::InvalidArgument(_)) + )); + assert!(matches!( + Sim::new(&[1.0], &[0.0]), + Err(SmolError::Smoldyn { .. }) + )); +} + +#[test] +fn build_and_step_to_completion() { + let _guard = lock(); + let mut sim = Sim::new(&[-1.0; 3], &[1.0; 3]).unwrap(); + sim.set_times(0.0, 0.01, 0.002).unwrap(); + sim.add_species("A", None).unwrap(); + sim.read_config("difc", "A 1").unwrap(); + sim.update().unwrap(); + + let mut steps = 0; + while sim.step().unwrap() == Progress::Running { + steps += 1; + assert!(steps < 100, "simulation never finished"); + } + assert!(steps >= 4, "only {steps} steps"); +} + +#[test] +fn bad_config_statement_is_error() { + let _guard = lock(); + let mut sim = Sim::new(&[-1.0; 3], &[1.0; 3]).unwrap(); + assert!(sim.read_config("no_such_statement", "1 2 3").is_err()); +} + +#[test] +fn getters_reflect_model() { + let _guard = lock(); + let mut sim = Sim::new(&[-1.0, -2.0], &[1.0, 3.0]).unwrap(); + sim.set_times(0.0, 0.01, 0.002).unwrap(); + sim.add_species("A", None).unwrap(); + sim.add_species("B", None).unwrap(); + sim.read_config("difc", "A 1").unwrap(); + sim.read_config("mol", "10 A u u").unwrap(); + sim.read_config("mol", "3 B u u").unwrap(); + sim.update().unwrap(); + + assert_eq!(sim.dim(), 2); + assert_eq!(sim.bounds(), (vec![-1.0, -2.0], vec![1.0, 3.0])); + assert_eq!(sim.time_start(), 0.0); + assert_eq!(sim.time_stop(), 0.01); + assert_eq!(sim.time_step(), 0.002); + assert_eq!(sim.species(), ["A", "B"]); + + assert_eq!(sim.molecule_count("A", MolecState::Soln).unwrap(), 10); + assert_eq!(sim.molecule_count("all", MolecState::All).unwrap(), 13); + assert!(sim.molecule_count("C", MolecState::All).is_err()); + + sim.run_until(0.006).unwrap(); + assert!(sim.time() >= 0.006 - 1e-12, "time is {}", sim.time()); + assert!(sim.time() < sim.time_stop()); +} + +#[test] +fn run_after_step_reaches_stop_time() { + let _guard = lock(); + let mut sim = Sim::new(&[-1.0; 3], &[1.0; 3]).unwrap(); + sim.set_times(0.0, 0.01, 0.002).unwrap(); + sim.update().unwrap(); + sim.step().unwrap(); + sim.step().unwrap(); + sim.run().unwrap(); + assert!( + (sim.time() - sim.time_stop()).abs() < 1e-9, + "stopped at {} instead of {}", + sim.time(), + sim.time_stop() + ); +} + +fn tiny_sim() -> Sim { + let mut sim = Sim::new(&[-1.0; 3], &[1.0; 3]).unwrap(); + sim.set_times(0.0, 0.01, 0.002).unwrap(); + sim.add_species("A", None).unwrap(); + sim.read_config("mol", "5 A u u u").unwrap(); + sim.update().unwrap(); + sim +} + +#[test] +fn run_with_calls_back_every_step() { + let _guard = lock(); + let mut sim = tiny_sim(); + let mut times = Vec::new(); + let progress = sim + .run_with(&AtomicBool::new(false), |s| { + times.push(s.time()); + assert_eq!(s.molecule_count("A", MolecState::All).unwrap(), 5); + }) + .unwrap(); + + assert_eq!(progress, Progress::Finished); + assert!(times.len() >= 5, "{times:?}"); + assert!(times.windows(2).all(|w| w[0] < w[1]), "{times:?}"); + assert!((sim.time() - sim.time_stop()).abs() < 1e-9); +} + +#[test] +fn run_with_stops_and_resumes() { + let _guard = lock(); + let mut sim = tiny_sim(); + let stop = AtomicBool::new(false); + + let mut steps = 0; + let progress = sim + .run_with(&stop, |_| { + steps += 1; + if steps == 2 { + stop.store(true, Ordering::Relaxed); + } + }) + .unwrap(); + assert_eq!(progress, Progress::Running); + assert_eq!(steps, 2); + assert!(sim.time() < sim.time_stop()); + + stop.store(false, Ordering::Relaxed); + assert_eq!(sim.run_with(&stop, |_| {}).unwrap(), Progress::Finished); + assert!((sim.time() - sim.time_stop()).abs() < 1e-9); +} + +#[test] +fn run_with_on_other_thread() { + let _guard = lock(); + let mut sim = tiny_sim(); + let (tx, rx) = std::sync::mpsc::channel(); + let worker = std::thread::spawn(move || { + sim.run_with(&AtomicBool::new(false), |s| tx.send(s.time()).unwrap()) + .unwrap() + }); + let times: Vec = rx.iter().collect(); + assert_eq!(worker.join().unwrap(), Progress::Finished); + assert!(!times.is_empty()); +} + +#[test] +fn run_with_resume_keeps_output() { + let _guard = lock(); + let out_dir = PathBuf::from(env!("CARGO_TARGET_TMPDIR")); + let out = out_dir.join("run_with_resume.txt"); + let _ = fs::remove_file(&out); + + let mut sim = tiny_sim(); + sim.read_config("output_root", &format!("{}/", out_dir.display())) + .unwrap(); + sim.read_config("output_files", "run_with_resume.txt") + .unwrap(); + sim.add_command("E molcount run_with_resume.txt").unwrap(); + sim.update().unwrap(); + + let stop = AtomicBool::new(false); + let mut steps = 0; + sim.run_with(&stop, |_| { + steps += 1; + if steps == 2 { + stop.store(true, Ordering::Relaxed); + } + }) + .unwrap(); + stop.store(false, Ordering::Relaxed); + sim.run_with(&stop, |_| {}).unwrap(); + + let lines = fs::read_to_string(&out).unwrap().lines().count(); + assert!(lines >= 5, "only {lines} lines in {out:?}"); +} diff --git a/rust/smoldyn-rs/Cargo.toml b/rust/smoldyn-rs/Cargo.toml new file mode 100644 index 00000000..4a0f4887 --- /dev/null +++ b/rust/smoldyn-rs/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "smoldyn" +version = "2.27.0" +edition = "2024" + +[dependencies] +anyhow = "1" +clap = { version = "4.6.6", features = ["derive"] } +strum = { version = "0.28.0", features = ["derive"] } +tracing = "0.1.44" +tracing-subscriber = { version = "0.3.23", features = ["env-filter"] } +libsmoldyn = { path = "../libsmoldyn/" } +ctrlc = "3.5.2" diff --git a/rust/smoldyn-rs/src/core/molecule.rs b/rust/smoldyn-rs/src/core/molecule.rs new file mode 100644 index 00000000..9f6eb1a2 --- /dev/null +++ b/rust/smoldyn-rs/src/core/molecule.rs @@ -0,0 +1,14 @@ + +#[derive(Debug, Eq, PartialEq, strum::Display, strum::EnumString)] +#[strum(serialize_all = "lowercase")] +pub enum MoleculeState { + Soln, + Front, + Back, + Down, + BSoln, + FSoln, // really Soln + All, + None, + Some, +} diff --git a/rust/smoldyn-rs/src/kinetics/mod.rs b/rust/smoldyn-rs/src/kinetics/mod.rs new file mode 100644 index 00000000..0e2bc051 --- /dev/null +++ b/rust/smoldyn-rs/src/kinetics/mod.rs @@ -0,0 +1,2 @@ +pub mod molecule; +pub use molecule::*; diff --git a/rust/smoldyn-rs/src/kinetics/molecule.rs b/rust/smoldyn-rs/src/kinetics/molecule.rs new file mode 100644 index 00000000..6b70c4c4 --- /dev/null +++ b/rust/smoldyn-rs/src/kinetics/molecule.rs @@ -0,0 +1,13 @@ +#[derive(Debug, Eq, PartialEq, strum::Display, strum::EnumString)] +#[strum(serialize_all = "lowercase")] +pub enum MoleculeState { + Soln, + Front, + Back, + Down, + BSoln, + FSoln, // really Soln + All, + None, + Some, +} diff --git a/rust/smoldyn-rs/src/lib.rs b/rust/smoldyn-rs/src/lib.rs new file mode 100644 index 00000000..30c3d791 --- /dev/null +++ b/rust/smoldyn-rs/src/lib.rs @@ -0,0 +1,7 @@ +pub mod kinetics; +pub use kinetics::*; + +pub mod simulation; +pub use simulation::*; + +pub use libsmoldyn::version; diff --git a/rust/smoldyn-rs/src/main.rs b/rust/smoldyn-rs/src/main.rs new file mode 100644 index 00000000..2f25240b --- /dev/null +++ b/rust/smoldyn-rs/src/main.rs @@ -0,0 +1,79 @@ +use smoldyn::{Progress, Simulation}; +use std::path::PathBuf; +use std::process::ExitCode; +use std::sync::Arc; +use std::sync::atomic::{AtomicBool, Ordering}; +use tracing_subscriber::{EnvFilter, fmt, prelude::*}; + +use clap::{Parser, Subcommand}; + +#[derive(Parser)] +#[command(about, long_about = None)] +struct Cli { + model: Option, + + #[arg(short, long, action = clap::ArgAction::Count)] + debug: u8, + + #[arg(long)] + version: bool, + + #[command(subcommand)] + command: Option, +} + +#[derive(Subcommand)] +enum Commands { + Simulate, +} + +fn main() -> anyhow::Result { + let cli = Cli::parse(); + + let log_level = match cli.debug { + 0 => "warn", + 1 => "info", + 2 => "debug", + _ => "trace", + }; + + unsafe { + std::env::set_var("SMOLDYN_LOG", log_level); + } + + tracing_subscriber::registry() + .with(fmt::layer()) + .with(EnvFilter::from_env("SMOLDYN_LOG")) + .init(); + + if cli.version { + show_version(); + return Ok(ExitCode::SUCCESS); + } + + if let Some(model) = cli.model { + let stop = stop_on_ctrlc()?; + let mut sim = Simulation::new().with_model_path(model)?; + if sim.run(&stop)? == Progress::Running { + return Ok(ExitCode::from(130)); + } + } + + Ok(ExitCode::SUCCESS) +} + +fn stop_on_ctrlc() -> anyhow::Result> { + let stop = Arc::new(AtomicBool::new(false)); + let handler_stop = stop.clone(); + ctrlc::set_handler(move || { + if handler_stop.swap(true, Ordering::Relaxed) { + std::process::exit(130); + } + eprintln!("stopping after the current time step (Ctrl-C again to quit now)"); + })?; + Ok(stop) +} + +fn show_version() { + println!("{}", libsmoldyn::version()); +} diff --git a/rust/smoldyn-rs/src/simulation/mod.rs b/rust/smoldyn-rs/src/simulation/mod.rs new file mode 100644 index 00000000..ffa81036 --- /dev/null +++ b/rust/smoldyn-rs/src/simulation/mod.rs @@ -0,0 +1,24 @@ +pub mod simulation; +pub use simulation::*; + +use std::path::Path; +use std::sync::atomic::AtomicBool; + +use anyhow::Context; +pub use libsmoldyn::Progress; +use libsmoldyn::Sim; + +pub fn run(model: &Path, flags: &str, stop: &AtomicBool) -> anyhow::Result { + tracing::info!("Running model {:?}", model); + + let mut sim = + Sim::from_file(model, flags).with_context(|| format!("failed to load model {model:?}"))?; + let progress = sim + .run_with(stop, |sim| tracing::trace!("t = {}", sim.time())) + .with_context(|| format!("failed to run model {model:?}"))?; + if progress == Progress::Running { + tracing::warn!("stopped at t = {} of {}", sim.time(), sim.time_stop()); + } + + Ok(progress) +} diff --git a/rust/smoldyn-rs/src/simulation/simulation.rs b/rust/smoldyn-rs/src/simulation/simulation.rs new file mode 100644 index 00000000..5f8bbe66 --- /dev/null +++ b/rust/smoldyn-rs/src/simulation/simulation.rs @@ -0,0 +1,78 @@ +use std::collections::HashMap; +use std::path::PathBuf; +use std::sync::atomic::AtomicBool; + +#[derive(Default, Debug)] +pub struct Simulation { + pub model_path: Option, + pub boundary: Boundary, + pub count: usize, + pub current_time: f64, +} + +impl Simulation { + pub fn new() -> Self { + Self::default() + } + + pub fn with_model_path(mut self, model_path: PathBuf) -> anyhow::Result { + self.model_path = Some(model_path); + Ok(self) + } + + pub fn initialize(&mut self) -> anyhow::Result<()> { + todo!() + } + + pub fn run(&mut self, stop: &AtomicBool) -> anyhow::Result { + let model_path = self + .model_path + .as_ref() + .ok_or_else(|| anyhow::anyhow!("no model to run"))?; + tracing::info!("running model {model_path:?}"); + super::run(model_path, "", stop) + } + + pub fn run_until( + &mut self, + _breaktime: f64, + _dt: f64, + _display: bool, + _overwrite: bool, + ) -> anyhow::Result<()> { + todo!() + } + + pub fn set_current_time(&mut self, current_time: f64) { + self.current_time = current_time + } + + pub fn get_current_time(&self) -> f64 { + self.current_time + } + + pub fn add_command( + &mut self, + _command: &str, + _cmd_type: char, + _kwargs: HashMap, + ) -> anyhow::Result<()> { + todo!() + } +} + +#[derive(Default, Debug, PartialEq)] +pub struct Boundary { + low: [f64; 3], + high: [f64; 3], +} + +#[cfg(test)] +mod test { + + #[test] + fn test_version() { + let version = crate::version(); + println!("{version}"); + } +} diff --git a/rust/smoldyn-rs/tests/e2e.rs b/rust/smoldyn-rs/tests/e2e.rs new file mode 100644 index 00000000..9e031b4e --- /dev/null +++ b/rust/smoldyn-rs/tests/e2e.rs @@ -0,0 +1,59 @@ +use std::fs; +use std::path::PathBuf; +use std::sync::Mutex; +use std::sync::atomic::AtomicBool; + +// The Smoldyn C library keeps global error/warning state, so tests that touch it +// must not run concurrently within the same process. +static LOCK: Mutex<()> = Mutex::new(()); + +const TINY_MODEL: &str = r#" +dim 3 +boundaries 0 -1 1 +boundaries 1 -1 1 +boundaries 2 -1 1 +time_start 0 +time_stop 0.01 +time_step 0.002 + +species A +difc A 1 +"#; + +fn tmp_dir() -> PathBuf { + let dir = PathBuf::from(env!("CARGO_TARGET_TMPDIR")); + fs::create_dir_all(&dir).unwrap(); + dir +} + +#[test] +fn runs_tiny_model() { + let _guard = LOCK.lock().unwrap(); + + let model = tmp_dir().join("tiny_e2e.txt"); + fs::write(&model, TINY_MODEL).unwrap(); + + println!("{}", smoldyn::version()); + let progress = + smoldyn::run(&model, "", &AtomicBool::new(false)).expect("tiny model should run"); + assert_eq!(progress, smoldyn::Progress::Finished); +} + +#[test] +fn missing_model_fails() { + let _guard = LOCK.lock().unwrap(); + + let model = tmp_dir().join("does_not_exist.txt"); + assert!(smoldyn::run(&model, "", &AtomicBool::new(false)).is_err()); +} + +#[test] +fn stop_flag_interrupts_run() { + let _guard = LOCK.lock().unwrap(); + + let model = tmp_dir().join("tiny_e2e_stop.txt"); + fs::write(&model, TINY_MODEL).unwrap(); + + let progress = smoldyn::run(&model, "", &AtomicBool::new(true)).unwrap(); + assert_eq!(progress, smoldyn::Progress::Running); +} diff --git a/source/Smoldyn/libsmoldyn.cpp b/source/Smoldyn/libsmoldyn.cpp index 9187b16a..e9bc4566 100644 --- a/source/Smoldyn/libsmoldyn.cpp +++ b/source/Smoldyn/libsmoldyn.cpp @@ -797,6 +797,7 @@ extern CSTRING enum ErrorCode smolRunCommand(simptr sim,const char *commandstrin strncpy(stringcopy,commandstring,STRCHARLONG-1); cmd=scmdalloc(); LCHECK(cmd,funcname,ECmemory,"failed to create a new command structure"); + cmd->cmds=sim->cmds; // commands that write output (e.g. molcount) need this strcpy(cmd->str,stringcopy); cmdcode=docommand((void*)sim,cmd,stringcopy); LCHECK(cmdcode==CMDok,funcname,ECwarning,cmd->erstr); @@ -2178,7 +2179,7 @@ extern CSTRING enum ErrorCode smolAddLattice(simptr sim,const char *lattice,cons simlattice = NULL; LCHECK(sim,funcname,ECmissing,"missing sim"); LCHECK(lattice,funcname,ECmissing,"missing lattice"); - p=stringfind(sim->latticess->latticenames,sim->latticess->nlattice,lattice); + p=sim->latticess?stringfind(sim->latticess->latticenames,sim->latticess->nlattice,lattice):-1; LCHECK(p<0,funcname,ECsame,"lattice name already exists"); er=latticeaddlattice(sim,&simlattice,lattice,min,max,dx,btype,LATTICEnsv); LCHECK(er==0,funcname,ECerror,"error adding lattice"); diff --git a/steve/src/main.rs b/steve/src/main.rs new file mode 100644 index 00000000..ae18da65 --- /dev/null +++ b/steve/src/main.rs @@ -0,0 +1,10 @@ +use autocxx::prelude::*; + +include_cpp! { + safety!(unsafe) + generate!("DeepThought") +} + +fn main() { + println!("main.."); +}