enh: vectorise spectral integrals and speed up spectrum construction

[tvare]

Summary

Four small internal speed-ups on the spectral-integration path. None of them changes any result: every patched routine is bit-for-bit identical to the current implementation (verified with np.array_equal, not allclose), and the added tests pin that.

1. DirectionalSpectrum._freq_spectrum / _SpectrumMixin._dir_spectrum — integrate all rows with one trapezoid(..., axis=-1) call instead of a Python list comprehension calling trapezoid once per row. A contiguous copy of the (strided) interpolated array keeps NumPy's summation order identical to the per-row loop, so the result is unchanged to the last bit. This is the hot spot when moment(), var(), std(), tz or spectrum1d() are evaluated for many sea states: for a 300 × 16 response spectrum the row loop is ~0.35 ms per call; the vectorised form is ~0.09 ms.
2. _SpectrumMixin.tz — integrate out the directional domain once and take both moments from it, instead of calling moment() twice. f**0 * s is exactly s, so tz equals sqrt(moment(0) / moment(2)) bit for bit (pinned by a test).
3. DirectionalSpectrum.from_spectrum1d — the direction offset theta - theta_p does not depend on frequency, so it is reduced once per direction instead of once per grid point. The spreading function is still called per grid point with the same arguments, so results are unchanged for any spread_fun, frequency-dependent or not.
4. _robust_modulus — scalar fast path (plain %, isfinite, conditional nextafter) that avoids asarray_chkfinite(...).copy() and the where=/out= machinery for 0-d inputs. Same arithmetic, same return type (0-d ndarray), same ValueError on non-finite input. This function is called once per grid point from from_spectrum1d and from BaseSpreading.__call__, so it dominated the cost of building spectra.

Measured effect

Response analysis of one vessel (84 × 16 RAO grid, 300 wave frequencies, 60 sea states, 12 response components, serial):

	before	after
build 60 WaveSpectrum via from_spectrum1d (cos² spreading)	2.78 s	1.07 s
calculate_response + std()/tz for 12 components × 60 sea states	3.46 s	0.78 s
cross-spectral moments, 21 RAO pairs × 20 sea states	1.99 s	0.42 s

All 288 output arrays from that run (wave spectra, std, tz, var, moment(1), moment(2, freq_hz=True), spectrum1d(axis=0/1), for both WaveSpectrum and WaveBinSpectrum input) are np.array_equal before and after.

Tests

• test_spectrum1d_axis1_equals_rowwise_integration / test_spectrum1d_axis0_equals_rowwise_integration (DirectionalSpectrum and DirectionalBinSpectrum): vectorised result assert_array_equal to the per-row trapezoid loop on a non-uniform direction grid.
• test_tz_equals_moment_ratio (both spectrum classes): tz == sqrt(moment(0) / moment(2)) exactly.
• test_from_spectrum1d_nonuniform_dirs_wrapping_dirp: frequency-dependent spreading, non-uniform directions, peak direction that wraps past 360°.
• _robust_modulus: 0-d return type, 0-d array input, non-finite scalar and array inputs raise ValueError.

Full suite: 690 passed; black --check and isort --check --profile black clean.

[vegard-solum-4ss]

Hi, thanks for the contribution! I will have a look at this when I get the time.