Skip to content

Commit

Permalink
Manually filter nans prior to fitting.
Browse files Browse the repository at this point in the history
The lmfit nan_policy needs all the variables to be aligned, but here we are
fitting with the "x" grid as the independent variable, but "z" has a different
shape.  But the nan_policy="omit" option in lmfit simply does a filtering prior
to doing any fit, so here we zero out any nans and ensure the weights are 0 for
these values.
  • Loading branch information
erykoff committed Jul 11, 2023
1 parent 55623ac commit fac7b07
Showing 1 changed file with 14 additions and 2 deletions.
16 changes: 14 additions & 2 deletions python/lsst/ip/diffim/dipoleFitTask.py
Original file line number Diff line number Diff line change
Expand Up @@ -580,6 +580,8 @@ def fitDipoleImpl(self, source, tol=1e-7, rel_weight=0.5,
subim = afwImage.MaskedImageF(self.diffim.getMaskedImage(), bbox=bbox, origin=afwImage.PARENT)

z = diArr = subim.image.array
# Make sure we don't overwrite buffers.
z = z.copy()
weights = 1. / subim.variance.array # get the weights (=1/variance)

if rel_weight > 0. and ((self.posImage is not None) or (self.negImage is not None)):
Expand Down Expand Up @@ -624,8 +626,8 @@ def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None

modelFunctor = dipoleModelFunctor # dipoleModel.makeModel does not work for now.
# Create the lmfit model (lmfit uses scipy 'leastsq' option by default - Levenberg-Marquardt)
# Note we can also tell it to drop missing values from the data.
gmod = lmfit.Model(modelFunctor, verbose=verbose, missing='drop')
# We have to (later) filter out the nans by hand in our input to gmod.fit().
gmod = lmfit.Model(modelFunctor, verbose=verbose)

# Add the constraints for centroids, fluxes.
# starting constraint - near centroid of footprint
Expand Down Expand Up @@ -748,6 +750,16 @@ def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None
if np.any(~mask):
weights[~mask] = 0.

# Filter out any nans, and make the weights 0.
nans = (np.isnan(z) | np.isnan(weights))
nNans = nans.sum()
if nNans > 0:
if nNans < len(z):
z[nans] = np.nanmedian(z)
else:
z[nans] = 0
weights[nans] = 0

# Note that although we can, we're not required to set initial values for params here,
# since we set their param_hint's above.
# Can add "method" param to not use 'leastsq' (==levenberg-marquardt), e.g. "method='nelder'"
Expand Down

0 comments on commit fac7b07

Please sign in to comment.