lsst::jointcal::PhotometryFit Class Reference

Class that handles the photometric least squares problem. More...

#include <PhotometryFit.h>

Inheritance diagram for lsst::jointcal::PhotometryFit:

Public Member Functions
	PhotometryFit (std::shared_ptr< Associations > associations, std::shared_ptr< PhotometryModel > photometryModel)
	Construct a photometry fitter.
	PhotometryFit (PhotometryFit const &)=delete
	No copy or move: there is only ever one fitter of a given type.
	PhotometryFit (PhotometryFit &&)=delete
PhotometryFit &	operator= (PhotometryFit const &)=delete
PhotometryFit &	operator= (PhotometryFit &&)=delete
void	assignIndices (std::string const &whatToFit) override
	Set parameters to fit and assign indices in the big matrix.
void	offsetParams (Eigen::VectorXd const &delta) override
	Offset the parameters by the requested quantities.
std::shared_ptr< PhotometryModel >	getModel () const
	Return the model being fit.
MinimizeResult	minimize (std::string const &whatToFit, double nSigmaCut=0, double sigmaRelativeTolerance=0, bool doRankUpdate=true, bool doLineSearch=false, std::string const &dumpMatrixFile="")
	Does a 1 step minimization, assuming a linear model.
Chi2Statistic	computeChi2 () const
	Returns the chi2 for the current state.
void	leastSquareDerivatives (TripletList &tripletList, Eigen::VectorXd &grad) const
	Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.
virtual void	saveChi2Contributions (std::string const &baseName) const
	Save the full chi2 term per star that was used in the minimization, for debugging.

Protected Member Functions
void	saveChi2MeasContributions (std::string const &filename) const override
	Save a CSV file containing residuals of measurement terms.
void	saveChi2RefContributions (std::string const &filename) const override
	Save a CSV file containing residuals of reference terms.
std::size_t	findOutliers (double nSigmaCut, MeasuredStarList &msOutliers, FittedStarList &fsOutliers, double &cut) const
	Find Measurements and references contributing more than a cut, computed as.
void	outliersContributions (MeasuredStarList &msOutliers, FittedStarList &fsOutliers, TripletList &tripletList, Eigen::VectorXd &grad)
	Contributions to derivatives from (presumably) outlier terms.
void	removeMeasOutliers (MeasuredStarList &outliers)
	Remove measuredStar outliers from the fit. No Refit done.
void	removeRefOutliers (FittedStarList &outliers)
	Remove refStar outliers from the fit. No Refit done.

Protected Attributes
std::shared_ptr< Associations >	_associations
std::string	_whatToFit
Eigen::Index	_lastNTrip
Eigen::Index	_nTotal
Eigen::Index	_nModelParams
Eigen::Index	_nStarParams
LOG_LOGGER	_log

Private Member Functions
void	accumulateStatImageList (CcdImageList const &ccdImageList, Chi2Accumulator &accum) const override
	Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.
void	accumulateStatRefStars (Chi2Accumulator &accum) const override
	Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.
void	getIndicesOfMeasuredStar (MeasuredStar const &measuredStar, IndexVector &indices) const override
	this routine is to be used only in the framework of outlier removal
void	leastSquareDerivativesMeasurement (CcdImage const &ccdImage, TripletList &tripletList, Eigen::VectorXd &grad, MeasuredStarList const *measuredStarList=nullptr) const override
	Compute the derivatives of the measured stars and model for one CcdImage.
void	leastSquareDerivativesReference (FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const override
	Compute the derivatives of the reference terms.

Detailed Description

Class that handles the photometric least squares problem.

Definition at line 46 of file PhotometryFit.h.

Constructor & Destructor Documentation

PhotometryFit() [1/3]

lsst::jointcal::PhotometryFit::PhotometryFit ( std::shared_ptr< Associations > associations, std::shared_ptr< PhotometryModel > photometryModel )

inline

Construct a photometry fitter.

Parameters

associations	The associations catalog to use in the fitter.
photometryModel	The model to build the fitter for.

Definition at line 54 of file PhotometryFit.h.

            : FitterBase(std::move(associations)),
              _fittingModel(false),
              _fittingFluxes(false),
              _photometryModel(std::move(photometryModel)) {
        _log = LOG_GET("lsst.jointcal.PhotometryFit");
    }

PhotometryFit() [2/3]

lsst::jointcal::PhotometryFit::PhotometryFit ( PhotometryFit const & )

delete

No copy or move: there is only ever one fitter of a given type.

PhotometryFit() [3/3]

lsst::jointcal::PhotometryFit::PhotometryFit ( PhotometryFit && )

delete

Member Function Documentation

accumulateStatImageList()

void lsst::jointcal::PhotometryFit::accumulateStatImageList ( CcdImageList const & ccdImageList, Chi2Accumulator & accum ) const

overrideprivatevirtual

Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.

Note

the math in this method and leastSquareDerivativesMeasurement() must be kept consistent, in terms of +/- convention, definition of model, etc.

Implements lsst::jointcal::FitterBase.

Definition at line 134 of file PhotometryFit.cc.

                                                                                                          {
    /**********************************************************************/
    /**********************************************************************/
    for (auto const &ccdImage : ccdImageList) {
        auto &catalog = ccdImage->getCatalogForFit();
 
        for (auto const &measuredStar : catalog) {
            if (!measuredStar->isValid()) continue;
            double sigma = _photometryModel->transformError(*ccdImage, *measuredStar);
            double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
 
            double chi2Val = std::pow(residual / sigma, 2);
            accum.addEntry(chi2Val, 1, measuredStar);
        }  // end loop on measurements
    }
}

accumulateStatRefStars()

void lsst::jointcal::PhotometryFit::accumulateStatRefStars ( Chi2Accumulator & accum ) const

overrideprivatevirtual

Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.

Note

the math in this method and leastSquareDerivativesReference() must be kept consistent, in terms of +/- convention, definition of model, etc.

Implements lsst::jointcal::FitterBase.

Definition at line 153 of file PhotometryFit.cc.

                                                                       {
    /**********************************************************************/
    /**********************************************************************/
 
    FittedStarList &fittedStarList = _associations->fittedStarList;
    for (auto const &fittedStar : fittedStarList) {
        auto refStar = fittedStar->getRefStar();
        if (refStar == nullptr) continue;
        double sigma = _photometryModel->getRefError(*refStar);
        double residual = _photometryModel->computeRefResidual(*fittedStar, *refStar);
        double chi2 = std::pow(residual / sigma, 2);
        accum.addEntry(chi2, 1, fittedStar);
    }
}

assignIndices()

void lsst::jointcal::PhotometryFit::assignIndices ( std::string const & whatToFit )

overridevirtual

Set parameters to fit and assign indices in the big matrix.

Parameters

[in]

whatToFit

Valid strings : "Model", "Fluxes", which define which parameter sets are going to be fitted. whatToFit="Model Fluxes" will set both parameter sets variable when computing derivatives. Provided it contains "Model", whatToFit is passed over to the PhotometryModel, and can hence be used to control more finely which subsets of the photometric model are being fitted, if the the actual PhotometryModel implements such a possibility.

Implements lsst::jointcal::FitterBase.

Definition at line 185 of file PhotometryFit.cc.

                                                            {
    _whatToFit = whatToFit;
    LOGLS_INFO(_log, "assignIndices: now fitting: " << whatToFit);
    _fittingModel = (_whatToFit.find("Model") != std::string::npos);
    _fittingFluxes = (_whatToFit.find("Fluxes") != std::string::npos);
    // When entering here, we assume that whatToFit has already been interpreted.
 
    _nModelParams = (_fittingModel) ? _photometryModel->assignIndices(whatToFit, 0) : 0;
    std::size_t ipar = _nModelParams;
 
    if (_fittingFluxes) {
        for (auto &fittedStar : _associations->fittedStarList) {
            // the parameter layout here is used also
            // - when filling the derivatives
            // - when updating (offsetParams())
            // - in getIndicesOfMeasuredStar
            fittedStar->setIndexInMatrix(ipar);
            ipar += 1;
        }
    }
    _nStarParams = ipar - _nModelParams;
    _nTotal = ipar;
    LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
                                            << " fluxes: " << _nStarParams);
}

computeChi2()

Chi2Statistic lsst::jointcal::FitterBase::computeChi2 ( ) const

inherited

Returns the chi2 for the current state.

Definition at line 43 of file FitterBase.cc.

                                            {
    Chi2Statistic chi2;
    accumulateStatImageList(_associations->getCcdImageList(), chi2);
    accumulateStatRefStars(chi2);
    // chi2.ndof contains the number of squares.
    // So subtract the number of parameters.
    chi2.ndof -= _nTotal;
    return chi2;
}

findOutliers()

std::size_t lsst::jointcal::FitterBase::findOutliers ( double nSigmaCut, MeasuredStarList & msOutliers, FittedStarList & fsOutliers, double & cut ) const

protectedinherited

Find Measurements and references contributing more than a cut, computed as.

\[ <chi2> + nSigmaCut + rms(chi2). \]

The outliers are NOT removed, and no refit is done.

After returning from here, there are still measurements that contribute above the cut, but their contribution should be evaluated after a refit before discarding them.

Parameters

[in]	nSigmaCut	Number of sigma to select on.
[out]	msOutliers	list of MeasuredStar outliers to populate
[out]	fsOutliers	list of FittedStar outliers to populate
[out]	cut	value of chi2 that defines which objects are outliers

Returns

Total number of outliers that were removed.

Definition at line 53 of file FitterBase.cc.

                                                                                    {
    // collect chi2 contributions
    Chi2List chi2List;
    chi2List.reserve(_associations->getMaxMeasuredStars() + _associations->refStarList.size());
    // contributions from measurement terms:
    accumulateStatImageList(_associations->ccdImageList, chi2List);
    // and from reference terms
    accumulateStatRefStars(chi2List);
 
    // compute some statistics
    size_t nval = chi2List.size();
    if (nval == 0) return 0;
    sort(chi2List.begin(), chi2List.end());
    double median = (nval & 1) ? chi2List[nval / 2].chi2
                               : 0.5 * (chi2List[nval / 2 - 1].chi2 + chi2List[nval / 2].chi2);
    auto averageAndSigma = chi2List.computeAverageAndSigma();
    LOGLS_DEBUG(_log, "findOutliers chi2 stat: mean/median/sigma " << averageAndSigma.first << '/' << median
                                                                   << '/' << averageAndSigma.second);
    cut = averageAndSigma.first + nSigmaCut * averageAndSigma.second;
    /* For each of the parameters, we will not remove more than 1
       measurement that contributes to constraining it. Keep track using
       of what we are touching using an integer vector. This is the
       trick that Marc Betoule came up to for outlier removals in "star
       flats" fits. */
    Eigen::VectorXi affectedParams(_nTotal);
    affectedParams.setZero();
 
    std::size_t nOutliers = 0;  // returned to the caller
    // start from the strongest outliers.
    for (auto chi2 = chi2List.rbegin(); chi2 != chi2List.rend(); ++chi2) {
        if (chi2->chi2 < cut) break;  // because the array is sorted.
        IndexVector indices;
        /* now, we want to get the indices of the parameters this chi2
           term depends on. We have to figure out which kind of term it
           is; we use for that the type of the star attached to the Chi2Star. */
        auto measuredStar = std::dynamic_pointer_cast<MeasuredStar>(chi2->star);
        std::shared_ptr<FittedStar> fittedStar;  // To add to fsOutliers if it is a reference outlier.
        if (measuredStar == nullptr) {
            // it is a reference outlier
            fittedStar = std::dynamic_pointer_cast<FittedStar>(chi2->star);
            if (fittedStar->getMeasurementCount() == 0) {
                LOGLS_WARN(_log, "FittedStar with no measuredStars found as an outlier: "
                                         << *fittedStar << " chi2: " << chi2->chi2);
                continue;
            }
            if (_nStarParams == 0) {
                LOGLS_TRACE(_log,
                            "RefStar is outlier but not removed when not fitting FittedStar-RefStar values: "
                                    << *(fittedStar->getRefStar()) << " chi2: " << chi2->chi2);
                continue;
            }
            // NOTE: Stars contribute twice to astrometry (x,y), but once to photometry (flux),
            // NOTE: but we only need to mark one index here because both will be removed with that star.
            indices.push_back(fittedStar->getIndexInMatrix());
            LOGLS_TRACE(_log, "Removing refStar " << *(fittedStar->getRefStar()) << " chi2: " << chi2->chi2);
            /* One might think it would be useful to account for PM
               parameters here, but it is just useless */
        } else {
            // it is a measurement outlier
            auto tempFittedStar = measuredStar->getFittedStar();
            if (tempFittedStar->getMeasurementCount() == 1 && tempFittedStar->getRefStar() == nullptr) {
                LOGLS_WARN(_log, "FittedStar with 1 measuredStar and no refStar found as an outlier: "
                                         << *tempFittedStar);
                continue;
            }
            getIndicesOfMeasuredStar(*measuredStar, indices);
            LOGLS_TRACE(_log, "Removing measStar " << *measuredStar << " chi2: " << chi2->chi2);
        }
 
        /* Find out if we already discarded a stronger outlier
        constraining some parameter this one constrains as well. If
         yes, we keep this one, because this stronger outlier could be
         causing the large chi2 we have in hand.  */
        bool drop_it = true;
        for (auto const &i : indices) {
            if (affectedParams(i) != 0) {
                drop_it = false;
            }
        }
 
        if (drop_it)  // store the outlier in one of the lists:
        {
            if (measuredStar == nullptr) {
                // reference term
                fsOutliers.push_back(fittedStar);
            } else {
                // measurement term
                msOutliers.push_back(measuredStar);
            }
            // mark the parameters as directly changed when we discard this chi2 term.
            for (auto const &i : indices) {
                affectedParams(i)++;
            }
            nOutliers++;
        }
    }  // end loop on measurements/references
    LOGLS_DEBUG(_log, "findOutliers: found " << msOutliers.size() << " meas outliers and "
                                             << fsOutliers.size() << " ref outliers ");
 
    return nOutliers;
}

getIndicesOfMeasuredStar()

void lsst::jointcal::PhotometryFit::getIndicesOfMeasuredStar ( MeasuredStar const & measuredStar, IndexVector & indices ) const

overrideprivatevirtual

this routine is to be used only in the framework of outlier removal

it fills the array of indices of parameters that a Measured star constrains. Not really all of them if you check.

Implements lsst::jointcal::FitterBase.

Definition at line 173 of file PhotometryFit.cc.

                                                                                                         {
    indices.clear();
    if (_fittingModel) {
        _photometryModel->getMappingIndices(measuredStar.getCcdImage(), indices);
    }
    if (_fittingFluxes) {
        std::shared_ptr<FittedStar const> const fs = measuredStar.getFittedStar();
        Eigen::Index fsIndex = fs->getIndexInMatrix();
        indices.push_back(fsIndex);
    }
}

getModel()

std::shared_ptr< PhotometryModel > lsst::jointcal::PhotometryFit::getModel ( ) const

inline

Return the model being fit.

Definition at line 87 of file PhotometryFit.h.

{ return _photometryModel; }

leastSquareDerivatives()

void lsst::jointcal::FitterBase::leastSquareDerivatives ( TripletList & tripletList, Eigen::VectorXd & grad ) const

inherited

Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.

The Jacobian is given as triplets in a sparse matrix, the gradient as a dense vector. The parameters which vary, and their indices, are to be set using assignIndices.

Parameters

tripletList	tripletList of (row,col,value) representing the Jacobian of the chi2.
grad	The gradient of the chi2.

Definition at line 347 of file FitterBase.cc.

                                                                                           {
    auto ccdImageList = _associations->getCcdImageList();
    for (auto const &ccdImage : ccdImageList) {
        leastSquareDerivativesMeasurement(*ccdImage, tripletList, grad);
    }
    leastSquareDerivativesReference(_associations->fittedStarList, tripletList, grad);
}

leastSquareDerivativesMeasurement()

void lsst::jointcal::PhotometryFit::leastSquareDerivativesMeasurement ( CcdImage const & ccdImage, TripletList & tripletList, Eigen::VectorXd & grad, MeasuredStarList const * measuredStarList = nullptr ) const

overrideprivatevirtual

Compute the derivatives of the measured stars and model for one CcdImage.

The last argument will process a sub-list for outlier removal.

Implements lsst::jointcal::FitterBase.

Definition at line 45 of file PhotometryFit.cc.

                                                                                                      {
    /**********************************************************************/
    /* @note the math in this method and accumulateStatImageList() must be kept consistent,
     * in terms of +/- convention, definition of model, etc. */
    /**********************************************************************/
 
    /* this routine works in two different ways: either providing the
       Ccd, of providing the MeasuredStarList. In the latter case, the
       Ccd should match the one(s) in the list. */
    if (measuredStarList) assert(&(measuredStarList->front()->getCcdImage()) == &ccdImage);
 
    std::size_t nparModel = (_fittingModel) ? _photometryModel->getNpar(ccdImage) : 0;
    std::size_t nparFlux = (_fittingFluxes) ? 1 : 0;
    std::size_t nparTotal = nparModel + nparFlux;
    IndexVector indices(nparModel, -1);
    if (_fittingModel) _photometryModel->getMappingIndices(ccdImage, indices);
 
    Eigen::VectorXd H(nparTotal);  // derivative matrix
    // current position in the Jacobian
    Eigen::Index kTriplets = tripletList.getNextFreeIndex();
    const MeasuredStarList &catalog = (measuredStarList) ? *measuredStarList : ccdImage.getCatalogForFit();
 
    for (auto const &measuredStar : catalog) {
        if (!measuredStar->isValid()) continue;
        H.setZero();  // we cannot be sure that all entries will be overwritten.
 
        double residual = _photometryModel->computeResidual(ccdImage, *measuredStar);
        double inverseSigma = 1.0 / _photometryModel->transformError(ccdImage, *measuredStar);
        double W = std::pow(inverseSigma, 2);
 
        if (_fittingModel) {
            _photometryModel->computeParameterDerivatives(*measuredStar, ccdImage, H);
            for (std::size_t k = 0; k < indices.size(); k++) {
                Eigen::Index l = indices[k];
                tripletList.addTriplet(l, kTriplets, H[k] * inverseSigma);
                grad[l] += H[k] * W * residual;
            }
        }
        if (_fittingFluxes) {
            Eigen::Index index = measuredStar->getFittedStar()->getIndexInMatrix();
            // Note: H = dR/dFittedStarFlux == -1
            tripletList.addTriplet(index, kTriplets, -1.0 * inverseSigma);
            grad[index] += -1.0 * W * residual;
        }
        kTriplets += 1;  // each measurement contributes 1 column in the Jacobian
    }
 
    tripletList.setNextFreeIndex(kTriplets);
}

leastSquareDerivativesReference()

void lsst::jointcal::PhotometryFit::leastSquareDerivativesReference ( FittedStarList const & fittedStarList, TripletList & tripletList, Eigen::VectorXd & grad ) const

overrideprivatevirtual

Compute the derivatives of the reference terms.

Note

the math in this method and accumulateStatReference() must be kept consistent, in terms of +/- convention, definition of model, etc.

Implements lsst::jointcal::FitterBase.

Definition at line 97 of file PhotometryFit.cc.

                                                                                                         {
    /**********************************************************************/
    /**********************************************************************/
 
    // Derivatives of terms involving fitted and refstars only contribute if we are fitting fluxes.
    if (!_fittingFluxes) return;
    // Can't compute anything if there are no refStars.
    if (_associations->refStarList.empty()) return;
 
    Eigen::Index kTriplets = tripletList.getNextFreeIndex();
 
    for (auto const &fittedStar : fittedStarList) {
        auto refStar = fittedStar->getRefStar();
        if (refStar == nullptr) continue;  // no contribution if no associated refstar
        // TODO: Can we actually work with multiple filters at a time in this scheme?
        // TODO: I feel we might only be able to fit one filter at a time, because these terms are
        // independent of the ccdImage, so we don't have a specific filter.
        // filter = ccdImage.getFilter();
 
        // W == inverseSigma^2
 
        double inverseSigma = 1.0 / _photometryModel->getRefError(*refStar);
        // Residual is fittedStar - refStar for consistency with measurement terms.
        double residual = _photometryModel->computeRefResidual(*fittedStar, *refStar);
 
        Eigen::Index index = fittedStar->getIndexInMatrix();
        // Note: H = dR/dFittedStar == 1
        tripletList.addTriplet(index, kTriplets, 1.0 * inverseSigma);
        grad(index) += 1.0 * std::pow(inverseSigma, 2) * residual;
        kTriplets += 1;
    }
    tripletList.setNextFreeIndex(kTriplets);
}

minimize()

MinimizeResult lsst::jointcal::FitterBase::minimize ( std::string const & whatToFit, double nSigmaCut = 0, double sigmaRelativeTolerance = 0, bool doRankUpdate = true, bool doLineSearch = false, std::string const & dumpMatrixFile = "" )

inherited

Does a 1 step minimization, assuming a linear model.

This is a complete Newton Raphson step. Compute first and second derivatives, solve for the step and apply it, with an optional line search.

It calls assignIndices, leastSquareDerivatives, solves the linear system and calls offsetParams, then removes outliers in a loop if requested. Relies on sparse linear algebra via Eigen's CholmodSupport package.

Parameters

[in]	whatToFit	See child method assignIndices for valid string values.
[in]	nSigmaCut	How many sigma to reject outliers at. Outlier rejection ignored for nSigmaCut=0.
[in]	sigmaRelativeTolerance	Percentage change in the chi2 cut for outliers tolerated for termination. If value is zero, minimization iterations will continue until there are no outliers.
[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.
[in]	doLineSearch	Use boost's brent_find_minima to perform a line search after the gradient solution is found, and apply the scale factor to the computed offsets. The line search is done in the domain [-1, 2], but if the scale factor is far from 1.0, then the problem is likely in a significantly non-linear regime.
[in]	dumpMatrixFile	Write the pre-fit Hessian matrix and gradient to the files with "-mat.txt" and "-grad.txt". Be aware, this requires a large increase in memory usage to create a dense matrix before writing it; the output file may be large. Writing the matrix can be helpful for debugging bad fits. Read it and compute the real eigenvalues (recall that the Hessian is symmetric by construction) with numpy: hessian = np.matrix(np.loadtxt("dumpMatrixFile-mat.txt")) values, vectors = np.linalg.eigh(hessian)

Returns

Return code describing success/failure of fit.

Note

When fitting one parameter set by itself (e.g. "Model"), the system is purely linear (assuming there are no cross-terms in the derivatives, e.g. the SimpleAstrometryModel), which should result in the optimal chi2 after a single step. This can be used to debug the fitter by fitting that parameter set twice in a row: the second run with the same "whatToFit" will produce no change in the fitted parameters, if the calculations and indices are defined correctly.

Definition at line 179 of file FitterBase.cc.

                                                                     {
    assignIndices(whatToFit);
 
    MinimizeResult returnCode = MinimizeResult::Converged;
 
    // For the initial vector size, use all measured stars + all fitted stars, which should give the
    // maximum possible number of triplets.
    std::size_t nTrip = (_lastNTrip)
                                ? _lastNTrip
                                : _associations->getMaxMeasuredStars() + _associations->fittedStarList.size();
    TripletList tripletList(nTrip);
    Eigen::VectorXd grad(_nTotal);
    grad.setZero();
    double scale = 1.0;
 
    // Fill the triplets
    leastSquareDerivatives(tripletList, grad);
    _lastNTrip = tripletList.size();
 
    LOGLS_DEBUG(_log, "End of triplet filling, ntrip = " << tripletList.size());
 
    SparseMatrixD hessian = createHessian(_nTotal, tripletList);
    tripletList.clear();  // we don't need it any more after we have the hessian.
 
    LOGLS_DEBUG(_log, "Starting factorization, hessian: dim="
                              << hessian.rows() << " non-zeros=" << hessian.nonZeros()
                              << " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
 
    if (dumpMatrixFile != "") {
        if (hessian.rows() * hessian.cols() > 2e8) {
            LOGLS_WARN(_log, "Hessian matrix is too big to dump to file, with rows, columns: "
                                     << hessian.rows() << ", " << hessian.cols());
        } else {
            dumpMatrixAndGradient(hessian, grad, dumpMatrixFile, _log);
        }
    }
 
    CholmodSimplicialLDLT2<SparseMatrixD> chol(hessian);
    if (chol.info() != Eigen::Success) {
        LOGLS_ERROR(_log, "minimize: factorization failed ");
        return MinimizeResult::Failed;
    }
 
    std::size_t totalMeasOutliers = 0;
    std::size_t totalRefOutliers = 0;
    double oldChi2 = computeChi2().chi2;
    double oldSigmaCut = 0.;
    double sigmaCut = 0.;
 
    while (true) {
        Eigen::VectorXd delta = chol.solve(grad);
        if (doLineSearch) {
            scale = _lineSearch(delta);
        }
        offsetParams(scale * delta);
        Chi2Statistic currentChi2(computeChi2());
        LOGLS_DEBUG(_log, currentChi2);
        if (!isfinite(currentChi2.chi2)) {
            LOGL_ERROR(_log, "chi2 is not finite. Aborting outlier rejection.");
            returnCode = MinimizeResult::NonFinite;
            break;
        }
        if (currentChi2.chi2 > oldChi2 && totalMeasOutliers + totalRefOutliers != 0) {
            LOGL_WARN(_log, "chi2 went up, skipping outlier rejection loop");
            returnCode = MinimizeResult::Chi2Increased;
            break;
        }
        oldChi2 = currentChi2.chi2;
 
        if (nSigmaCut == 0) break;  // no rejection step to perform
        MeasuredStarList msOutliers;
        FittedStarList fsOutliers;
        // keep nOutliers so we don't have to sum msOutliers.size()+fsOutliers.size() twice below.
        std::size_t nOutliers = findOutliers(nSigmaCut, msOutliers, fsOutliers, sigmaCut);
        double relChange = 0.;
        if(oldSigmaCut!=0.) relChange = (1 - sigmaCut / oldSigmaCut);
 
        LOGLS_DEBUG(_log, "findOutliers chi2 cut level: " << sigmaCut << ", relative change: " << relChange);
        // If sigmaRelativeTolerance is set and at least one iteration has been done, break loop when the
        // fractional change in sigmaCut levels is less than the sigmaRelativeTolerance parameter.
        if ((sigmaRelativeTolerance > 0) && (oldSigmaCut > 0 && relChange < sigmaRelativeTolerance)) {
            LOGLS_INFO(_log, "Iterations stopped with chi2 cut at " << sigmaCut << " and relative change of "
                                                                    << relChange);
            break;
        }
        totalMeasOutliers += msOutliers.size();
        totalRefOutliers += fsOutliers.size();
        oldSigmaCut = sigmaCut;
        if (nOutliers == 0) break;
        TripletList outlierTriplets(nOutliers);
        grad.setZero();  // recycle the gradient
        // compute the contributions of outliers to derivatives
        outliersContributions(msOutliers, fsOutliers, outlierTriplets, grad);
        // Remove significant outliers
        removeMeasOutliers(msOutliers);
        removeRefOutliers(fsOutliers);
        if (doRankUpdate) {
            // convert triplet list to eigen internal format
            SparseMatrixD H(_nTotal, outlierTriplets.getNextFreeIndex());
            H.setFromTriplets(outlierTriplets.begin(), outlierTriplets.end());
            chol.update(H, false /* means downdate */);
            // 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 {
            // don't reuse tripletList because we want a new nextFreeIndex.
            TripletList nextTripletList(_lastNTrip);
            grad.setZero();
            // Rebuild the matrix and gradient
            leastSquareDerivatives(nextTripletList, grad);
            _lastNTrip = nextTripletList.size();
            LOGLS_DEBUG(_log, "Triplets recomputed, ntrip = " << nextTripletList.size());
 
            hessian = createHessian(_nTotal, nextTripletList);
            nextTripletList.clear();  // we don't need it any more after we have the hessian.
 
            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;
            }
        }
    }
 
    if (totalMeasOutliers + totalRefOutliers > 0) {
        _associations->cleanFittedStars();
    }
 
    // only print the outlier summary if outlier rejection was turned on.
    if (nSigmaCut != 0) {
        LOGLS_INFO(_log, "Number of outliers (Measured + Reference = Total): "
                                 << totalMeasOutliers << " + " << totalRefOutliers << " = "
                                 << totalMeasOutliers + totalRefOutliers);
    }
    return returnCode;
}

offsetParams()

void lsst::jointcal::PhotometryFit::offsetParams ( Eigen::VectorXd const & delta )

overridevirtual

Offset the parameters by the requested quantities.

The used parameter layout is the one from the last call to assignIndices or minimize(). There is no easy way to check that the current setting of whatToFit and the provided Delta vector are compatible: we can only test the size.

Parameters

[in]

delta

vector of offsets to apply

Implements lsst::jointcal::FitterBase.

Definition at line 211 of file PhotometryFit.cc.

                                                           {
    if (delta.size() != _nTotal)
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          "PhotometryFit::offsetParams : the provided vector length is not compatible with "
                          "the current whatToFit setting");
    if (_fittingModel) _photometryModel->offsetParams(delta);
 
    if (_fittingFluxes) {
        for (auto &fittedStar : _associations->fittedStarList) {
            // the parameter layout here is used also
            // - when filling the derivatives
            // - when assigning indices (assignIndices())
            Eigen::Index index = fittedStar->getIndexInMatrix();
            _photometryModel->offsetFittedStar(*fittedStar, delta(index));
        }
    }
}

operator=() [1/2]

PhotometryFit & lsst::jointcal::PhotometryFit::operator= ( PhotometryFit && )

delete

operator=() [2/2]

PhotometryFit & lsst::jointcal::PhotometryFit::operator= ( PhotometryFit const & )

delete

outliersContributions()

void lsst::jointcal::FitterBase::outliersContributions ( MeasuredStarList & msOutliers, FittedStarList & fsOutliers, TripletList & tripletList, Eigen::VectorXd & grad )

protectedinherited

Contributions to derivatives from (presumably) outlier terms.

No discarding done.

Definition at line 322 of file FitterBase.cc.

                                                                                      {
    for (auto &outlier : msOutliers) {
        MeasuredStarList tmp;
        tmp.push_back(outlier);
        const CcdImage &ccdImage = outlier->getCcdImage();
        leastSquareDerivativesMeasurement(ccdImage, tripletList, grad, &tmp);
    }
    leastSquareDerivativesReference(fsOutliers, tripletList, grad);
}

removeMeasOutliers()

void lsst::jointcal::FitterBase::removeMeasOutliers ( MeasuredStarList & outliers )

protectedinherited

Remove measuredStar outliers from the fit. No Refit done.

Definition at line 333 of file FitterBase.cc.

                                                              {
    for (auto &measuredStar : outliers) {
        auto fittedStar = measuredStar->getFittedStar();
        measuredStar->setValid(false);
        fittedStar->getMeasurementCount()--;  // could be put in setValid
    }
}

removeRefOutliers()

void lsst::jointcal::FitterBase::removeRefOutliers ( FittedStarList & outliers )

protectedinherited

Remove refStar outliers from the fit. No Refit done.

Definition at line 341 of file FitterBase.cc.

                                                           {
    for (auto &fittedStar : outliers) {
        fittedStar->setRefStar(nullptr);
    }
}

saveChi2Contributions()

void lsst::jointcal::FitterBase::saveChi2Contributions ( std::string const & baseName ) const

virtualinherited

Save the full chi2 term per star that was used in the minimization, for debugging.

Saves results to text files "baseName-meas.csv" and "baseName-ref.csv" for the MeasuredStar and RefStar contributions, respectively. This method is mostly useful for debugging: we will probably want to create a better persistence system for jointcal's internal representations in the future (see DM-12446).

Definition at line 355 of file FitterBase.cc.

                                                                      {
    std::string replaceStr = "{type}";
    auto pos = baseName.find(replaceStr);
    std::string measFilename(baseName);
    measFilename.replace(pos, replaceStr.size(), "-meas.csv");
    std::string refFilename(baseName);
    refFilename.replace(pos, replaceStr.size(), "-ref.csv");
    saveChi2MeasContributions(measFilename);
    saveChi2RefContributions(refFilename);
}

saveChi2MeasContributions()

void lsst::jointcal::PhotometryFit::saveChi2MeasContributions ( std::string const & filename ) const

overrideprotectedvirtual

Save a CSV file containing residuals of measurement terms.

Implements lsst::jointcal::FitterBase.

Definition at line 229 of file PhotometryFit.cc.

                                                                             {
    std::ofstream ofile(filename.c_str());
    std::string separator = "\t";
 
    ofile << "id" << separator << "xccd" << separator << "yccd" << separator;
    ofile << "mag" << separator << "instMag" << separator << "instMagErr" << separator;
    ofile << "instFlux" << separator << "instFluxErr" << separator;
    ofile << "inputFlux" << separator << "inputFluxErr" << separator;
    ofile << "transformedFlux" << separator << "transformedFluxErr" << separator;
    ofile << "fittedFlux" << separator;
    ofile << "epoch" << separator;
    ofile << "fsindex" << separator;
    ofile << "ra" << separator << "dec" << separator;
    ofile << "chi2" << separator << "nm" << separator;
    ofile << "detector" << separator << "visit" << separator << std::endl;
 
    ofile << "id in source catalog" << separator << "coordinates in CCD" << separator << separator;
    ofile << "fitted magnitude" << separator << "measured magnitude" << separator
          << "measured magnitude error" << separator;
    ofile << "measured instrumental flux (ADU)" << separator << "measured instrument flux error" << separator;
    ofile << "measured flux (nJy)" << separator << "measured flux error" << separator;
    ofile << separator << separator;
    ofile << "fitted flux (nJy)" << separator;
    ofile << "Julian Epoch Year of the measurement" << separator;
    ofile << "unique index of the fittedStar" << separator;
    ofile << "on-sky position of fitted star" << separator << separator;
    ofile << "contribution to chi2 (1 dof)" << separator << "number of measurements of this FittedStar"
          << separator;
    ofile << "detector id" << separator << "visit id" << std::endl;
 
    const CcdImageList &ccdImageList = _associations->getCcdImageList();
    for (auto const &ccdImage : ccdImageList) {
        const MeasuredStarList &cat = ccdImage->getCatalogForFit();
        for (auto const &measuredStar : cat) {
            if (!measuredStar->isValid()) continue;
 
            double instFluxErr = _photometryModel->tweakFluxError(*measuredStar);
            double flux = _photometryModel->transform(*ccdImage, *measuredStar);
            double fluxErr = _photometryModel->transformError(*ccdImage, *measuredStar);
            std::shared_ptr<FittedStar const> const fittedStar = measuredStar->getFittedStar();
            double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
            double chi2Val = std::pow(residual / fluxErr, 2);
 
            ofile << std::setprecision(9);
            ofile << measuredStar->getId() << separator << measuredStar->x << separator << measuredStar->y
                  << separator;
            ofile << fittedStar->getMag() << separator << measuredStar->getInstMag() << separator
                  << measuredStar->getInstMagErr() << separator;
            ofile << measuredStar->getInstFlux() << separator << instFluxErr << separator;
            ofile << measuredStar->getFlux() << separator << measuredStar->getFluxErr() << separator;
            ofile << flux << separator << fluxErr << separator << fittedStar->getFlux() << separator;
            ofile << ccdImage->getEpoch() << separator;
            ofile << fittedStar->getIndexInMatrix() << separator;
            ofile << fittedStar->x << separator << fittedStar->y << separator;
            ofile << chi2Val << separator << fittedStar->getMeasurementCount() << separator;
            ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() << std::endl;
        }  // loop on measurements in image
    }      // loop on images
}

saveChi2RefContributions()

void lsst::jointcal::PhotometryFit::saveChi2RefContributions ( std::string const & filename ) const

overrideprotectedvirtual

Save a CSV file containing residuals of reference terms.

Implements lsst::jointcal::FitterBase.

Definition at line 289 of file PhotometryFit.cc.

                                                                            {
    std::ofstream ofile(filename.c_str());
    std::string separator = "\t";
 
    ofile << "ra" << separator << "dec " << separator;
    ofile << "mag" << separator;
    ofile << "refFlux" << separator << "refFluxErr" << separator;
    ofile << "fittedFlux" << separator << "fittedFluxErr" << separator;
    ofile << "fsindex" << separator << "chi2" << separator << "nm" << std::endl;
 
    ofile << "coordinates of fittedStar" << separator << separator;
    ofile << "magnitude" << separator;
    ofile << "refStar flux (nJy)" << separator << "refStar fluxErr" << separator;
    ofile << "fittedStar flux (nJy)" << separator << "fittedStar fluxErr" << separator;
    ofile << "unique index of the fittedStar" << separator << "refStar contribution to chi2 (1 dof)"
          << separator << "number of measurements of this FittedStar" << std::endl;
 
    // The following loop is heavily inspired from PhotometryFit::computeChi2()
    const FittedStarList &fittedStarList = _associations->fittedStarList;
    for (auto const &fittedStar : fittedStarList) {
        const RefStar *refStar = fittedStar->getRefStar();
        if (refStar == nullptr) continue;
 
        double chi2 = std::pow(((fittedStar->getFlux() - refStar->getFlux()) / refStar->getFluxErr()), 2);
 
        ofile << std::setprecision(9);
        ofile << fittedStar->x << separator << fittedStar->y << separator;
        ofile << fittedStar->getMag() << separator;
        ofile << refStar->getFlux() << separator << refStar->getFluxErr() << separator;
        ofile << fittedStar->getFlux() << separator << fittedStar->getFluxErr() << separator;
        ofile << fittedStar->getIndexInMatrix() << separator << chi2 << separator
              << fittedStar->getMeasurementCount() << std::endl;
    }  // loop on FittedStars
}

Member Data Documentation

_associations

std::shared_ptr<Associations> lsst::jointcal::FitterBase::_associations

protectedinherited

Definition at line 165 of file FitterBase.h.

_lastNTrip

Eigen::Index lsst::jointcal::FitterBase::_lastNTrip

protectedinherited

Definition at line 168 of file FitterBase.h.

_log

LOG_LOGGER lsst::jointcal::FitterBase::_log

protectedinherited

Definition at line 174 of file FitterBase.h.

_nModelParams

Eigen::Index lsst::jointcal::FitterBase::_nModelParams

protectedinherited

Definition at line 170 of file FitterBase.h.

_nStarParams

Eigen::Index lsst::jointcal::FitterBase::_nStarParams

protectedinherited

Definition at line 171 of file FitterBase.h.

_nTotal

Eigen::Index lsst::jointcal::FitterBase::_nTotal

protectedinherited

Definition at line 169 of file FitterBase.h.

_whatToFit

std::string lsst::jointcal::FitterBase::_whatToFit

protectedinherited

Definition at line 166 of file FitterBase.h.

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/jointcal/g82479be7b0+f568feb641/include/lsst/jointcal/PhotometryFit.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/jointcal/g82479be7b0+f568feb641/src/PhotometryFit.cc