LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
lsst::jointcal::PhotometryFit Class Reference

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

#include <PhotometryFit.h>

Inheritance diagram for lsst::jointcal::PhotometryFit:
lsst::jointcal::FitterBase

Public Member Functions

 PhotometryFit (std::shared_ptr< Associations > associations, std::shared_ptr< PhotometryModel > photometryModel)
 Construct a photometry fitter. More...
 
 PhotometryFit (PhotometryFit const &)=delete
 No copy or move: there is only ever one fitter of a given type. More...
 
 PhotometryFit (PhotometryFit &&)=delete
 
PhotometryFitoperator= (PhotometryFit const &)=delete
 
PhotometryFitoperator= (PhotometryFit &&)=delete
 
void assignIndices (std::string const &whatToFit) override
 Set parameters to fit and assign indices in the big matrix. More...
 
void offsetParams (Eigen::VectorXd const &delta) override
 Offset the parameters by the requested quantities. More...
 
std::shared_ptr< PhotometryModelgetModel () const
 Return the model being fit. More...
 
MinimizeResult minimize (std::string const &whatToFit, double const nSigmaCut=0, double sigmaRelativeTolerance=0, bool const doRankUpdate=true, bool const 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 photometric least squares problem.

Definition at line 45 of file PhotometryFit.h.

Constructor & Destructor Documentation

◆ PhotometryFit() [1/3]

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

Construct a photometry fitter.

Parameters
associationsThe associations catalog to use in the fitter.
photometryModelThe model to build the fitter for.

Definition at line 53 of file PhotometryFit.h.

55  : FitterBase(associations),
56  _fittingModel(false),
57  _fittingFluxes(false),
58  _photometryModel(photometryModel) {
59  _log = LOG_GET("lsst.jointcal.PhotometryFit");
60  }
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
FitterBase(std::shared_ptr< Associations > associations)
Definition: FitterBase.h:55

◆ PhotometryFit() [2/3]

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

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

◆ PhotometryFit() [3/3]

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

Member Function Documentation

◆ assignIndices()

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

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

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

Implements lsst::jointcal::FitterBase.

Definition at line 185 of file PhotometryFit.cc.

185  {
186  _whatToFit = whatToFit;
187  LOGLS_INFO(_log, "assignIndices: now fitting: " << whatToFit);
188  _fittingModel = (_whatToFit.find("Model") != std::string::npos);
189  _fittingFluxes = (_whatToFit.find("Fluxes") != std::string::npos);
190  // When entering here, we assume that whatToFit has already been interpreted.
191 
192  _nModelParams = (_fittingModel) ? _photometryModel->assignIndices(whatToFit, 0) : 0;
193  std::size_t ipar = _nModelParams;
194 
195  if (_fittingFluxes) {
196  for (auto &fittedStar : _associations->fittedStarList) {
197  // the parameter layout here is used also
198  // - when filling the derivatives
199  // - when updating (offsetParams())
200  // - in getIndicesOfMeasuredStar
201  fittedStar->setIndexInMatrix(ipar);
202  ipar += 1;
203  }
204  }
205  _nStarParams = ipar - _nModelParams;
206  _nTotal = ipar;
207  LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
208  << " fluxes: " << _nStarParams);
209 }
#define LOGLS_INFO(logger, message)
Log a info-level message using an iostream-based interface.
Definition: Log.h:639
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
Definition: Log.h:619
Eigen::Index _nModelParams
Definition: FitterBase.h:168
std::shared_ptr< Associations > _associations
Definition: FitterBase.h:163
T find(T... args)

◆ computeChi2()

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

Returns the chi2 for the current state.

Definition at line 42 of file FitterBase.cc.

42  {
43  Chi2Statistic chi2;
44  accumulateStatImageList(_associations->getCcdImageList(), chi2);
46  // chi2.ndof contains the number of squares.
47  // So subtract the number of parameters.
48  chi2.ndof -= _nTotal;
49  return chi2;
50 }
virtual void accumulateStatRefStars(Chi2Accumulator &accum) const =0
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.
virtual void accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const =0
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.

◆ 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]nSigmaCutNumber of sigma to select on.
[out]msOutlierslist of MeasuredStar outliers to populate
[out]fsOutlierslist of FittedStar outliers to populate
[out]cutvalue of chi2 that defines which objects are outliers
Returns
Total number of outliers that were removed.

Definition at line 52 of file FitterBase.cc.

53  {
54  // collect chi2 contributions
55  Chi2List chi2List;
56  chi2List.reserve(_associations->getMaxMeasuredStars() + _associations->refStarList.size());
57  // contributions from measurement terms:
58  accumulateStatImageList(_associations->ccdImageList, chi2List);
59  // and from reference terms
60  accumulateStatRefStars(chi2List);
61 
62  // compute some statistics
63  size_t nval = chi2List.size();
64  if (nval == 0) return 0;
65  sort(chi2List.begin(), chi2List.end());
66  double median = (nval & 1) ? chi2List[nval / 2].chi2
67  : 0.5 * (chi2List[nval / 2 - 1].chi2 + chi2List[nval / 2].chi2);
68  auto averageAndSigma = chi2List.computeAverageAndSigma();
69  LOGLS_DEBUG(_log, "findOutliers chi2 stat: mean/median/sigma " << averageAndSigma.first << '/' << median
70  << '/' << averageAndSigma.second);
71  cut = averageAndSigma.first + nSigmaCut * averageAndSigma.second;
72  /* For each of the parameters, we will not remove more than 1
73  measurement that contributes to constraining it. Keep track using
74  of what we are touching using an integer vector. This is the
75  trick that Marc Betoule came up to for outlier removals in "star
76  flats" fits. */
77  Eigen::VectorXi affectedParams(_nTotal);
78  affectedParams.setZero();
79 
80  std::size_t nOutliers = 0; // returned to the caller
81  // start from the strongest outliers.
82  for (auto chi2 = chi2List.rbegin(); chi2 != chi2List.rend(); ++chi2) {
83  if (chi2->chi2 < cut) break; // because the array is sorted.
84  IndexVector indices;
85  /* now, we want to get the indices of the parameters this chi2
86  term depends on. We have to figure out which kind of term it
87  is; we use for that the type of the star attached to the Chi2Star. */
88  auto measuredStar = std::dynamic_pointer_cast<MeasuredStar>(chi2->star);
89  std::shared_ptr<FittedStar> fittedStar; // To add to fsOutliers if it is a reference outlier.
90  if (measuredStar == nullptr) {
91  // it is a reference outlier
92  fittedStar = std::dynamic_pointer_cast<FittedStar>(chi2->star);
93  if (fittedStar->getMeasurementCount() == 0) {
94  LOGLS_WARN(_log, "FittedStar with no measuredStars found as an outlier: "
95  << *fittedStar << " chi2: " << chi2->chi2);
96  continue;
97  }
98  if (_nStarParams == 0) {
100  "RefStar is outlier but not removed when not fitting FittedStar-RefStar values: "
101  << *(fittedStar->getRefStar()) << " chi2: " << chi2->chi2);
102  continue;
103  }
104  // NOTE: Stars contribute twice to astrometry (x,y), but once to photometry (flux),
105  // NOTE: but we only need to mark one index here because both will be removed with that star.
106  indices.push_back(fittedStar->getIndexInMatrix());
107  LOGLS_TRACE(_log, "Removing refStar " << *(fittedStar->getRefStar()) << " chi2: " << chi2->chi2);
108  /* One might think it would be useful to account for PM
109  parameters here, but it is just useless */
110  } else {
111  // it is a measurement outlier
112  auto tempFittedStar = measuredStar->getFittedStar();
113  if (tempFittedStar->getMeasurementCount() == 1 && tempFittedStar->getRefStar() == nullptr) {
114  LOGLS_WARN(_log, "FittedStar with 1 measuredStar and no refStar found as an outlier: "
115  << *tempFittedStar);
116  continue;
117  }
118  getIndicesOfMeasuredStar(*measuredStar, indices);
119  LOGLS_TRACE(_log, "Removing measStar " << *measuredStar << " chi2: " << chi2->chi2);
120  }
121 
122  /* Find out if we already discarded a stronger outlier
123  constraining some parameter this one constrains as well. If
124  yes, we keep this one, because this stronger outlier could be
125  causing the large chi2 we have in hand. */
126  bool drop_it = true;
127  for (auto const &i : indices) {
128  if (affectedParams(i) != 0) {
129  drop_it = false;
130  }
131  }
132 
133  if (drop_it) // store the outlier in one of the lists:
134  {
135  if (measuredStar == nullptr) {
136  // reference term
137  fsOutliers.push_back(fittedStar);
138  } else {
139  // measurement term
140  msOutliers.push_back(measuredStar);
141  }
142  // mark the parameters as directly changed when we discard this chi2 term.
143  for (auto const &i : indices) {
144  affectedParams(i)++;
145  }
146  nOutliers++;
147  }
148  } // end loop on measurements/references
149  LOGLS_DEBUG(_log, "findOutliers: found " << msOutliers.size() << " meas outliers and "
150  << fsOutliers.size() << " ref outliers ");
151 
152  return nOutliers;
153 }
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:659
#define LOGLS_TRACE(logger, message)
Log a trace-level message using an iostream-based interface.
Definition: Log.h:599
virtual void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const =0
Set the indices of a measured star from the full matrix, for outlier removal.
T push_back(T... args)
T sort(T... args)

◆ getModel()

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

Return the model being fit.

Definition at line 86 of file PhotometryFit.h.

86 { return _photometryModel; }

◆ leastSquareDerivatives()

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

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

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

Parameters
tripletListtripletList of (row,col,value) representing the Jacobian of the chi2.
gradThe gradient of the chi2.

Definition at line 344 of file FitterBase.cc.

344  {
345  auto ccdImageList = _associations->getCcdImageList();
346  for (auto const &ccdImage : ccdImageList) {
347  leastSquareDerivativesMeasurement(*ccdImage, tripletList, grad);
348  }
349  leastSquareDerivativesReference(_associations->fittedStarList, tripletList, grad);
350 }
virtual void leastSquareDerivativesReference(FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const =0
Compute the derivatives of the reference terms.
virtual void leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList, Eigen::VectorXd &grad, MeasuredStarList const *measuredStarList=nullptr) const =0
Compute the derivatives of the measured stars and model for one CcdImage.

◆ minimize()

MinimizeResult lsst::jointcal::FitterBase::minimize ( std::string const &  whatToFit,
double const  nSigmaCut = 0,
double  sigmaRelativeTolerance = 0,
bool const  doRankUpdate = true,
bool const  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]whatToFitSee child method assignIndices for valid string values.
[in]nSigmaCutHow many sigma to reject outliers at. Outlier rejection ignored for nSigmaCut=0.
[in]sigmaRelativeTolerancePercentage change in the chi2 cut for outliers tolerated for termination. If value is zero, minimization iterations will continue until there are no outliers.
[in]doRankUpdateUse 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]doLineSearchUse 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]dumpMatrixFileWrite 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.

180  {
181  assignIndices(whatToFit);
182 
184 
185  // For the initial vector size, use all measured stars + all fitted stars, which should give the
186  // maximum possible number of triplets.
187  std::size_t nTrip = (_lastNTrip)
188  ? _lastNTrip
189  : _associations->getMaxMeasuredStars() + _associations->fittedStarList.size();
190  TripletList tripletList(nTrip);
191  Eigen::VectorXd grad(_nTotal);
192  grad.setZero();
193  double scale = 1.0;
194 
195  // Fill the triplets
196  leastSquareDerivatives(tripletList, grad);
197  _lastNTrip = tripletList.size();
198 
199  LOGLS_DEBUG(_log, "End of triplet filling, ntrip = " << tripletList.size());
200 
201  SparseMatrixD hessian = createHessian(_nTotal, tripletList);
202  tripletList.clear(); // we don't need it any more after we have the hessian.
203 
204  LOGLS_DEBUG(_log, "Starting factorization, hessian: dim="
205  << hessian.rows() << " non-zeros=" << hessian.nonZeros()
206  << " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
207 
208  if (dumpMatrixFile != "") {
209  if (hessian.rows() * hessian.cols() > 2e8) {
210  LOGLS_WARN(_log, "Hessian matrix is too big to dump to file, with rows, columns: "
211  << hessian.rows() << ", " << hessian.cols());
212  } else {
213  dumpMatrixAndGradient(hessian, grad, dumpMatrixFile, _log);
214  }
215  }
216 
218  if (chol.info() != Eigen::Success) {
219  LOGLS_ERROR(_log, "minimize: factorization failed ");
220  return MinimizeResult::Failed;
221  }
222 
223  std::size_t totalMeasOutliers = 0;
224  std::size_t totalRefOutliers = 0;
225  double oldChi2 = computeChi2().chi2;
226  double oldSigmaCut = 0;
227  double sigmaCut;
228 
229  while (true) {
230  Eigen::VectorXd delta = chol.solve(grad);
231  if (doLineSearch) {
232  scale = _lineSearch(delta);
233  }
234  offsetParams(scale * delta);
235  Chi2Statistic currentChi2(computeChi2());
236  LOGLS_DEBUG(_log, currentChi2);
237  if (!isfinite(currentChi2.chi2)) {
238  LOGL_ERROR(_log, "chi2 is not finite. Aborting outlier rejection.");
239  returnCode = MinimizeResult::NonFinite;
240  break;
241  }
242  if (currentChi2.chi2 > oldChi2 && totalMeasOutliers + totalRefOutliers != 0) {
243  LOGL_WARN(_log, "chi2 went up, skipping outlier rejection loop");
244  returnCode = MinimizeResult::Chi2Increased;
245  break;
246  }
247  oldChi2 = currentChi2.chi2;
248 
249  if (nSigmaCut == 0) break; // no rejection step to perform
250  MeasuredStarList msOutliers;
251  FittedStarList fsOutliers;
252  // keep nOutliers so we don't have to sum msOutliers.size()+fsOutliers.size() twice below.
253  std::size_t nOutliers = findOutliers(nSigmaCut, msOutliers, fsOutliers, sigmaCut);
254  double relChange = 1 - sigmaCut / oldSigmaCut;
255  LOGLS_DEBUG(_log, "findOutliers chi2 cut level: " << sigmaCut << ", relative change: " << relChange);
256  // If sigmaRelativeTolerance is set and at least one iteration has been done, break loop when the
257  // fractional change in sigmaCut levels is less than the sigmaRelativeTolerance parameter.
258  if ((sigmaRelativeTolerance > 0) && (oldSigmaCut > 0 && relChange < sigmaRelativeTolerance)) {
259  LOGLS_INFO(_log, "Iterations stopped with chi2 cut at " << sigmaCut << " and relative change of "
260  << relChange);
261  break;
262  }
263  totalMeasOutliers += msOutliers.size();
264  totalRefOutliers += fsOutliers.size();
265  oldSigmaCut = sigmaCut;
266  if (nOutliers == 0) break;
267  TripletList outlierTriplets(nOutliers);
268  grad.setZero(); // recycle the gradient
269  // compute the contributions of outliers to derivatives
270  outliersContributions(msOutliers, fsOutliers, outlierTriplets, grad);
271  // Remove significant outliers
272  removeMeasOutliers(msOutliers);
273  removeRefOutliers(fsOutliers);
274  if (doRankUpdate) {
275  // convert triplet list to eigen internal format
276  SparseMatrixD H(_nTotal, outlierTriplets.getNextFreeIndex());
277  H.setFromTriplets(outlierTriplets.begin(), outlierTriplets.end());
278  chol.update(H, false /* means downdate */);
279  // The contribution of outliers to the gradient is the opposite
280  // of the contribution of all other terms, because they add up to 0
281  grad *= -1;
282  } else {
283  // don't reuse tripletList because we want a new nextFreeIndex.
284  TripletList nextTripletList(_lastNTrip);
285  grad.setZero();
286  // Rebuild the matrix and gradient
287  leastSquareDerivatives(nextTripletList, grad);
288  _lastNTrip = nextTripletList.size();
289  LOGLS_DEBUG(_log, "Triplets recomputed, ntrip = " << nextTripletList.size());
290 
291  hessian = createHessian(_nTotal, nextTripletList);
292  nextTripletList.clear(); // we don't need it any more after we have the hessian.
293 
295  "Restarting factorization, hessian: dim="
296  << hessian.rows() << " non-zeros=" << hessian.nonZeros()
297  << " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
298  chol.compute(hessian);
299  if (chol.info() != Eigen::Success) {
300  LOGLS_ERROR(_log, "minimize: factorization failed ");
301  return MinimizeResult::Failed;
302  }
303  }
304  }
305 
306  if (totalMeasOutliers + totalRefOutliers > 0) {
307  _associations->cleanFittedStars();
308  }
309 
310  // only print the outlier summary if outlier rejection was turned on.
311  if (nSigmaCut != 0) {
312  LOGLS_INFO(_log, "Number of outliers (Measured + Reference = Total): "
313  << totalMeasOutliers << " + " << totalRefOutliers << " = "
314  << totalMeasOutliers + totalRefOutliers);
315  }
316  return returnCode;
317 }
Eigen::SparseMatrix< double, 0, Eigen::Index > SparseMatrixD
Definition: Eigenstuff.h:35
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition: Log.h:547
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
Definition: Log.h:679
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
Definition: Log.h:563
void leastSquareDerivatives(TripletList &tripletList, Eigen::VectorXd &grad) const
Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.
Definition: FitterBase.cc:344
void removeRefOutliers(FittedStarList &outliers)
Remove refStar outliers from the fit. No Refit done.
Definition: FitterBase.cc:338
Chi2Statistic computeChi2() const
Returns the chi2 for the current state.
Definition: FitterBase.cc:42
virtual void assignIndices(std::string const &whatToFit)=0
Set parameters to fit and assign indices in the big matrix.
virtual void offsetParams(Eigen::VectorXd const &delta)=0
Offset the parameters by the requested quantities.
void outliersContributions(MeasuredStarList &msOutliers, FittedStarList &fsOutliers, TripletList &tripletList, Eigen::VectorXd &grad)
Contributions to derivatives from (presumably) outlier terms.
Definition: FitterBase.cc:319
std::size_t findOutliers(double nSigmaCut, MeasuredStarList &msOutliers, FittedStarList &fsOutliers, double &cut) const
Find Measurements and references contributing more than a cut, computed as.
Definition: FitterBase.cc:52
void removeMeasOutliers(MeasuredStarList &outliers)
Remove measuredStar outliers from the fit. No Refit done.
Definition: FitterBase.cc:330
T isfinite(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
MinimizeResult
Return value of minimize()
Definition: FitterBase.h:40
T pow(T... args)

◆ offsetParams()

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

Offset the parameters by the requested quantities.

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

Parameters
[in]deltavector of offsets to apply

Implements lsst::jointcal::FitterBase.

Definition at line 211 of file PhotometryFit.cc.

211  {
212  if (delta.size() != _nTotal)
213  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
214  "PhotometryFit::offsetParams : the provided vector length is not compatible with "
215  "the current whatToFit setting");
216  if (_fittingModel) _photometryModel->offsetParams(delta);
217 
218  if (_fittingFluxes) {
219  for (auto &fittedStar : _associations->fittedStarList) {
220  // the parameter layout here is used also
221  // - when filling the derivatives
222  // - when assigning indices (assignIndices())
223  Eigen::Index index = fittedStar->getIndexInMatrix();
224  _photometryModel->offsetFittedStar(*fittedStar, delta(index));
225  }
226  }
227 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48

◆ operator=() [1/2]

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

◆ operator=() [2/2]

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

◆ outliersContributions()

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

Contributions to derivatives from (presumably) outlier terms.

No discarding done.

Definition at line 319 of file FitterBase.cc.

320  {
321  for (auto &outlier : msOutliers) {
322  MeasuredStarList tmp;
323  tmp.push_back(outlier);
324  const CcdImage &ccdImage = outlier->getCcdImage();
325  leastSquareDerivativesMeasurement(ccdImage, tripletList, grad, &tmp);
326  }
327  leastSquareDerivativesReference(fsOutliers, tripletList, grad);
328 }

◆ removeMeasOutliers()

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

Remove measuredStar outliers from the fit. No Refit done.

Definition at line 330 of file FitterBase.cc.

330  {
331  for (auto &measuredStar : outliers) {
332  auto fittedStar = measuredStar->getFittedStar();
333  measuredStar->setValid(false);
334  fittedStar->getMeasurementCount()--; // could be put in setValid
335  }
336 }

◆ removeRefOutliers()

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

Remove refStar outliers from the fit. No Refit done.

Definition at line 338 of file FitterBase.cc.

338  {
339  for (auto &fittedStar : outliers) {
340  fittedStar->setRefStar(nullptr);
341  }
342 }

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

352  {
353  std::string replaceStr = "{type}";
354  auto pos = baseName.find(replaceStr);
355  std::string measFilename(baseName);
356  measFilename.replace(pos, replaceStr.size(), "-meas.csv");
357  std::string refFilename(baseName);
358  refFilename.replace(pos, replaceStr.size(), "-ref.csv");
359  saveChi2MeasContributions(measFilename);
360  saveChi2RefContributions(refFilename);
361 }
virtual void saveChi2MeasContributions(std::string const &filename) const =0
Save a CSV file containing residuals of measurement terms.
virtual void saveChi2RefContributions(std::string const &filename) const =0
Save a CSV file containing residuals of reference terms.
T size(T... args)

◆ saveChi2MeasContributions()

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

Save a CSV file containing residuals of measurement terms.

Implements lsst::jointcal::FitterBase.

Definition at line 229 of file PhotometryFit.cc.

229  {
230  std::ofstream ofile(filename.c_str());
231  std::string separator = "\t";
232 
233  ofile << "id" << separator << "xccd" << separator << "yccd" << separator;
234  ofile << "mag" << separator << "instMag" << separator << "instMagErr" << separator;
235  ofile << "instFlux" << separator << "instFluxErr" << separator;
236  ofile << "inputFlux" << separator << "inputFluxErr" << separator;
237  ofile << "transformedFlux" << separator << "transformedFluxErr" << separator;
238  ofile << "fittedFlux" << separator;
239  ofile << "epoch" << separator;
240  ofile << "fsindex" << separator;
241  ofile << "ra" << separator << "dec" << separator;
242  ofile << "chi2" << separator << "nm" << separator;
243  ofile << "detector" << separator << "visit" << separator << std::endl;
244 
245  ofile << "id in source catalog" << separator << "coordinates in CCD" << separator << separator;
246  ofile << "fitted magnitude" << separator << "measured magnitude" << separator
247  << "measured magnitude error" << separator;
248  ofile << "measured instrumental flux (ADU)" << separator << "measured instrument flux error" << separator;
249  ofile << "measured flux (nJy)" << separator << "measured flux error" << separator;
250  ofile << separator << separator;
251  ofile << "fitted flux (nJy)" << separator;
252  ofile << "Julian Epoch Year of the measurement" << separator;
253  ofile << "unique index of the fittedStar" << separator;
254  ofile << "on-sky position of fitted star" << separator << separator;
255  ofile << "contribution to chi2 (1 dof)" << separator << "number of measurements of this FittedStar"
256  << separator;
257  ofile << "detector id" << separator << "visit id" << std::endl;
258 
259  const CcdImageList &ccdImageList = _associations->getCcdImageList();
260  for (auto const &ccdImage : ccdImageList) {
261  const MeasuredStarList &cat = ccdImage->getCatalogForFit();
262  for (auto const &measuredStar : cat) {
263  if (!measuredStar->isValid()) continue;
264 
265  double instFluxErr = _photometryModel->tweakFluxError(*measuredStar);
266  double flux = _photometryModel->transform(*ccdImage, *measuredStar);
267  double fluxErr = _photometryModel->transformError(*ccdImage, *measuredStar);
268  std::shared_ptr<FittedStar const> const fittedStar = measuredStar->getFittedStar();
269  double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
270  double chi2Val = std::pow(residual / fluxErr, 2);
271 
272  ofile << std::setprecision(9);
273  ofile << measuredStar->getId() << separator << measuredStar->x << separator << measuredStar->y
274  << separator;
275  ofile << fittedStar->getMag() << separator << measuredStar->getInstMag() << separator
276  << measuredStar->getInstMagErr() << separator;
277  ofile << measuredStar->getInstFlux() << separator << instFluxErr << separator;
278  ofile << measuredStar->getFlux() << separator << measuredStar->getFluxErr() << separator;
279  ofile << flux << separator << fluxErr << separator << fittedStar->getFlux() << separator;
280  ofile << ccdImage->getEpoch() << separator;
281  ofile << fittedStar->getIndexInMatrix() << separator;
282  ofile << fittedStar->x << separator << fittedStar->y << separator;
283  ofile << chi2Val << separator << fittedStar->getMeasurementCount() << separator;
284  ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() << std::endl;
285  } // loop on measurements in image
286  } // loop on images
287 }
T endl(T... args)
std::list< std::shared_ptr< CcdImage > > CcdImageList
Definition: CcdImage.h:46
T setprecision(T... args)

◆ saveChi2RefContributions()

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

Save a CSV file containing residuals of reference terms.

Implements lsst::jointcal::FitterBase.

Definition at line 289 of file PhotometryFit.cc.

289  {
290  std::ofstream ofile(filename.c_str());
291  std::string separator = "\t";
292 
293  ofile << "ra" << separator << "dec " << separator;
294  ofile << "mag" << separator;
295  ofile << "refFlux" << separator << "refFluxErr" << separator;
296  ofile << "fittedFlux" << separator << "fittedFluxErr" << separator;
297  ofile << "fsindex" << separator << "chi2" << separator << "nm" << std::endl;
298 
299  ofile << "coordinates of fittedStar" << separator << separator;
300  ofile << "magnitude" << separator;
301  ofile << "refStar flux (nJy)" << separator << "refStar fluxErr" << separator;
302  ofile << "fittedStar flux (nJy)" << separator << "fittedStar fluxErr" << separator;
303  ofile << "unique index of the fittedStar" << separator << "refStar contribution to chi2 (1 dof)"
304  << separator << "number of measurements of this FittedStar" << std::endl;
305 
306  // The following loop is heavily inspired from PhotometryFit::computeChi2()
307  const FittedStarList &fittedStarList = _associations->fittedStarList;
308  for (auto const &fittedStar : fittedStarList) {
309  const RefStar *refStar = fittedStar->getRefStar();
310  if (refStar == nullptr) continue;
311 
312  double chi2 = std::pow(((fittedStar->getFlux() - refStar->getFlux()) / refStar->getFluxErr()), 2);
313 
314  ofile << std::setprecision(9);
315  ofile << fittedStar->x << separator << fittedStar->y << separator;
316  ofile << fittedStar->getMag() << separator;
317  ofile << refStar->getFlux() << separator << refStar->getFluxErr() << separator;
318  ofile << fittedStar->getFlux() << separator << fittedStar->getFluxErr() << separator;
319  ofile << fittedStar->getIndexInMatrix() << separator << chi2 << separator
320  << fittedStar->getMeasurementCount() << std::endl;
321  } // loop on FittedStars
322 }

Member Data Documentation

◆ _associations

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

Definition at line 163 of file FitterBase.h.

◆ _lastNTrip

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

Definition at line 166 of file FitterBase.h.

◆ _log

LOG_LOGGER lsst::jointcal::FitterBase::_log
protectedinherited

Definition at line 172 of file FitterBase.h.

◆ _nModelParams

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

Definition at line 168 of file FitterBase.h.

◆ _nStarParams

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

Definition at line 169 of file FitterBase.h.

◆ _nTotal

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

Definition at line 167 of file FitterBase.h.

◆ _whatToFit

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

Definition at line 164 of file FitterBase.h.


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