Skip to content

Commit 80f4e5b

Browse files
pleasant-menlopleasant-menlo
committed
Harrison 3063/#36 - PMCL/LFAR - Maps spx: Add chi-squared variable
Currently not a "reduced" chi-squared. If we want it "reduced" we will need to divide by 2 because we have 2 degrees of freedom in the fit.
1 parent 45c7982 commit 80f4e5b

10 files changed

Lines changed: 211 additions & 26 deletions

File tree

‎imap_l3_processing/cdf/config/imap_hi_l3_variable_attrs.yaml‎

Lines changed: 26 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -427,4 +427,29 @@ ena_spectral_scalar_stat_uncert:
427427
LABL_PTR_2: longitude_label
428428
LABL_PTR_3: latitude_label
429429
SCALETYP: log
430-
VAR_NOTES: Calculated by propagating ENA intensity statistical uncertainty through Levenberg-Marquardt fitting
430+
VAR_NOTES: Calculated by propagating ENA intensity statistical uncertainty through Levenberg-Marquardt fitting
431+
ena_spectral_index_chisq:
432+
NAME: ena_spectral_index_chisq
433+
DATA_TYPE: CDF_REAL4
434+
CATDESC: ENA spectral index chi-squared
435+
DEPEND_0: epoch
436+
DEPEND_1: energy
437+
DEPEND_2: longitude
438+
DEPEND_3: latitude
439+
VAR_TYPE: data
440+
RECORD_VARYING: RV
441+
DISPLAY_TYPE: map_image
442+
VARIABLE_PURPOSE: ' '
443+
FIELDNAM: Spectral Index Chi-Squared
444+
FORMAT: F12.6
445+
LABLAXIS: spectral index chi-sq
446+
UNITS: ' '
447+
VALIDMIN: 0.0
448+
VALIDMAX: 1000000.0
449+
SCALEMIN: 0.0
450+
SCALEMAX: 10.0
451+
FILLVAL: -1.00E+31
452+
LABL_PTR_1: energy_label
453+
LABL_PTR_2: longitude_label
454+
LABL_PTR_3: latitude_label
455+
SCALETYP: linear

‎imap_l3_processing/cdf/config/imap_lo_l3_variable_attrs.yaml‎

Lines changed: 26 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -690,4 +690,29 @@ ena_spectral_scalar_stat_uncert:
690690
LABL_PTR_2: longitude_label
691691
LABL_PTR_3: latitude_label
692692
SCALETYP: log
693-
VAR_NOTES: Calculated by propagating ENA intensity statistical uncertainty through Levenberg-Marquardt fitting
693+
VAR_NOTES: Calculated by propagating ENA intensity statistical uncertainty through Levenberg-Marquardt fitting
694+
ena_spectral_index_chisq:
695+
NAME: ena_spectral_index_chisq
696+
DATA_TYPE: CDF_REAL4
697+
CATDESC: ENA spectral index chi-squared
698+
DEPEND_0: epoch
699+
DEPEND_1: energy
700+
DEPEND_2: longitude
701+
DEPEND_3: latitude
702+
VAR_TYPE: data
703+
RECORD_VARYING: RV
704+
DISPLAY_TYPE: map_image
705+
VARIABLE_PURPOSE: ' '
706+
FIELDNAM: Spectral Index Chi-Squared
707+
FORMAT: F12.6
708+
LABLAXIS: spectral index chi-sq
709+
UNITS: ' '
710+
VALIDMIN: 0.0
711+
VALIDMAX: 1000000.0
712+
SCALEMIN: 0.0
713+
SCALEMAX: 10.0
714+
FILLVAL: -1.00E+31
715+
LABL_PTR_1: energy_label
716+
LABL_PTR_2: longitude_label
717+
LABL_PTR_3: latitude_label
718+
SCALETYP: linear

‎imap_l3_processing/cdf/config/imap_ultra_l3_variable_attrs.yaml‎

Lines changed: 26 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -427,4 +427,29 @@ ena_spectral_scalar_stat_uncert:
427427
LABL_PTR_2: longitude_label
428428
LABL_PTR_3: latitude_label
429429
SCALETYP: log
430-
VAR_NOTES: Calculated by propagating ENA intensity statistical uncertainty through Levenberg-Marquardt fitting
430+
VAR_NOTES: Calculated by propagating ENA intensity statistical uncertainty through Levenberg-Marquardt fitting
431+
ena_spectral_index_chisq:
432+
NAME: ena_spectral_index_chisq
433+
DATA_TYPE: CDF_REAL4
434+
CATDESC: ENA spectral index chi-squared
435+
DEPEND_0: epoch
436+
DEPEND_1: energy
437+
DEPEND_2: longitude
438+
DEPEND_3: latitude
439+
VAR_TYPE: data
440+
RECORD_VARYING: RV
441+
DISPLAY_TYPE: map_image
442+
VARIABLE_PURPOSE: ' '
443+
FIELDNAM: Spectral Index Chi-Squared
444+
FORMAT: F12.6
445+
LABLAXIS: spectral index chi-sq
446+
UNITS: ' '
447+
VALIDMIN: 0.0
448+
VALIDMAX: 1000000.0
449+
SCALEMIN: 0.0
450+
SCALEMAX: 10.0
451+
FILLVAL: -1.00E+31
452+
LABL_PTR_1: energy_label
453+
LABL_PTR_2: longitude_label
454+
LABL_PTR_3: latitude_label
455+
SCALETYP: linear

‎imap_l3_processing/maps/map_models.py‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@
3939
ENA_SPECTRAL_INDEX_STAT_UNC_VAR_NAME = "ena_spectral_index_stat_uncert"
4040
ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME = "ena_spectral_scalar"
4141
ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME = "ena_spectral_scalar_stat_uncert"
42+
ENA_SPECTRAL_INDEX_CHISQ_VAR_NAME = "ena_spectral_index_chisq"
4243

4344
ENA_INTENSITY_VAR_NAME = "ena_intensity"
4445
ENA_INTENSITY_STAT_UNCERT_VAR_NAME = "ena_intensity_stat_uncert"
@@ -125,6 +126,7 @@ class SpectralIndexMapData(MapData):
125126
ena_spectral_index_stat_uncert: np.ndarray
126127
ena_spectral_index_scalar_coefficient: np.ndarray
127128
ena_spectral_index_scalar_coefficient_stat_uncert: np.ndarray
129+
ena_spectral_index_chisq: np.ndarray
128130

129131

130132
@dataclass
@@ -401,6 +403,9 @@ def to_healpix_skymap(self) -> HealpixSkyMap:
401403
"ena_spectral_index_scalar_coefficient": (full_shape, self.spectral_index_map_data.ena_spectral_index_scalar_coefficient),
402404
"ena_spectral_index_scalar_coefficient_stat_uncert": (full_shape,
403405
self.spectral_index_map_data.ena_spectral_index_scalar_coefficient_stat_uncert),
406+
"ena_spectral_index_chisq": (
407+
full_shape,
408+
self.spectral_index_map_data.ena_spectral_index_chisq),
404409
},
405410
coords={
406411
CoordNames.TIME.value: self.spectral_index_map_data.epoch,
@@ -511,6 +516,7 @@ def _spectral_index_data_variables(data: SpectralIndexMapData) -> list[DataProdu
511516
DataProductVariable(ENA_SPECTRAL_INDEX_STAT_UNC_VAR_NAME, data.ena_spectral_index_stat_uncert),
512517
DataProductVariable(ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME, data.ena_spectral_index_scalar_coefficient),
513518
DataProductVariable(ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME, data.ena_spectral_index_scalar_coefficient_stat_uncert),
519+
DataProductVariable(ENA_SPECTRAL_INDEX_CHISQ_VAR_NAME, data.ena_spectral_index_chisq),
514520
]
515521

516522

‎imap_l3_processing/maps/spectral_fit.py‎

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,8 @@ def calculate_spectral_index_for_multiple_ranges(intensity_data: IntensityMapDat
3232
ena_spectral_index=np.concat([m.ena_spectral_index for m in spectral_maps], axis=1),
3333
ena_spectral_index_stat_uncert=np.concat([m.ena_spectral_index_stat_uncert for m in spectral_maps], axis=1),
3434
ena_spectral_index_scalar_coefficient=np.concat([m.ena_spectral_index_scalar_coefficient for m in spectral_maps], axis=1),
35-
ena_spectral_index_scalar_coefficient_stat_uncert=np.concat([m.ena_spectral_index_scalar_coefficient_stat_uncert for m in spectral_maps], axis=1)
35+
ena_spectral_index_scalar_coefficient_stat_uncert=np.concat([m.ena_spectral_index_scalar_coefficient_stat_uncert for m in spectral_maps], axis=1),
36+
ena_spectral_index_chisq=np.concat([m.ena_spectral_index_chisq for m in spectral_maps], axis=1),
3637
)
3738

3839

@@ -82,7 +83,7 @@ def fit_spectral_index_map(intensity_data: IntensityMapData) -> SpectralIndexMap
8283
mean_energy = np.sqrt(min_energy * max_energy)
8384
new_energy_label = f"{min_energy} - {max_energy}"
8485

85-
output_scalar_coefficients, output_scalar_errors, output_gammas, output_gamma_errors = fit_arrays_to_power_law(fluxes, uncertainty, energy)
86+
output_scalar_coefficients, output_scalar_errors, output_gammas, output_gamma_errors, chisq = fit_arrays_to_power_law(fluxes, uncertainty, energy)
8687
mean_obs_date = calculate_datetime_weighted_average(intensity_data.obs_date,
8788
weights=intensity_data.exposure_factor,
8889
axis=1, keepdims=True)
@@ -94,7 +95,7 @@ def fit_spectral_index_map(intensity_data: IntensityMapData) -> SpectralIndexMap
9495
output_gammas[positive_gammas] = np.nan
9596
output_scalar_coefficients[positive_gammas] = np.nan
9697
output_gamma_errors[positive_gammas] = np.nan
97-
98+
chisq[positive_gammas] = np.nan
9899

99100
return SpectralIndexMapData(
100101
epoch=intensity_data.epoch,
@@ -113,11 +114,12 @@ def fit_spectral_index_map(intensity_data: IntensityMapData) -> SpectralIndexMap
113114
ena_spectral_index_stat_uncert=output_gamma_errors,
114115
ena_spectral_index_scalar_coefficient=output_scalar_coefficients,
115116
ena_spectral_index_scalar_coefficient_stat_uncert=output_scalar_errors,
117+
ena_spectral_index_chisq=chisq,
116118
)
117119

118120

119121
def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energy: np.ndarray) -> tuple[
120-
np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
122+
np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
121123
par_info = [
122124
{'limits': [0.0, 1000.0]},
123125
{'limits': [0.0, 1000.0]},
@@ -128,6 +130,7 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
128130
output_gamma_errors = np.full_like(output_gammas, np.nan)
129131
output_scalar_coefficients = np.full_like(output_gammas, np.nan)
130132
output_scalar_coefficients_errors = np.full_like(output_gammas, np.nan)
133+
output_chisqs = np.full_like(output_gammas, np.nan)
131134

132135
for epoch in range(fluxes.shape[0]):
133136
intensity = fluxes[epoch].reshape((fluxes[epoch].shape[0], -1))
@@ -137,6 +140,7 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
137140
gamma_errors = np.full_like(gammas, np.nan)
138141
scalar_coefficients = np.full(intensity.shape[-1], np.nan, dtype=float)
139142
scalar_coefficient_errors = np.full(intensity.shape[-1], np.nan, dtype=float)
143+
chisqs = np.full(intensity.shape[-1], np.nan, dtype=float)
140144

141145
for i in range(intensity.shape[-1]):
142146
flux = intensity[:, i]
@@ -163,11 +167,13 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
163167
scalar_coefficients[i] = a
164168
gamma_errors[i] = gamma_error
165169
scalar_coefficient_errors[i] = a_error
170+
chisqs[i] = fit.fnorm
166171
output_gammas[epoch, 0] = gammas.reshape(fluxes.shape[2:])
167172
output_gamma_errors[epoch, 0] = gamma_errors.reshape(fluxes.shape[2:])
168173
output_scalar_coefficients[epoch, 0] = scalar_coefficients.reshape(fluxes.shape[2:])
169174
output_scalar_coefficients_errors[epoch, 0] = scalar_coefficient_errors.reshape(fluxes.shape[2:])
170-
return output_scalar_coefficients, output_scalar_coefficients_errors, output_gammas, output_gamma_errors
175+
output_chisqs[epoch, 0] = chisqs.reshape(fluxes.shape[2:])
176+
return output_scalar_coefficients, output_scalar_coefficients_errors, output_gammas, output_gamma_errors, output_chisqs
171177

172178

173179
def power_law(params, **kwargs):

‎tests/hi/test_hi_processor.py‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -53,10 +53,12 @@ def test_process_spectral_fit(self, mock_save_data, mock_spectral_fit,
5353
descriptor="h90-spx-h-hf-sp-full-hae-4deg-6mo",
5454
)
5555

56-
expected_unc = np.full((1, 1, 3, 2), 1.0)
57-
expected_index = np.full_like(expected_unc, 2.0)
58-
expected_scalar = np.full_like(expected_unc, 3.0)
59-
mock_spectral_fit.return_value = expected_scalar, expected_index, expected_unc
56+
expected_index_unc = np.full((1, 1, 3, 2), 1.0)
57+
expected_index = np.full_like(expected_index_unc, 2.0)
58+
expected_scalar = np.full_like(expected_index_unc, 3.0)
59+
expected_scalar_unc = np.full_like(expected_index_unc, 4.0)
60+
expected_chisq = np.full_like(expected_index_unc, 5.0)
61+
mock_spectral_fit.return_value = expected_scalar, expected_scalar_unc, expected_index, expected_index_unc, expected_chisq
6062
processor = HiProcessor(upstream_dependencies, input_metadata)
6163
product = processor.process()
6264

@@ -79,8 +81,10 @@ def test_process_spectral_fit(self, mock_save_data, mock_spectral_fit,
7981
np.testing.assert_array_equal( spectral_index_data.energy_delta_minus,expected_energy_delta_minus)
8082
np.testing.assert_array_equal(spectral_index_data.energy_delta_plus, expected_energy_delta_plus)
8183
np.testing.assert_array_equal(spectral_index_data.ena_spectral_index,expected_index)
82-
np.testing.assert_array_equal(spectral_index_data.ena_spectral_index_stat_uncert, expected_unc)
84+
np.testing.assert_array_equal(spectral_index_data.ena_spectral_index_stat_uncert, expected_index_unc)
8385
np.testing.assert_array_equal(spectral_index_data.ena_spectral_index_scalar_coefficient, expected_scalar)
86+
np.testing.assert_array_equal(spectral_index_data.ena_spectral_index_scalar_coefficient_stat_uncert, expected_scalar_unc)
87+
np.testing.assert_array_equal(spectral_index_data.ena_spectral_index_chisq, expected_chisq)
8488
np.testing.assert_array_equal(spectral_index_data.energy, expected_energy)
8589
np.testing.assert_array_equal(spectral_index_data.latitude, intensity_data.latitude)
8690
np.testing.assert_array_equal(spectral_index_data.longitude, intensity_data.longitude)

‎tests/maps/test_builders.py‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,8 @@ def create_rectangular_spectral_index_map_data(epoch=None, epoch_delta=None, lon
157157
(len(epoch), len(energy), len(lon), len(lat)),
158158
fill_value=1)
159159
ena_spectral_index_scalar_coefficient = np.full_like(spectral_index, 3)
160+
ena_spectral_index_scalar_coefficient_stat_unc = np.full_like(spectral_index, 4)
161+
ena_spectral_index_chisq = np.full_like(spectral_index, 5)
160162

161163
if isinstance(spectral_index, np.ndarray):
162164
more_real_flux = spectral_index
@@ -179,7 +181,9 @@ def create_rectangular_spectral_index_map_data(epoch=None, epoch_delta=None, lon
179181
solid_angle=np.full_like(more_real_flux, 0),
180182
ena_spectral_index=spectral_index,
181183
ena_spectral_index_stat_uncert=spectral_index_stat_unc,
182-
ena_spectral_index_scalar_coefficient=ena_spectral_index_scalar_coefficient
184+
ena_spectral_index_scalar_coefficient=ena_spectral_index_scalar_coefficient,
185+
ena_spectral_index_scalar_coefficient_stat_uncert=ena_spectral_index_scalar_coefficient_stat_unc,
186+
ena_spectral_index_chisq=ena_spectral_index_chisq,
183187
),
184188
coords=RectangularCoords(
185189
latitude_delta=np.full_like(lat, 0),

‎tests/maps/test_map_models.py‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@ def test_rectangular_spectral_index_to_data_product_variables(self):
5353
ena_spectral_index_stat_uncert=sentinel.ena_spectral_index_stat_uncert,
5454
ena_spectral_index_scalar_coefficient=sentinel.ena_spectral_index_scalar_coefficient,
5555
ena_spectral_index_scalar_coefficient_stat_uncert=sentinel.ena_spectral_index_scalar_coefficient_stat_uncert,
56+
ena_spectral_index_chisq=sentinel.ena_spectral_index_chisq,
5657
),
5758
coords=RectangularCoords(
5859
latitude_delta=sentinel.latitude_delta,
@@ -84,6 +85,7 @@ def test_rectangular_spectral_index_to_data_product_variables(self):
8485
sentinel.ena_spectral_index_stat_uncert),
8586
DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient),
8687
DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient_stat_uncert),
88+
DataProductVariable(map_models.ENA_SPECTRAL_INDEX_CHISQ_VAR_NAME, sentinel.ena_spectral_index_chisq),
8789
DataProductVariable(map_models.LATITUDE_DELTA_VAR_NAME, sentinel.latitude_delta),
8890
DataProductVariable(map_models.LATITUDE_LABEL_VAR_NAME, sentinel.latitude_label),
8991
DataProductVariable(map_models.LONGITUDE_DELTA_VAR_NAME, sentinel.longitude_delta),
@@ -239,6 +241,7 @@ def test_healpix_spectral_index_to_data_product_variables(self):
239241
ena_spectral_index_stat_uncert=sentinel.ena_spectral_index_stat_uncert,
240242
ena_spectral_index_scalar_coefficient=sentinel.ena_spectral_index_scalar_coefficient,
241243
ena_spectral_index_scalar_coefficient_stat_uncert=sentinel.ena_spectral_index_scalar_coefficient_stat_uncert,
244+
ena_spectral_index_chisq=sentinel.ena_spectral_index_chisq,
242245
),
243246
coords=HealPixCoords(
244247
pixel_index=sentinel.pixel_index,
@@ -268,6 +271,7 @@ def test_healpix_spectral_index_to_data_product_variables(self):
268271
sentinel.ena_spectral_index_stat_uncert),
269272
DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient),
270273
DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient_stat_uncert),
274+
DataProductVariable(map_models.ENA_SPECTRAL_INDEX_CHISQ_VAR_NAME, sentinel.ena_spectral_index_chisq),
271275
DataProductVariable(map_models.PIXEL_INDEX_VAR_NAME, sentinel.pixel_index),
272276
DataProductVariable(map_models.PIXEL_INDEX_LABEL_VAR_NAME, sentinel.pixel_index_label),
273277
]
@@ -726,6 +730,7 @@ def test_healpix_spectral_index_map_data_to_skymap(self):
726730
ena_spectral_index_stat_uncert=fake_data_per_energy_per_pixel * 2.3,
727731
ena_spectral_index_scalar_coefficient=fake_data_per_energy_per_pixel * 2.4,
728732
ena_spectral_index_scalar_coefficient_stat_uncert=fake_data_per_energy_per_pixel * 2.5,
733+
ena_spectral_index_chisq=fake_data_per_energy_per_pixel * 2.6,
729734
)
730735

731736
healpix_spectral_index_map_data = HealPixSpectralIndexMapData(
@@ -761,9 +766,11 @@ def test_healpix_spectral_index_map_data_to_skymap(self):
761766

762767
np.testing.assert_array_equal(actual_dataset.data_vars["ena_spectral_index_scalar_coefficient"].values, spectral_index_map_data.ena_spectral_index_scalar_coefficient)
763768
np.testing.assert_array_equal(actual_dataset.data_vars["ena_spectral_index_scalar_coefficient_stat_uncert"].values, spectral_index_map_data.ena_spectral_index_scalar_coefficient_stat_uncert)
769+
np.testing.assert_array_equal(actual_dataset.data_vars["ena_spectral_index_chisq"].values, spectral_index_map_data.ena_spectral_index_chisq)
764770

765771
for key in [ "obs_date", "obs_date_range", "exposure_factor", "ena_spectral_index", "ena_spectral_index_stat_uncert",
766-
"ena_spectral_index_scalar_coefficient", "ena_spectral_index_scalar_coefficient_stat_uncert" ]:
772+
"ena_spectral_index_scalar_coefficient", "ena_spectral_index_scalar_coefficient_stat_uncert",
773+
"ena_spectral_index_chisq" ]:
767774
self.assertEqual((CoordNames.TIME.value, CoordNames.ENERGY_L2.value, CoordNames.GENERIC_PIXEL.value), actual_dataset.data_vars[key].dims)
768775
for key in [ "latitude", "longitude", "solid_angle" ]:
769776
self.assertEqual((CoordNames.GENERIC_PIXEL.value,), actual_dataset.data_vars[key].dims)

0 commit comments

Comments
 (0)