Skip to content

Commit

Permalink
Merge pull request #102 from lsst/tickets/DM-14786
Browse files Browse the repository at this point in the history 
DM-14786 implement ConstrainedMagnitudeModel
  • Loading branch information
@parejkoj
parejkoj committed Aug 23, 2018
2 parents 46f1b7d + b240722 commit 46b432b
Show file tree
Hide file tree
Showing 22 changed files with 1,098 additions and 421 deletions.
6 changes: 4 additions & 2 deletions bin.src/plot_photoCalib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,12 +120,14 @@ def makePhotoCalibImages(visit, butler, step=8, chips=[], tract=None,

calibScaling = 1
if singleCalib:
calib = butler.get('calexp_calib', dataId=dict(visit=int(visit), ccd=int(ccd), tract=0))
calib = butler.get('calexp_calib', dataId=dict(visit=int(visit), ccd=int(ccd), tract=tract))
calibScaling = calib.getFluxMag0()[0]
if verbose:
print('calib ccd %s: %s' % (ccd, calibScaling))

photoCalib = butler.get('jointcal_photoCalib', dataId=dict(visit=int(visit), ccd=int(ccd), tract=0))
photoCalib = butler.get('jointcal_photoCalib', dataId=dict(visit=int(visit),
ccd=int(ccd),
tract=tract))
if verbose:
print("photoCalib mean ccd %s: %s" % (ccd, photoCalib.getCalibrationMean()))
scaled = photoCalib.computeScaledCalibration()
Expand Down
10 changes: 7 additions & 3 deletions include/lsst/jointcal/AstrometryModel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,13 @@ class AstrometryModel {
//! assigned index.
virtual unsigned assignIndices(std::string const &whatToFit, unsigned firstIndex) = 0;

//! Offset the parameters by the provided amounts.
/*! The shifts are applied according to the indices given in
AssignIndices. */
/**
* Offset the parameters by the provided amounts (by -delta).
*
* The shifts are applied according to the indices given in assignIndices.
*
* @param[in] delta vector of offsets to apply
*/
virtual void offsetParams(Eigen::VectorXd const &delta) = 0;

//! The transformation used to project the positions of FittedStars.
Expand Down
130 changes: 101 additions & 29 deletions include/lsst/jointcal/ConstrainedPhotometryModel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,10 @@ class ConstrainedPhotometryModel : public PhotometryModel {
* @param[in] visitOrder The order of the visit polynomial.
*/
explicit ConstrainedPhotometryModel(CcdImageList const &ccdImageList,
afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7);
afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7)
: _fittingChips(false), _fittingVisits(false) {
_chipVisitMap.reserve(ccdImageList.size());
}

/// No copy or move: there is only ever one instance of a given model (i.e. per ccd+visit)
ConstrainedPhotometryModel(ConstrainedPhotometryModel const &) = delete;
Expand All @@ -47,18 +50,6 @@ class ConstrainedPhotometryModel : public PhotometryModel {
/// @copydoc PhotometryModel::offsetParams
void offsetParams(Eigen::VectorXd const &delta) override;

/// @copydoc PhotometryModel::offsetFittedStar
void offsetFittedStar(FittedStar &fittedStar, double delta) const { fittedStar.getFlux() -= delta; }

/// @copydoc PhotometryModel::computeResidual
double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::transform
double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::transformError
double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::freezeErrorTransform
void freezeErrorTransform() override;

Expand All @@ -72,27 +63,12 @@ class ConstrainedPhotometryModel : public PhotometryModel {
void computeParameterDerivatives(MeasuredStar const &measuredStar, CcdImage const &ccdImage,
Eigen::VectorXd &derivatives) const override;

/// @copydoc PhotometryModel::getRefError
double getRefError(RefStar const &refStar) const override { return refStar.getFluxErr(); }

/// @copydoc PhotometryModel::computeRefResidual
double computeRefResidual(FittedStar const &fittedStar, RefStar const &refStar) const override {
return fittedStar.getFlux() - refStar.getFlux();
};

/// @copydoc PhotometryModel::toPhotoCalib
std::shared_ptr<afw::image::PhotoCalib> toPhotoCalib(CcdImage const &ccdImage) const override;

/// @copydoc PhotometryModel::dump
void dump(std::ostream &stream = std::cout) const override;

private:
protected:
PhotometryMappingBase *findMapping(CcdImage const &ccdImage) const override;

// Which components of the model are we fitting currently?
bool _fittingChips;
bool _fittingVisits;

/* The per-ccdImage transforms, each of which is a composition of a chip and visit transform.
* Not all pairs of _visitMap[visit] and _chipMap[chip] are guaranteed to have an entry in
* _chipVisitMap (for example, one ccd in one visit might fail to produce a catalog).
Expand All @@ -106,6 +82,102 @@ class ConstrainedPhotometryModel : public PhotometryModel {
// The per-chip transforms that go into _chipVisitMap.
typedef std::map<CcdIdType, std::shared_ptr<PhotometryMapping>> ChipMapType;
ChipMapType _chipMap;

/**
* Initialize the chip, visit, and chipVisit mappings by creating appropriate transfos and mappings.
*/
template <class ChipTransfo, class VisitTransfo, class ChipVisitMapping>
void initialize(CcdImageList const &ccdImageList, afw::geom::Box2D const &focalPlaneBBox, int visitOrder);

/// Return the initial calibration to use from this photoCalib.
virtual double initialChipCalibration(std::shared_ptr<afw::image::PhotoCalib const> photoCalib) = 0;

private:
// Which components of the model are we fitting currently?
bool _fittingChips;
bool _fittingVisits;
};

class ConstrainedFluxModel : public ConstrainedPhotometryModel {
public:
explicit ConstrainedFluxModel(CcdImageList const &ccdImageList, afw::geom::Box2D const &focalPlaneBBox,
int visitOrder = 7)
: ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder) {
initialize<FluxTransfoSpatiallyInvariant, FluxTransfoChebyshev, ChipVisitFluxMapping>(
ccdImageList, focalPlaneBBox, visitOrder);
}

/// @copydoc PhotometryModel::offsetFittedStar
void offsetFittedStar(FittedStar &fittedStar, double delta) const override {
fittedStar.getFlux() -= delta;
}

/// @copydoc PhotometryModel::computeResidual
double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::computeRefResidual
double computeRefResidual(FittedStar const &fittedStar, RefStar const &refStar) const override {
return fittedStar.getFlux() - refStar.getFlux();
};

/// @copydoc PhotometryModel::getRefError
double getRefError(RefStar const &refStar) const override { return refStar.getFluxErr(); }

/// @copydoc PhotometryModel::transform
double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::transformError
double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::toPhotoCalib
std::shared_ptr<afw::image::PhotoCalib> toPhotoCalib(CcdImage const &ccdImage) const override;

protected:
/// @copydoc ConstrainedPhotometryModel::initialChipCalibration
double initialChipCalibration(std::shared_ptr<afw::image::PhotoCalib const> photoCalib) {
return photoCalib->getCalibrationMean();
}
};

class ConstrainedMagnitudeModel : public ConstrainedPhotometryModel {
public:
explicit ConstrainedMagnitudeModel(CcdImageList const &ccdImageList,
afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7)
: ConstrainedPhotometryModel(ccdImageList, focalPlaneBBox, visitOrder) {
initialize<MagnitudeTransfoSpatiallyInvariant, MagnitudeTransfoChebyshev, ChipVisitMagnitudeMapping>(
ccdImageList, focalPlaneBBox, visitOrder);
}

/// @copydoc PhotometryModel::offsetFittedStar
void offsetFittedStar(FittedStar &fittedStar, double delta) const override {
fittedStar.getMag() -= delta;
}

/// @copydoc PhotometryModel::computeResidual
double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::computeRefResidual
double computeRefResidual(FittedStar const &fittedStar, RefStar const &refStar) const override {
return fittedStar.getMag() - refStar.getMag();
};

/// @copydoc PhotometryModel::getRefError
double getRefError(RefStar const &refStar) const override { return refStar.getMagErr(); }

/// @copydoc PhotometryModel::transform
double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::transformError
double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override;

/// @copydoc PhotometryModel::toPhotoCalib
std::shared_ptr<afw::image::PhotoCalib> toPhotoCalib(CcdImage const &ccdImage) const override;

protected:
/// @copydoc ConstrainedPhotometryModel::initialChipCalibration
double initialChipCalibration(std::shared_ptr<afw::image::PhotoCalib const> photoCalib) {
return magFromFlux(photoCalib->getCalibrationMean());
}
};

} // namespace jointcal
Expand Down
3 changes: 2 additions & 1 deletion include/lsst/jointcal/FitterBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@ namespace jointcal {
enum class MinimizeResult {
Converged, // fit has converged - no more outliers
Chi2Increased, // still some ouliers but chi2 increases
Failed // factorization failed
Failed, // factorization failed
NonFinite // non-finite chi2 statistic
};

/**
Expand Down
76 changes: 52 additions & 24 deletions include/lsst/jointcal/PhotometryMapping.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,21 +36,21 @@ class PhotometryMappingBase {
* Return the on-sky transformed flux for measuredStar on ccdImage.
*
* @param[in] measuredStar The measured star position to transform.
* @param[in] instFlux The instrument flux to transform.
* @param[in] value The instrument flux or magnitude to transform.
*
* @return The on-sky flux transformed from instFlux at measuredStar's position.
* @return The on-sky value () transformed from value at measuredStar's position.
*/
virtual double transform(MeasuredStar const &measuredStar, double instFlux) const = 0;
virtual double transform(MeasuredStar const &measuredStar, double value) const = 0;

/**
* Return the on-sky transformed flux uncertainty for measuredStar on ccdImage.
* Matches the underlying PhotometryTransfo's `transformError()` until freezeErrorTransform() is called.
*
* @param[in] measuredStar The measured star position to transform.
* @param[in] value The value to transform.
* @param[in] valueErr The value uncertainty to transform.
* @param[in] value The flux or magnitude to transform.
* @param[in] valueErr The flux or magnitude uncertainty to transform.
*
* @return The on-sky flux transformed from instFlux at measuredStar's position.
* @return The on-sky value transformed from value at measuredStar's position.
*/
virtual double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const = 0;

Expand All @@ -67,10 +67,10 @@ class PhotometryMappingBase {
* Compute the derivatives with respect to the parameters (i.e. the coefficients).
*
* @param[in] measuredStar The measured star position to transform.
* @param[in] instFlux The instrument flux to compute the derivative at.
* @param[in] value The instrument flux or magnitude to compute the derivative at.
* @param[out] derivatives The computed derivatives, in the same order as the deltas in offsetParams.
*/
virtual void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
virtual void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
Eigen::Ref<Eigen::VectorXd> derivatives) const = 0;

/// Make this mapping's parameters fixed (i.e. not varied during fitting).
Expand Down Expand Up @@ -192,35 +192,27 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
ChipVisitPhotometryMapping(std::shared_ptr<PhotometryMapping> chipMapping,
std::shared_ptr<PhotometryMapping> visitMapping)
: PhotometryMappingBase(),
_nParChips(0),
_nParVisits(0),
_nParChip(0),
_nParVisit(0),
_chipMapping(std::move(chipMapping)),
_visitMapping(std::move(visitMapping)) {}

/// @copydoc PhotometryMappingBase::getNpar
unsigned getNpar() const override { return _nParChips + _nParVisits; }
unsigned getNpar() const override { return _nParChip + _nParVisit; }

/// @copydoc PhotometryMappingBase::transform
double transform(MeasuredStar const &measuredStar, double instFlux) const override {
double tempFlux = _chipMapping->getTransfo()->transform(measuredStar.x, measuredStar.y, instFlux);
double transform(MeasuredStar const &measuredStar, double value) const override {
double temp = _chipMapping->getTransfo()->transform(measuredStar.x, measuredStar.y, value);
return _visitMapping->getTransfo()->transform(measuredStar.getXFocal(), measuredStar.getYFocal(),
tempFlux);
temp);
}

/// @copydoc PhotometryMappingBase::transformError
double transformError(MeasuredStar const &measuredStar, double instFlux,
double instFluxErr) const override;

/// @copydoc PhotometryMappingBase::freezeErrorTransform
void freezeErrorTransform() override {
_chipMapping->freezeErrorTransform();
_visitMapping->freezeErrorTransform();
}

/// @copydoc PhotometryMappingBase::computeParameterDerivatives
void computeParameterDerivatives(MeasuredStar const &measuredStar, double instFlux,
Eigen::Ref<Eigen::VectorXd> derivatives) const override;

/// @copydoc PhotometryMappingBase::getParameters
Eigen::VectorXd getParameters() override {
Eigen::VectorXd joined(getNpar());
Expand Down Expand Up @@ -253,15 +245,51 @@ class ChipVisitPhotometryMapping : public PhotometryMappingBase {
std::shared_ptr<PhotometryMapping> getChipMapping() const { return _chipMapping; }
std::shared_ptr<PhotometryMapping> getVisitMapping() const { return _visitMapping; }

private:
unsigned getNParChip() const { return _nParChip; }
unsigned getNParVisit() const { return _nParVisit; }

protected:
// These are either transfo.getNpar() or 0, depending on whether we are fitting that component or not.
unsigned _nParChips, _nParVisits;
unsigned _nParChip, _nParVisit;

// the actual transformation to be fit
std::shared_ptr<PhotometryMapping> _chipMapping;
std::shared_ptr<PhotometryMapping> _visitMapping;
};

class ChipVisitFluxMapping : public ChipVisitPhotometryMapping {
public:
ChipVisitFluxMapping(std::shared_ptr<PhotometryMapping> chipMapping,
std::shared_ptr<PhotometryMapping> visitMapping)
: ChipVisitPhotometryMapping(chipMapping, visitMapping) {}

/// @copydoc PhotometryMappingBase::transformError
double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override;

/// @copydoc PhotometryMappingBase::computeParameterDerivatives
void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
Eigen::Ref<Eigen::VectorXd> derivatives) const override;
};

class ChipVisitMagnitudeMapping : public ChipVisitPhotometryMapping {
public:
ChipVisitMagnitudeMapping(std::shared_ptr<PhotometryMapping> chipMapping,
std::shared_ptr<PhotometryMapping> visitMapping)
: ChipVisitPhotometryMapping(chipMapping, visitMapping) {}

/**
* @copydoc PhotometryMappingBase::transformError
*
* @note This method takes instFlux and instFluxErr: the error calculation has to
* use fluxes to get the math right.
*/
double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override;

/// @copydoc PhotometryMappingBase::computeParameterDerivatives
void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
Eigen::Ref<Eigen::VectorXd> derivatives) const override;
};

} // namespace jointcal
} // namespace lsst

Expand Down
16 changes: 14 additions & 2 deletions include/lsst/jointcal/PhotometryModel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ class PhotometryModel {
virtual unsigned assignIndices(std::string const &whatToFit, unsigned firstIndex) = 0;

/**
* Offset the parameters by the provided amounts.
* Offset the parameters by the provided amounts (by -delta).
*
* The shifts are applied according to the indices given in assignIndices.
*
Expand All @@ -43,13 +43,25 @@ class PhotometryModel {
virtual void offsetParams(Eigen::VectorXd const &delta) = 0;

/**
* Offset the appropriate flux or magnitude.
* Offset the appropriate flux or magnitude (by -delta).
*
* @param fittedStar The star to update.
* @param delta The amount to update by.
*/
virtual void offsetFittedStar(FittedStar &fittedStar, double delta) const = 0;

/**
* Compute the residual between the model applied to a star and its associated fittedStar.
*
* @f[
* residual = Model(measuredStar) - fittedStar
* @f]
*
* @param ccdImage The ccdImage where measuredStar resides.
* @param measuredStar The measured star position to compute the residual of.
*
* @return The residual.
*/
virtual double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const = 0;

/**
Expand Down

0 comments on commit 46b432b

Please sign in to comment.