# Four tools for seismic data > Source: > Published: 2026-10-22 07:36:20+00:00 *Originally published at [amansinghnp.com](https://amansinghnp.com/writing/seismic-pytorch).* Over 2025 and 2026 I built four small open-source tools that sit at different points of a seismic data pipeline. Each solves a different problem, so start with the one you need: | Tool | What it does | Reach for it if you | Code | |---|---|---|---| | `torch-friction` | Rate-and-state friction as a differentiable PyTorch layer, so friction parameters can be fitted by gradient descent | want to recover *a* ,*b* and*Dc* from friction data, or put exact friction physics inside a larger network | [GitHub](https://github.com/AmanSinghNp/torch-friction) | | `Basalt` | C++ frequency-wavenumber (f-k) filtering, built around memory-mapped SEG-Y, aligned memory, tiled FFTs and AVX2/AVX-512 kernels | need to remove ground roll from active-source shot gathers | [GitHub](https://github.com/AmanSinghNp/Basalt) | | `fast-trigger` | Recursive STA/LTA event detection in C++23 with a NumPy-friendly Python API | need to flag events in single traces or batches of traces from Python | [GitHub](https://github.com/AmanSinghNp/fast-trigger) | | `seismic-linter` | A static analyser (on PyPI) that flags temporal data leakage in ML code before training runs | train models on earthquake catalogues or other time series and want leaky code blocked at commit time | [GitHub](https://github.com/AmanSinghNp/seismic-linter) | Each tool's section opens with a short entry: the clearest result, the main limitation and how it was tested. The physics, numerics and full test results follow. Every chart below was produced by running the published code (or, where noted, a line-for-line NumPy port of it) on the data listed in section 6. Some of the results are flattering and some are not. I've kept both. The earthquake-severity model that led to `seismic-linter` is on the Projects page as [South Asia Earthquake Severity](https://amansinghnp.com/projects#earthquake). **< 0.05%** friction parameter error after inversion **−20.7 dB** ground roll removed by the f-k mute **≈ 2 dB** less signal needed by x² STA/LTA at a 1% false-alarm rate Machine learning in seismology has moved fast. Deep networks now pick phases and detect events at scale (Mousavi & Beroza, 2022), and physics-informed neural networks (PINNs) put governing equations into the loss function (Raissi et al., 2019; Karniadakis et al., 2021). Two problems keep coming up anyway. The first is that the physics is often bolted on after the fact, instead of living inside the computational graph where gradients can reach it. The second is that evaluation is fragile: time-series leakage has inflated reported results across many fields (Kapoor & Narayanan, 2023), and earthquake forecasting has its own history of complex models that turn out no better than simple baselines (Mignan & Broccardo, 2019, 2020). The four tools line up against those two problems. `torch-friction` and `Basalt` handle the physics and signal-processing side. `fast-trigger` handles detection. `seismic-linter` guards the evaluation. **torch-friction in brief** `torch-friction` from PyPI in float64 on CPU. Laboratory rock friction isn't a constant. It depends on how fast the surfaces slide and on how long the contacts have had to age. Dieterich (1979) and Ruina (1983) captured this with a friction coefficient that depends on slip velocity *V* and a state variable *θ*, which you can read as the average lifetime of the asperity contacts: μ(V, θ) = μ0 + *a* ln(V / V0) + *b* ln(V0θ / Dc) The *a* term is the **direct effect**: friction jumps immediately when velocity changes. The *b* term is the **evolution effect**: as the contact population turns over across a characteristic slip distance *Dc*, friction relaxes to a new steady state. At steady state, *θ* = *Dc*/*V* and μss = μ0 + (*a* − *b*) ln(V/V0). When *a* − *b* < 0 the fault is *velocity-weakening*: sliding faster makes it weaker, which is the precondition for earthquake nucleation (Marone, 1998; Scholz, 1998). Whether that weakening actually produces an earthquake depends on the surroundings. In the classic spring-slider model, a block dragged through a spring of stiffness *k* slides stably if *k* is above a critical stiffness *kc* = σ(*b* − *a*)/*Dc*, and stick-slips below it (Ruina, 1983; Rice & Ruina, 1983). *Figure 1. The spring-slider model. The load point moves at a constant rate, stress builds in the spring while the block sticks, then drops when the block slips. With velocity-weakening friction and *k* < *kc*, this cycle repeats indefinitely: the laboratory analogue of the seismic cycle.* The friction equation says nothing about how *θ* changes over time. There are two standard choices, and `torch-friction` implements both as `law='aging'` and `law='slip'`: | Law | dθ/dt | Behaviour | Source | |---|---|---|---| | Aging (Dieterich) | 1 − Vθ/Dc | State grows with time even at rest, so contacts heal while stationary | Dieterich (1979) | | Slip (Ruina) | −(Vθ/Dc) ln(Vθ/Dc) | State only evolves with slip; symmetric response to velocity steps up and down | Ruina (1983) | The standard lab test is a velocity step. Figure 2 runs one through `FrictionLayer.forward_trajectory` with *a* = 0.015, *b* = 0.020 and *Dc* = 10 µm, stepping from 1 to 10 µm/s and back. *Figure 2. Velocity step up at t = 2 s and back down at t = 8 s, computed with torch-friction (Δt = 2 ms, RK4). The jump of *a*·ln10 ≈ 0.035 is the direct effect; the settling to (*a* − *b*)·ln10 ≈ −0.012 below the starting level is the evolution effect. The two laws agree on the peak and the steady state but differ on the path between them.* The two laws only differ in the transient, but the difference isn't small: up to Δμ = 0.012 about one second after the step up, which is a third of the whole direct effect. The asymmetry is easy to see from the equations. After a tenfold step up, Vθ/Dc = 10, so the aging law evolves at a rate of 1 − 10 = −9 while the slip law evolves at −10 ln 10 ≈ −23. The slip law weakens faster. After the step down (Vθ/Dc = 0.1) it's reversed: 0.9 for aging against 0.23 for slip, so the aging law heals faster. Six seconds after the step down, the aging curve has recovered to μ = 0.586 while the slip curve is still at 0.575. Lab experiments with large velocity steps tend to favour the slip law (Bhattacharya et al., 2015), which is why the library keeps both. Each call to the layer advances *θ* by one classical fourth-order Runge-Kutta (RK4) step, then evaluates the friction equation at the new state. The core is short: `torch_friction/layers.py` ``` python def rk4_step(self, t, theta, v, dt, Dc): k1 = self.physics(t, theta, v, Dc) k2 = self.physics(t + dt/2, theta + k1*dt/2, v, Dc) k3 = self.physics(t + dt/2, theta + k2*dt/2, v, Dc) k4 = self.physics(t + dt, theta + k3*dt, v, Dc) return theta + (dt/6.0)*(k1 + 2*k2 + 2*k3 + k4) ``` The state equation is stiff when *Dc* is small relative to *V*·Δt. Linearising the aging law about steady state gives an eigenvalue of −V/Dc, and RK4's stability interval on the negative real axis ends near −2.79 (Hairer & Wanner, 1996), so explicit integration needs Δt < 2.79 Dc/V. The layer warns once Δt goes past 0.5 Dc/Vmax, which leaves a margin of about 5.6× below the stability limit. That margin is sensible, because the accuracy limit arrives well before the stability limit does. To check the accuracy, Figure 3 relaxes the aging law from *θ*0 = 10 s at constant *V*, where the exact solution is *θ*(t) = Dc/V + (θ0 − Dc/V)e−Vt/Dc, and compares it with the library's RK4 at seven step sizes. *Figure 3. Error in *θ* after 2 s of relaxation against step size. Halving Δt cuts the error by 19.7× and then 16.8×, close to the 24 = 16 of a fourth-order method. Below Δt ≈ 0.06 s the error stops falling at about 3 × 10−7 s.* The floor isn't RK4's fault. The layer guards its logarithms and divisions with a 10−12 epsilon and a smooth `soft_clamp`, and with *Dc* = 10−5 m those guards shift the steady state by about one part in 107. That's a deliberate trade: gradients that never hit NaN in exchange for a relative bias far below lab measurement noise. If you're working at very small *Dc*, scale the epsilon with it. Because every operation in the RK4 step is a PyTorch tensor op, autograd records the whole integration. `forward_trajectory` unrolls the step T times into preallocated output tensors, and `loss.backward()` replays that graph in reverse to get ∂L/∂*a*, ∂L/∂*b*, ∂L/∂*Dc* and ∂L/∂μ0 (Paszke et al., 2019). This is "discretise-then-optimise": the gradients are exact for the discrete solver, at the cost of memory that grows linearly with T. The alternative is the adjoint method used by neural ODEs, which solves a second ODE backwards in time and uses constant memory, but its gradients are only as accurate as that backward solve (Chen et al., 2018). For sequences of a few thousand steps, exact gradients are the better deal. *Figure 4. One cell of the unrolled graph. Black arrows are the forward pass; green arrows are gradients flowing back through all four RK4 stages to the parameters. Parameters are stored in unconstrained form and passed through softplus, so *a*, *b* and *Dc* stay positive whatever the optimiser does.* Three other details matter for training. Parameters default to float64, because the evolution term takes logs of quantities that span many orders of magnitude. Negative velocities raise an error instead of being silently clamped. And `freeze_parameters(['a', 'Dc'])` lets you fix quantities you've measured independently, which matters a great deal for the next section. The test is to generate noisy friction data from a "teacher" layer with hidden parameters, then train a "student" layer, started from wrong values, to recover them. I ran two versions. **Setup A** is the demo in the repository's `examples/` folder: one step from 1 to 10 µm/s, *Dc* = 10 mm, Δt = 0.1 s, 1,000 samples, Adam at a learning rate of 10−3 for 150 epochs. **Setup B** uses a lab-scale *Dc* = 10 µm, a sequence of three steps (up ×10, down ×10, up ×3) over 12 s at Δt = 2 ms, and Adam at 10−2 for 600 epochs. Both add Gaussian noise with σ = 10−4. *Figure 5. Setup B: relative error of each recovered parameter against epoch, starting from guesses that were 33–100% wrong. All three are within 1% by epoch 166 and within 0.05% at the end. The loss levels off at 9.9 × 10−9, which is the noise variance (10−8): the student is fitting everything except the noise.* | Parameter | True | Start | Setup A (150 epochs) | Setup B (600 epochs) | |---|---|---|---|---| | *a* | 0.015 | 0.010 | 0.0110 | 0.01500 | | *b* | 0.020 | 0.010 | 0.0092 | 0.02000 | | *Dc* | 10 mm / 10 µm | 2× true | 21.7 mm | 9.996 µm | | μ0 | 0.600 | 0.600 | 0.607 | 0.60000 | Setup A doesn't recover the parameters, and *b* actually moves further from the truth. Part of the reason is optimisation (the loss is still 1,000× the noise floor after 150 epochs), but the main reason is **identifiability**. After the step, the block slides at 10 µm/s for 80 s, which is 0.8 mm, only 8% of *Dc*. The evolution effect has barely started, so the data mostly constrain the direct effect plus a linearised evolution term in which *b* and *Dc* appear together as roughly *b*/*Dc*. Many (*b*, *Dc*) pairs fit equally well. Setup B slides many times *Dc* after each step and steps in both directions, so the two parameters leave separate fingerprints. **Design rule** For rate-and-state inversion, design the experiment before the optimiser. Record at least several *Dc* of slip after each step, include steps up and down, and freeze any parameter you've measured some other way. A differentiable solver finds the minimum; it can't create information the data don't contain. `torch-friction` is deliberately narrow. It solves the zero-dimensional friction ODE, not the elastodynamic problem around a fault. Rucker and Erickson (2024) go further and train PINNs for the wave equation with a rate-and-state fault boundary in 1D and 2D, recovering depth-varying friction parameters. Borate et al. (2024) show a PINN can recover rate-and-state parameters from acoustic monitoring of lab stick-slip experiments. A differentiable friction layer is a building block for that kind of model: a component whose physics is exact, placed inside a larger network that learns everything else. **Basalt in brief** `PERFORMANCE.md` is marked TBD, and the C++ AVX2 kernels weren't run for this post. The kernel also reads only `mute_vel_max`, so the reject cone is always symmetric in ± In an active-source shot gather, the reflections geologists care about arrive with high apparent velocities across the receiver spread. Ground roll, the surface waves travelling along the ground, is slow (typically a few hundred m/s), low-frequency and much stronger. A 2D Fourier transform maps each linear event with apparent velocity *v* to a line *f* = *v*·*k* in the frequency-wavenumber plane, so the two separate by slope. Muting the slow wedge and transforming back is the "pie-slice" filter (Embree et al., 1963; Yilmaz, 2001). Basalt's mute is a smooth sigmoid rather than a hard cut, because a hard edge in the f-k plane rings in the time domain: w(f, k) = sigmoid( 10 · ( |f / (k + ε)| − vcut ) / taper ) *Figure 6. (a) A synthetic gather: 128 traces at 5 m spacing, 4 ms sampling, four reflectors (2,000–3,600 m/s, 30 Hz Ricker) and two ground-roll modes (320 m/s at 8 Hz, 450 m/s at 12 Hz). (b) Its f-k spectrum; the shaded region, |f/k| < 1,500 m/s, is muted. (c) The result, using a NumPy port of Basalt's mask with vcut = 1,500 m/s and taper = 100 m/s.* Measured on each component separately (the filter is linear), the mask keeps **99.97%** of the reflection energy and **0.85%** of the ground-roll energy, a 20.7 dB reduction. The panel also shows the two classic f-k artefacts. The small residue at the top centre is the ground-roll apex at zero offset, where the event is locally flat, so its energy sits near *k* = 0 inside the pass region. The faint smears at the edges near 1.0 s are wrap-around: the FFT treats the gather as periodic in offset, so energy that gets cut at one edge leaks back in at the other. Production workflows handle these with spatial tapering and padding before the transform. *Figure 7. Basalt's path from disk to filtered data. The lower strip shows the mask kernel sweeping one frequency row eight wavenumbers at a time, one AVX2 register per iteration.* Most of the engineering is about not moving data. SEG-Y files are memory-mapped rather than read into buffers, and the parser decodes the EBCDIC textual header and legacy IBM-360 floating-point samples in place (SEG Technical Standards Committee, 2017). All working memory comes from a `MemoryArena`, a bump allocator with 64-byte alignment, so one gather's scratch space is a single allocation and `reset()` frees it in O(1). Complex spectra are stored as a structure of arrays (separate `real` and `imag` planes), so that eight consecutive real parts load into one 256-bit register without shuffling. A 2D FFT is a set of 1D FFTs along rows, a transpose, then 1D FFTs along what used to be columns (Cooley & Tukey, 1965). The transpose is where naive code falls apart: striding down columns of a large matrix misses the cache on nearly every access. Basalt tiles the transpose so each block fits in cache, the approach cache-oblivious algorithms formalise (Frigo et al., 1999), and picks the 2D FFT schedule automatically by problem size, in the same spirit as FFTW's planner (Frigo & Johnson, 2005). The filter itself is embarrassingly parallel: every f-k bin gets multiplied by its own weight. The AVX2 path tiles the spectrum (32×64 bins for small gathers, 64×128 by default, 32×96 above 2,048), prefetches 8 to 24 elements ahead, and handles eight bins per instruction: `src/filter/fk_filter.cpp (AVX2 path, abridged)` ``` js const __m256 vk = _mm256_load_ps(k_vals.data()); const __m256 vel = _mm256_div_ps(vf, _mm256_add_ps(vk, v_eps)); __m256 x = _mm256_sub_ps(absolute_ps(vel), v_cut); x = _mm256_mul_ps(x, v_steep); const __m256 weight = fast_sigmoid_avx(x); real = _mm256_mul_ps(real, weight); imag = _mm256_mul_ps(imag, weight); ``` There are no branches in the inner loop: the "mute or keep" decision is a multiplication by a weight between 0 and 1, so no lane waits on a mispredicted branch. At runtime, the library uses the AVX-512 path if the CPU reports `avx512f`, `vl`, `bw` and `dq`, and otherwise falls back to AVX2. One detail worth knowing: the kernel reads only `mute_vel_max`, so the reject cone is always symmetric in ±*k*. `mute_vel_min` is currently unused. I'm not quoting throughput numbers here, because the repository hasn't published any yet. `PERFORMANCE.md` lays out the method: pinned CPU governors, the `basalt_bench` sweep across size, threads, SIMD mode, NUMA mode and memory policy, and PMU counters on self-hosted Linux. Every results cell is still marked TBD, because shared CI runners are fine for correctness testing but no good for performance claims. What *is* in CI is correctness: a golden checksum for single-thread AVX2 output, and a check that multi-threaded batch output matches single-threaded output bit for bit. **fast-trigger in brief** `fast-trigger` uses, roughly 2 dB less signal, and picks earlier (a median 0.26 s after onset at 3σ, against 0.39 s). See section 4.3.`recursive_sta_lta` on synthetic 100 Hz seismograms with known P and S onsets. Thresholds were calibrated on 300 noise-only traces, with 150 trials per P-wave amplitude and fixed seeds. The repository adds Hypothesis property-based tests and an AddressSanitizer and UBSan build. The short-term-average over long-term-average ratio is still the most widely used automatic event detector (Allen, 1978; Trnkoczy, 2012). A short window tracks the signal's recent energy, a long window tracks the background, and a trigger fires when their ratio crosses a threshold. The *recursive* form replaces the sliding windows with exponential moving averages: STAi = cs|xi| + (1 − cs)STAi−1, cs = 1/Nsta | | Classic (sliding window) | Recursive (exponential) | |---|---|---| | Work per sample | O(1) with running sums | O(1): two multiply-adds | | State per trace | Ring buffer of Nlta samples | Two floats | | Memory of the past | Hard cut-off at Nlta | Exponential decay, time constant Nlta | | Streaming | Needs the buffer carried between chunks | Needs only the two averages carried over | Withers et al. (1998) compare these and other trigger algorithms in detail. `fast-trigger` implements the recursive form in C++23, starts both averages at the first sample's magnitude, and exposes `compute_sta_lta` for a single trace and `compute_sta_lta_batch` for a 2D array of traces, the latter parallelised over traces with OpenMP. Inputs are validated at the boundary: empty arrays, zero windows, STA ≥ LTA and windows longer than the trace all raise `ValueError` in Python. Figure 8 applies it to a synthetic 100 Hz seismogram: Gaussian noise plus a 0.2 Hz microseism, a 6 Hz P arrival at 30 s and a stronger 3 Hz S arrival at 40 s. STA is 0.5 s and LTA is 10 s. For comparison, it also shows ObsPy's `recursive_sta_lta` (Beyreuther et al., 2010), which uses the same recursion on *x*2 instead of |*x*|. *Figure 8. Top: the seismogram (min/max envelope). Bottom: STA/LTA from both characteristic functions, with each one's threshold set for a 1% false-trigger rate on noise alone (dashed). Both trigger on P and S; the |*x*| ratio triggers at 30.23 s and the *x*2 ratio at 30.16 s. The green line sweeps the trace the way a real-time detector would.* Squaring the signal doubles its dynamic range in decibels. Where an arrival lifts the |*x*| ratio to about 3.4, it lifts the *x*2 ratio to about 7.7. The noise-only ratio gets peakier too, though, so the only fair comparison is at an equal false-alarm rate. I calibrated both thresholds on 300 noise-only traces (1.50 for |*x*|, 1.92 for *x*2), then measured how often each detects a P wave of known amplitude over 150 trials per level. *Figure 9. Detection probability against P-wave amplitude at a fixed 1% false-trigger rate. The *x*2 function reaches 50% detection at about 1.8σ; |*x*| needs about 2.25σ, roughly 2 dB more signal. The *x*2 function also picks earlier: a median 0.26 s after onset at 3σ, against 0.39 s.* For this signal and noise model, the energy function is simply the better detector, and no amount of C++ speed makes up for missing an event. Making the characteristic function a parameter (|*x*|, *x*2, or Allen's (1978) function, which adds a weighted derivative term) is the most useful next change to `fast-trigger`. It's also a reminder to benchmark detectors at a fixed false-alarm rate, not at a fixed threshold. The C++ core takes `std::span` inputs and uses an mdspan-style view for batches, with `__restrict` pointers so the compiler can vectorise. Release builds on Linux and macOS use `-ffast-math` by default, which is fine for a detector but not for bit-reproducible research, so `FAST_TRIGGER_STRICT_FLOAT=ON` restores strict IEEE behaviour. There are Hypothesis property-based tests, and an AddressSanitizer and UBSan build. One practical caveat: because the batch function already uses OpenMP internally, calling it from many Python threads at once oversubscribes the CPU. Use one batch call at a time, or set `OMP_NUM_THREADS`. **seismic-linter in brief** `train_df` gets through; T001 doesn't yet treat `ewm` or custom wrappers as safe; T003 only accepts a literal `shuffle=False`.`@verify_monotonicity` and `validate_split_integrity`, cover what static analysis can't see. Kaufman et al. (2012) define leakage as information in the training data that wouldn't be available at prediction time. In time series it usually comes from ordinary-looking code. Normalising with a global mean includes test-period statistics in every training row. Fitting a scaler before splitting does the same thing. And a shuffled train/test split lets the model train on 2023 and be tested on 2019. Standard k-fold cross-validation carries the same risk for time series unless it's done with care (Bergmeir & Benítez, 2012). Across 17 fields, Kapoor and Narayanan (2023) found leakage affecting 294 published papers. *Figure 10. Top: a shuffled split scatters test years (red) through the record, so the model trains on years that come after the ones it's tested on. Bottom: the temporal split used in the South Asia earthquake severity project, where every validation and test event is later than every training event.* *Figure 11. seismic-linter's pipeline. Notebooks are flattened to code cells, each file is parsed once and cached by content hash, and files are analysed in parallel processes.* The linter walks Python's abstract syntax tree and reasons about *context*. A call to `.mean()` is fine if it follows `groupby`, `rolling`, `expanding` or `resample` in the same chain, because those carry temporal boundaries; on a bare DataFrame, it gets flagged. | Rule | Severity | Flags | Fix | |---|---|---|---| | T001 | Warning | Global mean / std / min / max without a temporal predecessor | Rolling or expanding statistics, or statistics fitted on train only | | T002 | Info | `.fit()` on generic names such as`df` ,`data` ,`X` | Split first, fit on `X_train` | | T003 | Error | `train_test_split` without a literal`shuffle=False` | `shuffle=False` or a time-based split | Here's the real output of version 0.2.0 on an eight-line script that makes all three mistakes: ``` bash $ seismic-linter lintdemo/ File: lintdemo/train.py ❌ [Line 6] T001: Global mean() may cause temporal leakage. ❌ [Line 6] T001: Global std() may cause temporal leakage. ℹ️ [Line 7] T002: Model fit() called on potential non-training data 'df'. ❌ [Line 8] T003: train_test_split used with shuffle=True on time series Found 4 potential violation(s). # exit code 1 ``` A non-zero exit code is what turns this into a gate: it can run as a pre-commit hook or a GitHub Action, and `fail-on` / `ignore` in `pyproject.toml` control which rules block a merge. Deliberate exceptions get an inline `# seismic-linter: ignore T001`. Static analysis can't see data, so the package also ships two runtime checks: a `@verify_monotonicity(time_column="time")` decorator that rejects unsorted catalogues, and `validate_split_integrity(train, test, time_column="time")`, which raises `TemporalCausalityError` if any test timestamp is earlier than a training timestamp. The static rules are heuristics, and they have limits. T002 decides from variable names, so a leaky fit on a variable called `train_df` gets through. T001 doesn't yet recognise `ewm` or custom wrappers as safe, so those produce false positives. T003 only accepts a literal `shuffle=False`, so a variable holding `False` still gets reported. A proper fix would be data-flow analysis that follows which rows reach each `fit` call. For a linter whose job is to be run on every commit, false positives you can suppress are a better failure mode than silent misses. | Dataset | Provider | Used here for | Access | |---|---|---|---| | Synthetic velocity steps | Generated with torch-friction | Figures 2, 3, 5 (known ground truth, σ = 10 −4 noise) | This post | | Synthetic seismograms | Generated with NumPy | Figures 8–9 (known P and S onsets) | This post | | Synthetic shot gather | Generated with NumPy | Figure 6 (known reflection and ground-roll components) | This post | | Laboratory stick-slip (biaxial shear, Penn State) | Rouet-Leduc et al. (2017); LANL Earthquake Prediction on Kaggle | Not used yet: the natural real-data target for torch-friction inversion | [doi.org/10.1002/2017GL074677](https://doi.org/10.1002/2017GL074677) | | STanford EArthquake Dataset (STEAD) | Mousavi et al. (2019) | Not used yet: labelled P/S picks for benchmarking fast-trigger on real waveforms | [doi.org/10.1109/ACCESS.2019.2947848](https://doi.org/10.1109/ACCESS.2019.2947848) | I've used synthetics wherever the point was to test an algorithm, because only synthetic data come with exact ground truth. The two "not used yet" rows are where each tool should go next. `examples/inverse_problem_demo.py`, with the changes described in section 2.5. **You can backpropagate through physics. A friction law can be a layer in a neural network.** **Code & archives** • [torch-friction](https://github.com/AmanSinghNp/torch-friction) ([Zenodo DOI](https://doi.org/10.5281/zenodo.18366445)) • [Basalt](https://github.com/AmanSinghNp/Basalt) ([Zenodo DOI](https://doi.org/10.5281/zenodo.18907551)) • [fast-trigger](https://github.com/AmanSinghNp/fast-trigger) • [seismic-linter](https://github.com/AmanSinghNp/seismic-linter) ([PyPI](https://pypi.org/project/seismic-linter/)) *This article was first published at [amansinghnp.com](https://amansinghnp.com/writing/seismic-pytorch).*