Gallery
Worked examples, each a Python script that calls vipyr. Every page shows the code, what it printed and what it drew, from a run made when these docs were built.
Data
Check feature code for future leakage, get the same numbers in training and serving, and compile custom per-group logic to C.
Dataframes and features
Statistics
Probability distributions, summary statistics and regression, with the names and results of SciPy, statsmodels and R.
Distributions
Regression and GLMs
NIST's certified linear regressions, digit by digit
NIST's Statistical Reference Datasets include eleven linear least-squares problems with coefficients, standard errors, residual standard deviation and R^2 certified to 15 digits, from the trivial (Norris, a line) to the notorious (Filip, a degree-10 polynomial; Wampler4-5, degree 5 with large residuals). vp.stats.ols fits each with every solver of its slot, and the table gives the certified digits of the worst coefficient (NIST's log relative error, LRE) beside statsmodels 0.14 and R 4.3.3's lm.
A probit GLM, R's way
vp.stats.glm is R's glm.fit step for step. Here a probit model of mtcars' engine shape (vs: straight or V) against fuel economy, the data as R ships it. The coefficient table is R's summary.glm, and R 4.3.3 gives the same numbers within 2 × 10−15 relative. lib/stats/glm_check.json records the comparison for this fit and six more of R's datasets: at R's epsilon 10−14 the coefficients, standard errors and deviance agree within 1e-10. The fitting method and the information behind the standard errors are named slots: IRLS with the expected information is R's choice; Newton with the observed information (its second derivatives taken by Viv's deriv) is the other.
Machine learning
A small GPT trained in compiled Viv and returned to Python as named arrays, matching a NumPy program and PyTorch step for step; scikit-learn pipelines compiled to one C file.
Neural networks and transformers
A quantized model file, read in Viv
vp.gguf reads GGUF, the file format of ggml and llama.cpp, with lib/gguf: the header and metadata (header.viv), the tensor table (tensors.viv) and ggml's number formats (quants.viv), each a rewrite of ggml's code built one step at a time. The models are Karpathy's llama2.c models trained on TinyStories, in GGUF from huggingface.co/ggml-org/models (MIT). The smallest, stories260K, is in the repository (lib/gguf/tests/stories260K.gguf) and is read below as the page runs. The 15-million-parameter model comes in two files, one with 8-bit weights (Q8_0) and one with 4-bit weights (Q4_0); they are not in the repository, so gallery/ml/gguf_stories15m.py reads them with the same calls and writes the numbers this page prints to gallery/ml/gguf_stories15m.json. Every number equals what gguf-py's dequantize gives for the same bytes (tools/gguf/check.sh).
Models to C
Finance
Option pricing by closed form, binomial trees, Fourier inversion and finite differences, each checked against published values or SciPy.
Options and volatility
Implied volatility, binomial trees and Greeks of European options
Recover the implied volatility of 3,000 option quotes and plot the smile, watch a Cox-Ross-Rubinstein tree converge to the Black-Scholes price, and compare Greeks from Viv's reverse-mode grad with the closed forms. The pricing code is examples/options/*.viv, loaded as Python modules with vp.load.
The Heston model's volatility smile
Price European options under the Heston stochastic-volatility model with vp.finance.heston_cos, Fang and Oosterlee's COS method (a Fourier-cosine expansion of the density), convert the prices to Black-Scholes implied volatilities, and plot the smile at three maturities and for three correlations between the stock and its variance. The COS price converges to Fang and Oosterlee's published value as the number of terms grows; the convergence plot shows how fast.
Early exercise of an American put, priced two ways
Price an American put, which may be exercised at any time before expiry, with two independent methods: Crank-Nicolson finite differences on the Black-Scholes PDE with the early-exercise constraint (vp.finance.cn_american_put), and a Cox-Ross-Rubinstein binomial tree (vp.finance.crr_price). Plot the price against the stock price, the value of the right to exercise early, and the tree's convergence to the PDE.
Science
Classic research codes, callable from Python and checked against their originals: PySCF Hartree-Fock, BIGSTICK shell model, plane-wave DFT, Vlasov plasma, Monte Carlo.
Quantum chemistry
The Boys function, the integral behind Gaussian basis sets
Evaluate the Boys functions with vp.special.boys, which returns orders 0 through n at once, and measure their relative accuracy against mpmath at 50 digits over from 0 to 200. Every electron-repulsion and nuclear-attraction integral over Gaussian orbitals reduces to these functions, so a Hartree-Fock code calls them millions of times. This one is lib/hf/boys.viv, the routine Viv's Hartree-Fock library computes its integrals with.
Atoms and solids
Atoms from hydrogen to uranium, against NIST
Solve the Kohn-Sham equations of 23 atoms in the local density approximation with vp.dftatom.atom, the Viv port of Certik, Pask and Vackar's dftatom, and compare every orbital energy and total energy with NIST's atomic reference data (Standard Reference Database 141). Then plot radial densities of the noble gases, and redo gold and uranium with the Dirac equation to see what relativity does to their orbitals.
Nuclear
Level schemes of neon-20 and magnesium-24 in the sd shell
Diagonalize the USDB interaction for the valence nucleons above oxygen-16 with vp.nuclear.shell_model, a Lanczos solver in the M scheme (lib/nuclear/shell/shellmodel.viv), and draw the low-lying spectra with their angular momenta. Neon-20 has 4 valence nucleons and a basis of 640 states; magnesium-24 has 8 and a basis of 28,503. Every energy is compared with BIGSTICK 8.0.0, an independent shell-model code, run on the same Hamiltonian.
Lanczos and dense diagonalization in the sd shell
vp.nuclear.shell_model finds the lowest states with Lanczos. The program solves that step with the contract lowest_eigenpairs (lib/eigen/lowest), whose other methods form the Hamiltonian as a matrix, one column per basis state, and diagonalize all of it: householder_ql by Householder reduction and QL, jacobi by cyclic Jacobi. Lanczos and Householder + QL are run here on the same nuclei; their energies agree to within 10−12 MeV, and the cost of each is read from the contract's recorded sweep.
Passes to self-consistency in calcium-40 by density mixer
Solve calcium-40 with the Skyrme force SLy4 using vp.nuclear.skyrme_hf three times, once per density mixer of lib/scf/mix.viv: linear mixing, Pulay's mixing and Johnson's modified Broyden mixing (mixing=). Each pass of the self-consistent field (lib/scf/loop.viv) builds new densities from the last ones; the mixer decides what the next pass starts from. Count the passes each needs and check that all three reach the same nucleus.
Plasma
Landau damping at four wavelengths
Perturb a Maxwellian plasma with a small density wave and follow the electric field as it decays without any collisions. vp.plasma.landau runs a Vlasov-Poisson solver (plasma/landau.viv: Cheng and Knorr's split step, on a phase-space grid) and fits the decay rate and frequency to the peaks of the field energy. The exact values are the root of the Langmuir dispersion relation for a Maxwellian, which the program also computes, from the plasma dispersion function; the fitted rate matches it to better than one part in at every .
The Faddeeva and plasma dispersion functions on the complex plane
Map the Faddeeva function and the plasma dispersion function over a rectangle of the complex plane, and measure their accuracy against SciPy and against mpmath at 30 digits. vp.special.faddeeva_w and vp.special.plasma_z evaluate one point per call; for a whole grid this page loads a ten-line Viv file, gallery/science/zgrid.viv, that loops the same library function (lib/zfunc.viv, Weideman's rational expansion) over the grid in one call.
Statistical mechanics
The Ising model through its phase transition
Run Monte Carlo on the two-dimensional Ising model at 13 temperatures and three lattice sizes with vp.statmech.ising, and watch the magnetization vanish, the susceptibility peak and the Binder cumulants cross near Onsager's critical temperature T_c = 2/ln(1 + sqrt 2) = 2.269. The program is stat/ising/ising_mc.viv; at every run it also reports the exact answers it can compute for that lattice, which is how the numbers here are checked.
The Lennard-Jones fluid against its equation of state
Simulate 256 Lennard-Jones atoms by molecular dynamics at constant temperature with vp.liquids.md_nvt, the Viv port of Allen and Tildesley's md_nvt_lj, at several state points, and compare the energy and pressure with the equation of state Thol et al. fitted to simulations of the same model (vp.liquids.eos). Along the way, g(r) shows the fluid's structure change from a gas to a dense liquid. Every program here gives the bits of upstream's Fortran from the same start (ports/liquids/COMPARISON.md).
Energy conservation and diffusion in a Lennard-Jones liquid
Run constant-energy molecular dynamics of 256 Lennard-Jones atoms with vp.liquids.md_nve, the Viv port of Allen and Tildesley's md_nve_lj, and measure two things: how the fluctuation of the total energy shrinks with the time step, which is velocity Verlet's error made visible, and how the atoms diffuse, from the velocity autocorrelation function and from the mean-square displacement.
Math
Differential equations, linear algebra, special functions and randomized sketches written as readable algorithms, checked against SciPy and mpmath.
Differential equations
Robertson's stiff kinetics over sixteen decades
Three species, rate constants 0.04, 104 and 3 × 107: the fast reactions settle in microseconds and the slow one keeps going until t = 1e11. An explicit method would need steps of about 10−4 the whole way. vp.ode.solve with Radau IIA and with BDF crosses the range in under two thousand steps, lands on the published reference solution, and takes the steps SciPy's solve_ivp takes.
Ten thousand epidemics in one call
The SIR model run for 10,000 values of the basic reproduction number R0, each an adaptive Dormand-Prince solve of its own, in one vp.ode.batch call. The final size of each epidemic is checked against the closed form it must satisfy, and one row is checked against a single vp.ode.solve bit for bit.
Linear algebra
Fitting and optimization
Floating-point accuracy
Dynamics and control
Riccati and Lyapunov solvers on the Schur form, LQR design and closed-loop simulation; every solution is checked by the residual of its equation.
Linear control and stability
Balancing a cart-pole with LQR
Linearize a cart-pole about the upright, compute the optimal state feedback from the continuous algebraic Riccati equation, and run the nonlinear system under it from 25 degrees off vertical. The Riccati solver is Viv's Schur method, an algorithm after SLICOT (lib/are.viv, lib/schur.viv); its residual is checked here and its solution compared with SciPy's.
A stability map for the cart-pole's feedback gains
Scale the two halves of the cart-pole's LQR gain, the part that acts on the cart and the part that acts on the pole, and map which combinations keep the upright balanced. For every stable gain, solve the Lyapunov equation A X + X A^T = -I with vp.control.lyapunov (Bartels-Stewart, on Viv's Schur method, an algorithm after SLICOT). X is the steady-state covariance of the state when unit white noise drives every component, so its trace measures how far disturbances push each stable loop.
Hardware
Logic built from NAND gates up, simulated gate by gate; a C compiler, an assembler and an RV32IM simulator, each written in Viv.
Digital logic
Processors
Power systems
AC and DC power flow on MATPOWER cases, Newton-Raphson and fast-decoupled, with a sparse LU written in Viv, checked against MATPOWER's own solution of case9 and against an independent SciPy power flow on every case.
Power flow
Voltage profile of the IEEE 118-bus system
Solve the power flow of MATPOWER's case118 (the IEEE 118-bus test system, a model of part of the American Electric Power system in the Midwestern United States) with vp.power.runpf, and plot every bus's voltage magnitude and angle. The Newton-Raphson solution is compared with the linear DC model, which assumes flat 1.0 per-unit magnitudes and lossless branches, and with the solution when generators are held to their reactive limits.
Power flow convergence and timing from 9 to 9241 buses
Newton-Raphson against the fast-decoupled method on MATPOWER cases from 9 to 2383 buses: the mismatch left after each iteration, the iterations each needs, and the time per solve measured against an independent SciPy power flow, up to the 9241-bus PEGASE model of the European grid.
An interconnection study on the IEEE 118-bus system
A 200 MW generator asks to connect at bus 38 of MATPOWER's case118, a 345 kV bus. vp.power.interconnect solves the base case and all 240 single outages (every branch and generator) without the project and with it, lists the limits the project breaks, gives its distribution factors on every branch, and finds its transfer limit: the largest output with no new violation, by bisection on the full AC N-1.