Repository navigation
Final pass: kernels, pool, QC, vibroseis, stack and survey QC (0.3.0) - #4
Merged
Merged
Conversation
…ainer Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The four work arrays of a transform lay a whole number of pages apart at every power-of-two length, and so did the input streams of a radix-4 pass inside each: all the streams of a pass shared the same L1 (and for long transforms L2) cache sets. A gap of 17472 bytes between the arrays removes the conflict misses. The butterflies now run in fixed runs of the lane count, so the first stages are straight-line packed code with no remainder loop. 8 lanes of float32, ns per trace and transform, 3950X: n = 2048 2871 -> 1386, 32768 79484 -> 29616, 65536 292520 -> 65206, 131072 2219349 -> 150167. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Element 0 of each region is 64-byte aligned (the allocator gave 8), so a kernel that lays its buffers out in multiples of a line never has a SIMD load or store straddle two. The batched FFT is 15 % faster for it (n = 2048: 1330 -> 1130 ns per trace and transform). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…ked gain loop A full batch of rows of one length is packed and unpacked sample by sample across the lanes, so the lane buffer is written front to back (packing row by row with a stride missed L1 on every store). The trend and the spectrum's demean use the eight-accumulator sums of core/basic, now shared from core/sums. The gain multiply takes its lane count at run time: with a constant GCC vectorized across the bins through permutes, four times the instructions. The spectrum plan keeps the taper's sum. A band-pass of a 1000-sample trace: 6.3 -> 3.7 us. A row's result does not depend on the rows batched with it (new test, bit for bit, three build modes), which the pool relies on. The IIR module said its recursion vectorizes across traces; it does not (GCC keeps the eight lanes scalar). The doc and the vectorization test now say so. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…ny length The fixed gap put the arrays 1088 bytes apart modulo a page only when an array was a whole number of pages. The stride is now the least whole number of cache lines that is 17472 bytes more than a multiple of 64 KiB (review xE P6): the same speed at the lengths measured before, and no length whose arrays fall back into one set. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
One fixed range of rows per worker ran a block whose first half was dead at half speed, and a worker held up by another process held the whole job. The workers now take runs of about 256 KB of samples from a shared counter as they finish: measured by the review, 4 workers on a half-dead panel 2.06x -> 3.51x; on a uniform panel under load the median frame is 15-25 % shorter at 16-32 workers. Results stay bit-identical to the single-thread apply (new test with rows of uneven cost). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
… xE C2, xG C3) A kernel handed a workspace reserved for a smaller plan wrote past it in -d:danger (agc, the filters, the first-break picker, coherence SNR; ASan reproducers in both reviews): the check in `scratch` was compiled only with bound checks on. It is three compares per kernel call. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The Coppens ratio's numerator is the energy of the w samples ending at the sample; over the first w samples that window is cut by the trace's start, so the ratio climbs from 0 to about 1 whatever the data, and the rise picker took it: with noise a tenth of the onset 81 of 200 synthetic traces picked in the first 10 samples, and 3210 of the 3954 traces of the example SEG-D record in the first 12. The rise is now searched from the first whole window (200 of 200 within 3 samples; no pick of the example record in its first window). The burst test put its burst 4 ms into the trace, where only that warm-up could pick it; it now comes a window after the start, with the first breaks 100 ms later. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
A row with a zero head (a recording delay, a padded start) was picked on the step from zeros to data, with confidence 1: the pickers saw energy rise from nothing, and the confidence's noise window, all zeros, read as +60 dB (the stabilization alone). On 100 zero samples followed by noise and an onset, STA/LTA picked the edge on 55 of 100 traces and the energy ratio on 87, each with confidence 1.00. - A row is now picked over its liveSpan (the unfiltered row's first to last non-zero sample), so the edge is outside the search and the data a pick needs before it count from the first non-zero sample. The same case picks the onset on 100 of 100 traces with each method. - A noise window without any energy (zeros inside the trace, a dropout) gives a NaN snrDb and confidence 0 instead of the stabilization's ratio. Rows without leading or trailing zeros pick as before. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
… edge is in `nearest` took floor(u + 0.5). One ulp below a bin's upper edge the sum rounds up to the next integer (pred(0.5) + 0.5 == 1.0), so the position got the next bin: on an axis of one bin that is bin `len`. A lineBinning or cdpBinning of one bin returned 1 of 1, a 3 x 1 inlineCrosslineBinning gave the next inline's bin (and 3 of 3 on the last inline), and StackAccumulator.add's doAssert then ended the process from a prefetch worker of bruteStack. Bin 0 also excluded its lower edge (u > -0.5) while every other bin included its own. The half-way test now uses u - floor(u), which is exact, on the range [-0.5, n - 0.5), with a min() as a guard. tbinning's bin_edges block failed before (binAt(pred(0.5)) == 1 with count = 1) and passes in debug, release and danger. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…d (review xE C1, C3, P3) Edge handling. Padding by odd reflection pins the end samples of a low-cut or band-pass output to H(0)·x0 = 0, whatever the data: against the same samples inside a longer trace the first L/2 samples were off by 41 % (low-cut 2-4 Hz) and 92 % (4-8-60-80 Hz) of the signal there, and a zoom window at margin L by 2.3 % of the peak. When the response rejects 0 Hz and the line was removed, the FFT filter now pads with zeros and the IIR passes start from rest: 4 % and 15 % at the ends, 0.3 % in the zoom window (IIR 4-60 Hz: 48 % -> 7 %, 1.6 % -> 0.2 %). High-cuts, notches and trMean/trNone keep the reflection, which is the better choice for them (measured). With zero pads 2L between the trace's ends is enough (wrap leak -79 dB), so nfft is samples + 2L there: 2500 -> 2250 for 4-8-60-80 on 2001 samples, 5760 -> 3840 for a 0.5-1 Hz low-cut. IIR kernel. The lane loop of 8 was unrolled by GCC and left scalar (48 scalar operations per step). With 24 lanes it is kept and packed: a zero-phase band-pass of a 2001-sample trace 13.7 -> 6.3 us. A state left to decay after the input stops is flushed to 0 every 256 steps instead of running on in subnormal numbers (33639 of 40000 output samples were subnormal after a spike at 0.25 ms). tvectorized now checks the IIR kernel by its share of packed arithmetic: 'some packed instruction' had passed the scalar kernel. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
They indexed the flat storage with bin * samples + sample, so only the product was bound-checked: acc.foldAt(0, 9) of a 2 x 8 stack returned bin 1's sample 1 instead of the documented IndexDefect. Both now check each index on its own when bound checks are on (as RowBlock.rowPtr does), and the message names the bin, the sample and the stack's shape. tstack's refused_traces block: foldAt(0, 5) of a 2 x 4 stack raised nothing before, raises IndexDefect now (debug and release; danger has no bound checks). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
… statics staticsInto(meta, d, statics, staticsMs) on a SEG-Y file with unset datum fields gave -93.0 ms where the keys form gave 6.96 ms: parseSEG's traceStatics returned 0, not NaN, for bytes 53-60 left zero, and a 0 here is a datum at sea level. That is fixed in parseSEG (unset datums become NaN); this side already takes d.datumZ for a NaN datum and must not guess that 0 means unset. The new tstatics block pins the contract the parseSEG fix relies on: columns with NaN datums, columns holding d.datumZ and the form without columns give the same statics (bit for bit), each end falls back on its own, and a datum of 0 stays a datum. The module documentation says so. No code change. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
`plan.rowsNeeded` holds rows of the planning view, but the documented compact read was `rowsOf(f, plan.rowsNeeded)`, which takes trace numbers. A plan made on `chunkMeta(f, 3 ..< 12)` therefore read traces 0..2 for rows 0..2 and cut the wrong segment (sample 500 where the stream at the shot time holds 1500), without an error. The plan now also stores `tracesNeeded`, `traceOf` of each needed row, and the documentation reads from the file with it. `rowsNeeded` stays for blocks already in memory and for `SlicePiece.row`. Test: traces_needed_of_a_view_inside_the_file (a SEG-Y of 4 segments x 3 receivers, planned on traces 3 ..< 12, compact read against the full block). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…ocities Silent misreads of velocity files: - `HANDVEL nan` gave a pick with a NaN location and `HANDVEL 1e400` (or a VELF or CSV location of 1e400 or inf) an infinite one; only CSV refused NaN. Every format now refuses a NaN or infinite location on its line. - A CSV with decimal commas and no header, `100,0,1500,5`, read as a 3D pick of 5 m/s at 1500 ms. A velocity below MinPickVelocity (100 distance units per second, exported and documented) now raises ValueError naming the line, the pick and the likely causes (in CSV: the decimal comma). rmsVelocity and picks built in code keep no bound. - A VELF continuation card that lost a number reads as a new function. The module documentation claimed the CDP order catches that; it does only when a function with a smaller CDP follows. No rule tells the last card of a file from a real function (a first time above the previous function's last is legal), so the documentation now says what is and is not caught, and a test pins both cases. tvelocity's refused_locations_and_velocities block failed before (it did not compile without MinPickVelocity; the HANDVEL and CSV cases parsed). The fuzz target also checks finite locations and the bound; libFuzzer's runtime is not installed here, so the same checks ran over 200 000 random mutations of the seeds under ASan (gcc): 24 721 parsed, the rest ValueError. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
With a fixed or guided search window, pickInto and pick still copied, band-passed and scanned every trace to its end: on 2000 traces of 3000 samples with a 150 ms guided window the IIR over all 3000 samples was 45 of 54 ms (review xG P1). STA/LTA and the energy ratio now work on the stretch their picks depend on: the attribute's and the confidence's windows before the search window plus the causal filter's requiredMargin, and the longer of windowMs and snrWindowMs after it. The stretches are band-passed in the IIR lanes as whole rows were (a meta of the lanes gives the filter their lengths). Coppens still reads the whole trace: its ratio is cumulative from the trace's start and its beta relative to the whole trace's energy. Measured, -d:danger --opt:speed -march=x86-64-v3, one thread, minima of 9 on the reviewer's case (2000 x 3000, 150 ms window): energy ratio 8.2 -> 1.0 ms + 8-60 Hz band-pass 40.7 -> 4.0 ms STA/LTA 30.1 -> 1.3 ms + 8-60 Hz band-pass 62.7 -> 4.3 ms Coppens 15.4 -> 8.1 ms + 8-60 Hz band-pass 47.9 -> 40.7 ms bench/bprestack (2000 x 2001, 150 ms window, medians): STA/LTA 20.1 -> 1.36, energy ratio 5.7 -> 1.04, Coppens 10.3 -> 5.55 ms; with a 5-60 Hz band-pass 42.9 -> 5.84, 28.4 -> 5.21, 33.1 -> 28.1 ms. Results. Unwindowed plans and Coppens: unchanged, bit for bit. Windowed, no band-pass: the same column on every one of 33 600 compared picks, the SNR within 0.002 dB (the 1e-6 stabilization is now relative to the stretch, not to the trace). Windowed with a band-pass: the filter starts at the stretch and is not preceded by a detrend of the whole trace (a line fitted to a short stretch is set by the break itself, so a stretch is filtered as it is), which moves 5 of 12 000 picks of the reviewer's case (four by one sample, one by four) and leaves the mean pick error the same to 0.01 ms. Where the stretch is cut by the trace's start within the filter's margin (a 3 Hz low cut settles over 261 samples) the start-up differs more: up to 7 % of such picks move, mostly by one sample. A windowed pick is now a function of its stretch alone, which the new test checks bit for bit. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…and no pow per sample off the table (review xE P1) A plan compiled for a zoom window's width, as the integration doc advised, tabulated its gain curves over the first samples of the record; the window lies further down, so every sample fell off the table and cost a pow: 38 ms for 2000 x 1101 samples against 0.6 ms on the table. - compile(T, dtMs, samples, recordSamples = 0): the trace length when samples is a window's width; gain-curve tables then cover the record (plan.recordSamples). - applyGainCurve(b, meta, plan, ws), used by the pipeline: the factors off the table are computed once per run of rows with one time axis. The same case without the record length: 2.6 ms. Traces with a recording delay gain the same way. Docs and the viewer simulator pass the record length. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
A viewer picks on the window it decoded, and had no way to know how much of each trace before the search window a pick depends on (review xG A6: windows that start too close to the break change the picks), nor a pooled pickInto. - FirstBreakPlan.requiredMargin: the reach of a pick's stretch before its search window (the attribute's and the confidence's windows plus the band-pass's margin; the reach after it is shorter). A block holding that many samples on each side of every row's search window, as decodeWindow(searchColumns, plan.requiredMargin, traceLen) gives, picks as the whole traces do, bit for bit (tested with STA/LTA and the energy ratio, with and without a band-pass). A plan that reads whole traces (no window, or Coppens) returns plan.samples. - ProcessPool.pickInto(plan, block | tensor, meta, picks): the rows in the pool's runs, each worker reserving plan.workspaceSize; identical picks for 1, 3 and 8 workers on rows of uneven cost, no allocation once warm, clean under ASan. The checks are checkPick, new and exported for it, on the caller's thread. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
A bin size typed as 0.01 for 10 gave a grid of 1.06e10 bins, and `initFoldMap` or `preplotFold` then asked for 42 GB of counts: the process ran out of memory, which no `except` can catch. `surveyGrid`, `initFoldMap` and `preplotFold` take `maxBins` (default 1 shl 26: 256 MB of fold counts) and raise `ValueError` beyond it, before anything is allocated. `surveyGrid` counts the bins as floats first, so a count beyond an int is refused too, and its message names the bin sizes and the midpoint extent; `initFoldMap` names the counts and, for a grid, the bin sizes and the extent measured from the bin centres. The `strays` test asked for the 2.9e9-bin grid that covers a stray station: that is what the budget refuses, so it now checks the refusal and passes `maxBins = high(int)` to get the grid. Test: bin_budget. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
flow[0].agc.windowMs = 300 on a filter op is a FieldDefect in release, but -d:danger has no field checks: the write landed in the filter's spec and compile crashed. A viewer that keeps 'the index of the AGC op' and lets the user reorder the flow is one index away from it. The variant's fields are now private; op.kind, op.agc, op.filter, op.dB, ... are accessors with the same syntax that raise ValueError for an op of another kind, in every build mode. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The prefix and suffix sums of a block are two serial add chains, and that loop was over half of an AGC's time. Two blocks per iteration give the core four independent chains: 1.77 -> 1.31 ns per sample (2001 samples, 500 ms window), 1.65 -> 1.15 at 100 ms, with the same sums in the same order (identical output). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The first-break module said the band-pass delays picks "by a few ms, the same on every trace". It does not: the causal filter's break grows out of the noise, and a picker finds it where it clears the noise. Measured on 400 traces with a 30 Hz onset at 2 ms in a guided window, median pick minus onset in samples at noise of 2.6 / 6.6 / 13 / 26 % of the peak: pmStaLta 1.9 1.9 1.9 1.8 + 8-60 Hz 7.7 8.7 9.5 10.3 pmEnergyRatio 0.7 0.9 1.2 1.6 2.9 3.6 4.1 4.8 pmCoppens 0.6 1.0 1.4 1.7 2.8 3.6 9.7 10.7 The module documentation and the QC guide now give the table and what it means for a spread whose SNR falls with offset. Also corrected: - qc/firstbreak: pmCoppens was called the modified Coppens method of Sabbione and Velis (2010); it takes that method's stabilization term but not its edge-preserving smoothing. - vibroseis/correlate: the example's output of a 10001-sample record and an 8 s sweep was annotated "2001: 4 s listen"; it is 6001, 12 s. - the QC guide called FirstBreakPick.col a fractional column; it is always a whole one. - The price of not picking the step from zeros to data is now stated: a noise-free synthetic is not picked on its first non-zero sample. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
`generateRelations` with `ppNumbers` and a NaN source line raised `IndexDefect` (`fractionalIndex`: `upperBound` of NaN is 0, then `xs[-1]`), an out-of-bounds read in a `-d:danger` build. An X record names its source, so a shot whose line or point is not finite now gets no records in any placement and edge mode, as a shot without coordinates gets none by coordinates. Two guards beside it: `fractionalIndex` returns NaN for a NaN position instead of indexing, and `roundUp` stops at +-1e15, so a finite number map that sends the source to 1e303 no longer converts that to an int (undefined in C; it gave no records by luck). Output for numbered shots is unchanged: 6215 template, box and static cases plus five Morgantown runs give the same records as before. Test: shots_without_a_station_number. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
- ws.reserve(plan) for any plan with a workspaceSize, not the pipeline's alone; - dtMs on IirPlan and SpectrumPlan, as the other plans have; - an AGC window is at most 2^26 samples each side (its tables were allocated from the window alone, up to int32.high entries); - filter/apply's example called applyFilter with the wrong argument order; - the pool model cites procs, not line numbers that every edit moved. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…ne stack
Each prefetch worker owned a full StackAccumulator, merged at the end. On
a stack of 40 000 bins x 1501 samples (60 000 traces of 1501 samples,
danger, x86-64-v3) more workers were slower and the job held workers + 1
stacks: 583 / 523 / 846 ms and 1.3 / 2.9 / 5.7 GB peak RSS for 1 / 4 / 8
workers. The result also depended on which worker took which chunk: runs
agreed to rounding only.
Now the transform does statics, the flow and NMO on the worker and writes
the corrected trace back into its chunk row (rows not stacked are zeroed);
the calling thread takes the chunks in order and adds each row to its bin
in one accumulator, whose storage, normalized in place, becomes the
result (StackAccumulator.finishInto). Chunks are as wide as the longer of
the input trace and the stack; for a longer stack the rows are moved to
the front first, so the flow sees the block it was compiled for.
Measured with the same setup (min of 5, peak RSS includes the 389 MB file
mapping):
1396 bins: 502 / 149 / 101 ms, 425 / 487 / 550 MB
-> 432 / 119 / 91 ms, 408 / 413 / 421 MB
40 000 bins: 583 / 523 / 846 ms, 1313 / 2932 / 5711 MB
-> 496 / 167 / 142 ms, 855 / 861 / 867 MB
The stack is bit-identical for 1, 2, 4 and 8 workers, and to the old
1-worker result in all six benchmark modes (plain, statics, flow, all,
field, cubic: same hashes), since the order of the sums is the old
1-worker order.
A bin from the keys is checked against the stack before it reaches
StackAccumulator.add on either thread: a Defect there ends the process.
StackAccumulator keeps its sums in tensors so that finishInto can hand
them over; its =copy and =dup clone them, so a copy stays a full copy.
stackInto divides by max(fold, 1) (a sum of fold 0 is 0), which makes the
normalization loops straight-line; merge runs on a pointer kernel.
Tests: tbatch same_bits_for_any_worker_count (failed before, with 2
workers), equals_the_stack_by_hand (also onto a longer and a shorter time
axis), memory_bound (it encoded workers + 1 stacks: now one); tstack
finish_without_a_copy. tbatch also runs clean under ASan and TSan.
bruteStackMemory and the memory and reproducibility statements of the
module, the guide and integration.md are updated.
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…annels Four cases where `recordCompleteness` reported a record wrongly: - A header source station that no X record lists gave "no X record" under `rmAuto`, although the FFID named one. `rmAuto` now falls back to the FFID when the station has no candidate at all (`recordMatch` says `rmFieldRecord`, `sourceMismatch` reports the station); an explicit `rmSourceStation` does not fall back, and several candidates stay ambiguous. - A station of line 0, point 0 (a recorder that left the fields empty) counted as a station: source 0/0 found no record, receivers 0/0 made one station of the whole record (0 present, 9 duplicates). The module treats NaN as "no station", so both zero is now no station either, for the source and the receivers; the traces then match by channel. - An auxiliary row that failed to decode had the status "failed", so it took part in channel matching (a duplicate of seismic channel 1). The header status decides: auxiliary rows never match; a failed one still counts in `failed`. - A header channel of NaN counted as a channel mismatch on every station. It has no number to compare. Tests: auto_falls_back_to_the_ffid, zeroed_stations_are_no_stations, unknown_channels, failed_auxiliary_row. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…rgins as measured now Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
With `rcvX = Inf` the source-receiver distance is infinite, so the tolerance `max(tolerance, relTolerance * distance)` was infinite too and `|offset - distance| <= tolerance` held: a trace with a header offset of 100 and a receiver at infinity got no flag. The relative part now applies to finite distances only, so the trace is flagged `gfOffsetMismatch` (distance Inf, offsetError -Inf). An infinite coordinate stays a known, wrong value: it is not reported as `gfNoCoordinates`. Test: infinite_coordinates. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
… a ValueError `sliceKeys` copied the receiver's keys and replaced the source line and point by the shot's, but kept the segment's `srcIndex`: a node file whose general header holds source 77/88 index 3 gave records "9/4 index 3". `SliceShot` gains `srcIndex` (0: unknown, the default of `initSliceShot`), and the output keys take it. A row whose sample count is near `high(int)` (a caller's array: no file holds such a trace) overflowed `position + length` on the sample grid: an `OverflowDefect` in checked builds, a wrapped position in `-d:danger`. `planSlices` now raises `ValueError` for a row of more than 2^48 samples. Tests: slice_keys_take_the_source_from_the_shot_alone, a_sample_count_beyond_any_file. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
fitTwoSegment took a crossover whenever it lowered the L1 misfit by a part in 10^6. Picks that one line fits exactly have a misfit of rounding size (1e-13 ms), which any split lowers by more than that: 24 flat picks came back as two segments of 2e19 and 8e19 units/s meeting at 425 (review xG A8; in release and danger builds, where the fit's sums are vectorized). A second segment must now also lower the misfit by more than 1e-9 ms a pick, far below any pick's precision. Flat picks and picks on one sloping line give crossover Inf and v2 NaN; tested in debug, release and danger. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Four names or types of the QC modules said something other than what they hold (review xG A4). Breaking, and wanted before a viewer depends on them: - FirstBreakPick.col is an int, -1 without a pick (it was a float64 that only ever held a whole column, NaN without a pick). timeMs keeps the float. The fields now default to "no pick" (NaN, -1, NaN, 0), so FirstBreakPick() equals noPick(). - AmplitudeParams.spikeWindowMs -> spikeHalfWindowMs: the neighbours are those within that many ms on each side, while every other window of the library is a full or a one-sided length. - ShotFitParams.residualMs -> maxResidualMs: a threshold, next to ShotFit.residualMs, which holds residuals. - ShotFit.lagMs: its sign is documented (the neighbour's arrival minus the trace's, each from its model time: negative for a trace that arrives late against its neighbour) and tested with a trace shifted 10 ms each way. Callers updated: tests/tfirstbreak, tproperty, tamplitude, tshotfit, docs/api.md and the QC guide. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The block forms compared a block's interval with their plan's each in its own way (review xG A5): pickInto and RecordStack.merge exactly (!=), correlation and the IIR and FFT filters within a relative 1e-9, addRecord not at all until the previous fix. An interval read from a header (2000 us / 1000) and one computed (1000 / 500 Hz) can differ in the last bit, and then one plan took the block and the next refused it. checkPick (pickInto, pooled or not) and RecordStack.merge now use the filters' test, abs(a - b) > 1e-9 * b, as checkCorrelate and addRecord do. The expression is written out at each of the four sites, as it is in filter/apply, filter/iir and gain; a shared helper belongs in core/types. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
- recordAmplitude(traces, params = AmplitudeParams(), ws): the default on `params` could never be used, because the required `ws` follows it. `params` is now required, as it is in traceAmplitudeInto. No call changes: every caller already passed all three. - recordStackMemory(T, rows, samples, mode) took a sample type that it ignored (the stack's sums are float64 whatever T). Breaking: it is now recordStackMemory(rows, samples, mode = rsVertical). The rest of A2 (where `ws` and `dst` stand in pickInto and coherenceSnrInto, recordSamples before dtMs in initCorrelationPlan) is left as it is and reported as a recommendation. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Interim report of worktree `stack` (review xF), after the correctness findings, P1, P6 and the signature-changing API items. Still to come, in this order: P4 (shiftSamples), P2 (NMO kernel), P3 (velocity field columns), API 6 (bin midpoints, trace counts), API 3 (binning steps), API 2 (NMO preview by binning), bench/README. A final report follows. FINDINGS C3 fixed, 52e6950 (binning) + f428aba (batch guard). `nearest` now tests u - floor(u) >= 0.5 (exact) on [-0.5, n - 0.5) with a min() guard, instead of floor(u + 0.5): no position gives bin `len`, bin 0's lower edge is in. Exact form chosen over the sketched min(floor(u+0.5), n-1), which still put pred(0.5) in bin 1 of 4. Test: tbinning `bin_edges` (failed before: binAt(pred(0.5)) == 1 with count = 1). Decision on stackChunk: no Defect is caught; a bin from the keys is checked against the stack on both threads before StackAccumulator.add (worker: `stackable`, consumer: `stackChunk`), so no data value reaches the doAssert. Reviewer's rbin1/rdefect no longer reproduce. C4 fixed, a8051b4. foldAt/sumAt check bin and sample separately (IndexDefect with bound checks on). Test: tstack `refused_traces`. C1 (processSEG side) done, 41a6b03: no code change; tstatics `unset_datums_take_the_requested_datum` pins that TraceStatics with NaN datums, with datums == d.datumZ and the keys-only form agree bit for bit, and that a datum of 0 stays a datum. With the parse-stable snapshot the reviewer's rstat1 still shows -93 ms: it needs the parseSEG fix (unset SEG-Y datums -> NaN). C5 fixed, f848bd0. HANDVEL, VELF and CSV refuse NaN and infinite locations ("a NaN location" / "an infinite location", with the line); a velocity below MinPickVelocity = 100 (exported const) raises ValueError naming line, pick and likely causes (CSV: decimal comma). VELF last-card continuation: not detectable safely (a first time above the previous function's last time is legal), so the doc is corrected and a test pins both cases. Test: tvelocity `refused_locations_and_velocities`. Fuzz target fvelocity checks the new invariants; libFuzzer's runtime (libclang_rt) is not installed on this machine, so I ran the same checks over 200 000 random mutations under ASan/gcc instead: no finding. C2 (processSEG side) done, e612a95: thread contract in batch.nim's module doc, the guide and integration.md "Batch results" (+ one sentence in the viewer checklist item 9, outside my sections). Placeholder "parseSEG 0.10.0" in those three places, to confirm at release. P1 fixed, f428aba. Workers do statics, flow and NMO into the chunk row; the calling thread stacks chunks in order into one accumulator whose storage, normalized in place, becomes the result (StackAccumulator.finishInto). Chunk width max(inSamples, samples); for a longer stack the rows are compacted first so the flow sees the block it was compiled for. Measured (danger, x86-64-v3, 60 000 traces x 1501, min of 5, 1/4/8 workers; peak RSS includes the 389 MB mapping): 1396 bins: 502/149/101 ms, 425/487/550 MB -> 432/119/91 ms, 408/413/421 MB 40 000 bins: 583/523/846 ms, 1313/2932/5711 MB -> 496/167/142 ms, 855/861/867 MB Bit-identical for 1, 2, 4, 8 workers, and identical to the OLD 1-worker result in all six benchmark modes (plain, statics, flow, all, field, cubic: same hashes): the order of the sums is the old 1-worker order. Tests: tbatch `same_bits_for_any_worker_count` (failed before, at 2 workers), `equals_the_stack_by_hand` (also on a longer and a shorter time axis), `memory_bound`; tstack `finish_without_a_copy`. tbatch clean under ASan and TSan (gcc). P6 fixed, 6b3a603. allFinite as an exponent-bit OR (kernel 1.07 -> 0.07 us per 1501 f32; addLive of 2000 x 2001 into 200 bins 3.48 -> 1.13 ms), maxFold on the fold's buffer (release: 9.2 -> 0.3 ms for 1396 bins), once per bin in writeStack (release 32 -> 14 ms; danger unchanged at 10.7 ms). liveSpan once per trace comes with P1. Results unchanged (same hashes). Tests: tstack `non_finite_samples_refused_anywhere`, tbatch `max_fold_of_any_layout`, tvectorized checks `allFinite`. API 5 done, 923a093 (BREAKING). API 4 done, e443142. API 8 done, b1c0fad. pC C8 done, 123dc4b: StackResult.intervalUs, written by stackKeys (1001 us gave 1000.9999999999999). Test: tbatch `interval_kept_exact`. API 12: nothing in code. PLAN.md's architecture listing still says "batch.nim bruteStack, first-break pass over a file": there is no first-break pass in batch.nim. P5 skipped as briefed (the flow's mutes and gains would see shifted times). API 7, 9, 10, 11: no change (7 needs pipeline.nim; 9 and 10 are renames across the whole API; 11's one cost, maxFold twice per bin, is in P6). PUBLIC API CHANGES SO FAR - BREAKING batch.BruteStackOptions: `useStatics` removed; `statics` is `Option[DatumStatics]` (some(datumStatics(...)); none = no statics). batch re-exports std/options. - batch.bruteStack: two new overloads taking `traces: openArray[int]` (with VelocityField / VelocityFunction). - batch.bruteStack: raises ValueError when a VelocityField ends before the stack's last sample (was accepted; velocity froze). Result is now bit-identical for any worker count. - batch.bruteStackMemory: one stack instead of workers + 1, chunks of max(inSamples, samples); raises ValueError for bad arguments or overflow (was OverflowDefect / wrapped in danger). - batch.StackResult: new field `intervalUs`. - stack.StackAccumulator: new `finishInto(data, fold, norm = fnFold)`; exported `=copy` / `=dup` (sums now live in tensors; a copy is still a full copy). `stackMemory` raises ValueError for negative or overflowing sizes. `foldAt` / `sumAt`: IndexDefect per index. - binning.binAt / binOf: bin 0's lower edge inclusive; one ulp below a bin's upper edge stays in that bin. - velocity/formats: new const `MinPickVelocity`; parseHandvel, parseVelf, parseVelocityCsv, parseVelocities, readVelocities raise ValueError for NaN / infinite locations and velocities < 100. - velocity/model: new `VelocityField.tMaxMs`. FOR THE LEAD WHEN MERGING - src/processSEG.nim (not mine) still describes batch as "(prefetch transform, per-worker stacks)": now "stacked on the calling thread". - "parseSEG 0.10.0" placeholders: src/processSEG/batch.nim module doc, docs/guide/velocity-and-stack.md, docs/integration.md. - docs/integration.md checklist item 9: one sentence added. - Callers of opts.useStatics elsewhere: none in this repo. - The shared machine was memory-starved at times by another worktree's test (index-tmp/bin_rel/tviewerread at 23 GB RSS). SUITE `nim test` (debug) at e612a95: all 50 test files passed (48 tests + 2 reject tests). The tests of the touched modules (tbinning, tstack, tstatics, tvelocity, tbatch) also pass in -d:release and -d:danger; tstack and tbatch under -d:addressSanitizer; tdocs and tvectorized pass. `nim testRelease` is still to run, after the kernel changes. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…lper each `preplotFold`, `describePatches` (and `compareRelations`) and `recordCompleteness` walk a record's relations channel by channel with parseSEG's `relations`, which scans the record's earlier runs for every channel. parseSEG is about to export a run-level iterator; to let it be swapped in with a few lines, each module now does that walk in one private iterator with the yield type of `relations` itself: qc/fold.nim channelRelations(g, record, domain) survey/patch.nim channelRelations(g, record) qc/record.nim channelRelations(g, record) Replacing the body of each (a loop over runs and their stations) is the whole change on this side; the callers stay as they are. The remaining `continue` statements of the two walks became conditions. No behaviour change: on Morgantown the fold map, the patch description and the completeness report hash the same as before, and the times are within noise (release, min of 5: preplotFold 1105 -> 1065 ms, describePatches 2554 -> 2395 ms, 100 x recordCompleteness 63.6 -> 62.6 ms). A helper that yields a re-packed tuple was 6 % slower in preplotFold, so the helpers pass the relation through unchanged. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
shiftSamples ran one scalar loop with per-sample bounds and "saved sample" branches: 1.46 ns per sample linear, 3.7-3.9 cubic. Now: - a whole number of samples is a move (memmove and a zero fill); - the few outputs with a tap outside the row (at most 1 per end linear, 3 cubic) are computed first, from the untouched row, by the general formula; - the interior, where every tap is in the row, runs branch-free in blocks of 256 into a stack buffer (so the row and the output cannot alias and the loop vectorizes), copied back with each block's last sample held back until the next block is computed (the one sample the cubic reads next to the write position); ascending for shift > 0, descending otherwise, as before. The sums stay in float64 whatever T, so float32 rows keep their bits: 0 of 28 million samples differed from the old kernel over 120 000 random shifts, and bruteStack with statics gives the same hash as before. float64 rows differ by at most 1 ulp (4.4e-16) in 0.03 % of samples (which product the compiler folds into the multiply-add). Measured (danger, x86-64-v3): one row of 2001 float32 samples, linear 2.93 -> 0.84 us, cubic 7.3-7.8 -> 1.25 us; float64 2.7 -> 0.53-0.72 and 7.0 -> 0.64-0.79 us. Panel of 2000 x 2001 float32: applyStatics linear 6.95 -> 2.69 ms, cubic 15.3 -> 3.48 ms; applyLmo + removeLmo 11.1 -> 1.7-2.9 ms. In bruteStack (60 000 traces, 1 worker) datum statics now add 50 ms to the plain stack (425 -> 475 ms); at the start of this branch they added 150 ms (502 -> 653 ms). Interpolating float32 rows in float32 would halve the row times again (0.34-0.46 us) but change 18 % of the samples by up to 4.8e-7: not done. tshift random_shifts_equal_reference (6000 random lengths and shifts around the block size, both types, against the copy-based reference) and whole_shifts_move_any_sample pass before and after; ASan clean. Both kernels hold packed ymm instructions in release and danger builds. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Interim report of the qc engineer (review xG), branch fix/qc on bd27552. Correctness and API items are done; P4, P2 and the rest of P5 follow and will be covered by a later "REPORT fix/qc (final)" commit. FINDINGS C1 fixed 2cf5eda. pmCoppens rise search starts at the first whole numerator window; doc "one window". Tests (tfirstbreak): coppens_does_not_pick_its_own_warm_up (200 traces, noise 0.1 of the onset: no col < w, >= 95 % within 3 samples) and coppens_on_the_example_record (0 picks in the first window; was 3210 of 3954 in the first 12 samples). The burst test moved its burst from 4 ms to 24 ms (a window after the start), with the reason in the commit. C2 fixed d4b0c7d. (a) noise window without energy -> snrDb NaN, confidence 0. (b) a row is picked over its liveSpan (taken on the unfiltered row), so the zero-to-data edge is outside the search. Test zeros_before_the_data_are_not_noise: edge picks 55/100 (STA/LTA) and 87/100 (energy ratio) with confidence 1.00 -> 0/100, onset 100/100 for all three methods, with and without band-pass. CONSEQUENCE TO KNOW: a noise-free synthetic, whose first break is the first non-zero sample, is no longer picked on it (median error 0.5 -> 5.5 samples for the energy ratio, 16.8 for STA/LTA, Coppens none; confidence 0). This follows from the requested fix in either form (the reviewer's `lo = first + w` does the same); documented in the module and the guide. All tests and guide examples use noise. C3 not mine (closed by the lead on the base commit). C4 fixed b426927. sweepFromHeaders validates at the trace's interval (maxNyquistFraction 1) and, with the new `samples` parameter, refuses a sweep longer than the trace (not "of another length": an aux pilot trace is often as long as the record). pilotFromTrace also drops a description whose upper half-band holds < 5 % of the trace's energy: without that the 0.1 Hz header mistake passes at <= 0.5 ms (tested with a written file at 0.25 ms). resamplePilot judges the samples (energyAbove, new) and keeps a description only if it still validates; a pilot without a description is now checked too, and `bandLimit = true` (new) is the way to cut a pilot on purpose. sweepOverrides RAISES ValueError for a fractional Hz or ms instead of rounding: the fields are read back as Pilot.sweep, and a 7.5 Hz start written as 8 is a wrong sweep nothing can detect. Tests (tpilot): headers_that_do_not_fit_the_trace, the_example_pilot (pilotFromTrace(Sweep0.sgy).resamplePilot(2).len == 6001, skipped when absent), sweep_overrides, resample_down_and_up. C5 fixed 0e37697. Passband edge = cutoffFraction - 2.736 / halfWidth of the lower Nyquist; energy above it is refused (or cut with bandLimit) and the kept description gets the edge as maxNyquistFraction; taps beyond the pilot's ends are made up for by rescaling (0.61 -> 1.0000). DEFAULT halfWidth 24 -> 32 so that the default passband (0.8145 Nyquist) holds a default sweep (0.8); at 24 such a sweep would be refused. Test resampling_keeps_its_passband_flat: within 0.001 dB at 150/190/199 Hz, 10-210 and 10-220 Hz raise. C6 fixed 437c939. Both slopes go to butterBandPass. Test in instrument_filter_causality (3 Hz at 6 dB/oct with 200 Hz at 72). C7 fixed 1352427. sweepSamples validates first; sweepLen returns 0 for a spec with an invalid interval or length. Test specs_built_as_objects (debug, release, danger). C8 fixed 8a3b05a. addRecord(meta form) raises for another dtMs (relative 1e-9) or firstSample != 0. Test in bad_arguments (tvstack). C9 fixed a1d14ee. Each polarity has its own clip level when both reach the maximum magnitude within 1 % (int24 sine: 500 -> 1000 clipped samples counted); afClipped needs two runs or one of 2*clipRun (`clipped` still counts every run). Tests clips_of_both_polarities, one_short_run_at_the_maximum_is_not_a_clip, the_example_record_is_not_clipped (4 flagged -> 0). Two existing tests that encoded "one run of clipRun is a clip" were changed (said in the commit). C10 documented dfcb935, measured on the current code in a guided window (400 traces, 30 Hz onset, 2 ms; median pick - onset in samples at noise 2.6/6.6/13/26 % of the peak): pmStaLta 1.9 1.9 1.9 1.8, with 8-60 Hz 7.7 8.7 9.5 10.3; pmEnergyRatio 0.7 0.9 1.2 1.6, with filter 2.9 3.6 4.1 4.8; pmCoppens 0.6 1.0 1.4 1.7, with filter 2.8 3.6 9.7 10.7. In the module doc and docs/guide/qc.md. C11 fixed 247b9b7. validate bounds the total boost (300 dB for dB/octave, 150 dB for dB/Hz; the dB/Hz limit was 600); frequencyAt returns f1Hz/f2Hz exactly at the ends and stays in band; the dB/octave law is computed from f2/f1 (no overflow at high bands; a rounding-level change of its samples). Test steep_boosts_stay_finite. C12 fixed dfcb935 (correlate.nim "6001: 12 s listen"; qc.md column; the pmCoppens description). P1 fixed 417b0bb. Windowed STA/LTA and energy-ratio picks read, copy and band-pass only their stretch. -d:danger --opt:speed -march=x86-64-v3, one thread, minima of 9, the reviewer's case (2000 x 3000, 150 ms window): energy ratio 8.2 -> 1.0 ms, with 8-60 Hz 40.7 -> 4.0 ms; STA/LTA 30.1 -> 1.3 and 62.7 -> 4.3 ms; Coppens (whole trace by its definition) 15.4 -> 8.1 and 47.9 -> 40.7 ms. bench/bprestack medians: 20.1/5.7/10.3 -> 1.36/1.04/5.55 ms, with 5-60 Hz 42.9/28.4/33.1 -> 5.84/5.21/28.1 ms (bench/README.md). PICKS: unwindowed plans and Coppens bit-identical; windowed without band-pass the same column on all 33 600 compared picks, SNR within 0.002 dB; windowed WITH band-pass NOT identical: on fb6's case 5 of 12 000 picks move (four by 1 sample, one by 4), mean pick error equal to 0.01 ms. Cause: the filter starts at the stretch and the stretch is not detrended (a line fitted to a short stretch is set by the break). Where the stretch is cut by the trace start inside the filter's margin (3 Hz low cut: 261 samples) up to 7 % of picks move, mostly by one sample. Exact identity would need the whole trace filtered. Test a_windowed_pick_reads_only_its_stretch (bit for bit). A6 fixed 5a45efc. FirstBreakPlan.requiredMargin; a block holding it on each side of the search window picks bit-identically to whole traces (test a_zoom_window_with_the_required_margin_picks_the_same); plan.samples for whole-trace plans (no window, or Coppens). ProcessPool.pickInto(plan, block | tensor, meta, picks), test the_pool_picks_the_same (1, 3, 8 workers; ASan clean; no allocation once warm). New exported checkPick. P5 part: resamplePilot fixed 2b46151 (kernels tabulated per phase): 30 s at 0.25 ms -> 4 ms 358 -> 5.5 ms (danger, v3, min of 3); test resampling_by_tabulated_kernels (1e-12 against a tap-by-tap reference, ten ratios). coherenceSnrInto and fitShot: pending. A3 fixed a476d41 (breaking). A7 fixed 223f9c5. A8 fixed 4bddbd0 (the defect shows in release and danger builds only). A4 fixed fd8012d (breaking). A5 fixed 07a13be. A2: the two clear cases fixed 5c9c39e; the rest is a recommendation (below). P2, P4, P5 (coherence, fitShot): pending, next. P3: deferred (ideas of the review: split windowPowers' second loop into head / interior / tail so the interior vectorizes; median of every h-th window power). PUBLIC API CHANGES (B = breaking) qc/firstbreak B FirstBreakPick.col: float64 (NaN: none) -> int (-1: none); the fields default to "no pick", FirstBreakPick() == noPick(). + requiredMargin(plan), checkPick(plan, b, meta, picks). Behaviour: liveSpan rule (C2), windowed stretch (P1), interval within a relative 1e-9 (A5). core/pool + pickInto(pool, plan, b | samples, meta, dst). qc/amplitude B AmplitudeParams.spikeWindowMs -> spikeHalfWindowMs. B recordAmplitude(traces, params, ws): `params` lost its (unusable) default; no caller changes. afClipped semantics (C9). qc/shotfit B ShotFitParams.residualMs -> maxResidualMs. vibroseis/pilot B InstrumentPhase, ipMinimumPhase, ipZeroPhase removed: InstrumentFilter.phase is filter/iir's IirPhase (ipMinimum, ipZero), default ipMinimum. B sweepFromHeaders is a proc with a new last parameter `samples = 0`; returns none for descriptions that do not validate. B sweepOverrides raises ValueError for fractional Hz / ms. B resamplePilot(p, dtMs, cutoffFraction = 0.9, halfWidth = 32, bandLimit = false): default halfWidth was 24; refuses by the samples and at the passband edge; halfWidth must leave a passband. + energyAbove(pilot, fHz). applyInstrumentFilter takes unequal slopes. vibroseis/sweep B validate: boost bounds (300 dB per-octave total, 150 dB per-Hz total; was 600 for dB/Hz). sweepLen returns 0 on an invalid spec. vibroseis/taper SweepTaper.param defaults to NaN (the shape's default). vibroseis/vstack B recordStackMemory(rows, samples, mode = rsVertical): the ignored `T` is gone. B addRecord (meta form) raises for another interval or a window. adapters pilotFromTrace: sweep only when the trace can be it (C4). A2, RECOMMENDATION NOT APPLIED. Rule read from the Into procs outside my area (extractInto, amplitudeInto, averageAmplitudeInto, FFT forward / inverse, PipelinePlan.apply): a compiled plan first; then the input and its meta; then parameters; then the output; `ws` the last required parameter; defaults after it; constructors (T, spec, dtMs, samples, ...). Against it, still as they were: plan.pickInto(b, meta, ws, dst) -> (b, meta, dst, ws) coherenceSnrInto(b, meta, windowMs, ws, dst, floorDb = ...) -> (..., dst, ws, floorDb) initCorrelationPlan(T, pilot, recordSamples, dtMs, ...) -> (T, pilot, dtMs, recordSamples) pool.pickInto follows plan.pickInto without `ws`, so it would not change. Conforming already: windowSnrInto, traceAmplitudeInto, fitShot, correlateInto (it mirrors PipelinePlan.apply: src, dst, meta, ws), initFirstBreakPlan, initRecordStack (as initStackAccumulator). FOR THE MERGE - Other modules' callers of the renamed names: none left in this tree (grep clean for spikeWindowMs, ShotFitParams.residualMs, InstrumentPhase, recordStackMemory(T, ...)). tests/tproperty.nim (not mine) was adapted to the int column. - tests/tdocs.nim resolves `export a.B` of a module to a sibling file, so pilot.nim does not re-export iir.IirPhase; the umbrella exports both. - A5: the interval test abs(a - b) > 1e-9 * b is written out at four sites of mine (checkPick, checkCorrelate, addRecord, RecordStack.merge) as in filter/apply, filter/iir and gain; a helper in core/types would replace all of them. - The lead's final-pass changes that touch my code at merge: the IIR's 24 lanes (pickInto batches fftLanes stretches per applyFilter call through a lanes meta of fftLanes keys: correct with any lane count, but the batch size should follow the IIR's), and zero padding for low-cut / band-pass filters (changes the start-up of a stretch cut by the trace start, so the P1 pick differences near the trace start will change). - bench/README.md has a new subsection "First-break picking in a search window"; bench/bqc.nim (new, QC and vibroseis) comes with the performance commits. SUITE - nim test (debug) at 5c9c39e: "All 50 test files passed" (48 tests, 2 reject tests). - -d:release and -d:danger at 5c9c39e, the tests of every module touched (tfirstbreak, tshotfit, tsnr, tamplitude, tsweep, ttaper, tpilot, tcorrelate, tvstack, tpool, tproperty, tadapters): 24 of 24 pass. The whole `nim testRelease` is left for the final report. - tfirstbreak also under -d:addressSanitizer (after the pooled pickInto): clean. - nimpretty --indent:2 --maxLineLen:100 run on every changed .nim file. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…eriors nmoInto computed each input column with two divisions (x / V, / dt) and a sqrt, computed the columns a second time, one by one, for the stretch scan, and interpolated every sample through a general path with bounds tests, in float64. Now: - the kernel works on the sloth 1/V2: x2 * sloth and * (1/dt), no division. nmoInto(vrms) inverts its velocities per block of 256 (one vectorized division); the new nmoSlothInto takes sloths, for a caller that corrects many traces with one column (vrmsToSloth converts it once, in place). applyNmo and bruteStack do so. Both forms give the same bits. - the stretch scan reads the vectorized columns, block by block, from the first output time >= 0, with the old test (d > 0 and 1/d - 1 <= maxStretch) on the same differences: the mute starts where it did. - a block whose columns all have their taps inside the input (counted by a vectorized pass) is interpolated branch-free: whole parts, fractions and cubic weights in vectorized loops (floor, hence -fno-trapping-math for this module's C code, next to -fno-math-errno), then a loop of plain loads. float32 rows run it in float32. Blocks that touch the input's ends or hold a NaN column keep the general path. Measured (danger, x86-64-v3, one 1501-sample float32 trace, offsets 200/1500/2500): linear with a 0.3 stretch mute 5.56/6.48/7.07 us -> 3.49/3.58/3.42 (velocities), 2.99/3.07/2.91 (sloths); without a mute 5.2 -> 3.2 / 2.8; cubic 9.5/9.1/8.6 -> 4.8/4.5/4.0 / 4.3/4.0/3.5. float64 rows: linear 5.2-6.9 -> 3.2-3.5 / 2.5-3.0, cubic 8.0-8.3 -> 4.1-4.9 / 3.6-4.5. bruteStack of 60 000 x 1501, 1 worker: plain 427 -> 210 ms, cubic 266 ms (695 at the start of the branch). Rounding: over 20 000 random cases per type the mute start is identical in all; float32 samples differ from the old kernel by at most 7.7e-7 relative to max(1, |sample|) on unit-variance data, float64 by 2.3e-12. On the stacks of the benchmark line: at most 2.4e-7 (plain), 3.0e-7 (cubic) for a peak of 1.08, and 5.6e-7 relative to max(|sample|, rms). A float32 sample's last bit can depend on whether its block took the fast path (e.g. on where a mute starts). Tests: tnmo nmo_matches_a_reference now checks linear and cubic, both types, the samples and the mute start against a direct evaluation, and the sloth form bit for bit against the velocity form. Its float64 tolerance goes from 1e-12 to 1e-9: the reference divides where the kernel multiplies by a reciprocal, so the columns differ by a few ulps; that is the rounding this change is allowed, not a looser kernel. New block unknown_velocities_and_input_ends. tvectorized checks moveoutCols, slothCols, insideCount and splitCols in release and danger. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Interim report: correctness and API items done; performance P4 and P7 still to come (P8 is in). Branch fix/surveyqc, 14 commits on bd27552. Full suite `nim test` (debug) at 5127cd0: all 51 test files passed (48 tests, 3 reject tests). The tests of the touched modules also pass in -d:release and -d:danger; tslice also under -d:addressSanitizer. FINDINGS C6 fixed, 6333dfc. `plan.tracesNeeded` (traceOf of each needed row); `rowsNeeded` kept (rows of the planning view: SlicePiece.row, blocks in memory). Test tslice traces_needed_of_a_view_inside_the_file (SEG-Y of 4 segments x 3 receivers, plan on chunkMeta(f, 3 ..< 12), compact read against the full block). C9 fixed, 6d5d99a. `maxBins = 1 shl 26` on surveyGrid, initFoldMap and preplotFold; ValueError names bin sizes and extent, nothing allocated. Test tfold bin_budget. The existing `strays` test asked for a 2.9e9-bin grid: it now checks the refusal and passes maxBins = high(int). Not budgeted: offset classes and azimuth sectors multiply the memory (bins x layers); a mistyped azimuthSectors can still exhaust memory. C13 fixed, 3f0c613. Shots whose source line or point is not finite get no X records in any placement and edge mode; fractionalIndex is NaN-safe; roundUp stops at +-1e15 (a finite number map of 1e300 converted 1e303 to int: undefined in C). Test tpatchgen shots_without_a_station_number. Output for numbered shots identical (6215 cases + 5 Morgantown runs). C7 (processSEG side) fixed, 1c0373f. rmAuto falls back to the FFID when the header source station has no X record at all (recordMatch then says rmFieldRecord, sourceMismatch reports the station); rmSourceStation asked for explicitly does not fall back; several candidates stay ambiguous. CHOICE on (0, 0): the existing code treats a NaN station as absent (hasSource false, no sourceMismatch), so a station whose line and point are both 0 is treated as absent too, for the source and for the receivers (traces then match by channel). Tests trecord auto_falls_back_to_the_ffid, zeroed_stations_are_no_stations. C14 fixed: - geometry Inf coordinate, 25cdba3: the relative tolerance applies to finite distances only; the trace is flagged gfOffsetMismatch (an Inf coordinate is a known wrong value, not gfNoCoordinates). Test tgeometry infinite_coordinates. - record, 1c0373f: an auxiliary row never takes part in matching, failed or not (it still counts in `failed`); a NaN header channel is no channelMismatch. Tests trecord failed_auxiliary_row, unknown_channels. - slice, 5acd2ba: sliceKeys takes srcIndex from the shot (new SliceShot.srcIndex, 0 = unknown); a row of more than 2^48 samples is a ValueError, not an OverflowDefect. Tests tslice slice_keys_take_the_source_from_the_shot_alone, a_sample_count_beyond_any_file. - keys.nim SEG-D Rev 3 timestamps: parseSEG, not mine. P8 fixed, adad4ca. copyMem/zeroMem per run. Min of 5, 2000 x 3000 float32, --opt:speed -march=x86-64-v3: -d:release saNearest 3.76 -> 2.09 ms, saLinear(frac 0) 5.66 -> 2.37, saLinear 14.35 -> 10.89, saCubic 26.55 -> 23.16 -d:danger saNearest 2.13 -> 2.19, saLinear(0) 2.56 -> 2.53, saLinear 11.01 -> 11.06, saCubic 24.60 -> 24.35 The reviewer's 1.9-2.3x holds for checked builds; a danger build was already at memcpy speed. Bit-identical (7200 random plans). API 7 fixed, e310311. A row of another sample interval is left out and listed in the new `plan.otherInterval`; no ValueError. The bad_arguments assertion of the ValueError is replaced by a_row_of_another_interval_is_reported. API 8 fixed, bf7ea40. FoldMap: `=copy` and `=dup` are compile errors, new `clone`. Tests tfold clone_is_independent, tests/reject/rfoldcopy.nim. No other type in my files holds a tensor. API 9 fixed, 6444094. RecordCompleteness: `missing` is now the count, `missingRanges: seq[StationRange]` is the former `missing`, new `missingStations: seq[StationKey]` (X stations in X order; also with channel matching); the func missingStations(c) is removed. API 10 fixed, 0a7bb30. One rule: a receiver is stationKey(rcvLine, rcvPoint, rcvIndex) (1/100 rounding, with the index). planSlices groups by it, then channel and component; both checkGeometryInto forms leave non-live rows unchecked. Stated in docs/guide/survey-qc.md ("One rule for a receiver") and vibroseis.md. Tests tslice a_receiver_is_its_station_with_the_index, tgeometry both_forms_skip_rows_that_are_not_live. The tgeometry assertion that the view form flags a dead row encoded the disagreement and is replaced. API 11/13 fixed, 796ca94. patch: one clamp, stationIndex (private), where patchShots narrows SpsPoint.index (the only such place in my files); record: `$`(StationKey) is private now (stationText). Test tpatchgen point_index_out_of_range. pC C8 fixed, 2e50e7c. sliceKeys keeps the receiver's intervalUs bit for bit when its keys hold one (1e6/31 us does not survive / 1000 * 1000). Test tslice slice_keys_keep_the_source_interval. P3 consumer side done, 5127cd0. One private iterator per module, with the yield type of parseSEG's `relations` itself: src/processSEG/qc/fold.nim channelRelations(g, record, domain) src/processSEG/survey/patch.nim channelRelations(g, record) src/processSEG/qc/record.nim channelRelations(g, record) Replace each body by a loop over parseSEG's runs. No behaviour change: Morgantown fold, patch description and completeness hash the same; release times within noise (preplotFold 1105 -> 1065 ms, describePatches 2554 -> 2395 ms). P4, P7: not yet committed (next). Deferred as briefed: the multi-file slicing path, surveyGrid taking a ReceiverLayout, the duplicated station-index lookup rule (patch.lookup, fold.find) which needs a parseSEG export. PUBLIC API CHANGES (B = breaking) continuous/slice + plan.tracesNeeded + plan.otherInterval + SliceShot.srcIndex; initSliceShot(..., srcIndex = 0'i16) B planSlices no longer raises for a row of another interval (reported) B default receivers: station with the point index and 1/100 rounding (was exact line, point, channel, component without the index) ~ planSlices raises ValueError for a row of more than 2^48 samples ~ sliceKeys: srcIndex from the shot; intervalUs as the keys hold it ~ rowsNeeded documented as rows of the planning view (read files with tracesNeeded) qc/fold + maxBins parameter on surveyGrid, initFoldMap, preplotFold B more than 2^26 bins raises ValueError unless maxBins is raised B FoldMap cannot be copied (`=copy`, `=dup` errors); + clone(m) qc/record B RecordCompleteness.missing: seq[StationRange] -> int + RecordCompleteness.missingRanges (the former missing) + RecordCompleteness.missingStations: seq[StationKey] B func missingStations(c): int removed (use c.missing) B `$`(StationKey) no longer exported ~ rmAuto FFID fallback; (0, 0) stations absent; failed auxiliary rows; NaN channels (see C7, C14) qc/geometry B checkGeometryInto(rows: TraceRows, ...) leaves dead and auxiliary rows unchecked (as the ChunkMeta form) ~ an infinite coordinate is flagged gfOffsetMismatch survey/patch ~ generateRelations: no records for a shot without a finite source line and point; patchShots clamps the point index FOR THE LEAD WHEN MERGING - parseSEG has no `$` for StationKey: processSEG's is private now (qc/record.nim stationText: "line/point", ".index" above 1). Move it to parseSEG if wanted. - parseSEG's index clamp (`idx` in survey/join.nim) is private; patch.nim's stationIndex mirrors it. Export one and drop the other. - When parseSEG exports the run-level iterator: replace the bodies of the three channelRelations helpers above. - docs/api.md, docs/guide/survey-qc.md and the slicing section of docs/guide/vibroseis.md are updated. Not mine and untouched: docs/integration.md and docs/guide/qc.md mention checkGeometryInto for a meta only, which is still right. - New test file: tests/reject/rfoldcopy.nim. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Around its two FFTs, correlateInto spent half as much time again (review
xG P4), and its one-record form ran the eight-lane path for one trace.
- The multiplication of the spectra by the pilot's was vectorized across
bins, with the eight lanes shuffled in and out of each register. The
per-bin kernel is no longer inlined, so its loop is one vector operation
(measured alone, 250 batches of 8101 bins: 17.3 -> 4.1 ms).
- The lane buffer is gathered sample by sample (a read stream per record,
one pass over the buffer) and zeroed only past the longest record, where
the whole buffer was zeroed and then written lane by lane; the output is
scattered the same way, with the lags' three spans found once instead of
a test per sample (2000 records: gather 9.8 -> 3.3 ms, scatter 6.9 ->
3.1 ms).
- correlateInto(x, y, ws) uses the one-trace FFT.
Measured, -d:danger --opt:speed -march=x86-64-v3, one thread:
2000 records of 10001, 12 s sweep (nfft 16200), minima of 5, three runs:
77.8-80.7 -> 64.4-65.9 ms (the FFT pair alone: 53 ms)
bench/bqc (new), medians: 200 records 7.90 -> 6.18 ms; one record
254 -> 86 us.
Results: equal to the last bit or one from it (the fused products round
in another order): of 105 064 output values in four cases (two modes,
pilot start times before, at and after the time break) 3727 differ, by at
most 1.9e-8 at a peak of 0.2. The one-record form now agrees with its row
of the block form to the same rounding instead of exactly. New test
ragged_batches_and_the_one_record_form (unequal live lengths within a
batch, a batch that does not fill the lanes, dead rows, both forms against
the direct sum; no allocation); tcorrelate passes in debug, release,
danger and under ASan.
bench/bqc.nim is new (QC and vibroseis, picked up by `nim bench`);
bench/README.md gets its section, with the resampling numbers of 2b46151.
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
`planSlices` rebuilt the whole stream layout on every call: it sorted the
rows into receivers, placed each segment on its receiver's sample grid and
found the overlaps, and only then looked at the shots. A viewer that plans
one shot at a time over a million segment rows paid 128 to 408 ms per
shot for 0.05 ms of work.
The first half is now `initSliceStreams(meta, [times, [receivers]])`,
returning a `SliceStreams`: a value with no pointer into the meta, kept
for the dataset. `planSlices(streams, shots, window, order, align)` lays
out the shots from it. The forms with today's signature build the streams
and plan in one call.
Per-plan cost no longer depends on the rows either: the rows a plan
needs are a bit per row (250 KB of zeroes at a million rows, with the
count per word; the flag and the index per row were 9 MB), from which
`rowsNeeded` comes out
ascending without a sort and a row's place among them without a search.
The receiver ids from the keys are made once per row instead of inside
the sort's comparisons.
1e6 segment rows (1000 receivers x 1000 segments), 6 s windows,
-d:danger --opt:speed -march=x86-64-v3, minimum of 5 runs:
before after
planSlices(meta, 1 shot), from keys 407.5 ms 248.9 ms
planSlices(meta, 1 shot), caller ids 128.0 ms 128.2 ms
initSliceStreams, from keys / ids - 250.0 / 127.3 ms
planSlices(streams, 1 shot) - 0.063 ms
planSlices(streams, 1001 shots) - 46.9 ms
per extra shot 0.05 ms 0.047 ms
Results are unchanged: 7200 random plans (three forms, both orders, three
alignments, gaps, overlaps, failed, dead and untimed rows) and their
extracted samples and reports are identical to the old code, from the
one-call forms and from streams reused across plans; a compact read of
`rowsNeeded` extracts what the full block does.
A plan copies the streams' receiver keys and their lists of data problems
(`overlaps`, `timingIssues`, `untimed`, `otherInterval`), so `plan.overlaps`
still means every overlap of the dataset: free for clean data, a copy per
plan for data with many overlaps.
Tests: streams_made_once_serve_every_plan, streams_arguments; tslice in
debug, release, danger and under AddressSanitizer.
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…h a workspace traceAmplitudeInto cost 6.4 ns a sample (review xG P2): one scalar loop for the sums, the extremes and a finite test per sample, a second for the clip runs, a third for a sliding sum of squares that was rebuilt whenever a large term left it. Per finite trace now: - sum and sum of squares in eight partial sums; - the largest positive and the largest negative magnitude, the finite test and the constant test as three integer reductions on the samples' bits (float min/max reductions do not vectorize without reordering); - a vector count of the samples at the clip levels: only traces with clipRun of them get the run by run pass; - the spike windows from per-block suffix and exclusive prefix sums of the squared deviations (no subtraction, as vstack's window powers), so that the loop over the windows inside the trace is straight-line code. A trace with a NaN or an Inf takes the former path, sample by sample. BREAKING: the block sums need 2 * cols + 1 float64, so traceAmplitudeInto(b | samples, meta, params, dst, ws) takes a Workspace holding traceAmplitudeWorkspace(cols) (new; ValueError when it holds less, as recordAmplitude raises for its own). The pooled form keeps its signature: each worker reserves its workspace. The one-trace traceAmplitude allocates its scratch. Callers updated: tests, the survey QC guide, docs/api.md. Measured, -d:danger --opt:speed -march=x86-64-v3, one thread, bench/bqc medians, 2000 x 3000 float32 (the reviewer's panel): 38.6 -> 14.2 ms; the reviewer's own kernel measured 22.8 ms. The Morgantown record through the pool of 4 in a debug build: 790 -> 197 ms. Results: counts and flags equal on 8800 traces compared with the previous kernel (spikes, clips, DC, quantized, short, constant, dead, with windows of 1 to 100 samples); rmsDb and maxDb to the printed 1e-9 dB; the mean differs in the 13th digit on three traces (the sums' order). New test the_kernels_against_the_definitions: 360 traces of every length from 0 to 90 and five window lengths against a sample by sample reference, float32, float64 and the one-trace form. tamplitude passes in debug, release, danger and under ASan. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
coherenceSnrInto summed each column's window in a loop, cols x window operations: 21 ms for a 12 s window on a 24 x 6000 gather (review xG P5). The four window sums now cost O(cols). The residual and the coherent energy come from block sums (per block of one window, a suffix and an exclusive prefix; a window is one or two of them added), not from differences of running sums: a quiet window after a loud stretch keeps its precision, and a NaN spoils the windows that hold it and no others, as before. The two counts are whole numbers and are exact as differences of running sums. Measured, -d:danger --opt:speed -march=x86-64-v3, one thread, bench/bqc medians, 24 x 6000 float32: window 0.1 s 0.375 -> 0.153 ms, 1 s 2.42 -> 0.162 ms, 12 s 21.0 -> 0.185 ms. Results within rounding: the sums are formed in another order. New test coherence_window_sums_against_the_definition compares with a loop over every window, for windows of 1 sample to more than the block, with mutes, a 10^4 : 1 loud and quiet stretch and a NaN: within 1e-9 dB, the same NaN columns. The existing exact-equality tests pass unchanged. coherenceSnrWorkspace(cols) grows from 3 * cols to 6 * cols + 2 float64; a workspace that holds less now raises ValueError (it was the scratch bound check's IndexDefect), as the amplitude QC does. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…field fill A VelocityField column cost more than the NMO it served, and bruteStack computed one per trace on shot-sorted input: - VelocityField.vrmsInto blended the nodes twice for every output sample. The output times increase, so the grid cell only moves forward: the nodes are now blended once per grid time, and the grid position takes a multiplication by 1/dt (the point query vrmsAt takes it the same way, so column and point query agree to an ulp: 4e-16 over the test's cases). 1501 samples of a 2D field between two nodes: 9.1 -> 6.8 us (danger), 12.3 -> 7.5 us (release). - bruteStack kept one column per worker (the last bin's), so shot-sorted input read the field about once per trace. Each worker now keeps the sloth columns of the bins it met last, by bin number, up to BruteStackOptions.velocityCacheBytes (new, 8 MiB: 698 columns of 1501 samples; bruteStackMemory counts it). A column depends on its bin only, so the cache never shows in the result (tested with caches of 1, 3 and 7 columns, and the default). 60 000 shot-sorted traces x 1501, a 2D field, 1 / 4 / 8 workers: 825 / - / - ms after the NMO change (1058 / 287 / 174 at the start of the branch) -> 241 / 101 / 87 ms, the plain stack's time; a 3D field on 40 000 bins: 782 / 563 -> 284 / 130 ms (1 / 4 workers). - initVelocityField evaluated every function at every grid time through vrmsAt (a knot search, a square root, then its square). A node is now filled by one forward sweep over its knots, storing W/t as it comes: 200 x 250 nodes of 12 knots x 1501 times: 740 -> 217 ms (danger). V2 of a node may differ in its last bit from before (no sqrt-then-square). vrmsAt and the sweep share the two formulas (knotsV2, gradientV2). Tests: tvelocity field_columns_equal_point_queries (2D and 3D, on and off the nodes, intervals finer and coarser than the grid, columns past tMax) and field_nodes_equal_their_functions (knots on and off the grid, a knot at zero, one knot, gradients, a constant); tbatch velocity_cache_never_shows and memory_bound. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
fitShot allocated about once per row with the two-segment model and twice with the binned one (review xG P5): 1793 and 4111 allocations for a second call on 2000 rows. The median of each trace's neighbours sorted a seq (the sort's buffer), the binned model copied its bins for every row (`let bs = fit.bins`) and built its lists anew, as did the symmetry test. - The medians select in place (a quickselect; the same value as from a sort: the middle one, or the mean of the two middle ones). - The index, key and side lists are scratch of the ShotFit, kept from call to call; the bins and the spread order are read in place. A second call on 2000 rows now allocates 2 times (two-segment) or 3 (binned): the buffers of its sorts. Outputs equal to the last bit on five shots of 7 to 2000 rows, both models, with and without samples (9488 lines of %.17g compared). Time: bench/bqc 22.6 -> 22.3 ms and 22.8 -> 22.4 ms; the L1 fit over the crossover candidates is the cost, not the allocator. New test a_fit_allocates_a_few_times_whatever_the_rows. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The sentence on the lag's sign had been appended to a line already at the limit. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
A stack by header numbers (cdpBinning, inlineCrosslineBinning) had no coordinates: stackKeys wrote NaN as cmpX/cmpY, so a viewer could not place it on its map; and a fold map meant scanning the bins x samples fold tensor. StackAccumulator now counts, per bin, the traces added (tracesAt: the adds that added a sample) and keeps the mean of the midpoints given for them (addMidpoint; cmpXAt, cmpYAt, NaN when none): 24 bytes per bin, which stackMemory counts (bins x (samples x (sizeof(T) + 4) + 24)). merge, clear and copies carry them. bruteStack gives each stacked trace's midpoint, and StackResult has binTraces (traces per bin), cmpX and cmpY (their mean midpoint). stackKeys, and so writeStack, give the bin centre for line and grid binnings as before, and the mean midpoint for bins of header numbers. All of it is the same for any number of workers. Tests: tstack traces_and_midpoints_per_bin, validation_memory_allocation (the formula); tbatch fold_map_and_midpoints (counts and means against the keys, an empty bin, the written file read back, a line binning keeps its centres). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
lineBinning numbered its bins firstCdp + i and gridBinning its inlines and crosslines in ones, while cdpBinning and inlineCrosslineBinning take steps: a line whose CDPs go in twos, or a grid numbered every 2 inlines and 5 crosslines, could not be described, and a VelocityField picked at such numbers was read at the wrong places. Both take steps now, as last parameters with default 1 (cdpStep; inlineStep, crosslineStep), finite and positive. Additive: no existing call changes its meaning. The parameter order of gridBinning stays as it is (spacings crossline first, counts and first numbers inline first). Reordering would silently swap the meaning of positional calls, and an options object would only add a second way to build the same grid next to the one whose order stays inconsistent; Nim's named arguments already give the order-free form, so the doc comment and api.md now say to pass them by name. tbinning numbering_steps: the numbers with and without steps, the bin of a midpoint unchanged, and the refused steps. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
For every shot `generateRelations` looked at every station of every template line (to keep those in the patch's slots), and with distance limits at every station of the survey. - Template: the slots of a line ascend with its points, so the stations in slots s0 .. s1 are a run found by bisection. Whether a line's slots ascend is checked on its first use (a layout built by hand may have them in any order: such a line is still scanned whole). Short lines are listed once instead of being looked for among all lines per shot. - Distance limits: the limits lie inside a rectangle along the receiver lines (the box turned into the layout frame, cut by the radius and the maximum offset). The lines that can reach it are a range, found from the running extremes of the lines' crossline positions; on a line whose inline positions run one way, so are its stations. Every candidate still passes the exact test, so the stations picked are the same. - A radius or maximum offset of Inf no longer costs a `hypot` per station, and the picked stations are sorted only when they are not in order already. -d:danger --opt:speed -march=x86-64-v3, bench_survey, minimum of 5: receivers (lines x stations), shots, patch template distance box 100 x 1000, 2000 shots, 10 x 100 74.9 -> 31.6 226.4 -> 30.1 ms 20 x 200, 2000 shots, 10 x 100 30.9 -> 25.3 53.3 -> 20.1 ms 300 x 2000, 500 shots, 20 x 200 76.5 -> 32.3 802.3 -> 49.0 ms (per shot at 1e5 receivers: 37 -> 16 µs and 113 -> 15 µs; what is left is mostly building the X records.) Identical output: 6215 synthetic cases (three survey sizes, five azimuths, three numberings, jitter, holes, short lines, 46 templates) and five Morgantown runs give the same records as the code before; so do 840 cases on layouts altered by hand (reversed and shuffled stations, odd point steps, NaN positions, reversed lines, all lines short). Tests: tpatchgen stations_out_of_order; tpatchdist irregular_layout_against_brute_force (jittered stations, a line that doubles back, rotated boxes, a strip, radius, maximum offset, against a test of every station) and the_order_of_a_layout_does_not_matter. All three pass on the old code too, and fail when the scan fallbacks are taken out. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
bruteStack reads a VelocityField at each bin's numbers (binning.location(bin)), while applyNmo and moveoutTimesInto read it at the traces' header numbers (ensemble, or inline and crossline). The two agree for cdpBinning and inlineCrosslineBinning only: with lineBinning or gridBinning a viewer could not preview the NMO, or draw the hyperbolas, with the velocities the stack would use. New, all additive: - Binning.locationOf(keys): the location of a trace's bin (all NaN when it is in none); never raises. - moveoutTimesInto(rows | meta, field, binning, t0Ms, times) and moveoutColumnsInto(meta, field, binning, t0Ms, cols): the hyperbolas with the field's velocity at the row's bin. - applyNmo(b | samples, meta, field, binning, p, ws): NMO with the field read at the row's bin; rows in no bin are left alone. The header form and this one share one loop. - checkSuits(field, binning): the 2D/3D check bruteStack made privately (same messages), now also made by these forms. Overloads taking the Binning fit the existing set better than a locationsInto: the locations forms are 2D only, and applyNmo has none. Tests: tmoveout field_times_by_binning (a line and a grid, a row in no bin, against the header form), tnmo field_by_binning (the same samples as the header form on keys that hold the bin's CDP). The guide's overlay example and integration.md "Overlays" show the binning forms. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…file The review found no test of bruteStack running while another thread uses the file. tbatch other_threads_read_the_file runs the job on a thread of its own while the main thread keeps reading the same SeismicFile (readTracesInto, rowsOf, traceKeys, binOf), then compares the result bit for bit with a stack made on the main thread. It is the part of the documented thread contract that holds with any parseSEG; creating prefetchers meanwhile needs the parseSEG release that locks the map's pin count, and is not exercised here. Clean under ThreadSanitizer (gcc, 0 reports for the whole of tbatch) and AddressSanitizer. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Final report of the qc engineer (review xG), branch fix/qc on bd27552: 24 commits of work, head 52036f3, plus the interim report ba8999b, which this replaces (the correctness and API items are as there; the performance items P2, P4, P5 are new). FINDINGS C1 fixed 2cf5eda. pmCoppens' rise search starts at the first whole numerator window; doc "one window". Tests (tfirstbreak): coppens_does_not_pick_its_own_warm_up (200 traces, noise 0.1 of the onset: no col < w, >= 95 % within 3 samples), coppens_on_the_example_record (0 picks in the first window; was 3210 of 3954 in the first 12 samples). The existing burst test moved its burst from 4 ms to 24 ms, reason in the commit. C2 fixed d4b0c7d. (a) a noise window without energy gives snrDb NaN, confidence 0; (b) a row is picked over its liveSpan (of the unfiltered row), so the zero-to-data edge is outside the search. Test zeros_before_the_data_are_not_noise: edge picks 55/100 (STA/LTA) and 87/100 (energy ratio) at confidence 1.00 -> 0/100; the onset 100/100 with all three methods, with and without band-pass. CONSEQUENCE: a noise-free synthetic, whose first break is its first non-zero sample, is no longer picked on it (median error 0.5 -> 5.5 samples for the energy ratio, 16.8 for STA/LTA, none for Coppens; confidence 0). Either form of the requested fix does this (the reviewer's `lo = first + w` too). Documented in the module and the guide; every test and guide example has noise. C3 not mine (closed by the lead on the base commit). C4 fixed b426927. sweepFromHeaders validates at the trace's interval (maxNyquistFraction 1) and, with the new `samples` parameter, refuses a sweep longer than the trace (not "of another length": an aux pilot trace is often as long as the record). pilotFromTrace also drops a description whose upper half-band holds < 5 % of the trace's energy: without it the 0.1 Hz header mistake passes at <= 0.5 ms (tested with a file written at 0.25 ms). resamplePilot judges the samples (energyAbove, new), keeps a description only if it still validates, checks pilots without a description too; `bandLimit = true` (new) cuts a pilot on purpose. sweepOverrides RAISES ValueError for a fractional Hz or ms instead of rounding: the fields are read back as Pilot.sweep, and 7.5 Hz written as 8 is a wrong sweep nothing can detect. Tests (tpilot): headers_that_do_not_fit_the_trace, the_example_pilot (pilotFromTrace(Sweep0.sgy).resamplePilot(2).len == 6001; skipped when absent), sweep_overrides, resample_down_and_up. C5 fixed 0e37697. Passband edge = cutoffFraction - 2.736 / halfWidth of the lower Nyquist; energy above it is refused (or cut with bandLimit); the kept description gets the edge as maxNyquistFraction; taps beyond the pilot's ends are made up for by rescaling (0.61 -> 1.0000). DEFAULT halfWidth 24 -> 32, so that the default passband (0.8145 Nyquist) holds a default sweep (0.8); at 24 it would be refused. Test resampling_keeps_its_passband_flat (within 0.001 dB at 150/190/199 Hz; 10-210 and 10-220 Hz raise). C6 fixed 437c939. Both slopes go to butterBandPass; test in instrument_filter_causality (3 Hz at 6 dB/oct with 200 Hz at 72). C7 fixed 1352427. sweepSamples validates first; sweepLen is 0 for a spec with an invalid interval or length. Test specs_built_as_objects. C8 fixed 8a3b05a. addRecord (meta form) raises for another dtMs (relative 1e-9) or firstSample != 0. Test in bad_arguments (tvstack). C9 fixed a1d14ee. Each polarity has its own clip level when both reach the maximum magnitude within 1 % (int24 sine: 500 -> 1000 clipped samples counted); afClipped needs two runs or one of 2 * clipRun (`clipped` still counts every run). Tests clips_of_both_polarities, one_short_run_at_the_maximum_is_not_a_clip, the_example_record_is_not_clipped (4 flagged -> 0). Two existing tests that encoded "one run of clipRun is a clip" were changed. C10 documented dfcb935; measured on the current code, guided window, 400 traces, 30 Hz onset, 2 ms; median pick - onset in samples at noise 2.6/6.6/13/26 % of the peak: pmStaLta 1.9 1.9 1.9 1.8, with 8-60 Hz 7.7 8.7 9.5 10.3; pmEnergyRatio 0.7 0.9 1.2 1.6, with filter 2.9 3.6 4.1 4.8; pmCoppens 0.6 1.0 1.4 1.7, with filter 2.8 3.6 9.7 10.7. In the module doc and docs/guide/qc.md. C11 fixed 247b9b7. validate bounds the total boost (300 dB dB/octave, 150 dB dB/Hz; the dB/Hz limit was 600); frequencyAt returns f1Hz/f2Hz exactly at the ends and stays in band; the dB/octave law is computed from f2/f1 (rounding-level change of its samples). Test steep_boosts_stay_finite. C12 fixed dfcb935 (correlate.nim "6001: 12 s listen"; the guide's column; the pmCoppens description). P1 fixed 417b0bb. Windowed STA/LTA and energy-ratio picks read, copy and band-pass only their stretch. -d:danger --opt:speed -march=x86-64-v3, one thread, minima of 9, the reviewer's fb6 case (2000 x 3000, 150 ms window): energy ratio 8.2 -> 1.0 ms, with 8-60 Hz 40.7 -> 4.0 ms; STA/LTA 30.1 -> 1.3 and 62.7 -> 4.3 ms; Coppens (whole trace by its definition) 15.4 -> 8.1 and 47.9 -> 40.7 ms. bench/bprestack medians: 20.1/5.7/10.3 -> 1.36/1.04/5.55 ms, with 5-60 Hz 42.9/28.4/33.1 -> 5.84/5.21/28.1 ms. PICKS: unwindowed plans and Coppens bit-identical. Windowed, no band-pass: the same column on all 33 600 compared picks, SNR within 0.002 dB. Windowed WITH band-pass: NOT identical. On fb6's case 5 of 12 000 picks move (four by 1 sample, one by 4); the mean pick error is equal to 0.01 ms. Cause: the filter starts at the stretch and the stretch is not detrended (a line fitted to a short stretch is set by the break). Where the trace start cuts the stretch inside the filter's margin (3 Hz low cut: 261 samples) up to 7 % of picks move, mostly by one sample. Exact identity needs the whole trace filtered. Test a_windowed_pick_reads_only_its_stretch (bit for bit). P2 fixed 263856a. bench/bqc medians, 2000 x 3000 float32 (the reviewer's panel): 38.6 -> 14.2 ms (the reviewer's kernel: 22.8 ms). Counts and flags equal to the previous kernel on 8800 traces; rmsDb/maxDb to 1e-9 dB; the mean differs in the 13th digit on three (sum order). Test the_kernels_against_the_definitions. BREAKING signature (below). P3 DEFERRED (diversity stack, 10.1 ns/sample). Not attempted. Ideas of the review: split windowPowers' second loop into head / interior / tail so that the interior vectorizes (the exclusive-prefix block sums of amplitude.spikeCount and snr.blockSums do exactly that and could be shared); take the median of every h-th window power. P4 fixed 51c4f05. 2000 records of 10001 samples, 12 s sweep (nfft 16200), minima of 5, three runs: 77.8-80.7 -> 64.4-65.9 ms (the FFT pair alone: 53 ms). bench/bqc: 200 records 7.90 -> 6.18 ms; one record 254 -> 86 us. Beyond the brief's gather: the spectral multiply was vectorized across bins (17.3 -> 4.1 ms alone) and the scatter is sample-major. Outputs equal to the last bit or one from it (3727 of 105 064 values differ, by <= 1.9e-8 at a peak of 0.2); the one-record form agrees with its block row to the same rounding, no longer exactly. Test ragged_batches_and_the_one_record_form. P5 fixed in three commits. coherenceSnrInto 6209fcf: O(cols); 24 x 6000, window 0.1/1/12 s: 0.375/2.42/21.0 -> 0.153/0.162/0.185 ms. TOLERANCE: within 1e-9 dB of a loop over every window (test coherence_window_sums_against_the_definition, with a 1e4 : 1 loud/quiet stretch and a NaN, which stays in its own windows). The sums are block sums, not differences of running sums, to keep that. fitShot 39f3003: a second call on 2000 rows allocates 2 (two-segment) or 3 (binned) times, was 1793 and 4111; outputs bit-identical; time 22.6 -> 22.3 ms (the L1 fit is the cost). Test a_fit_allocates_a_few_times_whatever_the_rows. resamplePilot 2b46151: kernels tabulated per phase; 30 s at 0.25 ms to 4 ms 358 -> 5.5 ms (same halfWidth 32); within 1e-12 of a tap-by-tap reference (test resampling_by_tabulated_kernels). A2 the two clear cases fixed 5c9c39e; the rest NOT applied, as told (recommendation below). A3 fixed a476d41 (breaking). A4 fixed fd8012d (breaking); the lag's sign tested (the_sign_of_a_lag). A5 fixed 07a13be. A6 fixed 5a45efc. FirstBreakPlan.requiredMargin: a block holding it on each side of the search window picks bit-identically to whole traces (test a_zoom_window_with_the_required_margin_picks_the_same); plan.samples for whole-trace plans (no window, or Coppens). ProcessPool.pickInto(plan, block | tensor, meta, picks): test the_pool_picks_the_same (1, 3, 8 workers; no allocation once warm; ASan clean). New exported checkPick. A7 fixed 223f9c5. A8 fixed 4bddbd0 (shows in release and danger builds). PUBLIC API CHANGES (B = breaking, + = new) qc/firstbreak B FirstBreakPick.col: float64 (NaN: none) -> int (-1: none); the fields default to "no pick", so FirstBreakPick() == noPick(). + requiredMargin(plan); checkPick(plan, b, meta, picks). Behaviour: liveSpan rule (C2); windowed stretch (P1); the interval is compared within a relative 1e-9 (A5). core/pool + pickInto(pool, plan, b | samples, meta, dst). traceAmplitudeInto(pool, ...): same signature; each worker reserves traceAmplitudeWorkspace. qc/amplitude B traceAmplitudeInto(b | samples, meta, params, dst, ws): new last parameter `ws: var Workspace`; ValueError when it holds less than + traceAmplitudeWorkspace(cols) (2 * cols + 1 float64). B AmplitudeParams.spikeWindowMs -> spikeHalfWindowMs. B recordAmplitude(traces, params, ws): `params` lost its (unusable) default; no caller changes. traceAmplitude (one trace) allocates its scratch. afClipped needs two runs or a long one (C9). qc/snr B coherenceSnrWorkspace(cols) is 6 * cols + 2 (was 3 * cols); coherenceSnrInto raises ValueError for a smaller workspace. qc/shotfit B ShotFitParams.residualMs -> maxResidualMs. vibroseis/pilot B InstrumentPhase, ipMinimumPhase, ipZeroPhase removed: InstrumentFilter.phase is filter/iir's IirPhase (ipMinimum, ipZero), default ipMinimum. B sweepFromHeaders: a proc, with a new last parameter `samples = 0`; none for a description that does not validate. B sweepOverrides raises ValueError for fractional Hz / ms. B resamplePilot(p, dtMs, cutoffFraction = 0.9, halfWidth = 32, bandLimit = false): default halfWidth was 24; refuses by the samples and at the passband edge; halfWidth must leave a passband. + energyAbove(pilot, fHz). applyInstrumentFilter takes unequal slopes. vibroseis/sweep B validate: boost bounds (300 dB total for dB/octave, 150 dB for dB/Hz; was 600 for dB/Hz). sweepLen returns 0 on an invalid spec. vibroseis/taper SweepTaper.param defaults to NaN (the shape's default). vibroseis/vstack B recordStackMemory(rows, samples, mode = rsVertical): no `T`. B addRecord (meta form) raises for another interval or a window. vibroseis/correlate No signature change; rounding-level change of the outputs (P4). adapters pilotFromTrace: `sweep` only when the trace can be it (C4). A2, RECOMMENDATION NOT APPLIED. The rule read from the Into procs outside my area (extractInto, amplitudeInto, averageAmplitudeInto, FFT forward / inverse, PipelinePlan.apply): a compiled plan first; then the input and its meta; then parameters; then the output; `ws` the last required parameter; defaults after it; constructors (T, spec, dtMs, samples, ...). Still against it: plan.pickInto(b, meta, ws, dst) -> (b, meta, dst, ws) coherenceSnrInto(b, meta, windowMs, ws, dst, floorDb = ...) -> (..., dst, ws, floorDb) initCorrelationPlan(T, pilot, recordSamples, dtMs, ...) -> (T, pilot, dtMs, recordSamples) pool.pickInto follows plan.pickInto without `ws`, so it would not change. Conforming: windowSnrInto, traceAmplitudeInto (dst, then ws), fitShot, correlateInto (it mirrors PipelinePlan.apply: src, dst, meta, ws), initFirstBreakPlan, initRecordStack (as initStackAccumulator). FOR THE MERGE - Callers of the renamed names outside my files: none left in this tree (grep clean for spikeWindowMs, ShotFitParams.residualMs, InstrumentPhase, recordStackMemory(T, ...), the four-argument traceAmplitudeInto). I adapted tests/tproperty.nim (int column) and tests/tcorrelate.nim (ipMinimum), which are not on my list. - tests/tdocs.nim resolves `export a.B` to a sibling file of the exporting module, so pilot.nim does not re-export iir.IirPhase; the umbrella exports both. - A5: abs(a - b) > 1e-9 * b is written out at four sites of mine (checkPick, checkCorrelate, addRecord, RecordStack.merge), as in filter/apply, filter/iir and gain. A `sameInterval` in core/types would replace all of them. - Block sums without subtraction exist three times now (vstack.windowPowers, amplitude.spikeCount, snr.blockSums); core/sums could hold one. - The final-pass changes that meet my code: the IIR's 24 lanes (pickInto batches fftLanes stretches per applyFilter call through a lanes meta with fftLanes keys: correct with any lane count, but the batch should follow the IIR's), and zero padding for low-cut / band-pass filters (it changes the start-up of a stretch cut by the trace start, so the P1 differences near the trace start will change; the zoom-window and the pooled identity tests do not depend on it). - bench/bqc.nim is new and runs under `nim bench`; bench/README.md has two new sections ("First-break picking in a search window", "QC and vibroseis"). - Not done, worth doing: tests/tvectorized.nim (not mine) does not know the new vector kernels (amplitude.sumAndSquares / extremes / countBeyond / spikeCount's interior loop, correlate.mulKernel, pilot.dot). mulKernel in particular is fast only while it is not inlined; a check there would keep it so. SUITE (at 52036f3, clean tree) - debug: All 50 test files passed (48 tests, 2 reject tests). - release: All 50 test files passed. - danger: All 50 test files passed. (tests/tester.nim run three times side by side, the release and danger runs with their own --outdir, Nim cache and TMPDIR; tvectorized runs in the debug suite only, by its own design. The example files were present: no sample-file test was skipped.) - Under -d:addressSanitizer: tfirstbreak (pooled pickInto), tcorrelate, tamplitude, tsnr: clean. - nimpretty --indent:2 --maxLineLen:100 is stable on every .nim file changed since bd27552; nim doc builds the eleven changed modules. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The lock names parseSEG's 0.10.0 commit, and its copy of the package file is the current one (it still said 0.1.0). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The final pass over processSEG before viewer development starts: a review of every module with reproducers, then fixes with a regression test each. The release notes below are the summary. The lock moves to parseSEG 0.10.0 (xrfez/parseSEG#3).
processSEG 0.3.0 is a correctness and speed pass over the whole library before viewer development starts: the FFT and filter kernels, the pool, gains, QC, vibroseis, velocity, NMO and stacking, and survey QC. Every fix has a regression test. It needs parseSEG 0.10.0.
Install
Or in your
.nimblefile:requires "https://github.com/xrfez/processSEG >= 0.3.0".Faster
Measured on a Ryzen 9 3950X, gcc 12, x86-64-v3, single thread unless said.
traceAmplitudeIntoProcessPoolhands rows out in runs from a shared counter instead of one fixed share per worker: 4 workers on a panel with half its traces dead go from 2.06x to 3.51x.powper sample off the table: givecompilethe record length.shiftSamples, velocity field columns blended once with a column cache, a coherence SNR whose time no longer grows with its window, slice extraction that copies whole runs.Fixed
bruteStackstacks on the calling thread, in file order, into one stack: the result is bit-identical for any worker count and takes one stack of memory. Size overflow is aValueError; a velocity field must cover the stack.maxBinsbins raisesValueError.New
compile(p, T, dtMs, samples, recordSamples): the record length a zoom window is cut from.applyGainCurve(..., ws)andworkspaceSize(GainCurvePlan);reserve(ws, plan)for any plan with aworkspaceSize.requiredMargin(FirstBreakPlan)for zoom windows, andpickIntoon the pool.bruteStacktakes a list of traces; a stack has a trace count and a mean midpoint per bin; numbering steps for line and grid bins.SliceStreams, made once per dataset; a slice plan lists the traces to read.FoldMap.clone.Breaking changes
Opparameters are read and set through accessors that check the kind (ValueErrorfor the wrong one); the variant fields are private.FilterPlanpads less under a low-cut.BruteStackOptions.staticsis anOption[DatumStatics].FirstBreakPick.colis anint, -1 for no pick.traceAmplitudeIntotakes a workspace.coherenceSnrWorkspaceis 6 x cols + 2.spikeHalfWindowMs,maxResidualMs;missingis a count everywhere.InstrumentPhaseis removed:InstrumentFilter.phaseis anIirPhase.sweepFromHeaderstakes the trace's sample count;sweepOverridesraises for fractional values;resamplePilotdefaults tohalfWidth = 32;recordStackMemory(rows, samples, mode).FoldMapis moved or cloned, never copied.Verification
nim test,nim testRelease(release and danger),nim asanandnim docspass; the formatting check is clean.tests/tvectorized.nimchecks from the disassembly that the FFT, filter, IIR and NMO kernels are packed.🤖 Generated with Claude Code