Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-42766: Update linearity fitting code to fit temperature dependence and other updates recommended by Pierre Astier #239

Merged
merged 21 commits into from
May 13, 2024

Conversation

erykoff
Copy link
Contributor

@erykoff erykoff commented Apr 18, 2024

No description provided.

@erykoff erykoff requested a review from Alex-Broughton April 18, 2024 15:17
Copy link
Contributor

@czwa czwa left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mostly docstring comments.

raise ValueError("Config requests fitting temperature coefficient for "
f"{self.config.splineFitTemperatureColumn} but this column "
"is not available in inputPtc.auxValues.")
tempValue = inputPtc.auxValues[self.config.splineFitTemperatureColumn]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I'm going to suggest a different variable name for this. Am I correct that this is a vector of temperatures?

# here, for proper residual computation.
# We must modify the inputOrdinate according to the
# nuisance terms in the linearity fit for the residual
# computationcode to work properly.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo

def estimate_p0(self):
"""Estimate initial fit parameters.

Returns
-------
p0 : `np.ndarray`
Parameter array, with spline values (one for each node) followed
by proportionality constants (one for each group).
by proportionality constants (one for each group), and one extra
for the offset O (if fitOffset was set to True).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And +2 if fit_weights=True, and +1 if fit_temperature=True, right?

pars : `np.ndarray`
Parameter array, with spline values (one for each node) followed
by proportionality constants (one for each group.)
by proportionality constants (one for each group.), followed by
(optionally) one offset O.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar comment here: I don't think these docstrings were updated after the fit parameters were extended.

@@ -1075,7 +1185,7 @@ def fit(self, p0, min_iter=3, max_iter=20, max_rejection_per_iteration=5, n_sigm
----------
p0 : `np.ndarray`
Initial fit parameters (one for each knot, followed by one for
each grouping).
each grouping, and optionally one for the offset.)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And here.

@@ -324,6 +373,11 @@ def _check_linearity_spline(self, do_pd_offsets=False, n_points=200):
self.input_dims,
).outputLinearizer

if do_weight_fit:
# These checks currently fail, and weight fitting is not
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a ticket to investigate/fix the weight fitting that could be added here?

splineFitWeightParsStart = pexConfig.ListField(
dtype=float,
doc="Starting parameters for weight fit, if doSplineFitWeights=True. "
"Parameters are such that sigma = sqrt(par[0]**2. + par[1]**2./mu).",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have questions about this, but after reading the rest of the changes (and the test comment), it seems like these aren't being fit as part of the standard pipeline code, so I think my questions aren't relevant at the moment.

Copy link
Contributor

@czwa czwa left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few comments as I worked my way through the fitting process, with only docstring changes I think.

"If doSplineFitWeights=False then these are used as-is; otherwise "
"they are used as the initial values for fitting these parameters.",
length=2,
default=[1.0, 0.0],
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Am I correct that these are essentially a fixed-error and a signal-to-noise based term?

@@ -365,12 +420,19 @@ def run(self, inputPtc, dummy, camera, inputDims,
else:
mask = inputPtc.expIdMask[ampName].copy()

if self.config.linearityType == "Spline" and temperatureValues is not None:
mask &= np.isfinite(temperatureValues)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these NaNs? Do we understand what's going on in these exposures?

# Sum[(S(mu_i) + mu_i)/(k_j * D_i) - 1]**2, where S(mu_i) is
# an Akima Spline function of mu_i, the observed flat-pair
# Sum[(S(mu_i) + mu_i - O)/(k_j * D_i) - 1]**2, where S(mu_i)
# is an Akima Spline function of mu_i, the observed flat-pair
# mean; D_j is the photo-diode measurement corresponding to
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

D_i?

# mean; D_j is the photo-diode measurement corresponding to
# that flat-pair; and k_j is a constant of proportionality
# which is over index j as it is allowed to
# be different based on different photodiode settings (e.g.
# CCOBCURR).
# CCOBCURR); and O is a constant offset to allow for light
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this O a single value for all pairs? I'll find out when I get to the fitting code, but I'm just checking if this needs any index as well.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And I agree now that this is O, with no indexing.

sigmad = np.nan
else:
sigmad = median_abs_deviation(residuals[finite]/inputOrdinate[finite], scale="normal")
linearizer.fitResidualsSigmaMad[ampName] = sigmad
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not clear on the logic here. Is this checking for NaN (in which case any NaN value will result in sigmad=NaN), or is this checking for infinite values (which I see scipy.stats.median_abs_deviation does handle).

else:
mu_corr = mu

numerator = mu_corr - spl.interpolate(mu_corr)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this line contradict the docstrings that list S(mu_i) + mu_i - O? It's also entirely possible I've confused myself.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I've convinced myself that I am not confused in this case.

@czwa
Copy link
Contributor

czwa commented May 13, 2024

I continue to be happy with approval, including the two most recent commits (on failed amp/length matching),

@erykoff erykoff merged commit 87dad42 into main May 13, 2024
2 checks passed
@erykoff erykoff deleted the tickets/DM-42766 branch May 13, 2024 22:45
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants