Abins interpolated broadening: tolerate direction changes in sigma

The initial implementation assumed that the input series of sigma
values would monotonically increase, which is reasonable for TOSCA but
not for other instruments.

The logic is modified so that instead of
the input data, the list of sigma values used for convolution is
required to be sorted. As this list was generated within the same
function and ranges from (min(sigma), max(sigma)), the condition is
always satisfied.

The unit tests are not affected, but there may be a small numerical
difference as another implementation detail has been fixed in the
process; the interpolation x-values ('sigma factors') should be the
ratio between sigma and sigma at the lower bound of the current
window. In the initial implementation this was mistakenly identified
as sigma_chunk[0], the left-most value in the block _out of the actual
sigma values_; it should be the exact convolved width from `sigma_samples`.

scripts/AbinsModules/Instruments/Broadening.py

@@ -558,16 +558,18 @@ def interpolated_broadening(sigma=None, points=None, bins=None,
     spectra = np.array([convolve(hist, kernel, mode='same') for kernel in kernels])
 
     # Interpolate with parametrised relationship
-    sigma_locations = np.searchsorted(sigma, sigma_samples) # locations in full sigma of sampled values
+    sigma_locations = np.searchsorted(sigma_samples, sigma) # locations in sampled values of points from sigma
     spectrum = np.zeros_like(points)
-    for sample_i, left_sigma_i, right_sigma_i in zip(range(len(sigma_samples)),
-                                                     sigma_locations[:-1],
-                                                     sigma_locations[1:]):
-        sigma_chunk = sigma[left_sigma_i:right_sigma_i]
-
-        sigma_factors = sigma_chunk / sigma_chunk[0]
-        left_mix = np.polyval(mix_functions[function][spacing]['lower'], sigma_factors)
-        right_mix = np.polyval(mix_functions[function][spacing]['upper'], sigma_factors)
-        spectrum[left_sigma_i:right_sigma_i] = (left_mix * spectra[sample_i, left_sigma_i:right_sigma_i]
-                                                + right_mix * spectra[sample_i + 1, left_sigma_i:right_sigma_i])
+    # Samples with sigma == min(sigma) are a special case: copy directly from spectrum
+    spectrum[sigma_locations==0] = spectra[0, sigma_locations==0]
+
+    for i in range(1, len(sigma_samples)):
+        masked_block = (sigma_locations == i)
+        sigma_factors = sigma[masked_block] / sigma_samples[i - 1]
+        lower_mix = np.polyval(mix_functions[function][spacing]['lower'], sigma_factors)
+        upper_mix = np.polyval(mix_functions[function][spacing]['upper'], sigma_factors)
+
+        spectrum[masked_block] = (lower_mix * spectra[i-1, masked_block]
+                                  + upper_mix * spectra[i, masked_block])
+
     return points, spectrum