Harrison - PMCL/EYEA - SWE: Test in new breakpoint fitting algorithm and incorporate Quality Flags

Author: Menlo Innovations - CAVA Project

imap_l3_processing/cdf/config/imap_swe_l3_variable_attrs.yaml

@@ -1671,18 +1671,18 @@ inst_az:
     bin.
   DATA_TYPE: CDF_REAL4
   RECORD_VARYING: NRV
-swp_flags:
-  NAME: swp_flags
+swe_flags:
+  NAME: swe_flags
   DATA_TYPE: CDF_UINT2
-  CATDESC: Quality flags for SWAPI L3a
+  CATDESC: Quality flags for SWE L3
   DEPEND_0: epoch
   VAR_TYPE: support_data
   RECORD_VARYING: RV
-  FIELDNAM: Quality Flags for SWAPI L3a
+  FIELDNAM: Quality Flags for SWE L3
   FORMAT: I5
   UNITS: ' '
   VALIDMIN: 0
   VALIDMAX: 65534
   FILLVAL: 65535
-  LABLAXIS: Quality flags for SWAPI L3a
+  LABLAXIS: Quality flags for SWE L3
 

imap_l3_processing/swe/l3/models.py

@@ -89,7 +89,7 @@
 CORE_TEMPERATURE_TENSOR_INTEGRATED_CDF_VAR_NAME = "core_temperature_tensor_integrated"
 HALO_TEMPERATURE_TENSOR_INTEGRATED_CDF_VAR_NAME = "halo_temperature_tensor_integrated"
 TOTAL_TEMPERATURE_TENSOR_INTEGRATED_CDF_VAR_NAME = "total_temperature_tensor_integrated"
-SWP_FLAGS_VAR_NAME = "swp_flags"
+SWE_FLAGS_VAR_NAME = "swe_flags"
 INTEGRATED_LABEL = "integrated_label"
 TEMPERATURE_TENSOR_LABEL = "temperature_tensor_label"
 RTN_LABEL = "rtn_label"
@@ -221,7 +221,7 @@ class SweL3Data(DataProduct):
     raw_psd_by_phi_rebinned: np.ndarray
     raw_psd_by_theta_rebinned: np.ndarray
     # other
-    swp_flags: np.ndarray
+    swe_flags: np.ndarray
 
     def to_data_product_variables(self) -> list[DataProductVariable]:
         return [
@@ -381,7 +381,7 @@ def to_data_product_variables(self) -> list[DataProductVariable]:
                                 value=self.moment_data.halo_temperature_tensor_integrated),
             DataProductVariable(TOTAL_TEMPERATURE_TENSOR_INTEGRATED_CDF_VAR_NAME,
                                 value=self.moment_data.total_temperature_tensor_integrated),
-            DataProductVariable(SWP_FLAGS_VAR_NAME, value=self.swp_flags),
+            DataProductVariable(SWE_FLAGS_VAR_NAME, value=self.swe_flags),
             DataProductVariable("raw_1d_psd_rebinned", value=self.raw_1d_psd_rebinned),
             DataProductVariable("raw_psd_by_theta_rebinned", value=self.raw_psd_by_theta_rebinned),
             DataProductVariable("raw_psd_by_phi_rebinned", value=self.raw_psd_by_phi_rebinned),

imap_l3_processing/swe/l3/science/pitch_calculations.py

@@ -10,15 +10,6 @@
 from imap_l3_processing.swe.l3.models import SweConfiguration
 from imap_l3_processing.swe.quality_flags import SweL3Flags
 
-def piece_wise_model(x: np.ndarray, b0: float, b1: float,
-                     b2: float, b3: float, b4: float, b5: float) -> np.ndarray:
-    return np.log(np.piecewise(x, [x <= b2, (x > b2) & (x <= b4), x > b4],
-                               [
-                                   lambda x: b0 * np.exp(-b1 * x),
-                                   lambda x: b0 * np.exp(b2 * (b3 - b1)) * np.exp(-b3 * x),
-                                   lambda x: b0 * np.exp(b2 * (b3 - b1)) * np.exp(b4 * (b5 - b3)) * np.exp(-b5 * x),
-                               ]))
-
 def mec_breakpoint_finder(energies: np.ndarray, averaged_psd: np.ndarray) -> tuple[float, float, SweL3Flags]:
     """
     Input:
@@ -165,101 +156,6 @@ def refine_breakpoint_value(energy, psd, breakpoint_value, num_points):
     return sc_pot_output, ch_break_output, FALLBACK_POTENTIAL_ESTIMATE | BACKUP_SPLINE_UNRESOLVED | POTENTIAL_FIT_UNCONVERGED | BREAKPOINT_FIT_UNCONVERGED
 
 
-def find_breakpoints(energies: np.ndarray, averaged_psd: np.ndarray, latest_spacecraft_potentials: list[float],
-                     latest_core_halo_break_points: list[float],
-                     config: SweConfiguration) -> tuple[
-    float, float]:
-    log_psd = np.log(averaged_psd)
-    slopes = -np.diff(log_psd) / np.diff(energies)
-    slope_ratios = slopes[1:] / slopes[:-1]
-    numb = np.max(np.nonzero(slope_ratios > config['slope_ratio_cutoff_for_potential_calc']), initial=0)
-
-    if not numb:
-        return latest_spacecraft_potentials[-1], latest_core_halo_break_points[-1]
-
-    energies = energies[:numb]
-    log_psd = log_psd[:numb]
-    b1: float = slopes[0]
-    core_index = np.searchsorted(energies, config["core_energy_for_slope_guess"]) - 1
-    halo_index = np.searchsorted(energies, config["halo_energy_for_slope_guess"]) - 1
-    b3: float = slopes[core_index]
-    b5: float = slopes[halo_index]
-    b0: float = np.exp(log_psd[0] + b1 * energies[0])
-    initial_guesses = (
-        b0, b1, np.average(latest_spacecraft_potentials), b3, np.average(latest_core_halo_break_points), b5)
-
-    first_min_index = 0
-    for i in range(1, len(slope_ratios) - 1):
-        if slope_ratios[i - 1] > slope_ratios[i] < slope_ratios[i + 1]:
-            first_min_index = i
-            break
-
-    last_max_index = 0
-    for i in reversed(range(1 + first_min_index, len(slopes) - 1)):
-        if slopes[i - 1] < slopes[i] > slopes[i + 1]:
-            last_max_index = i
-            break
-
-    if last_max_index < config["refit_core_halo_breakpoint_index"]:
-        delta_b2 = -1.5
-        delta_b4 = -10
-    else:
-        delta_b2 = -1.0
-        delta_b4 = 10
-    return try_curve_fit_until_valid(energies, log_psd, initial_guesses, latest_spacecraft_potentials[-1],
-                                     latest_core_halo_break_points[-1], delta_b2, delta_b4)
-
-
-def ls_fit(x, y, initial_b):
-    b = np.copy(initial_b)
-    flamda = 0.01
-    deltaa = 0.05 * np.abs(b)
-    nloops = 14
-    prev_chisqr = 0.0
-    sigmay = None
-    nterms = len(b)
-    yfit = np.zeros(len(x))
-    sigmaa = np.zeros(nterms)
-
-    for k in range(nloops):
-        chisqr = curve_fit(6, len(x), x, y, sigmay, b, deltaa, flamda, sigmaa, yfit)
-        delta_chisqr = np.abs(prev_chisqr - chisqr)
-        if delta_chisqr < 0.001:
-            break
-        prev_chisqr = chisqr
-    return b
-
-
-def try_curve_fit_until_valid(energies: np.ndarray, log_psd: np.ndarray, initial_guesses: tuple[float, ...],
-                              latest_spacecraft_potential: float, latest_core_halo_breakpoint: float,
-                              delta_b2: float, delta_b4: float) -> tuple[float, float]:
-    b = ls_fit(energies, log_psd, initial_guesses)
-
-    def bad_fit(b):
-        return (b[1] <= 0 or
-                b[3] <= 0 or
-                b[5] <= 0 or
-                b[2] >= b[4] or
-                b[4] <= 15 or
-                b[2] <= energies[0] or
-                b[2] >= 20 or
-                b[2] >= 2 * latest_spacecraft_potential)
-
-    attempt_count = 0
-
-    modified_guesses = list(initial_guesses)
-    while bad_fit(b):
-        if attempt_count < 3:
-            modified_guesses[2] += delta_b2
-            modified_guesses[4] += delta_b4
-            b = ls_fit(energies, log_psd, modified_guesses)
-            attempt_count += 1
-        else:
-            return latest_spacecraft_potential, latest_core_halo_breakpoint
-
-    return b[2], b[4]
-
-
 def average_over_look_directions(phase_space_density: np.ndarray, geometric_weights: np.ndarray,
                                  minimum_psd_value: float) -> np.ndarray:
     enforced_minimum_psd = np.maximum(phase_space_density, minimum_psd_value)
@@ -530,101 +426,3 @@ def swe_rebin_intensity_by_pitch_angle_and_gyrophase(intensity_data: np.ndarray[
         averaged_rebinned_intensity_by_pa_and_gyro, averaged_rebinned_intensity_by_pa_only,
         uncertainties_by_pa_and_gyro,
         uncertainties_by_pa_only)
-
-
-def curve_fit(nterms: int, npts: int, x, y, sigmay, b, deltaa, flamda, sigmaa, yfit):
-    nfree = npts - nterms
-    ngoes = 1
-    alpha = np.zeros((nterms, nterms))
-    weights = np.ones(npts)
-    beta = np.zeros(nterms)
-    b2 = np.zeros(nterms)
-    array = np.full_like(alpha, np.nan)
-    invarray = np.full_like(array, np.nan)
-
-    if nfree <= 0:
-        chisqr = 0.0
-        return chisqr
-
-    for i in range(npts):
-
-        deriv = fderiv(nterms, x[i], b, deltaa)
-
-        for j in range(nterms):
-            beta[j] += weights[i] * (y[i] - functn(nterms, x[i], b)) * deriv[j]
-
-            for k in range(nterms):
-                alpha[j][k] += weights[i] * deriv[j] * deriv[k]
-
-    for i in range(npts):
-        yfit[i] = functn(nterms, x[i], b)
-
-    chisq1 = fchisqr(npts, y, weights, nfree, yfit)
-    chisqr = np.inf
-
-    while ((chisq1 - chisqr < 0.0) and (ngoes <= 5)):
-        for j in range(nterms):
-            for k in range(nterms):
-                array[j][k] = alpha[j][k] / np.sqrt(alpha[j][j] * alpha[k][k])
-            array[j][j] = 1 + flamda
-
-        invarray = np.linalg.inv(array)
-
-        for j in range(nterms):
-            b2[j] = b[j]
-
-            for k in range(nterms):
-                b2[j] += beta[k] * invarray[j][k] / np.sqrt(alpha[j][j] * alpha[k][k])
-
-        for i in range(npts):
-            yfit[i] = functn(nterms, x[i], b2)
-
-        chisqr = fchisqr(npts, y, weights, nfree, yfit)
-        ngoes += 1
-        flamda *= 10.0
-
-    for j in range(nterms):
-        b[j] = b2[j]
-        sigmaa[j] = np.sqrt(invarray[j][j] / alpha[j][j])
-
-    flamda /= 10.0
-    return chisqr
-
-
-def fderiv(nterms, x, b, deltaa):
-    b2 = np.copy(b)
-    deriv = np.full(nterms, np.nan)
-
-    for j in range(nterms):
-        b2[j] = b[j]
-
-    for j in range(nterms):
-        bj = b2[j]
-        delta = deltaa[j]
-
-        b2[j] = bj + delta
-        yfit1 = functn(nterms, x, b2)
-
-        b2[j] = bj - delta
-        yfit2 = functn(nterms, x, b2)
-
-        deriv[j] = (yfit1 - yfit2) / (2 * delta)
-
-        b2[j] = bj
-
-    return deriv
-
-
-def functn(nterms, x, b):
-    assert nterms == len(b)
-    return piece_wise_model(x, *b)
-
-
-def fchisqr(npts, y, weights, free, yfit):
-    chisq = 0
-    if free <= 0:
-        return 0.0
-    for i in range(npts):
-        chisq += weights[i] * (y[i] - yfit[i]) * (y[i] - yfit[i])
-
-    return chisq / free