23167106439

Committed 16 Mar 2026 09:37PM UTC coverage: 99.962%. Remained the same

Build # 23167106439

Build Type

push

github

Committed by

aburrell

Commit Message

DOC: update changelog

Added a description of the changes to the changelog.

Run Details

2652 of 2653 relevant lines covered (99.96%)

13.89 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

99.22

/pyIntensityFeatures/proc/fitting.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# DOI: 10.5281/zenodo.15102100
# Full license can be found in License.md
#
# DISTRIBUTION STATEMENT A: Approved for public release. Distribution is
# unlimited.
# -----------------------------------------------------------------------------
"""Functions used for fitting functions to intensity data."""

import numpy as np
from scipy.optimize import leastsq
from scipy import stats

from pyIntensityFeatures.utils import checks
from pyIntensityFeatures.utils import distributions


def estimate_peak_widths(data, data_loc, peak_inds, peak_mags):
    """Estimate the full width at half max of the peaks in a curve.

    Parameters
    ----------
    data : array-like
        Input data (y-axis)
    data_loc : array-like
        Input data location (x-axis)
    peak_inds : array-like
        Indexes of the peaks
    peak_mags : array-like
        Magnitude of the peaks

    Returns
    -------
    peak_sigmas : list-like
        FWHM values for each peak, or the closest estimate possible

    Notes
    -----
    Differs from Longden, et al. (2010) function by not assuming a fixed data
    location increment of 1.0 degrees.

    """
    # Ensure the input is array-like
    data = np.asarray(data)
    data_loc = np.asarray(data_loc)
    peak_inds = np.asarray(peak_inds)
    peak_mags = np.asarray(peak_mags)

    # Get the half-max for each peak
    half_max = 0.5 * peak_mags
    peak_sigmas = list()

    for i, ipeak in enumerate(peak_inds):
        # Step through all values of the curve data (`data`) with lower indices
        # than this peak until either the half-max threshold or another peak is
        # reached
        reached = False
        good_down = False
        j = ipeak
        while not reached and j > 0:
            if (data[j] >= half_max[i]) and (data[j - 1] <= half_max[i]):
                reached = True
                fwhm_loc = data_loc[j - 1] + (
                    data_loc[j] - data_loc[j - 1]) * (
                        half_max[i] - data[j - 1]) / (data[j] - data[j - 1])
                width_down = data_loc[ipeak] - fwhm_loc
                good_down = True
            else:
                if data[j] < data[j - 1]:
                    # Data is no longer decreasing. Since the correct width was
                    # not reached before finding the desired half-max value,
                    # set an approximate value
                    reached = True
                    width_down = data_loc[ipeak] - data_loc[j - 1]
            j -= 1

        if not good_down and not reached:
            # Reached edge of data before defining the width
            width_down = data_loc[ipeak] - data_loc[j]

        # Step through all values of the curve data (`data`) with higher
        # indices than this peak until either the half-max threshold or another
        # peak is reached
        reached = False
        good_up = False
        j = ipeak
        while not reached and j < len(data) - 1:
            if (data[j] >= half_max[i]) and (data[j + 1] <= half_max[i]):
                reached = True
                fwhm_loc = data_loc[j + 1] + (
                    data_loc[j] - data_loc[j + 1]) * (
                        half_max[i] - data[j + 1]) / (data[j] - data[j + 1])
                width_up = fwhm_loc - data_loc[ipeak]
                good_up = True
            else:
                if data[j] < data[j + 1]:
                    # Data is no longer decreasing. Since the correct width was
                    # not reached before finding the desired half-max value,
                    # set an approximate value
                    reached = True
                    width_up = data_loc[j + 1] - data_loc[ipeak]
            j += 1

        if not good_up and not reached:
            # Reached edge of data before defining the width
            width_up = data_loc[j] - data_loc[ipeak]

        # Calculate the peak sigmas using logic based on whether or not the
        # HWHM was found on each leg
        if good_up and good_down:
            peak_sigmas.append((width_up + width_down) / checks.fwhm_const)
        elif good_up:
            peak_sigmas.append((2.0 * width_up) / checks.fwhm_const)
        elif good_down:
            peak_sigmas.append((2.0 * width_down) / checks.fwhm_const)
        else:
            peak_sigmas.append(2.0 * min([width_up, width_down])
                               / checks.fwhm_const)

    return peak_sigmas


def gauss_quad_err(params, xvals, yvals, weights):
    """Calculate the error of a quadriatic Gaussian fit.

    Parameters
    ----------
    params : set
        Set of parameters used by `mult_gauss_quad`
    xvals : float or array-like
        Location values
    yvals : float or array-like
        Intensity values
    weights : float or array-like
        Weights for each data value

    Returns
    -------
    err : float or array-like
        Weighted error of the fit defined by `params`

    See Also
    --------
    utils.distributions.mult_gauss_quad

    """
    # Define the error function to fit
    err = (distributions.mult_gauss_quad(xvals, params) - yvals) * weights

    return err


def get_fitting_params(mlat_bins, ilats, ilt, mean_intensity, num_gauss=3):
    """Find the parameters for a Gaussian fit with a quadratic background.

    Parameters
    ----------
    mlat_bins : np.array
        Magnetic latitude of output bin centres
    ilats : np.array
        Indices for magnetic latitudes with valid intensity values
    ilt : int
        Index for the magnetic local time
    mean_intensity : np.array
       2D array of mean intensity values
    num_gauss : int
        Maximum number of Gaussians to fit (default=3)

    Returns
    -------
    params : list
        List of the parameters needed to form a mult-Gaussian fit with a
        quadratic background.
    ipeaks : list
        List containing indexes of the normalized intensity peaks for the
        2D `mean_intesity` array down-selected by [`ilats`, `ilt`].

    """
    # Perform a first-order polynomial fit to the data to obtain estimates for
    # the background levels in the intensity profile. The first term is the
    # slope and the second term is the intercept.
    bg = np.polyfit(mlat_bins[ilats], mean_intensity[ilats, ilt], 1)

    # Normalize the intensity curve
    norm_intensity = (mean_intensity[ilats, ilt]
                      - mean_intensity[ilats, ilt].min()) / (
                          mean_intensity[ilats, ilt].max()
                          - mean_intensity[ilats, ilt].min())

    # Find the main peak and its characteristics
    ipeaks = [norm_intensity.argmax()]
    peak_sigmas = estimate_peak_widths(norm_intensity, mlat_bins[ilats], ipeaks,
                                       [norm_intensity[ipeaks[0]]])

    # Locate any additional peaks, up to a desired maximum
    mlt_gauss = num_gauss
    while len(ipeaks) < mlt_gauss:
        # Get a Gaussian for the prior peak
        prior_gauss = distributions.gauss(
            mlat_bins[ilats], norm_intensity[ipeaks[-1]],
            mlat_bins[ilats][ipeaks[-1]], peak_sigmas[-1], 0.0)

        # Remove the gaussian fit from the data
        norm_intensity -= prior_gauss

        # Get the next peak
        inew = norm_intensity.argmax()

        # Test to see if the new peak is significant
        if not np.any((mlat_bins[ilats][inew] > mlat_bins[ilats][ipeaks]
                       - np.asarray(peak_sigmas) * 2.0)
                      & (mlat_bins[ilats][inew] < mlat_bins[ilats][ipeaks]
                         + np.asarray(peak_sigmas) * 2.0)):
            ipeaks.append(inew)
            peak_sigmas.extend(estimate_peak_widths(
                norm_intensity, mlat_bins[ilats], [ipeaks[-1]],
                [norm_intensity[ipeaks[-1]]]))
        else:
            mlt_gauss = len(ipeaks)

    # Use the peak information to set the Gaussian fitting parameters
    params = [bg[1], bg[0], 0.0]
    for i, ipeak in enumerate(ipeaks):
        params.extend([mean_intensity[ilats, ilt][ipeak],
                       mlat_bins[ilats][ipeak], peak_sigmas[i]])

    return params, ipeaks


def get_gaussian_func_fit(mlat_bins, mlt_bins, mean_intensity, std_intensity,
                          num_intensity, num_gauss=3, min_num=3,
                          min_intensity=0.0, min_lat_perc=70.0):
    """Fit intensity data using least-squares minimization in MLT bins.

    Parameters
    ----------
    mlat_bins : np.array
        Magnetic latitude of output bin centres
    mlt_bins : np.array
        Magnetic local time of output bin centres
    mean_intensity : np.array
       2D array of mean intensity values
    std_intensity : np.array
       2D array of intensity standard deviations
    num_intensity : np.array
       2D array of intensity counts per bin
    num_gauss : int
        Maximum number of Gaussians to fit (default=3)
    min_num : int
        Minimum number of samples contributing to the intensity mean
        (default=3)
    min_intensity : float
        Minimum intensity magnitude in Rayleighs (default=0.0)
    min_lat_perc : int
        Minimum percentage of latitude bins needed to define an intensity
        function (default=70.0)

    Returns
    -------
    func_params : list-like
        Set of parameters used by `mult_gauss_quad` for each MLT bin that
        defines the mean intensity as a function of magnetic latitude or set
        to a string output with a reason if no successful fit was performed.
    fit_cov : list-like
        Covariance matrix for each fit
    fit_pearsonr : list-like
        Pearson r-value for correlation between the observations and the fit
    fit_pearsonp : list-like
        Pearson p-value for correlation between the observations and the fit
    num_peaks : list-like
        Number of peaks used in the fit

    See Also
    --------
    utils.distributions.mult_gauss_quad

    """
    # Initialize the output
    func_params = list()
    fit_cov = list()
    fit_pearsonr = list()
    fit_pearsonp = list()
    num_peaks = list()

    # Get the percentage multiplier
    perc_mult = 100.0 / len(mlat_bins)

    # Loop through each MLT bin
    for ilt, mlt in enumerate(mlt_bins):
        ilats = np.where(np.isfinite(mean_intensity[:, ilt])
                         & (num_intensity[:, ilt] >= min_num)
                         & (mean_intensity[:, ilt] >= min_intensity))[0]

        if len(ilats) * perc_mult >= min_lat_perc:
            # Get the peak information to set the Gaussian fitting parameters
            params, ipeaks = get_fitting_params(mlat_bins, ilats, ilt,
                                                mean_intensity, num_gauss)

            # Find the desired Gaussian + Quadratic fit using least-squares
            # if there are enough latitudes to provide a fit
            lsq_args = (mlat_bins[ilats], mean_intensity[ilats, ilt],
                        1.0 / std_intensity[ilats, ilt])
            num_peaks.append(len(ipeaks))
            if len(params) <= len(ilats):
                lsq_result = leastsq(gauss_quad_err, params, args=lsq_args,
                                     full_output=True)

                # Evaluate the least squares output and save the results
                if lsq_result[-1] in [1, 2, 3, 4]:
                    gauss_out = distributions.mult_gauss_quad(
                        mlat_bins[ilats], lsq_result[0])
                    fmask = np.isfinite(gauss_out)
                    if fmask.sum() > 3:
                        pres = stats.pearsonr(
                            mean_intensity[ilats, ilt][fmask],
                            gauss_out[fmask])
                        func_params.append(lsq_result[0])

                        # SciPy versions below 1.9.0 will not work
                        fit_pearsonr.append(pres.statistic)
                        fit_pearsonp.append(pres.pvalue)

                        fit_cov.append(lsq_result[1])
                        found_fit = True
                    else:
                        found_fit = False
                else:
                    found_fit = False
            else:
                found_fit = False
                lsq_result = ["insufficient mlat coverage for expected peaks",
                              -1]

            if not found_fit:
                while not found_fit and len(ipeaks) > 1:
                    # Try reducing the number of peaks to obtain a valid fit
                    ipeaks.pop()
                    params = list(np.array(params)[:-3])

                    if len(params) <= len(ilats):
                        lsq_result = leastsq(gauss_quad_err, params,
                                             args=lsq_args, full_output=True)
                        num_peaks[-1] = len(ipeaks)
                        if lsq_result[-1] in [1, 2, 3, 4]:
                            gauss_out = distributions.mult_gauss_quad(
                                mlat_bins[ilats], lsq_result[0])
                            gmask = np.isfinite(gauss_out)

                            if sum(gmask) > 3:
                                pres = stats.pearsonr(
                                    mean_intensity[ilats, ilt][gmask],
                                    gauss_out[gmask])
                                found_fit = True
                                func_params.append(lsq_result[0])
                                fit_pearsonr.append(pres.statistic)
                                fit_pearsonp.append(pres.pvalue)
                                fit_cov.append(lsq_result[1])
                if not found_fit:
                    func_params.append(lsq_result[-2])
                    fit_cov.append(None)
                    fit_pearsonr.append(None)
                    fit_pearsonp.append(None)
        else:
            func_params.append("insufficient mlat coverage ({:} < {:})".format(
                len(ilats), min_lat_perc / perc_mult))
            fit_cov.append(None)
            fit_pearsonr.append(None)
            fit_pearsonp.append(None)
            num_peaks.append(0)

    return func_params, fit_cov, fit_pearsonr, fit_pearsonp, num_peaks

1	#!/usr/bin/env python
2	# -- coding: utf-8 --
3	# DOI: 10.5281/zenodo.15102100
4	# Full license can be found in License.md
5	#
6	# DISTRIBUTION STATEMENT A: Approved for public release. Distribution is
7	# unlimited.
8	# -----------------------------------------------------------------------------
9	"""Functions used for fitting functions to intensity data."""
10
11	import numpy as np	14✔
12	from scipy.optimize import leastsq	14✔
13	from scipy import stats	14✔
14
15	from pyIntensityFeatures.utils import checks	14✔
16	from pyIntensityFeatures.utils import distributions	14✔
17
18
19	def estimate_peak_widths(data, data_loc, peak_inds, peak_mags):	14✔
20	"""Estimate the full width at half max of the peaks in a curve.
21
22	Parameters
23	----------
24	data : array-like
25	Input data (y-axis)
26	data_loc : array-like
27	Input data location (x-axis)
28	peak_inds : array-like
29	Indexes of the peaks
30	peak_mags : array-like
31	Magnitude of the peaks
32
33	Returns
34	-------
35	peak_sigmas : list-like
36	FWHM values for each peak, or the closest estimate possible
37
38	Notes
39	-----
40	Differs from Longden, et al. (2010) function by not assuming a fixed data
41	location increment of 1.0 degrees.
42
43	"""
44	# Ensure the input is array-like
45	data = np.asarray(data)	14✔
46	data_loc = np.asarray(data_loc)	14✔
47	peak_inds = np.asarray(peak_inds)	14✔
48	peak_mags = np.asarray(peak_mags)	14✔
49
50	# Get the half-max for each peak
51	half_max = 0.5 * peak_mags	14✔
52	peak_sigmas = list()	14✔
53
54	for i, ipeak in enumerate(peak_inds):	14✔
55	# Step through all values of the curve data (`data`) with lower indices
56	# than this peak until either the half-max threshold or another peak is
57	# reached
58	reached = False	14✔
59	good_down = False	14✔
60	j = ipeak	14✔
61	while not reached and j > 0:	14✔
62	if (data[j] >= half_max[i]) and (data[j - 1] <= half_max[i]):	14✔
63	reached = True	14✔
64	fwhm_loc = data_loc[j - 1] + (	14✔
65	data_loc[j] - data_loc[j - 1]) * (
66	half_max[i] - data[j - 1]) / (data[j] - data[j - 1])
67	width_down = data_loc[ipeak] - fwhm_loc	14✔
68	good_down = True	14✔
69	else:
70	if data[j] < data[j - 1]:	14✔
71	# Data is no longer decreasing. Since the correct width was
72	# not reached before finding the desired half-max value,
73	# set an approximate value
74	reached = True	14✔
75	width_down = data_loc[ipeak] - data_loc[j - 1]	14✔
76	j -= 1	14✔
77
78	if not good_down and not reached:	14✔
79	# Reached edge of data before defining the width
80	width_down = data_loc[ipeak] - data_loc[j]	14✔
81
82	# Step through all values of the curve data (`data`) with higher
83	# indices than this peak until either the half-max threshold or another
84	# peak is reached
85	reached = False	14✔
86	good_up = False	14✔
87	j = ipeak	14✔
88	while not reached and j < len(data) - 1:	14✔
89	if (data[j] >= half_max[i]) and (data[j + 1] <= half_max[i]):	14✔
90	reached = True	14✔
91	fwhm_loc = data_loc[j + 1] + (	14✔
92	data_loc[j] - data_loc[j + 1]) * (
93	half_max[i] - data[j + 1]) / (data[j] - data[j + 1])
94	width_up = fwhm_loc - data_loc[ipeak]	14✔
95	good_up = True	14✔
96	else:
97	if data[j] < data[j + 1]:	14✔
98	# Data is no longer decreasing. Since the correct width was
99	# not reached before finding the desired half-max value,
100	# set an approximate value
101	reached = True	14✔
102	width_up = data_loc[j + 1] - data_loc[ipeak]	14✔
103	j += 1	14✔
104
105	if not good_up and not reached:	14✔
106	# Reached edge of data before defining the width
107	width_up = data_loc[j] - data_loc[ipeak]	14✔
108
109	# Calculate the peak sigmas using logic based on whether or not the
110	# HWHM was found on each leg
111	if good_up and good_down:	14✔
112	peak_sigmas.append((width_up + width_down) / checks.fwhm_const)	14✔
113	elif good_up:	14✔
114	peak_sigmas.append((2.0 * width_up) / checks.fwhm_const)	14✔
115	elif good_down:	14✔
116	peak_sigmas.append((2.0 * width_down) / checks.fwhm_const)	14✔
117	else:
118	peak_sigmas.append(2.0 * min([width_up, width_down])	14✔
119	/ checks.fwhm_const)
120
121	return peak_sigmas	14✔
122
123
124	def gauss_quad_err(params, xvals, yvals, weights):	14✔
125	"""Calculate the error of a quadriatic Gaussian fit.
126
127	Parameters
128	----------
129	params : set
130	Set of parameters used by `mult_gauss_quad`
131	xvals : float or array-like
132	Location values
133	yvals : float or array-like
134	Intensity values
135	weights : float or array-like
136	Weights for each data value
137
138	Returns
139	-------
140	err : float or array-like
141	Weighted error of the fit defined by `params`
142
143	See Also
144	--------
145	utils.distributions.mult_gauss_quad
146
147	"""
148	# Define the error function to fit
149	err = (distributions.mult_gauss_quad(xvals, params) - yvals) * weights	14✔
150
151	return err	14✔
152
153
154	def get_fitting_params(mlat_bins, ilats, ilt, mean_intensity, num_gauss=3):	14✔
155	"""Find the parameters for a Gaussian fit with a quadratic background.
156
157	Parameters
158	----------
159	mlat_bins : np.array
160	Magnetic latitude of output bin centres
161	ilats : np.array
162	Indices for magnetic latitudes with valid intensity values
163	ilt : int
164	Index for the magnetic local time
165	mean_intensity : np.array
166	2D array of mean intensity values
167	num_gauss : int
168	Maximum number of Gaussians to fit (default=3)
169
170	Returns
171	-------
172	params : list
173	List of the parameters needed to form a mult-Gaussian fit with a
174	quadratic background.
175	ipeaks : list
176	List containing indexes of the normalized intensity peaks for the
177	2D `mean_intesity` array down-selected by [`ilats`, `ilt`].
178
179	"""
180	# Perform a first-order polynomial fit to the data to obtain estimates for
181	# the background levels in the intensity profile. The first term is the
182	# slope and the second term is the intercept.
183	bg = np.polyfit(mlat_bins[ilats], mean_intensity[ilats, ilt], 1)	14✔
184
185	# Normalize the intensity curve
186	norm_intensity = (mean_intensity[ilats, ilt]	14✔
187	- mean_intensity[ilats, ilt].min()) / (
188	mean_intensity[ilats, ilt].max()
189	- mean_intensity[ilats, ilt].min())
190
191	# Find the main peak and its characteristics
192	ipeaks = [norm_intensity.argmax()]	14✔
193	peak_sigmas = estimate_peak_widths(norm_intensity, mlat_bins[ilats], ipeaks,	14✔
194	[norm_intensity[ipeaks[0]]])
195
196	# Locate any additional peaks, up to a desired maximum
197	mlt_gauss = num_gauss	14✔
198	while len(ipeaks) < mlt_gauss:	14✔
199	# Get a Gaussian for the prior peak
200	prior_gauss = distributions.gauss(	14✔
201	mlat_bins[ilats], norm_intensity[ipeaks[-1]],
202	mlat_bins[ilats][ipeaks[-1]], peak_sigmas[-1], 0.0)
203
204	# Remove the gaussian fit from the data
205	norm_intensity -= prior_gauss	14✔
206
207	# Get the next peak
208	inew = norm_intensity.argmax()	14✔
209
210	# Test to see if the new peak is significant
211	if not np.any((mlat_bins[ilats][inew] > mlat_bins[ilats][ipeaks]	14✔
212	- np.asarray(peak_sigmas) * 2.0)
213	& (mlat_bins[ilats][inew] < mlat_bins[ilats][ipeaks]
214	+ np.asarray(peak_sigmas) * 2.0)):
215	ipeaks.append(inew)	14✔
216	peak_sigmas.extend(estimate_peak_widths(	14✔
217	norm_intensity, mlat_bins[ilats], [ipeaks[-1]],
218	[norm_intensity[ipeaks[-1]]]))
219	else:
220	mlt_gauss = len(ipeaks)	14✔
221
222	# Use the peak information to set the Gaussian fitting parameters
223	params = [bg[1], bg[0], 0.0]	14✔
224	for i, ipeak in enumerate(ipeaks):	14✔
225	params.extend([mean_intensity[ilats, ilt][ipeak],	14✔
226	mlat_bins[ilats][ipeak], peak_sigmas[i]])
227
228	return params, ipeaks	14✔
229
230
231	def get_gaussian_func_fit(mlat_bins, mlt_bins, mean_intensity, std_intensity,	14✔
232	num_intensity, num_gauss=3, min_num=3,
233	min_intensity=0.0, min_lat_perc=70.0):
234	"""Fit intensity data using least-squares minimization in MLT bins.
235
236	Parameters
237	----------
238	mlat_bins : np.array
239	Magnetic latitude of output bin centres
240	mlt_bins : np.array
241	Magnetic local time of output bin centres
242	mean_intensity : np.array
243	2D array of mean intensity values
244	std_intensity : np.array
245	2D array of intensity standard deviations
246	num_intensity : np.array
247	2D array of intensity counts per bin
248	num_gauss : int
249	Maximum number of Gaussians to fit (default=3)
250	min_num : int
251	Minimum number of samples contributing to the intensity mean
252	(default=3)
253	min_intensity : float
254	Minimum intensity magnitude in Rayleighs (default=0.0)
255	min_lat_perc : int
256	Minimum percentage of latitude bins needed to define an intensity
257	function (default=70.0)
258
259	Returns
260	-------
261	func_params : list-like
262	Set of parameters used by `mult_gauss_quad` for each MLT bin that
263	defines the mean intensity as a function of magnetic latitude or set
264	to a string output with a reason if no successful fit was performed.
265	fit_cov : list-like
266	Covariance matrix for each fit
267	fit_pearsonr : list-like
268	Pearson r-value for correlation between the observations and the fit
269	fit_pearsonp : list-like
270	Pearson p-value for correlation between the observations and the fit
271	num_peaks : list-like
272	Number of peaks used in the fit
273
274	See Also
275	--------
276	utils.distributions.mult_gauss_quad
277
278	"""
279	# Initialize the output
280	func_params = list()	14✔
281	fit_cov = list()	14✔
282	fit_pearsonr = list()	14✔
283	fit_pearsonp = list()	14✔
284	num_peaks = list()	14✔
285
286	# Get the percentage multiplier
287	perc_mult = 100.0 / len(mlat_bins)	14✔
288
289	# Loop through each MLT bin
290	for ilt, mlt in enumerate(mlt_bins):	14✔
291	ilats = np.where(np.isfinite(mean_intensity[:, ilt])	14✔
292	& (num_intensity[:, ilt] >= min_num)
293	& (mean_intensity[:, ilt] >= min_intensity))[0]
294
295	if len(ilats) * perc_mult >= min_lat_perc:	14✔
296	# Get the peak information to set the Gaussian fitting parameters
297	params, ipeaks = get_fitting_params(mlat_bins, ilats, ilt,	14✔
298	mean_intensity, num_gauss)
299
300	# Find the desired Gaussian + Quadratic fit using least-squares
301	# if there are enough latitudes to provide a fit
302	lsq_args = (mlat_bins[ilats], mean_intensity[ilats, ilt],	14✔
303	1.0 / std_intensity[ilats, ilt])
304	num_peaks.append(len(ipeaks))	14✔
305	if len(params) <= len(ilats):	14✔
306	lsq_result = leastsq(gauss_quad_err, params, args=lsq_args,	14✔
307	full_output=True)
308
309	# Evaluate the least squares output and save the results
310	if lsq_result[-1] in [1, 2, 3, 4]:	14✔
311	gauss_out = distributions.mult_gauss_quad(	14✔
312	mlat_bins[ilats], lsq_result[0])
313	fmask = np.isfinite(gauss_out)	14✔
314	if fmask.sum() > 3:	14✔
315	pres = stats.pearsonr(	14✔
316	mean_intensity[ilats, ilt][fmask],
317	gauss_out[fmask])
318	func_params.append(lsq_result[0])	14✔
319
320	# SciPy versions below 1.9.0 will not work
321	fit_pearsonr.append(pres.statistic)	14✔
322	fit_pearsonp.append(pres.pvalue)	14✔
323
324	fit_cov.append(lsq_result[1])	14✔
325	found_fit = True	14✔
326	else:
327	found_fit = False	×
328	else:
329	found_fit = False	14✔
330	else:
331	found_fit = False	14✔
332	lsq_result = ["insufficient mlat coverage for expected peaks",	14✔
333	-1]
334
335	if not found_fit:	14✔
336	while not found_fit and len(ipeaks) > 1:	14✔
337	# Try reducing the number of peaks to obtain a valid fit
338	ipeaks.pop()	14✔
339	params = list(np.array(params)[:-3])	14✔
340
341	if len(params) <= len(ilats):	14✔
342	lsq_result = leastsq(gauss_quad_err, params,	14✔
343	args=lsq_args, full_output=True)
344	num_peaks[-1] = len(ipeaks)	14✔
345	if lsq_result[-1] in [1, 2, 3, 4]:	14✔
346	gauss_out = distributions.mult_gauss_quad(	14✔
347	mlat_bins[ilats], lsq_result[0])
348	gmask = np.isfinite(gauss_out)	14✔
349
350	if sum(gmask) > 3:	14✔
351	pres = stats.pearsonr(	14✔
352	mean_intensity[ilats, ilt][gmask],
353	gauss_out[gmask])
354	found_fit = True	14✔
355	func_params.append(lsq_result[0])	14✔
356	fit_pearsonr.append(pres.statistic)	14✔
357	fit_pearsonp.append(pres.pvalue)	14✔
358	fit_cov.append(lsq_result[1])	14✔
359	if not found_fit:	14✔
360	func_params.append(lsq_result[-2])	14✔
361	fit_cov.append(None)	14✔
362	fit_pearsonr.append(None)	14✔
363	fit_pearsonp.append(None)	14✔
364	else:
365	func_params.append("insufficient mlat coverage ({:} < {:})".format(	14✔
366	len(ilats), min_lat_perc / perc_mult))
367	fit_cov.append(None)	14✔
368	fit_pearsonr.append(None)	14✔
369	fit_pearsonp.append(None)	14✔
370	num_peaks.append(0)	14✔
371
372	return func_params, fit_cov, fit_pearsonr, fit_pearsonp, num_peaks	14✔

aburrell / pyIntensityFeatures / 23167106439

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous