Docs 2.4.4 · 7 min read

Science, engineering and medicine

Apsi ships the numerical tools of a research lab as standard modules. They need nothing beyond NumPy, accept plain numbers, tensors and physical quantities, and are checked in the test suite against SciPy, statsmodels, lifelines, scikit-learn, pydicom and exact analytic solutions (tests/test_science.py).

The existing use sci (FFT, integration, roots, optimisation, curve fitting with uncertainties, stiff ODEs), use linalg, units, measurement uncertainty and automatic differentiation work together with all of them.

Special functions

use special
print(gamma(4.5), lgamma(200.0), beta(2.0, 3.0), factorial(25), comb(60, 30))
print(erf(1.0), erfinv(0.5), erfcinv(1e-300))
print(gammainc(2.5, 1.0), betainc(2.0, 3.0, 0.4), gammaincinv(2.5, 0.9))
print(besselj(0, 2.404825557695773), besseli(1, 3.0), ellipk(0.5), ellipe(0.5))
print(legendre(5, 0.3), hermite(4, 1.2), laguerre(3, 0.5, 1.0), sph_harm(2, 1, 0.7, 1.1))

Every function works element-wise on tensors. Accuracy is close to double precision (relative errors of 1e-16 to 1e-13 against SciPy). gamma, lgamma, erf, erfc, erfinv, factorial, comb, beta, besselj and legendre are also in use math; erf is a differentiable built-in.

Statistics

use stats
let a = randn(40) + 0.5
let b = randn(35)
let t = ttest_ind(a, b)                       # Welch by default
print(t.statistic, t.p, t.df, t.ci, t.cohen_d)

let d = student_t(12)                          # pdf, cdf, sf, ppf, sample, mean, var, interval
print(d.ppf(0.975), d.interval(0.95))

let fit = ols(X, y, names=["dose", "age"])     # coefficients, SE, t, p, CI, R², adjusted R², F, AIC
print(fit.coef.dose, fit.p, fit.r2)
let lr = logit(X, outcome)                     # odds ratios with CI
print(lr.odds_ratio, lr.odds_ratio_ci)

print(sample_size_ttest(0.5, power=0.8).n)     # 64 per group, exact noncentral t
print(adjust_pvalues([0.01, 0.04, 0.03, 0.2], "bh"))

Every test takes alternative="two-sided" | "less" | "greater" where it applies, and returns a record with the statistic, the p-value and what a paper reports (degrees of freedom, confidence intervals, effect sizes). The results agree with SciPy and statsmodels to 1e-10 or better.

Signal processing

use signal
let sos = butter(4, [0.5, 40.0], "bandpass", 500.0)     # second-order sections
let clean = filtfilt(sos, ecg)                          # zero phase, SciPy-compatible edges
let y = notch(clean, 50.0, 500.0)                        # mains hum
let pk = find_peaks(y, height=0.5, distance=100, prominence=0.3)
let psd = welch(y, 500.0, nperseg=1024)                 # record(f, psd)
let sg = spectrogram(y, 500.0, 256, 128)                 # record(f, t, power)

Also: sosfilt, freqz, lowpass, highpass, bandpass, bandstop, hann, hamming, blackman, kaiser, convolve, correlate, savgol, moving_average, detrend, resample, hilbert, envelope. Filter design, filtering, peaks, Welch spectra, Savitzky–Golay, resampling and the Hilbert transform agree with scipy.signal to 1e-8 or better. When SciPy is installed its compiled filter loop is used for speed; otherwise a built-in loop gives the same result.

Mechanics

use mechanics

# any force law: x'' = accel(x, v, t)
fn spring(x) = -x
let r = simulate(spring, tensor([1.0]), tensor([0.0]), 0.01, 10000, "yoshida4")

# Lagrangian mechanics: the equations of motion are derived by automatic differentiation
let g = 9.81
fn L(q, qd) = qd[0]^2 + 0.5*qd[1]^2 + qd[0]*qd[1]*cos(q[0]-q[1]) + 2*g*cos(q[0]) + g*cos(q[1])
let pend = lagrangian_flow(L, tensor([1.5, -0.5]), tensor([0.0, 0.0]), 0.001, 5000)
print(max(pend.energy) - min(pend.energy))            # energy is conserved to ~1e-8

# gravity
let solar = nbody(masses, positions, velocities, 3600.0, 24 * 365, softening=1e3)
print(solar.energy_drift)

Partial differential equations

use pde
let r = heat1d(u0, dx, dt, 1000, alpha, bc="neumann", every=100)   # record(t, x, u[frame, i])
let w = wave1d(u0, 0.0, c, dx, dt, 2000, bc="fixed")
let phi = poisson2d(rho / epsilon, dx)                             # exact solution of the 5-point scheme
let flow = navier_stokes2d(omega0, nu=1e-3, dt=1e-3, steps=5000, every=500)
let s = cg(A, b)                                                   # A: matrix or function x -> A·x

heat2d (ADI), advection1d (upwind, Lax–Wendroff), bicgstab and tridiag_solve (batched Thomas) are also available. The tests compare every solver with exact solutions: decaying sine modes, standing waves, the discrete Laplacian and the Taylor–Green vortex (to 1e-15).

Quantum mechanics

use quantum
let x = linspace(-10.0, 10.0, 2001)
let s = schrodinger1d(fn(x) => 0.5 * x * x, x, mass=1.0, n=5, hbar=1.0)   # 0.5, 1.5, 2.5, … (1e-8)
let wave = evolve1d(gaussian_packet(x, -3.0, 0.5, 4.0), fn(x) => 0.0, x, 0.001, 2000, mass=1.0, hbar=1.0)

var psi = ket("00")
psi = apply(psi, gates.H, 0)
psi = apply(psi, gates.CNOT, [0, 1])
print(measure(psi, shots=1000, seed=1), concurrence(psi), entropy(partial_trace(density(psi), [0])))

SI units are the default (mass = m_e, hbar = ħ); pass hbar=1.0, mass=1.0 for natural units. Bound states use fourth-order finite differences with ψ = 0 beyond the ends of the grid. Also: expect_x, expect_p, uncertainty, probability, dagger, commutator, anticommutator, kron, expm, expm_hermitian, time_evolution, pauli, spin_ops(s), rx, ry, rz, phase, controlled, probabilities, purity, fidelity, bloch, ghz, bell, hydrogen_energy, hydrogen_radial.

Physics formulas

use physics now also contains formulas that take and return quantities with units:

use physics
print(wien_peak(5778K))                      # 5.0e-07 m
print(planck(500e-9[m], 5778K))               # spectral radiance
print(lorentz_factor(0.9 * c), time_dilation(1[s], 0.9 * c))
print(coulomb_field([1e-9], [[0.0, 0.0, 0.0]], [[1.0, 0.0, 0.0]]))
print(biot_savart(loop_points, 2.0, [[0.0, 0.0, 0.0]]))

Thermal radiation (planck, wien_peak, blackbody_power), relativity (lorentz_factor, relativistic_energy, relativistic_momentum, lorentz_boost, time_dilation, length_contraction, velocity_add), waves and optics (doppler, relativistic_doppler, snell, critical_angle, brewster_angle, fresnel, thin_lens, single_slit, double_slit, photon_energy, de_broglie), electromagnetism (coulomb_field, coulomb_potential, biot_savart (exact finite segments), lorentz_force, cyclotron_frequency, skin_depth), thermodynamics (ideal_gas_pressure, carnot_efficiency, maxwell_speed_pdf, mean_speed, rms_speed, heat_conduction) and decay (decay, activity, half_life_from_rate).

Medical and health research

These functions are for research and education. They are not a medical device and must not be used on their own to make decisions about the care of a patient. The text is also available in programs as DISCLAIMER.

use medical

# pharmacokinetics
let c = pk_one_compartment(500, 40, times, CL=5.0, route="oral", ka=1.1, tau=12, doses=6)
let r = nca(times, measured, dose=500)        # Cmax, Tmax, AUC (linear-up/log-down), λz, t½, CL/F, Vz/F, MRT

# survival analysis
let km = kaplan_meier(time, event)            # survival, Greenwood SE, log-log CI, median
let lr = logrank(time, event, group)
let cox = cox_ph(time, event, X, names=["age", "treatment"])
print(cox.hazard_ratio, cox.hr_ci, cox.p, cox.concordance)

# diagnostic accuracy
print(diagnostic_metrics(90, 20, 10, 180))    # sensitivity, specificity, PPV, NPV, LR±, DOR with CIs
let roc_ = roc(scores, labels)                 # AUC with DeLong CI, Youden threshold
print(roc_compare(model_a, model_b, labels).p)

# signals and images
let beats = ecg_rpeaks(ecg, 360.0)
print(hrv(beats.rr))                           # SDNN, RMSSD, pNN50, LF, HF, LF/HF
let ct = read_dicom("slice.dcm")               # pixels in Hounsfield units, spacing in mm
let view = window_level(ct.pixels, 40, 400)    # soft-tissue window

Kaplan–Meier, the log-rank test and Cox regression agree with lifelines to 1e-9 or better. ROC AUC and kappa agree with scikit-learn, and the 2×2 and Mantel–Haenszel measures with statsmodels. DICOM pixels are bit-identical to pydicom for explicit/implicit little-endian and big-endian files.

read_dicom reads uncompressed DICOM files and leaves patient-identifying fields out unless you pass include_identifiers=true. Compressed transfer syntaxes (JPEG, JPEG 2000, RLE) give a clear message asking you to decompress the file first. Reading DICOM needs --allow-read, like any other file access.

Reproducible research

use lab
set_seed(42)                                   # Apsi, numpy and Python's random module
let run = experiment("dropout-sweep", record(dropout=0.1, lr=3e-4))   # needs --allow-write
for epoch in 0..20:
    run.log(epoch, record(loss=loss, val_auc=auc))
run.result(record(test_auc=0.912))
run.finish()                                   # experiments/dropout-sweep.jsonl

print(format_measurement(9.8123, 0.0312, "m/s^2"))         # 9.81(3) m/s^2
print(table_latex([record(model="A", auc=0.91), record(model="B", auc=0.87)], caption="Results"))
let runs = load_runs("experiments")                          # every run with params, logs and provenance

Every run records its provenance: the Apsi, Python and NumPy versions, OS, CPU, GPUs, seed and time (UTC). User names, host names and paths are not recorded. Tables come as Markdown (table_markdown), LaTeX with booktabs and escaping (table_latex) or CSV (table_csv). sigfig rounds to significant figures.

Found a mistake? Report it · Markdown source