Speed up ApplyReliabilityCalibration (#1989)

* Use np.interp instead of scipy interp1d

* Revert "Use np.interp instead of scipy interp1d"

This reverts commit 3ea6cd6.

* Use np.interp instead of scipy interp1d

* update comments

* update comments

---------

Co-authored-by: Belinda Trotta <btrotta-bom@users.noreply.github.com>

improver/calibration/reliability_calibration.py

@@ -1258,10 +1258,26 @@ def _interpolate(
 
         forecast_probabilities = np.ma.getdata(forecast_threshold).flatten()
 
+        # Interpolate using scipy first to get extrapolated values at endpoints
+        # since np.interp does not allow extrapolation. We would need to change back
+        # to scipy.interpolate if we want non-linear interpolation in future.
         interpolation_function = scipy.interpolate.interp1d(
             reliability_probabilities, observation_frequencies, fill_value="extrapolate"
         )
-        interpolated = interpolation_function(forecast_probabilities.data)
+        y_0, y_1 = interpolation_function([0, 1])
+        xp = np.copy(reliability_probabilities)
+        # Extrapolation preserves the slope of the first and last segments of the piecewise
+        # linear function. Thus the slope betweeen [0, y_0] and [xp[0], fp[0]] is the same as that
+        # between [xp[0], fp[0]] and [xp[1], fp[1]], so we can replace
+        # [xp[0], fp[0]] with [0, y_0] to extend the width of the first segment of the
+        # piecewise linear function. A similar argument applies for the last segment.
+        xp[0] = 0
+        xp[-1] = 1
+        fp = np.copy(observation_frequencies)
+        fp[0] = y_0
+        fp[-1] = y_1
+
+        interpolated = np.interp(forecast_probabilities.data, xp, fp)
 
         interpolated = interpolated.reshape(shape).astype(np.float32)