Repository navigation
Fancyvoigtfit - #76
Merged
Merged
Conversation
get_mean_flux_scale and its _rescale_mean_flux wrapper are gone from
py_module.cpp; fluxstatistics.mean_flux now does the same Newton-Raphson
iteration in numpy, splitting the two reductions into one chunk per core
and summing them in a ThreadPoolExecutor. The numpy ufuncs release the
GIL, so the chunks really do run in parallel. The thread pool is a module
global created on first use, and the scratch buffers are allocated once
per call, so the iteration itself does not allocate.
numpy's SIMD exp beats scalar libm, and the C loop was memory bandwidth
bound on tau (only ~3x on 8 threads), so the python is close to par.
Wall seconds on 8 cores, tol=1e-5:
n C 8thr C 1thr numpy ratio to C 8thr
1e6 0.008 0.018 0.017 2.1
5e6 0.033 0.105 0.049 1.5
2e7 0.135 0.396 0.148 1.1
8e7 0.601 1.74 0.580 0.96
Agrees with the C result to 1e-14 relative for mean fluxes from 0.05 to
0.95 and n from 5e3 to 2e6. Partial sums are combined with math.fsum so
the answer does not depend on the thread count. mean_flux takes an
optional nthreads, and stays serial below 250k elements where the
handoff costs more than it saves.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
When mean_flux_desired was not set we computed np.mean(np.exp(-tau)) over
the whole array purely to normalise the flux, and then computed
exp(-tau) a second time in the batch loop. But the k=0 mode of each
sightline is just the flux summed over pixels, so the mean flux falls out
of the FFT we already do and the extra pass can go.
The rest follows from the FFT being linear. d_F = F/mean(F) - 1: dividing
by the mean flux rescales every mode and subtracting one shifts k=0
alone, so both can be applied to the summed power at the end rather than
to the flux of every batch. That drops two full passes over each batch,
and folding vmax and the npix**2 FFT normalisation into the same final
multiply drops two more.
Timings, best of 9 alternating runs on 8 cores:
shape case old new speedup
32000x2048 mean from data 1.859 1.348 1.38
32000x2048 rescaled 1.827 1.920 0.95
8000x2048 mean from data 0.363 0.184 1.97
8000x2048 rescaled 0.378 0.359 1.06
1000x512 mean from data 0.006 0.003 1.86
1000x512 rescaled 0.012 0.010 1.17
The rescaled branch never did the redundant pass, so it gains only the
folded passes: 1.2-1.3x on the loop body at 8000x2048 and below, but 4%
slower at 32000x2048, which reproduces across processes and looks like
allocator behaviour at a 524MB working set rather than arithmetic.
Agrees with the previous output to 2.3e-15 relative (float32 input: 4e-7,
which is that dtype's epsilon) over seven shapes, both mean flux
branches and window on and off. The one exception is nspec=1, where the
k=0 bin is an exact cancellation and both versions return ~1e-29.
_powerspectrum is now unused by flux_power but is kept: the tests call it
directly.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
np.fft.rfft is single threaded; scipy.fft.rfft takes a workers argument
and threads over the transforms. scipy is already a dependency.
The worker threads spin for a while after each transform, and in this
loop they spend that time competing with the (serial) np.exp of the next
batch, which made small batches slower rather than faster: measured
0.62-0.81x for batches between 2.5e3 and 2e5 elements. So ask for
threads only once a batch is worth it, at 1<<20 elements. overwrite_x
lets the transform reuse the flux buffer, which is a temporary here.
Best of 9 alternating runs, mean flux taken from the data:
shape batch elts np.fft scipy speedup
32000x2048 6553600 1.156 0.892 1.30
8000x2048 1638400 0.167 0.121 1.38
2000x1024 204800 0.014 0.014 1.04
1000x512 51200 0.003 0.003 1.00
200x128 2560 0.000 0.000 0.94
ie, 1.3-1.4x once the batches are big enough to thread, and unchanged
below the threshold. Same for the rescaled branch: 1.35x at 32000x2048
measured with the Newton solve pinned.
Output is unchanged to 6.5e-14 relative in float64 against the reference
captured before this series of changes, and 5.0e-07 for float32 input,
which is that dtype's epsilon.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
testPowerspectrum and testFluxPower were built around _powerspectrum, the
private FFT helper, which flux_power stopped calling when the mean flux
moved to the k=0 Fourier mode. Testing the helper no longer said much
about the function that matters, and the old testFluxPower compared
flux_power against _powerspectrum of the same data divided by a magic
12.5569 (the square of the mean flux of that particular input), so it
could not survive any change to how the normalisation is applied.
The replacements go at flux_power itself, with analytic expectations
rather than a reference implementation:
testFluxPowerSine optical depths whose flux is exactly 1+a*sin put
a^2/4 in one bin and nothing anywhere else, and
vmax only sets the units
testFluxPowerParseval summed over modes the power is the variance of
the flux variation, which pins the npix**2 and
the mean flux normalisation
testFluxPowerMean the power is the mean over sightlines, checked
against single sightlines and over more
sightlines than the batch loop has batches
testFluxPowerWindow the window function is divided out
testFluxPowerRescaled the rescaled branch hits the desired mean flux
and normalises to it
All keep the old tests' even and odd npix coverage. Mutation tested:
dropping vmax, npix for npix**2, dropping the mean over spectra, mean
flux instead of its square, dropping or flipping the sign of the k=0
shift, multiplying by the window instead of dividing, perturbing the mean
flux, and truncating the last batch are all caught. The one deliberate
mutation not caught is the FFT threading threshold, which by design does
not change the output.
_powerspectrum then has no callers and is removed. flux_power output is
unchanged: 6.5e-14 relative in float64 against the reference from before
this series, 5.0e-07 for float32 input.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
_flux_power_bins built its bins with np.fft.rfftfreq while flux_power does the transform with scipy.fft.rfft. The two agree today, and the 0..n/2 ordering with a positive Nyquist is documented API in both, so this is consistency rather than a fix: if the bin convention ever moves, the bins now move with the transform that produced the power. Nothing tied the returned bins to the power ordering, so add that to testFluxPowerSine: the sine peak must land at bins[freq] == 2 pi freq in the convention flux_power returns, and bins[0] must be zero. Mutation tested with shifted and reversed bins, both caught. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
mean_flux had a pool of its own; make it a module level pool and put the
flux_power batch loop through the same one. They never overlap, since
flux_power calls mean_flux before the batch loop rather than from inside
it, and the pool grows if a later call asks for more threads than an
earlier one set it up with. Anything submitted to the pool must not
itself submit to the pool, or a fixed size pool waits on itself: noted at
_get_pool.
The batch loop body moves to _batch_power. Its exp, transform and abs
all drop the GIL, so the batches now run at the same time rather than
just the transform being threaded internally; the transform gets
workers=1 in that case so the two levels do not oversubscribe. When we
are not using the threads ourselves the transform still gets all of them.
Batches go out in waves of nthreads so the new nthreads argument really
does cap the threads used, and inputs below 1<<20 optical depths stay on
the serial path, where the handoff costs more than it saves.
Best of 9 alternating runs against the previous commit:
shape elts mean flux before after speedup
32000x2048 6.6e+07 from data 1.036 0.456 2.27
32000x2048 6.6e+07 rescaled 1.567 1.004 1.56
8000x2048 1.6e+07 from data 0.139 0.068 2.05
8000x2048 1.6e+07 rescaled 0.259 0.203 1.27
2000x1024 2.0e+06 from data 0.019 0.010 1.80
2000x1024 2.0e+06 rescaled 0.037 0.025 1.45
1000x512 5.1e+05 from data 0.006 0.007 0.91
200x128 2.6e+04 from data 0.001 0.001 0.97
The last two are below the threshold and so run the same code as before;
they are milliseconds and the ratios are noise.
flux_power is bit identical across thread counts. mean_flux is not, since
its chunk boundaries follow nthreads and so the partial sums round
differently: nthreads=3 moves the scale by 6e-16, and the rescaled power
by 1.1e-15. Unchanged against the reference from before this series:
6.5e-14 relative in float64, 5.0e-07 for float32 input.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Same treatment as flux_power: the flux is exponentiated and histogrammed
in batches, which go to the shared thread pool. np.exp releases the GIL
and so does most of np.histogram, which for uniform bins is mostly
ufuncs: measured 4.0x and 2.4x on 8 threads respectively.
Counts add, so the batches can simply be summed, and the normalisation
to a density is applied once at the end instead of by np.histogram. It
is written in numpy's order, counts/np.diff(bins)/np.sum(counts), so the
result is bit identical rather than merely close: multiplying by nbins
instead of dividing by the bin width was off by up to 1e-14.
Note this normalises by the number of samples which landed in a bin, not
by the size of the input, which is what np.histogram(density=True) does
and so what this always did. Optical depths below zero give a flux above
one, which falls outside the bins and is not counted.
flux_pdf also takes nthreads now. Inputs below 1<<20 optical depths stay
on the old path, a single batch computed in this thread.
Best of 9 alternating runs:
shape elts mean flux before after speedup
32000x2048 6.6e+07 from data 1.1187 0.4491 2.49
32000x2048 6.6e+07 rescaled 1.6438 0.9740 1.69
8000x2048 1.6e+07 from data 0.2605 0.0754 3.45
8000x2048 1.6e+07 rescaled 0.3969 0.2141 1.85
2000x1024 2.0e+06 from data 0.0329 0.0114 2.88
2000x1024 2.0e+06 rescaled 0.0461 0.0279 1.65
1000x512 5.1e+05 from data 0.0076 0.0075 1.00
200x128 2.6e+04 from data 0.0002 0.0002 0.90
Bit identical to the previous flux_pdf over 52 cases: six bin counts from
2 to 100, 1D and 2D, float32, non-contiguous, empty, a single value,
optical depths of exactly zero so the flux lands on the last bin edge,
and negative optical depths which fall outside the bins. Thread count
does not change the answer.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
If the scale is NaN, assert. With a test thanks to Claude Opus 5
_filter_single_tau_complex replaced the optical depths in and around each strong absorber with the mean optical depth over the pre-mask spectrum, then let fluxstatistics rescale the whole array to the desired mean flux. This means that the DLA pixels are not necessarily at exactly the mean flux, but a few percent off. Now we first find a mask (using a threshold that is scaled by an estimate of the mean flux scaling factor pre-mask), then rescale only the non-masked regions, then set the masked regions to the desired mean tau. This is implemented by moving rescaling into the filtering function. The region-finding threshold is a simple fixed value, rather than a function of the desired mean flux. _filter_tau_rescale returns a new array instead of modifying its argument. The flux power spectra no longer modify the stored tau array.
_filter_tau_rescale called _mask_single_tau once per sightline, and that
function returns an all-false mask whenever the sightline has no pixel
above tau_thresh. So for the great majority of sightlines the whole cost
was two allocations and a np.max which could only ever come back empty
handed. Take the maximum over the whole array in one vectorised pass
instead and walk only the rows that can contribute, as _filter_tau did
before the masking rewrite.
The two threshold divisions move out of the loop, where they were redone
per sightline, and the emptiness assertion becomes mask.all(), which
short circuits and does not build a temporary. The compressed copy for
the second mean flux fit is skipped when nothing is masked.
The masking itself is 3.2x faster (0.057s to 0.018s at 8000x2048), but
it is only 14% of the function: the two Newton solves in mean_flux are
73%, and they already run threaded and converge in four iterations.
Best of five, seconds:
shape absorbers tau_thresh before after speedup
8000x2048 2% 1e6 0.4093 0.3692 1.11
8000x2048 2% 1e4 0.5122 0.4879 1.05
8000x2048 0% 1e6 0.3998 0.3477 1.15
8000x2048 0% 1e4 0.5303 0.4414 1.20
4000x1024 10% 1e6 0.1172 0.1046 1.12
2000x2048 2% 1e6 0.0744 0.0661 1.12
2000x2048 2% None 0.0303 0.0301 1.01
A lower tau_thresh qualifies more sightlines and so gains less. With
tau_thresh unset the prefilter is not reached and nothing changes.
Bit identical to the previous version over 25 cases: five shapes with
absorber fractions from 0 to 50%, crossed with tau_thresh of 1e6, 1e4
and None and mean_flux_desired of 0.65, 0.4 and None. Skipping a
sightline is exactly equivalent because _mask_single_tau's own loop
tests the same condition before doing anything.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This tries to fit a clipped envelope, masks the DLA center and then subtracts the wings. Currently the fit amplitude is biased high, because it uses least squares
Absorption blended with an HCD - the forest, neighbouring systems, the velocity structure of the HCD itself - is one-sided: it can only add to the optical depth, never subtract from it. A least squares fit can therefore only account for it by deepening the profile, which biased the amplitude high by a factor of a few even on noiseless data. An 8e6 HCD in a forest of mean tau 0.4 came out 3.3x too deep, inflating the mask from 237 to 433 bins and biasing the mean flux in the surviving region by 12%. Since tau/voigt_shape is an upper bound on the amplitude at every pixel, take a low quantile of it over the shoulder of the profile instead. The saturated core carries no amplitude information and the far wings are dominated by whatever else is in the spectrum, so both are excluded, with a fallback to the least saturated pixels for an absorber saturated over the whole spectrum. The mask is now honoured when choosing those pixels, as an already subtracted core would otherwise drag the quantile down. The same case now recovers the amplitude to 1%, masks 239 bins, and leaves the mean flux outside the mask within 0.1% of the truth. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
np.trapz was removed in numpy 2.0 and renamed to np.trapezoid, so equivalent_width, eq_width_hist, line_density_eq_w, omega_abs_cddf and emission.calculator all raised AttributeError. This drops support for numpy 1.x. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
_rho_abs defaulted to upthresh=None and then compared the column
densities against it, so rho_DLA() raised a TypeError on every call.
Default to np.inf instead (still accepting None), and test the upper
limit with np.isfinite.
replace_not_DLA sized H1_DLA and cofm_DLA from the number of sightlines
rather than the number of absorbers requested, so asking for more
absorbers than there are sightlines silently truncated the output and
left uninitialised rows in the saved column densities. Allocate both
with ndla rows, drop the stale [:top] slice when storing the result, and
raise if more absorbers are requested than can be found, which was never
supportable because we cannot invent new axes. Also guard the
discarded-fraction division against found == 0 and store the result
under (elem, ion) rather than a hardcoded ("H", 1).
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
_densweightdens read the particle data twice, once for the ion and once with ion=-1 for the element as a whole. Only the first of these applies the elem_den > 0 filter, so for every metal the two arrays covered different particles and multiplying them raised a ValueError. get_dens_weighted_density only worked for HI, which takes the unfiltered branch. Add an optional get_species_den argument to _read_particle_data, which returns the density in all ionisation states alongside the ionic density, sliced by the same filter. As well as being correct this halves the reads for this quantity. Existing callers are unaffected. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The check for None cofm or axis came after the astype calls that use them, so it never fired: callers got an AttributeError instead of the message telling them to use reload_file=False. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
add_noise and add_cont_error seeded the generator with spec_num, so the default of -1 raised a ValueError from np.random.seed for a single spectrum. Fall back to a fixed seed instead, keeping spec_num itself for the snr and CE lookups. The continuum error is now seeded with an odd number and the noise with the sightline index, so a spectrum no longer gets the same seed for both. Previously spectrum 0 drew its continuum error from the same stream as its noise. _read_particle_data only loaded the temperature when get_tau or (ion != -1 and elem != 'H'), but then indexed it with the cloudy filter for every ion except HI. Hydrogen ions other than HI therefore indexed a length one dummy and raised an IndexError. Move the condition into _need_temp so the load site and the concatenation in _interpolate_single_file cannot drift apart again. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
_read_particle_data converted a whole snapshot block to float32 and only then reduced it to the particles near a sightline, so every field was copied in full for the two percent of particles we keep. astype copies even when the dtype already matches, so this happened even for data that was float32 on disc. get_mass_frac was the worst of it, converting the whole nine component GFM_Metals array to take a single column. Index first and convert the result instead. On four million particles with two percent near a sightline this is about 120ms of copying per segment per species rather than 5ms, and avoids roughly 200MB of transient allocation. Position and the smoothing length are still needed in full to find the particles, so they only get copy=False. The primordial abundance fallback now builds the array it returns directly, which also drops a get_blocklen call. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
When loading all the data first, which is what we do for Arepo, each new segment was concatenated onto the arrays read so far, recopying all of them every time. Collect the segments in a list and join them once: for 32 segments of 125000 particles this is 24ms rather than 267ms. Peak memory is unchanged, as the join still needs both copies. This also removes the special case for an empty first segment, since a segment with no particles near a sightline is now simply not collected. find_all_particles did the same thing, and additionally started from a float64 empty array, which upcast the float32 positions it was given. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
get_observer_tau held a line x sightline x pixel block in double, even though compute_spectra is float32, and then handed the whole thing to res_corr, which returned a second copy of it. Store the block as float32 and convolve one line at a time inside the loop that fills it. For 29 lines of 20000 sightlines this is 4.6GB rather than 17.7GB. The line selection is unchanged: convolving in single precision does not move any sightline across the tau < 3 or tau > 0.1 cuts in testing. add_noise grew its noise array with np.append once per sightline. Fill a preallocated array instead, and return the noise from the single spectrum branch, which previously returned an empty array. The values are unchanged, but the noise is now returned with the same shape as the flux rather than flattened. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
res_corr computed sigma = 0 for a spec_res of zero, which is the default, so gaussian_filter1d divided by zero and returned a spectrum which was all NaN. get_observer_tau then found no line above its thresholds and failed to assign. Return the spectrum unsmoothed instead, which is how fluxstatistics already treats spec_res = 0. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Sometimes we have two HCDs close together and the fit fails, leaving distorting wings in the spectrum. In this case the peak absorption is not necessarily the peak of the profile. These are tests for cases of blended profiles. Co-Authored-by: Claude Opus 5
mean_flux converted the whole array to double and allocated scratch space as large as the array again, so it needed four times the memory of the optical depths it was given. Compute the flux in the precision of the data instead, promoting only the sums: it is the accumulation of many terms that needs the extra bits, not the flux of a single pixel, which is no more accurate than the optical depth it came from. On 115200 sightlines this is 2.4 times faster and the peak memory drops from 1.46 to 0.37 GB. The scale changes by 2e-8, against a convergence tolerance of 1e-6 and single precision input. The flux is now computed one fixed size block at a time rather than one block per thread, so that the scratch space does not scale with the input either. This is the same speed, and saves an array sized buffer. mean_flux also takes a mask of pixels to leave out of the mean, so that filtering DLAs no longer has to compress the array, which was a full copy. The masked mean agrees with the compressed one to 1e-15. Note integer optical depths are no longer accepted. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The module level ThreadPoolExecutor in fluxstatistics is replaced by a pool argument to mean_flux, flux_pdf and flux_power, defaulting to None for serial execution. Spectra owns a pool as a member variable and passes it in. A caller supplied pool caps its own concurrency, so the batches no longer need to be submitted in waves. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
particles_near_lines only needs to know whether a particle is within its smoothing length of some sightline, not which sightlines those are, so it is a nearest-neighbour query rather than the range query the C++ IndexTable performed. scipy's cKDTree answers it exactly, periodicity included, and threads the query itself. Querying the tree for every particle is wasteful when almost none are near a sightline, so particles first go through a chaining mesh of cells within reach of a line, binned by smoothing length in powers of two. Smoothing lengths beyond 2**maxbin cells skip the mesh, so that a few large low-density particles cannot force a very wide dilation. On 1e7 particles with h spanning 0.1 to 100, this is 1.54s against 1.03s for the C++ at 5000 sightlines, and 3.31s against 3.69s at 32000, where the omp critical insert into the std::set starts to serialise. _near_lines is left in the extension module: the new test uses it as the reference implementation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Comparing against _near_lines tied the test to the extension module, which should be free to go now that nothing else calls it. Comparing every particle to every sightline is slow but obviously correct, and the test sizes are small enough that it does not matter. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Nothing calls it now that particles_near_lines is done with numpy and a KD-tree, and the test for that compares against a brute force search rather than this. IndexTable stays: get_near_particles, and so get_near_lines, are still needed by the SPH interpolation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…articles The C++ neighbour search demanded float32, so a double precision snapshot had to be converted a whole block at a time, holding both copies at once. The numpy search has no such requirement, so pos and h can now be indexed before they are converted, like every other block. Peak allocation for the search drops from 150MB to 86MB on 4e6 double precision particles. The mesh lookup is still done in single precision whatever the input dtype, so that a double precision snapshot does not simply move the temporary somewhere else. The spare cell of dilation radius covers the rounding. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The mesh lookup was the largest single cost of the neighbour search, and was the only part still running in serial: the tree query already threads itself through scipy. The numpy calls it is made of, including the mesh gather, release the GIL, so splitting the particles into blocks and handing the blocks to the thread pool speeds the lookup up by 4.5x on 8 cores, and the whole search by 1.5x at 5000 sightlines. Spectra already keeps a pool for the flux statistics, so pass that one in. Blocking also bounds the lookup temporaries at a block rather than a snapshot, which is most of the difference in peak allocation. What is left is dominated by the tree query, which grows with the number of sightlines: at 32000 sightlines the mesh only rejects three quarters of the particles, so blocking the lookup gains proportionally less there. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
get_near_particles built a std::map for every sightline, filling them under one omp critical. That serialised the whole search: with three million particle-line pairs it took 11.5s on one core and 6.6s on eight, getting slower above four threads. Instead each thread searches a contiguous chunk of the particles into its own buffer, counting as it goes how many it finds for each line. Those counts give every thread somewhere to write, so the results can be gathered in one pass with no lock, into a compressed layout of three vectors rather than a std::map per line. Chunks are assigned to threads in order and searched in order, so the particles of a line still come out sorted by particle index, and the interpolation is unchanged to the bit. The search is now the whole cost: the gather is free, and it scales 4.4x over eight cores rather than 1.7x. A column density interpolation of two million particles takes 1.3s rather than 6.7s at 5000 sightlines with h ~ 200, and 8.7s rather than 62.5s at 32000 sightlines. Below about a million particle-line pairs the critical section was rarely contended and this changes nothing: the common case of h ~ 20 and 5000 sightlines is 0.30s either way. The traversal is now a template taking a callback, so that the search does not have to build a container per particle, and there is still only one copy of it: get_near_lines is a wrapper that collects the callback into the map it always returned. assign_cells took its neighbour list by value, which copied every list for every line. It now takes a reference. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…face Nothing outside the tests has called get_near_lines since the neighbour search moved to python, and the interpolation wants the transpose of it anyway. The traversal it wrapped is now private, leaving get_near_particles as the way in. The geometry it was used to check is still worth checking, and some of it was checked nowhere else: get_near_particles is asserted on counts and particle indices, never on a distance, so the squared distances and the periodic wrapping that produces them were only covered through get_near_lines. Those checks now run against a single particle put through get_near_particles, in a test of their own. Added a case for the distance cut itself: a particle within h of a line in each coordinate separately but further than h away. The old checks all had the line outside the search bounds as well, so loosening the cut to 4h^2 left them passing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The C _interpolate_2d rebuilt a gsl_interp2d_bicubic interpolator on every call. Build a RectBivariateSpline once in build_interp instead and evaluate it with ev(). Accuracy against the exact rate network is unchanged (3.16e-6 vs 3.11e-6 max error over the test sample); the two interpolants differ from each other by at most 1.3e-5 because GSL uses finite-difference derivatives where FITPACK fits a not-a-knot spline. This loses the OpenMP parallelism, so wall-clock for 2M points goes from 0.27s on 8 cores to 1.86s on one, though total CPU time is slightly lower. The spline extrapolates silently below the density grid where GSL raised GSL_EDOM, so check that bound explicitly. GSL was used only by this function, so drop it from the build entirely. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
RectBivariateSpline.ev holds the GIL, so handing chunks of it to the Spectra thread pool gained nothing (1.01x on 4 threads). The grids are uniform, so map_coordinates interpolates the same bicubic spline, and it does release the GIL. Prefilter the coefficients once in build_interp and evaluate with prefilter=False. For 4M particles on 8 cores: 1357ms serial with RectBivariateSpline, 241ms on the pool, against 270ms for the OpenMP GSL code we started from. Below MINPARALLEL particles the dispatch is not worth it and we interpolate in serial. map_coordinates takes the spline outside the table from a boundary condition rather than a not-a-knot end condition, which costs ~1e-3 in the outermost cell. Pad the table with real rate network values instead: the boundary error decays by ~3.7 per cell, so NPAD=8 leaves it at 3e-7, below the ~2e-6 error of the interpolation itself. Tested at the edges against the exact rate network, since a spline over the unpadded table is the less accurate reference there. The pool is set as an attribute on gasprop rather than passed to its constructor, so that gasprop classes which do not use one need not accept it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
os.sched_getaffinity is Linux-only, so constructing a Spectra raised AttributeError on macOS and Windows. Prefer os.process_cpu_count, which is the portable successor and respects the affinity mask where the platform has one, falling back to the affinity mask on python < 3.13 and to os.cpu_count where neither exists. Neither fallback can run on linux with a recent python, so the test fakes their absence. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
assign_cells wrote to arr2[0] and arr2[1] on a zero-length allocation whenever a sightline had no cells near it: Ncells is then 0, but the loop over the grid points still ran with min_ind left at its initialiser. This corrupted the heap from several threads at once, as the callers in part_int run it under an omp parallel for. Return the empty allocation before the grid loop; the particle loop in the callers was already empty for such a line. In the python module, PyArray_GETCONTIGUOUS and PyArray_SimpleNew were dereferenced through PyArray_DATA before the checks for whether they had failed, so an allocation failure segfaulted rather than raising MemoryError. The error returns also dropped the references that PyArray_GETCONTIGUOUS had taken. Hold those references in a small RAII class, which releases them on every path out, and check for NULL before using the data. Also release the descriptor that PyArray_DescrFromType hands out in check_type. The diagnostic print at the end of load_hdf5_snapshot indexed length-1, which is not a valid index when nothing was read: length is unsigned, and we reach that print with length == 0 whenever reading the file failed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
NearParticles.part was a vector<int>, so the particle indices found by get_near_particles were silently truncated for a snapshot with more than 2^31 particles, even though the interface takes the particle count as a long long. part_int then computed Pos[3*ipart+axis-1] in int arithmetic, which is signed overflow above around 716 million particles and reads from an arbitrary offset rather than failing. Widen part to long long, and with it the index variables that are derived from it in part_int and assign_cells. The loop over the lines becomes long long as well, so that the offset into the output, i*nbins, no longer wraps when the number of lines times the number of bins exceeds 2^32. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Only dens and h were checked against the length of pos, but vel and temp are indexed with the same particle indices, so a short array was an out of bounds read rather than an exception. Check them when computing tau, which is when they are read: the caller passes a length one dummy for both when it only wants the column density. The second dimension of pos was never checked either, and cofm's was checked with PyArray_DIM, which reads past the shape when the array has too small a rank. Check the rank of each array before asking for a dimension beyond the first. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Almost all of the time in compute_tau was spent in Faddeeva::w, called seven times per bin per particle from SingleAbsorber::tau_kern_inner. The damping parameter a = gamma lambda / (4 pi b) is small: 5e-4 for Lyman alpha at 10^4 K. For such a the expansion of the Voigt function to first order in a is accurate to a few times 1e-6, and the two functions of u it needs can be tabulated. Past u = 6 the Gaussian has underflowed and only the asymptotic series for the damping term is left, which is where the wings of a strong absorber are, so those bins become cheap. Larger a still goes through Faddeeva, so nothing is silently degraded. Also hoist the sample points of the inner integral and their kernel weights out of the bin loop and into the constructor: neither depends on the bin, and evaluating the quintic kernel for every bin is expensive. This part is exact. compute_tau is 6-8x faster on eight threads for a test problem with two million particles and 512 sightlines, and the absorption itself 13-16x; the rest is now mostly the neighbour search. Optical depths agree with the old code to one part in 1e10. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The index tables were sorted on one transverse coordinate only, so a search had to look at every line within a smoothing length in that coordinate and reject nearly all of them on the second: work that grows with the number of lines, and that reads the axis and cofm tables at random as it goes. Walking the two red-black trees cost a further 30ns per particle whether or not anything was near. Bucket the lines onto a uniform grid in both transverse coordinates instead, one grid per sightline axis so that every line in a grid is indexed on the same two coordinates. The lines of a cell are contiguous and carry their own coordinates, so a search reads only the cells it needs, in order, and the number of lines it looks at is set by how densely the lines are spread rather than by how many there are. The search is 1.5x faster for 512 sightlines and 7.8x for 20000, and its cost is now flat in the number of sightlines: 0.026s against 0.045s for a million particles, where the map went from 0.038s to 0.35s. compute_colden is 1.1x to 3.1x faster over the same range. The neighbour lists are identical to the last digit, which needs the line positions kept in double as they are in cofm. The index table now buckets a line by its axis, so an axis outside 1 to 3 would be an out of bounds access rather than a bad answer. Check for it in the python module, which did not check before. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Remove the casts in spectra.py
The new C++ searching is so fast that the python pre-filtering dominates the time cost. We are a factor of 2 faster on a 3072^3 simulation snapshot without it. Claude confirms that we are faster even on a small box. The cost is a higher peak memory.
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
1 participant
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.
This introduces a better masking of HCDs that corrects for the damping wings. It is also the result of running claude over the codebase. It is ~10x faster.
In a small number of cases there are significant changes to the output tau. Confirmed that this is due to the move to double-precision: there are a very small number of particles where the floating point roundoff previously took them in or out of the kernel (about 1particle in 10^7). The new version is more accurate.