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
	AstrometryFit (AstrometryFit const &)=delete
	No copy or move: there is only ever one fitter of a given type.
	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.
void	freezeErrorTransform ()
	The transformations used to propagate errors are freezed to the current state.
void	offsetParams (Eigen::VectorXd const &delta) override
	Offset the parameters by the requested quantities.
std::shared_ptr< AstrometryModel >	getModel () const
	Return the model being fit.
void	checkStuff ()
	DEBUGGING routine.
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	leastSquareDerivativesMeasurement (CcdImage const &ccdImage, TripletList &tripletList, Eigen::VectorXd &grad, MeasuredStarList const *msList=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.
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

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

accumulateStatImageList()

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

overrideprivatevirtual

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

Implements lsst::jointcal::FitterBase.

Definition at line 372 of file AstrometryFit.cc.

                                                                                                          {
    for (auto const &ccdImage : ccdImageList) {
        accumulateStatImage(*ccdImage, accum);
    }
}

accumulateStatRefStars()

void lsst::jointcal::AstrometryFit::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 378 of file AstrometryFit.cc.

                                                                       {
    /**********************************************************************/
    /**********************************************************************/
 
    /* If you wonder why we project here, read comments in
       AstrometryFit::leastSquareDerivativesReference(TripletList &TList, Eigen::VectorXd &Rhs) */
    FittedStarList &fittedStarList = _associations->fittedStarList;
    TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
    for (auto const &fs : fittedStarList) {
        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);
        // TO DO : account for proper motions.
        double chi2 = computeProjectedRefStarChi2(rsProj);
        accum.addEntry(chi2, 2, fs);
    }
}

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 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;
}

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(); }

getIndicesOfMeasuredStar()

void lsst::jointcal::AstrometryFit::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 404 of file AstrometryFit.cc.

                                                                                                         {
    if (_fittingDistortions) {
        const AstrometryMapping *mapping = _astrometryModel->getMapping(measuredStar.getCcdImage());
        mapping->getMappingIndices(indices);
    }
    std::shared_ptr<FittedStar const> const fs = measuredStar.getFittedStar();
    Eigen::Index fsIndex = fs->getIndexInMatrix();
    if (_fittingPos) {
        indices.push_back(fsIndex);
        indices.push_back(fsIndex + 1);
    }
    // For securing the outlier removal, the next block is just useless
    if (_fittingPM) {
        for (std::size_t k = 0; k < 2; ++k) indices.push_back(fsIndex + 2 + k);
    }
    /* Should not put the index of refaction stuff or we will not be
       able to remove more than 1 star at a time. */
}

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

                                                                                            {
    /**********************************************************************/
    /* @note the math in this method and accumulateStatImage() must be kept consistent,
     * in terms of +/- convention, definition of model, etc. */
    /**********************************************************************/
 
    /* Setup */
    /* this routine works in two different ways: either providing the
       ccdImage, of providing the MeasuredStarList. In the latter case, the
       ccdImage should match the one(s) in the list. */
    if (msList) assert(&(msList->front()->getCcdImage()) == &ccdImage);
 
    // get the Mapping
    const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
    // count parameters
    std::size_t npar_mapping = (_fittingDistortions) ? mapping->getNpar() : 0;
    std::size_t npar_pos = (_fittingPos) ? 2 : 0;
    std::size_t npar_pm = (_fittingPM) ? 2 : 0;
    std::size_t npar_tot = npar_mapping + npar_pos + npar_pm;
    // if (npar_tot == 0) this CcdImage does not contribute
    // any constraint to the fit, so :
    if (npar_tot == 0) return;
    IndexVector indices(npar_tot, -1);
    if (_fittingDistortions) mapping->getMappingIndices(indices);
 
    // FittedStar is "observed" epoch, MeasuredStar is "baseline"
    double deltaYears = _epoch - ccdImage.getEpoch();
    // transformation from sky to TP
    auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
    // reserve matrices once for all measurements
    AstrometryTransformLinear dypdy;
    // the shape of H (et al) is required this way in order to be able to
    // separate derivatives along x and y as vectors.
    Eigen::MatrixX2d H(npar_tot, 2), halpha(npar_tot, 2), HW(npar_tot, 2);
    Eigen::Matrix2d transW(2, 2);
    Eigen::Matrix2d alpha(2, 2);
    Eigen::VectorXd grad(npar_tot);
    // current position in the Jacobian
    Eigen::Index kTriplets = tripletList.getNextFreeIndex();
    const MeasuredStarList &catalog = (msList) ? *msList : ccdImage.getCatalogForFit();
 
    for (auto &i : catalog) {
        const MeasuredStar &ms = *i;
        if (!ms.isValid()) continue;
        // tweak the measurement errors
        FatPoint inPos = ms;
        tweakAstromMeasurementErrors(inPos, ms, _posError);
        H.setZero();  // we cannot be sure that all entries will be overwritten.
        FatPoint outPos;
        // should *not* fill H if whatToFit excludes mapping parameters.
        if (_fittingDistortions)
            mapping->computeTransformAndDerivatives(inPos, outPos, H);
        else
            mapping->transformPosAndErrors(inPos, outPos);
 
        std::size_t ipar = npar_mapping;
        double det = outPos.vx * outPos.vy - std::pow(outPos.vxy, 2);
        if (det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
            LOGLS_WARN(_log, "Inconsistent measurement errors: drop measurement at "
                                     << Point(ms) << " in image " << ccdImage.getName());
            continue;
        }
        transW(0, 0) = outPos.vy / det;
        transW(1, 1) = outPos.vx / det;
        transW(0, 1) = transW(1, 0) = -outPos.vxy / det;
        // compute alpha, a triangular square root
        // of transW (i.e. a Cholesky factor)
        alpha(0, 0) = sqrt(transW(0, 0));
        // checked that  alpha*alphaT = transW
        alpha(1, 0) = transW(0, 1) / alpha(0, 0);
        // DB - I think that the next line is equivalent to : alpha(1,1) = 1./sqrt(outPos.vy)
        // PA - seems correct !
        alpha(1, 1) = 1. / sqrt(det * transW(0, 0));
        alpha(0, 1) = 0;
 
        std::shared_ptr<FittedStar const> const fs = ms.getFittedStar();
 
        Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
 
        // compute derivative of TP position w.r.t sky position ....
        if (npar_pos > 0)  // ... if actually fitting FittedStar position
        {
            sky2TP->computeDerivative(*fs, dypdy, 1e-3);
            // sign checked
            // TODO Still have to check with non trivial non-diagonal terms
            H(npar_mapping, 0) = -dypdy.A11();
            H(npar_mapping + 1, 0) = -dypdy.A12();
            H(npar_mapping, 1) = -dypdy.A21();
            H(npar_mapping + 1, 1) = -dypdy.A22();
            indices[npar_mapping] = fs->getIndexInMatrix();
            indices.at(npar_mapping + 1) = fs->getIndexInMatrix() + 1;
            ipar += npar_pos;
        }
        /* only consider proper motions of objects allowed to move,
        unless the fit is going to be degenerate */
        // TODO: left as reference for when we implement PM fitting
        // TODO: mjd here would become either deltaYears or maybe associations.epoch? Check the math first!
        // if (_fittingPM && fs->mightMove) {
        //     H(ipar, 0) = -mjd;  // Sign unchecked but consistent with above
        //     H(ipar + 1, 1) = -mjd;
        //     indices[ipar] = fs->getIndexInMatrix() + 2;
        //     indices[ipar + 1] = fs->getIndexInMatrix() + 3;
        //     ipar += npar_pm;
        // }
 
        // We can now compute the residual
        Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
 
        // do not write grad = H*transW*res to avoid
        // dynamic allocation of a temporary
        halpha = H * alpha;
        HW = H * transW;
        grad = HW * res;
        // now feed in triplets and fullGrad
        for (std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
            for (std::size_t ic = 0; ic < 2; ++ic) {
                double val = halpha(ipar, ic);
                if (val == 0) continue;
                tripletList.addTriplet(indices[ipar], kTriplets + ic, val);
            }
            fullGrad(indices[ipar]) += grad(ipar);
        }
        kTriplets += 2;  // each measurement contributes 2 columns in the Jacobian
    }                    // end loop on measurements
    tripletList.setNextFreeIndex(kTriplets);
}

leastSquareDerivativesReference()

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

overrideprivatevirtual

Compute the derivatives of the reference terms.

Implements lsst::jointcal::FitterBase.

Definition at line 236 of file AstrometryFit.cc.

                                                                                   {
    /**********************************************************************/
    /* @note the math in this method and accumulateStatRefStars() must be kept consistent,
     * in terms of +/- convention, definition of model, etc. */
    /**********************************************************************/
 
    /* We compute here the derivatives of the terms involving fitted
       stars and reference stars. They only provide contributions if we
       are fitting positions: */
    if (!_fittingPos) return;
    /* the other case where the accumulation of derivatives stops
       here is when there are no RefStars */
    if (_associations->refStarList.size() == 0) return;
    Eigen::Matrix2d W(2, 2);
    Eigen::Matrix2d alpha(2, 2);
    Eigen::Matrix2d H(2, 2), halpha(2, 2), HW(2, 2);
    AstrometryTransformLinear der;
    Eigen::Vector2d res, grad;
    Eigen::Index indices[2];
    Eigen::Index kTriplets = tripletList.getNextFreeIndex();
    /* We cannot use the spherical coordinates directly to evaluate
       Euclidean distances, we have to use a projector on some plane in
       order to express least squares. Not projecting could lead to a
       disaster around the poles or across alpha=0.  So we need a
       projector. We construct a projector and will change its
       projection point at every object */
    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);
        // Compute the derivative of the projector to incorporate its effects on the errors.
        proj.computeDerivative(fs, der, 1e-4);
        // sign checked. TODO check that the off-diagonal terms are OK.
        H(0, 0) = -der.A11();
        H(1, 0) = -der.A12();
        H(0, 1) = -der.A21();
        H(1, 1) = -der.A22();
        // TO DO : account for proper motions.
        double det = rsProj.vx * rsProj.vy - std::pow(rsProj.vxy, 2);
        if (rsProj.vx <= 0 || rsProj.vy <= 0 || det <= 0) {
            LOGLS_WARN(_log, "RefStar error matrix not positive definite for:  " << *rs);
            continue;
        }
        W(0, 0) = rsProj.vy / det;
        W(0, 1) = W(1, 0) = -rsProj.vxy / det;
        W(1, 1) = rsProj.vx / det;
        // compute alpha, a triangular square root
        // of W (i.e. a Cholesky factor)
        alpha(0, 0) = sqrt(W(0, 0));
        // checked that  alpha*alphaT = transW
        alpha(1, 0) = W(0, 1) / alpha(0, 0);
        alpha(1, 1) = 1. / sqrt(det * W(0, 0));
        alpha(0, 1) = 0;
        indices[0] = fs.getIndexInMatrix();
        indices[1] = fs.getIndexInMatrix() + 1;
        unsigned npar_tot = 2;
        /* TODO: account here for proper motions in the reference
        catalog. We can code the effect and set the value to 0. Most
        (all?)  catalogs do not even come with a reference epoch. Gaia
        will change that. When refraction enters into the game, one should
        pay attention to the orientation of the frame */
 
        /* The residual should be Proj(fs)-Proj(*rs) in order to be consistent
        with the measurement terms. Since P(fs) = 0, we have: */
        res[0] = -rsProj.x;
        res[1] = -rsProj.y;
        halpha = H * alpha;
        // grad = H*W*res
        HW = H * W;
        grad = HW * res;
        // now feed in triplets and fullGrad
        for (std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
            for (unsigned ic = 0; ic < 2; ++ic) {
                double val = halpha(ipar, ic);
                if (val == 0) continue;
                tripletList.addTriplet(indices[ipar], kTriplets + ic, val);
            }
            fullGrad(indices[ipar]) += grad(ipar);
        }
        kTriplets += 2;  // each measurement contributes 2 columns in the Jacobian
    }
    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::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 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::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-9.0.0/Linux64/jointcal/g82479be7b0+e2bd23ab8b/include/lsst/jointcal/AstrometryFit.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/jointcal/g82479be7b0+e2bd23ab8b/src/AstrometryFit.cc