Skip to content

Commit

Permalink
Add option to skip rank update when minimizing.
Browse files Browse the repository at this point in the history 
Add jointcalConfig option to skip photometry rank update.
Add FitterBase support for skipping rank update.
Add decam photometry tests with and without rank update.

Lift out Hessian creation into anonymous function.
  • Loading branch information
@parejkoj
parejkoj committed Apr 5, 2018
1 parent 0df1c64 commit 4e396a9
Show file tree
Hide file tree
Showing 5 changed files with 114 additions and 30 deletions.
6 changes: 5 additions & 1 deletion include/lsst/jointcal/FitterBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,9 @@ class FitterBase {
* @param[in] whatToFit See child method assignIndices for valid string values.
* @param[in] nSigmaCut How many sigma to reject outliers at. Outlier
* rejection ignored for nSigmaCut=0.
* @param[in] doRankUpdate Use CholmodSimplicialLDLT2.update() to do a fast rank update after outlier
* removal; otherwise do a slower full recomputation of the matrix.
* Only matters if nSigmaCut != 0.
*
* @return Return code describing success/failure of fit.
*
Expand All @@ -59,7 +62,8 @@ class FitterBase {
* the second run with the same "whatToFit" will produce no change in
* the fitted parameters, if the calculations and indices are defined correctly.
*/
MinimizeResult minimize(std::string const &whatToFit, double nSigmaCut = 0);
MinimizeResult minimize(std::string const &whatToFit, double const nSigmaCut = 0,
bool const doRankUpdate = true);

/**
* Returns the chi2 for the current state.
Expand Down
2 changes: 1 addition & 1 deletion python/lsst/jointcal/fitter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ namespace {
void declareFitterBase(py::module &mod) {
py::class_<FitterBase, std::shared_ptr<FitterBase>> cls(mod, "FitterBase");

cls.def("minimize", &FitterBase::minimize, "whatToFit"_a, "nSigRejCut"_a = 0);
cls.def("minimize", &FitterBase::minimize, "whatToFit"_a, "nSigRejCut"_a = 0, "doRankUpdate"_a = true);
cls.def("computeChi2", &FitterBase::computeChi2);
cls.def("saveChi2Contributions", &FitterBase::saveChi2Contributions);
}
Expand Down
36 changes: 28 additions & 8 deletions python/lsst/jointcal/jointcal.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,22 @@ class JointcalConfig(pexConfig.Config):
dtype=int,
default=7,
)
photometryDoRankUpdate = pexConfig.Field(
doc="Do the rank update step during minimization. "
"Skipping this can help deal with models that are too non-linear.",
dtype=bool,
default=True,
)
maxPhotometrySteps = pexConfig.Field(
doc="Maximum number of minimize iterations to take when fitting photometry.",
dtype=int,
default=20,
)
maxAstrometrySteps = pexConfig.Field(
doc="Maximum number of minimize iterations to take when fitting photometry.",
dtype=int,
default=20,
)
astrometryRefObjLoader = pexConfig.ConfigurableField(
target=LoadIndexedReferenceObjectsTask,
doc="Reference object loader for astrometric fit",
Expand Down Expand Up @@ -574,7 +590,8 @@ def _fit_photometry(self, associations, dataName=None):
model.freezeErrorTransform()
self.log.debug("Photometry error scales are frozen.")

chi2 = self._iterate_fit(fit, model, 20, "photometry", "Model Fluxes")
chi2 = self._iterate_fit(fit, model, self.config.maxPhotometrySteps, "photometry", "Model Fluxes",
doRankUpdate=self.config.photometryDoRankUpdate)

add_measurement(self.job, 'jointcal.photometry_final_chi2', chi2.chi2)
add_measurement(self.job, 'jointcal.photometry_final_ndof', chi2.ndof)
Expand Down Expand Up @@ -646,27 +663,30 @@ def _fit_astrometry(self, associations, dataName=None):
raise FloatingPointError('Pre-iteration chi2 is invalid: %s'%chi2)
self.log.info("Fit prepared with %s", str(chi2))

chi2 = self._iterate_fit(fit, model, 20, "astrometry", "Distortions Positions")
chi2 = self._iterate_fit(fit, model, self.config.maxAstrometrySteps,
"astrometry", "Distortions Positions")

add_measurement(self.job, 'jointcal.astrometry_final_chi2', chi2.chi2)
add_measurement(self.job, 'jointcal.astrometry_final_ndof', chi2.ndof)

return Astrometry(fit, model, sky_to_tan_projection)

def _iterate_fit(self, fit, model, max_steps, name, whatToFit):
def _iterate_fit(self, fit, model, max_steps, name, whatToFit, doRankUpdate=True):
"""Run fit.minimize up to max_steps times, returning the final chi2."""

for i in range(max_steps):
r = fit.minimize(whatToFit, 5) # outlier removal at 5 sigma.
# outlier removal at 5 sigma.
r = fit.minimize(whatToFit, 5, doRankUpdate=doRankUpdate)
chi2 = fit.computeChi2()
if not np.isfinite(chi2.chi2):
raise FloatingPointError('Fit iteration chi2 is invalid: %s'%chi2)
self.log.info(str(chi2))
if r == MinimizeResult.Converged:
self.log.debug("fit has converged - no more outliers - redo minimization "
"one more time in case we have lost accuracy in rank update.")
# Redo minimization one more time in case we have lost accuracy in rank update
r = fit.minimize(whatToFit, 5) # outliers removal at 5 sigma.
if doRankUpdate:
self.log.debug("fit has converged - no more outliers - redo minimization "
"one more time in case we have lost accuracy in rank update.")
# Redo minimization one more time in case we have lost accuracy in rank update
r = fit.minimize(whatToFit, 5) # outliers removal at 5 sigma.
chi2 = fit.computeChi2()
self.log.info("Fit completed with: %s", str(chi2))
break
Expand Down
56 changes: 37 additions & 19 deletions src/FitterBase.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,15 @@ unsigned FitterBase::findOutliers(double nSigmaCut, MeasuredStarList &msOutliers
return nOutliers;
}

MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaCut) {
namespace {
SparseMatrixD createHessian(int nParTot, TripletList const &tripletList) {
SparseMatrixD jacobian(nParTot, tripletList.getNextFreeIndex());
jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
return jacobian * jacobian.transpose();
}
} // namespace

MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaCut, bool doRankUpdate) {
assignIndices(whatToFit);

MinimizeResult returnCode = MinimizeResult::Converged;
Expand All @@ -133,14 +141,7 @@ MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaC

LOGLS_DEBUG(_log, "End of triplet filling, ntrip = " << tripletList.size());

SparseMatrixD hessian;
{
SparseMatrixD jacobian(_nParTot, tripletList.getNextFreeIndex());
jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
// release memory shrink_to_fit is C++11
tripletList.clear(); // tripletList.shrink_to_fit();
hessian = jacobian * jacobian.transpose();
} // release the Jacobian
SparseMatrixD hessian = createHessian(_nParTot, tripletList);

LOGLS_DEBUG(_log, "Starting factorization, hessian: dim="
<< hessian.rows() << " non-zeros=" << hessian.nonZeros()
Expand Down Expand Up @@ -176,21 +177,38 @@ MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaC
totalMeasOutliers += msOutliers.size();
totalRefOutliers += fsOutliers.size();
if (nOutliers == 0) break;
TripletList tripletList(nOutliers);
TripletList outlierTriplets(nOutliers);
grad.setZero(); // recycle the gradient
// compute the contributions of outliers to derivatives
outliersContributions(msOutliers, fsOutliers, tripletList, grad);
outliersContributions(msOutliers, fsOutliers, outlierTriplets, grad);
// Remove significant outliers
removeMeasOutliers(msOutliers);
removeRefOutliers(fsOutliers);
// convert triplet list to eigen internal format
SparseMatrixD H(_nParTot, tripletList.getNextFreeIndex());
H.setFromTriplets(tripletList.begin(), tripletList.end());
int update_status = chol.update(H, false /* means downdate */);
LOGLS_DEBUG(_log, "cholmod update_status " << update_status);
// The contribution of outliers to the gradient is the opposite
// of the contribution of all other terms, because they add up to 0
grad *= -1;
if (doRankUpdate) {
// convert triplet list to eigen internal format
SparseMatrixD H(_nParTot, outlierTriplets.getNextFreeIndex());
H.setFromTriplets(outlierTriplets.begin(), outlierTriplets.end());
int update_status = chol.update(H, false /* means downdate */);
LOGLS_DEBUG(_log, "cholmod update_status " << update_status);
// The contribution of outliers to the gradient is the opposite
// of the contribution of all other terms, because they add up to 0
grad *= -1;
} else {
// Rebuild the matrix and gradient
leastSquareDerivatives(tripletList, grad);
_lastNTrip = tripletList.size();
LOGLS_DEBUG(_log, "Triplets recomputed, ntrip = " << tripletList.size());
hessian = createHessian(_nParTot, tripletList);
LOGLS_DEBUG(_log,
"Restarting factorization, hessian: dim="
<< hessian.rows() << " non-zeros=" << hessian.nonZeros()
<< " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
chol.compute(hessian);
if (chol.info() != Eigen::Success) {
LOGLS_ERROR(_log, "minimize: factorization failed ");
return MinimizeResult::Failed;
}
}
}

// only print the outlier summary if outlier rejection was turned on.
Expand Down
44 changes: 43 additions & 1 deletion tests/test_jointcal_decam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,8 @@ def setUp(self):
input_dir=input_dir,
all_visits=all_visits,
other_args=other_args,
do_plot=do_plot)
do_plot=do_plot,
log_level="DEBUG")

def test_jointcalTask_2_visits(self):
self.config = lsst.jointcal.jointcal.JointcalConfig()
Expand Down Expand Up @@ -107,6 +108,47 @@ def test_jointcalTask_2_visits_constrainedPoly(self):

self._testJointcalTask(2, relative_error, self.dist_rms_absolute, pa1, metrics=metrics)

def test_jointcalTask_2_visits_constrainedPhotometry(self):
self.config = lsst.jointcal.jointcal.JointcalConfig()
self.config.astrometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
self.config.photometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
self.config.photometryModel = "constrained"
self.config.doAstrometry = False
self.jointcalStatistics.do_astrometry = False

pa1 = 0.15
metrics = {'collected_photometry_refStars': 8189,
'selected_photometry_refStars': 773,
'associated_photometry_fittedStars': 8241,
'selected_photometry_fittedStars': 2261,
'selected_photometry_ccdImages': 17,
'photometry_final_chi2': 3487.39,
'photometry_final_ndof': 2190,
}

self._testJointcalTask(2, None, None, pa1, metrics=metrics)

def test_jointcalTask_2_visits_constrainedPhotometry_no_rank_update(self):
self.config = lsst.jointcal.jointcal.JointcalConfig()
self.config.astrometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
self.config.photometryRefObjLoader.retarget(LoadAstrometryNetObjectsTask)
self.config.photometryModel = "constrained"
self.config.photometryDoRankUpdate = False
self.config.doAstrometry = False
self.jointcalStatistics.do_astrometry = False

pa1 = 0.15
metrics = {'collected_photometry_refStars': 8189,
'selected_photometry_refStars': 773,
'associated_photometry_fittedStars': 8241,
'selected_photometry_fittedStars': 2261,
'selected_photometry_ccdImages': 17,
'photometry_final_chi2': 3522.10,
'photometry_final_ndof': 2197,
}

self._testJointcalTask(2, None, None, pa1, metrics=metrics)


class MemoryTester(lsst.utils.tests.MemoryTestCase):
pass
Expand Down

0 comments on commit 4e396a9

Please sign in to comment.