Crate multicalc

Expand description

§multicalc

Downloads

Scientific computing that fits on a microcontroller, built and tested from scratch in one integrated package. Estimation, control, kinematics, Lie groups, calculus, autodiff and linear algebra in stable no_std Rust with no heap, no panics and no unsafe.

§Why use it

1 kHz loop rates: No heap, fixed-size types, bounded work per call. Results in a full robotics control loop at 1 kHz.
Tested on six embedded targets: Every commit is built and tested on six targets: the x86_64 and aarch64 Linux hosts and on four bare-metal ABIs (thumbv7em soft-float, thumbv7em hardware-FPU, thumbv6m, and riscv32imc), running the real math under QEMU. no_std, no-alloc, and no-panic rules hold on each target.
Measured against external references: Each module’s results are verified against established libraries like numpy, scipy, and filterpy fixtures within ~1 ulp, thus validating the rust implementation. See the benchmarks.
Pure safe and panic-free. #![forbid(unsafe_code)], no C dependencies, and unwrap/ panic denied on library paths; every fallible call returns a typed error. Types are fixed-size and stack-allocated, and iteration counts are bounded.

§What it does

§Robotics and control

Estimation: linear, extended, and unscented KalmanFilters (autodiff Jacobians, no hand-derived ones; the unscented one needs no derivatives at all), an ErrorStateKalmanFilter that fuses an IMU with position and heading fixes, MahonyFilter and MadgwickFilter for attitude estimation, and a ParticleFilter for nonlinear, non-Gaussian problems (alloc only), with a Monte Carlo Localization built on top of it.
Control: Pid control, infinite horizon Lqr,; GeometricAttitudeController for drones, the pure pursuit path-following law; and FollowTheGap reactive obstacle avoidance.
Spatial math: Quaternion, the SO2/SE2/SO3/SE3 Lie groups for 2D and 3D rotations and rigid-body transforms with left and right Jacobians and their inverses on all four, and Twist/Wrench screw-theory types.
Rigid-body dynamics: RigidBody computes the motion of a single rigid body, from a SpatialInertia saying how its mass is spread out and a FreeJointState for a body free to move in all six directions — loadable straight from MuJoCo model files with multicalc-mjcf.
Plant: What sits between a command and the force a body actually feels — MultirotorMixer shares a wanted lift and turn out across the rotors, and RotorLag models the moment a rotor takes to catch up to what it was asked for.
Kinematics: differential-drive and unicycle maps between wheel and body motion, with exact SE(2) odometry.
Motion: PolylinePath for waypoint paths with arc-length, closest-point, and lookahead queries, and MinimumSnapPlanner for the smoothest trajectory through them.
Mapping: 2D OccupancyGrid and ScanGeometry

§Core math

Automatic differentiation: Exact autodiff of any order (total and partial), plus Jacobian and Hessian matrices.
Linear algebra: fixed-size, stack-allocated Matrix and Vector with LU, Cholesky, column-pivoted QR, SVD, symmetric eigendecomposition, and the matrix exponential expm. General N×N determinant and inverse, pseudo-inverse, eigenvalue clamping, solve_discrete_riccati and solve_discrete_lyapunov.
Least-squares optimization: LevenbergMarquardt and GaussNewton solvers for nonlinear curve fitting.
Root finding: bracketed bisection and Newton solvers for scalar equations and square systems, with an optional damped line search.
Polynomials: Polynomial for evaluation with any number of derivatives in one pass, arithmetic, calculus, fitting and real roots; PiecewisePolynomial for curves made of pieces; and MultivariatePolynomial for several variables with symbolic partial derivatives.
Integration: iterative Newton-Cotes rules (Boole, Simpson, Trapezoidal) and Gaussian quadrature (Legendre, Hermite, Laguerre) over finite, semi-infinite, and infinite limits.
ODE integrators: fixed-step Rk4 and adaptive Rk45 (Dormand-Prince 5(4)) with PI step control and dense output, plus ExponentialMap, which is a purely orientation integrator.
Discretization: zero-order hold, Van Loan, and discrete white-noise models for continuous-time linear systems.
Signal processing: Biquad low-pass, high-pass, band-pass, and notch filters; with cascades, motor-harmonic notches, and per-channel filtering. Plus MovingAverage, RunningMedian, SavitzkyGolay smoothing, Deadband, Hysteresis and SlewRateLimiter conditioning.
Vector calculus: curl, divergence, and line and flux integrals.
Approximation: linear and quadratic Taylor models with goodness-of-fit metrics.
Random: Pcg32 and the RandomSource trait, a seedable no_std generator for the particle filter and for stochastic models.

§Quick start

Two formulas, written once, carried through six modules, each step feeding the next:

use multicalc::prelude::*;
use multicalc::{Hessian, Jacobian, KalmanFilter, KalmanModel, Matrix, Newton, SE3, SO3, Vector, c};
use multicalc::{scalar_fn, scalar_fn_vec};

fn main() -> Result<(), CalcError> {
    // Written once, evaluated at f64 here and at an autodiff number wherever a derivative is asked
    // for — the formula text never changes.
    let f = scalar_fn!(|x| x * x * x - c(2.0) * x);                     // f(x)    = x³ - 2x
    let g = scalar_fn!(|v: &[f64; 2]| v[0] * v[0] * v[1] + v[0].sin()); // g(x, y) = x²y + sin x

    // Derivatives — exact, by forward-mode autodiff. No step size, no truncation error.
    let single_point = 2.0_f64;
    let slope = derivative(&f, single_point);                // f'(2)  = 10
    let bend = second_derivative(&f, single_point);          // f''(2) = 12

    let point = [1.0_f64, 2.0];
    let x_index = 0;
    let dg_dx = partial(&g, x_index, &point)?;

    // The derivative matrices of those same two formulas.
    let hessian = Hessian::new().evaluate(&g, &point)?;      // 2x2 second derivatives
    let both = scalar_fn_vec!(|v: &[f64; 2]| [
        v[0] * v[0] * v[1] + v[0].sin(),
        v[0] * v[0] * v[0] - c(2.0) * v[0],
    ]);
    let jacobian = Jacobian::new().evaluate(&both, &point)?; // 2x2 first derivatives

    // Integration — f again, this time over an interval.
    let limits = [0.0, 2.0];
    let area = integral(&|x: f64| f.eval(x), limits)?;       // ∫₀² f = 0

    // Linear algebra — solve H·x = b with the Hessian computed three lines up.
    let b = Vector::new([1.0, 2.0]);
    let x = hessian.solve(b)?;

    // Root finding — Newton on the same f, its derivative supplied by autodiff.
    let initial_guess = 2.0;
    let root = Newton::new().solve(&f, initial_guess)?.root; // √2 ≈ 1.41421356

    // Rigid-body motion — SO(3)/SE(3), generic over the scalar like everything above.
    let quarter_turn_about_z = Vector::new([0.0, 0.0, core::f64::consts::FRAC_PI_2]);
    let translation = Vector::new([1.0, 2.0, 3.0]);
    let start = Vector::new([1.0, 0.0, 0.0]);

    let pose = SE3::from_parts(SO3::exp(quarter_turn_about_z), translation);
    let moved = pose.act(start);                  // rotate, then translate → (1, 3, 3)

    // Estimation — a Kalman filter recovering the velocity it never measures.
    let initial_state = Vector::new([0.0, 0.0]);  // [position, velocity]
    let initial_covariance = Matrix::new([[1.0, 0.0], [0.0, 1.0]]);
    let model = KalmanModel {
        state_transition: Matrix::new([[1.0, 1.0], [0.0, 1.0]]),
        measurement_model: Matrix::new([[1.0, 0.0]]),        // position only
        process_noise: Matrix::new([[0.01, 0.0], [0.0, 0.01]]),
        measurement_noise: Matrix::new([[0.1]]),
    };

    let mut filter = KalmanFilter::new(initial_state, initial_covariance, model);
    filter.predict();

    let measurement = Vector::new([1.0]);         // the target moved about 1 m
    filter.update(measurement)?;
    let velocity = filter.state()[1];             // recovered, though never measured

    Ok(())
}

Every fallible call propagates with ?: each module has its own error enum, and all of them convert into the CalcError umbrella, so one return type covers a program that mixes modules.

§Full tutorial

Refer to the tutorials for a comprehensive tutorial for each module. They show the full imports, expected outputs in comments, error-path notes, and pointers to runnable demos. Start there when you need the complete picture of a feature.

§Accuracy

Verified against external-library fixtures (mpmath, numpy, scipy, filterpy) in the multicalc-qa crate, with per-module tables generated from those fixtures. See benchmarks/README.md for the index, or go straight to calculus, linear_algebra, optimization, ode, estimation, or root_finding.

§Runnable demos

Runnable, self-contained programs for each module live in the demos/ crate. See demos/README.md. Run one with:

cargo run -p multicalc-demos --example <name>

§Feature flags

alloc (off by default): enables the heap-based methods for inputs too large for the stack. See Heap allocation.

§Heap allocation

The library allocates nothing by default: every type is fixed-size and lives on the stack. Turning on alloc pulls in extern crate alloc and unlocks exactly two things:

estimation::ParticleFilter, whose cloud of samples is sized at runtime and so cannot be a fixed-size stack type.
numerical_derivative::jacobian::Jacobian::get_on_heap, which returns a Vec<Vec<_>> for Jacobians too large to sit on the stack. The stack-allocated get is always available.

Nothing else changes: no_std, forbid(unsafe_code), and the no-panic rules hold either way, and the feature never pulls in std.

§MSRV and edition

Edition 2024, minimum supported Rust version 1.85.

§Contributing

See CONTRIBUTING.md.

§Acknowledgements

The least-squares solvers and QR factorization port the public-domain MINPACK routines lmder, lmpar, qrfac, and qrsolv (Moré, Garbow, Hillstrom; netlib), following Moré (1978), “The Levenberg-Marquardt algorithm: Implementation and theory”, and Nocedal & Wright, Numerical Optimization (chapters 4 and 10).

§License

multicalc is licensed under the MIT license.

§Contact

anmolkathail@gmail.com

Re-exports§

pub use scalar::Numeric;
pub use scalar::Dual;
pub use scalar::HyperDual;
pub use scalar::Jet;
pub use scalar::ScalarFn;
pub use scalar::ScalarFnN;
pub use scalar::VectorFn;
pub use scalar::Const;
pub use scalar::Primal;
pub use scalar::c;
pub use numerical_derivative::AutoDiffMulti;
pub use numerical_derivative::AutoDiffSingle;
pub use numerical_derivative::DerivatorMultiVariable;
pub use numerical_derivative::DerivatorSingleVariable;
pub use numerical_derivative::FiniteDifferenceConfig;
pub use numerical_derivative::FiniteDifferenceMode;
pub use numerical_derivative::FiniteDifferenceMulti;
pub use numerical_derivative::FiniteDifferenceSingle;
pub use numerical_derivative::Hessian;
pub use numerical_derivative::Jacobian;
pub use numerical_derivative::derivative;
pub use numerical_derivative::partial;
pub use numerical_derivative::second_derivative;
pub use numerical_integration::GaussianConfig;
pub use numerical_integration::GaussianMulti;
pub use numerical_integration::GaussianQuadratureMethod;
pub use numerical_integration::GaussianSingle;
pub use numerical_integration::IntegratorMultiVariable;
pub use numerical_integration::IntegratorSingleVariable;
pub use numerical_integration::IterativeConfig;
pub use numerical_integration::IterativeMethod;
pub use numerical_integration::IterativeMulti;
pub use numerical_integration::IterativeSingle;
pub use numerical_integration::SummationMethod;
pub use numerical_integration::integral;
pub use approximation::LinearApproximation;
pub use approximation::LinearApproximationPredictionMetrics;
pub use approximation::LinearApproximator;
pub use approximation::QuadraticApproximation;
pub use approximation::QuadraticApproximationPredictionMetrics;
pub use approximation::QuadraticApproximator;
pub use linear_algebra::Matrix;
pub use linear_algebra::Vector;
pub use linear_algebra::Matrix2D;
pub use linear_algebra::Matrix3D;
pub use linear_algebra::Matrix4D;
pub use linear_algebra::Matrix6D;
pub use linear_algebra::Vector2D;
pub use linear_algebra::Vector3D;
pub use linear_algebra::Vector6D;
pub use linear_algebra::solve_discrete_lyapunov;
pub use linear_algebra::solve_discrete_riccati;
pub use discretization::q_discrete_white_noise;
pub use discretization::van_loan;
pub use discretization::zoh;
pub use ode::ExponentialMap;
pub use ode::Rk4;
pub use ode::Rk45;
pub use signal_processing::Biquad;
pub use signal_processing::BiquadCascade;
pub use signal_processing::BiquadCoefficients;
pub use signal_processing::Deadband;
pub use signal_processing::Hysteresis;
pub use signal_processing::MovingAverage;
pub use signal_processing::MultiChannelBiquad;
pub use signal_processing::OnePoleLowPass;
pub use signal_processing::RunningMedian;
pub use signal_processing::SavitzkyGolay;
pub use signal_processing::SlewRateLimiter;
pub use signal_processing::harmonic_notch_coefficients;
pub use spatial::Quaternion;
pub use spatial::SE2;
pub use spatial::SE3;
pub use spatial::SO2;
pub use spatial::SO3;
pub use spatial::Twist;
pub use spatial::Wrench;
pub use spatial::SpatialInertia;
pub use spatial::FreeJointState;
pub use kinematics::BodyArc;
pub use kinematics::BodyTwist;
pub use kinematics::DifferentialDrive;
pub use kinematics::WheelRotations;
pub use kinematics::WheelVelocities;
pub use mapping::MutableOccupancyMap;
pub use mapping::OccupancyMap;
pub use mapping::ScanGeometry;
pub use mapping::DynamicOccupancyGrid;alloc
pub use estimation::CovarianceUpdate;
pub use estimation::ExtendedKalmanFilter;
pub use estimation::KalmanFilter;
pub use estimation::KalmanModel;
pub use estimation::UnscentedKalmanFilter;
pub use estimation::ErrorStateKalmanFilter;
pub use estimation::ImuNoise;
pub use estimation::NominalState;
pub use estimation::NominalStateFn;
pub use estimation::MadgwickFilter;
pub use estimation::MahonyFilter;
pub use estimation::ConstantTurnAndSpeed;
pub use estimation::DirectMeasurement;
pub use estimation::residual_with_wrapped_angles;
pub use estimation::GaussianLikelihood;alloc
pub use estimation::Likelihood;alloc
pub use estimation::ParticleFilter;alloc
pub use estimation::ResamplingScheme;alloc
pub use estimation::BeamModel;alloc
pub use estimation::InitialParticleCloud;alloc
pub use estimation::MonteCarloLocalizer;alloc
pub use random::Pcg32;
pub use random::RandomSource;
pub use optimization::GaussNewton;
pub use optimization::LevenbergMarquardt;
pub use optimization::MinimizationReport;
pub use optimization::TerminationReason;
pub use root_finding::Bisection;
pub use root_finding::Newton;
pub use root_finding::NewtonSystem;
pub use root_finding::RootReport;
pub use root_finding::RootReportN;
pub use root_finding::RootTermination;
pub use polynomial::MultivariatePolynomial;
pub use polynomial::MultivariateTerm;
pub use polynomial::PiecewisePolynomial;
pub use polynomial::Polynomial;
pub use polynomial::RealRoots;
pub use control::Curvature;
pub use control::FollowTheGap;
pub use control::FollowTheGapOutput;
pub use control::GeometricAttitudeController;
pub use control::Lqr;
pub use control::Pid;
pub use control::ThrustCommand;
pub use control::pure_pursuit_curvature;
pub use control::thrust_command_from_acceleration;
pub use motion::BoundaryDerivatives;
pub use motion::EndOfPath;
pub use motion::MinimumSnapPlanner;
pub use motion::PathProjection;
pub use motion::PolylinePath;
pub use motion::durations_from_average_speed;
pub use dynamics::RigidBody;
pub use dynamics::RigidBodyAcceleration;
pub use plant::MultirotorMixer;
pub use plant::RotorCommands;
pub use plant::RotorLag;
pub use plant::RotorSpin;
pub use error::CalcError;
pub use error::ControlError;
pub use error::DiffError;
pub use error::DynamicsError;
pub use error::EstimationError;
pub use error::IntegrateError;
pub use error::KinematicsError;
pub use error::LinalgError;
pub use error::MappingError;
pub use error::MotionError;
pub use error::PlantError;
pub use error::PolynomialError;
pub use error::SignalError;
pub use error::SolveError;
pub use error::SpatialError;
pub use libm;

Modules§

approximation: Least-squares function approximation with goodness-of-fit metrics.
control: Control: feedback controllers, signal filters, and path-following laws.
discretization: Discretization of continuous-time linear systems. All build on Matrix::expm.
dynamics: Dynamics: how a rigid body moves under the forces put on it.
error: Error types for the crate. Each module family has its own enum; CalcError is the umbrella they all convert into.
estimation: State estimation from noisy measurements.
gaussian_tables: Gaussian quadrature node and weight tables.
kinematics: Kinematics: maps between actuator motion and body motion, and pose integration.
linear_algebra: Fixed-size, stack-allocated linear algebra.
mapping: Maps a robot can measure against: a grid of free and blocked cells, and the range a sensor reads across it.
motion: Motion: waypoint paths, planned trajectories, and the geometric queries a path-following controller consumes.
numerical_derivative: Differentiation: exact automatic differentiation and finite differences.
numerical_integration: Numerical integration.
ode: Ordinary differential equation integrators.
optimization: Nonlinear least-squares optimization.
plant: Plant: the machinery between a command and the force a body actually feels.
polynomial: Polynomials: as coefficients, as pieces, and in several variables.
prelude: The traits and one-call functions worth importing together.
random: Credits: Pcg32 follows the permuted congruential generator designed and published by Melissa O’Neill, with the constants from her reference implementation (see https://www.pcg-random.org/). The normal draw uses the polar form of the Box–Muller transform. Thanks to both for putting the method and code in the open.
root_finding: Root finding for scalar equations and square systems.
scalar: The scalar number system the calculus modules are generic over.
signal_processing: Signal processing: filters, smoothers, and signal conditioning.
spatial: Spatial math: rotations, Lie groups, and spatial-algebra types.
vector_field: Vector calculus on 2D/3D fields.

Macros§

matrix: Builds a Matrix from bracketed row literals.
multivariate_polynomial: Builds a MultivariatePolynomial from pairs of a number and the power each variable is raised to.
polynomial: Builds a Polynomial from its coefficients, lowest power first.
scalar_fn: Builds a ScalarFn (one variable) or ScalarFnN (N variables) from closure-style syntax.
scalar_fn_vec: Builds a VectorFn from a closure returning a fixed-size array, for Jacobians.
vector: Builds a Vector from a comma-separated list of components.

Crate multicalc

Crate multicalc Copy item path

§multicalc

§Why use it

§What it does

§Robotics and control

§Core math

§Quick start

§Full tutorial

§Accuracy

§Runnable demos

§Feature flags

§Heap allocation

§MSRV and edition

§Contributing

§Acknowledgements

§License

§Contact

Re-exports§

Modules§

Macros§

Crate multicalc