lsst::jointcal::AstrometryFit Class Reference

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

#include <AstrometryFit.h>

Inheritance diagram for lsst::jointcal::AstrometryFit:

Public Member Functions
	AstrometryFit (std::shared_ptr< Associations > associations, std::shared_ptr< AstrometryModel > astrometryModel, double posError)
	this is the only constructor More...
	AstrometryFit (AstrometryFit const &)=delete
	No copy or move: there is only ever one fitter of a given type. More...
	AstrometryFit (AstrometryFit &&)=delete
AstrometryFit &	operator= (AstrometryFit const &)=delete
AstrometryFit &	operator= (AstrometryFit &&)=delete
void	assignIndices (std::string const &whatToFit) override
	Set parameters to fit and assign indices in the big matrix. More...
void	freezeErrorTransform ()
	The transformations used to propagate errors are freezed to the current state. More...
void	offsetParams (Eigen::VectorXd const &delta) override
	Offset the parameters by the requested quantities. More...
std::shared_ptr< AstrometryModel >	getModel () const
	Return the model being fit. More...
void	checkStuff ()
	DEBUGGING routine. More...
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. More...
Chi2Statistic	computeChi2 () const
	Returns the chi2 for the current state. More...
void	leastSquareDerivatives (TripletList &tripletList, Eigen::VectorXd &grad) const
	Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting. More...
virtual void	saveChi2Contributions (std::string const &baseName) const
	Save the full chi2 term per star that was used in the minimization, for debugging. More...

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

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

Detailed Description

Class that handles the astrometric least squares problem.

This is the class that actually computes the quantities required to carry out a LS astrometric fit wrt distortion mappings and coordinates of common objects. Namely it computes the Jacobian and gradient of the chi2 (w.r.t. parameters), and the Chi2 itself. It interfaces with the actual modelling of distortions via a mimimum virtual interface AstrometryModel, and the actual mappings via an other virtual interface : Mapping.

In short AstrometryFit aims at computing derivatives of least quares. The terms of the chi2 are of two kinds:

kind 1 -> (T(X_M) - p(F))^T W (T(X_M) - p(F))

with X_M is a measured (2d) position in CCD coordinates, F refers to the position of the object in some space, defined in practise by p. There is one such term per measurement. The default setup would be that p is the projection from sky to some tangent plane and hence T maps the CCD coordinates onto this TP. p is obtained via the DistorsionModel and can be different for all CcdImage's. Depending on what is beeing fitted, one could imagine cases where the projector p is the same for all CcdImages.

Kind 2 -> (p'(F)-p'(R))^T W_R (p'(F)-p'(R)) R refers to some externally-provided reference object position, and p' to some projector from sky to some plane. The reference objects define the overall coordinate frame, which is required when all T and all F are fitted simultaneously. There is one such term per external reference object. There can be more F (fitted) objects than R (reference) objects.

In the same framework, one can fit relative transforms between images by setting p = Identity for all input CcdImages and not fitting T for one of the CcdImage's. One does not need reference object and would then naturally not have any Kind 2 terms.

Definition at line 78 of file AstrometryFit.h.

Constructor & Destructor Documentation

AstrometryFit() [1/3]

lsst::jointcal::AstrometryFit::AstrometryFit	(	std::shared_ptr< Associations >	associations,
		std::shared_ptr< AstrometryModel >	astrometryModel,
		double	posError
	)

this is the only constructor

Definition at line 66 of file AstrometryFit.cc.

        : FitterBase(std::move(associations)),
          _astrometryModel(std::move(astrometryModel)),
          _epoch(_associations->getEpoch()),
          _posError(posError) {
    _log = LOG_GET("lsst.jointcal.AstrometryFit");
}

AstrometryFit() [2/3]

lsst::jointcal::AstrometryFit::AstrometryFit ( AstrometryFit const &  )

delete

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

AstrometryFit() [3/3]

lsst::jointcal::AstrometryFit::AstrometryFit ( AstrometryFit &&  )

delete

Member Function Documentation

assignIndices()

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

overridevirtual

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

Parameters

[in]

whatToFit

Valid strings: zero or more of "Distortions", "Positions", "PM" which define which parameter set is going to be variable when computing derivatives (leastSquareDerivatives) and minimizing (minimize()). whatToFit="Positions Distortions" will minimize w.r.t mappings and objects positions, and not w.r.t proper motions and refraction modeling. However if proper motions and/or refraction parameters have already been set, then they are accounted for when computing residuals. The string is forwarded to the AstrometryModel, and it can then be used to turn subsets of distortion parameter on or off, if the AstrometryModel implements such a thing.

Implements lsst::jointcal::FitterBase.

Definition at line 423 of file AstrometryFit.cc.

                                                            {
    _whatToFit = whatToFit;
    LOGLS_INFO(_log, "assignIndices: Now fitting " << whatToFit);
    _fittingDistortions = (_whatToFit.find("Distortions") != std::string::npos);
    _fittingPos = (_whatToFit.find("Positions") != std::string::npos);
    _fittingPM = (_whatToFit.find("PM") != std::string::npos);
    // When entering here, we assume that whatToFit has already been interpreted.
 
    _nModelParams = 0;
    if (_fittingDistortions) _nModelParams = _astrometryModel->assignIndices(_whatToFit, 0);
    std::size_t ipar = _nModelParams;
 
    if (_fittingPos) {
        FittedStarList &fittedStarList = _associations->fittedStarList;
        for (auto &fittedStar : fittedStarList) {
            // the parameter layout here is used also
            // - when filling the derivatives
            // - when updating (offsetParams())
            // - in GetMeasuredStarIndices
            fittedStar->setIndexInMatrix(ipar);
            ipar += 2;
            // TODO: left as reference for when we implement PM fitting
            // if ((_fittingPM)&fittedStar->mightMove) ipar += 2;
        }
    }
    _nStarParams = ipar - _nModelParams;
    _nTotal = ipar;
    LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
                                            << " values: " << _nStarParams);
}

checkStuff()

void lsst::jointcal::AstrometryFit::checkStuff

(

)

DEBUGGING routine.

Definition at line 480 of file AstrometryFit.cc.

                               {
    const char *what2fit[] = {"Positions", "Distortions", "Positions Distortions"};
    // DEBUG
    for (unsigned k = 0; k < sizeof(what2fit) / sizeof(what2fit[0]); ++k) {
        assignIndices(what2fit[k]);
        TripletList tripletList(10000);
        Eigen::VectorXd grad(_nTotal);
        grad.setZero();
        leastSquareDerivatives(tripletList, grad);
        SparseMatrixD jacobian(_nTotal, tripletList.getNextFreeIndex());
        jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
        SparseMatrixD hessian = jacobian * jacobian.transpose();
        LOGLS_DEBUG(_log, "npar : " << _nTotal << ' ' << _nModelParams);
    }
}

computeChi2()

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

inherited

Returns the chi2 for the current state.

Definition at line 42 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 52 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;
}

freezeErrorTransform()

void lsst::jointcal::AstrometryFit::freezeErrorTransform ( )

inline

The transformations used to propagate errors are freezed to the current state.

The routine can be called when the mappings are roughly in place. After the call, the transformations used to propage errors are no longer affected when updating the mappings. This allows to have an exactly linear fit, which can be useful.

Definition at line 117 of file AstrometryFit.h.

{ _astrometryModel->freezeErrorTransform(); }

getModel()

std::shared_ptr< AstrometryModel > lsst::jointcal::AstrometryFit::getModel ( ) const

inline

Return the model being fit.

Definition at line 122 of file AstrometryFit.h.

{ return _astrometryModel; }

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 346 of file FitterBase.cc.

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

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 178 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::AstrometryFit::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 454 of file AstrometryFit.cc.

                                                           {
    if (delta.size() != _nTotal)
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          "AstrometryFit::offsetParams : the provided vector length is not compatible with "
                          "the current whatToFit setting");
    if (_fittingDistortions) _astrometryModel->offsetParams(delta);
 
    if (_fittingPos) {
        FittedStarList &fittedStarList = _associations->fittedStarList;
        for (auto const &i : fittedStarList) {
            FittedStar &fs = *i;
            // the parameter layout here is used also
            // - when filling the derivatives
            // - when assigning indices (assignIndices())
            Eigen::Index index = fs.getIndexInMatrix();
            fs.x += delta(index);
            fs.y += delta(index + 1);
            // TODO: left as reference for when we implement PM fitting
            // if ((_fittingPM)&fs.mightMove) {
            //     fs.pmx += delta(index + 2);
            //     fs.pmy += delta(index + 3);
            // }
        }
    }
}

operator=() [1/2]

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

delete

operator=() [2/2]

AstrometryFit & lsst::jointcal::AstrometryFit::operator= ( AstrometryFit 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 321 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 332 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 340 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 354 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::AstrometryFit::saveChi2MeasContributions ( std::string const &  filename ) const

overrideprotectedvirtual

Save a CSV file containing residuals of measurement terms.

Implements lsst::jointcal::FitterBase.

Definition at line 496 of file AstrometryFit.cc.

                                                                             {
    std::ofstream ofile(filename.c_str());
    std::string separator = "\t";
 
    ofile << "id" << separator << "xccd" << separator << "yccd " << separator;
    ofile << "rx" << separator << "ry" << separator;
    ofile << "xtp" << separator << "ytp" << separator;
    ofile << "mag" << separator << "deltaYears" << separator;
    ofile << "xErr" << separator << "yErr" << separator << "xyCov" << separator;
    ofile << "xtpi" << separator << "ytpi" << separator;
    ofile << "rxi" << separator << "ryi" << separator;
    ofile << "fsindex" << separator;
    ofile << "ra" << separator << "dec" << separator;
    ofile << "chi2" << separator << "nm" << separator;
    ofile << "detector" << separator << "visit" << std::endl;
 
    ofile << "id in source catalog" << separator << "coordinates in CCD (pixels)" << separator << separator;
    ofile << "residual on TP (degrees)" << separator << separator;
    ofile << "transformed coordinate in TP (degrees)" << separator << separator;
    ofile << "rough magnitude" << separator << "Julian epoch year delta from fit epoch" << separator;
    ofile << "transformed measurement uncertainty (degrees)" << separator << separator << separator;
    ofile << "as-read position on TP (degrees)" << separator << separator;
    ofile << "as-read residual on TP (degrees)" << separator << separator;
    ofile << "unique index of the fittedStar" << separator;
    ofile << "on sky position of fittedStar" << separator << separator;
    ofile << "contribution to chi2 (2D dofs)" << 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();
        const AstrometryMapping *mapping = _astrometryModel->getMapping(*ccdImage);
        const auto readTransform = ccdImage->getReadWcs();
        // FittedStar is "observed" epoch, MeasuredStar is "baseline"
        double deltaYears = _epoch - ccdImage->getEpoch();
        for (auto const &ms : cat) {
            if (!ms->isValid()) continue;
            FatPoint tpPos;
            FatPoint inPos = *ms;
            tweakAstromMeasurementErrors(inPos, *ms, _posError);
            mapping->transformPosAndErrors(inPos, tpPos);
            auto sky2TP = _astrometryModel->getSkyToTangentPlane(*ccdImage);
            const std::unique_ptr<AstrometryTransform> readPixToTangentPlane =
                    compose(*sky2TP, *readTransform);
            FatPoint inputTpPos = readPixToTangentPlane->apply(inPos);
            std::shared_ptr<FittedStar const> const fs = ms->getFittedStar();
 
            Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
            Point res = tpPos - fittedStarInTP;
            Point inputRes = inputTpPos - fittedStarInTP;
            double det = tpPos.vx * tpPos.vy - std::pow(tpPos.vxy, 2);
            double wxx = tpPos.vy / det;
            double wyy = tpPos.vx / det;
            double wxy = -tpPos.vxy / det;
            double chi2 = wxx * res.x * res.x + wyy * res.y * res.y + 2 * wxy * res.x * res.y;
 
            ofile << std::setprecision(9);
            ofile << ms->getId() << separator << ms->x << separator << ms->y << separator;
            ofile << res.x << separator << res.y << separator;
            ofile << tpPos.x << separator << tpPos.y << separator;
            ofile << fs->getMag() << separator << deltaYears << separator;
            ofile << tpPos.vx << separator << tpPos.vy << separator << tpPos.vxy << separator;
            ofile << inputTpPos.x << separator << inputTpPos.y << separator;
            ofile << inputRes.x << separator << inputRes.y << separator;
            ofile << fs->getIndexInMatrix() << separator;
            ofile << fs->x << separator << fs->y << separator;
            ofile << chi2 << separator << fs->getMeasurementCount() << separator;
            ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() << std::endl;
        }  // loop on measurements in image
    }      // loop on images
}

saveChi2RefContributions()

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

overrideprotectedvirtual

Save a CSV file containing residuals of reference terms.

Implements lsst::jointcal::FitterBase.

Definition at line 569 of file AstrometryFit.cc.

                                                                            {
    std::ofstream ofile(filename.c_str());
    std::string separator = "\t";
 
    ofile << "ra" << separator << "dec " << separator;
    ofile << "rx" << separator << "ry" << separator;
    ofile << "mag" << separator;
    ofile << "xErr" << separator << "yErr" << separator << "xyCov" << separator;
    ofile << "fsindex" << separator;
    ofile << "chi2" << separator << "nm" << std::endl;
 
    ofile << "coordinates of fittedStar (degrees)" << separator << separator;
    ofile << "residual on TP (degrees)" << separator << separator;
    ofile << "magnitude" << separator;
    ofile << "refStar transformed measurement uncertainty (degrees)" << separator << separator << separator;
    ofile << "unique index of the fittedStar" << separator;
    ofile << "refStar contribution to chi2 (2D dofs)" << separator
          << "number of measurements of this FittedStar" << std::endl;
 
    // The following loop is heavily inspired from AstrometryFit::computeChi2()
    const FittedStarList &fittedStarList = _associations->fittedStarList;
    TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
    for (auto const &i : fittedStarList) {
        const FittedStar &fs = *i;
        auto rs = fs.getRefStar();
        if (rs == nullptr) continue;
        proj.setTangentPoint(fs);
        // fs projects to (0,0), no need to compute its transform.
        FatPoint rsProj;
        proj.transformPosAndErrors(*rs, rsProj);
        double chi2 = computeProjectedRefStarChi2(rsProj);
 
        ofile << std::setprecision(9);
        ofile << fs.x << separator << fs.y << separator;
        ofile << rsProj.x << separator << rsProj.y << separator;
        ofile << fs.getMag() << separator;
        ofile << rsProj.vx << separator << rsProj.vy << separator << rsProj.vxy << separator;
        ofile << fs.getIndexInMatrix() << separator;
        ofile << chi2 << separator << fs.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-4.0.1/Linux64/jointcal/ge2d134c3d5+cc4dbb2e3f/include/lsst/jointcal/AstrometryFit.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-4.0.1/Linux64/jointcal/ge2d134c3d5+cc4dbb2e3f/src/AstrometryFit.cc