Skip to content

Commit

Permalink
Remove refraction and color code
Browse files Browse the repository at this point in the history 
The refraction code here has not been used in ages, and we are not
supplying relevant color information as input anyway.
  • Loading branch information
@parejkoj
parejkoj committed Jun 14, 2021
1 parent 5624b44 commit 102ebf1
Show file tree
Hide file tree
Showing 6 changed files with 19 additions and 143 deletions.
7 changes: 0 additions & 7 deletions include/lsst/jointcal/Associations.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,13 +161,6 @@ class Associations {
//! track of that.
void deprojectFittedStars();

//! Set the color field of FittedStar 's from a colored catalog.
/* If Color is "g-i", then the color is assigned from columns "g" and "i" of the colored catalog. */
#ifdef TODO
void setFittedStarColors(std::string const &dicStarListName, std::string const &color,
double matchCutArcSec);
#endif

/**
* Prepare the fittedStar list by making quality cuts and normalizing measurements.
*
Expand Down
14 changes: 3 additions & 11 deletions include/lsst/jointcal/AstrometryFit.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ class AstrometryFit : public FitterBase {
* Set parameters to fit and assign indices in the big matrix.
*
* @param[in] whatToFit Valid strings: zero or more of "Distortions", "Positions",
* "Refrac", "PM" which define which parameter set
* "PM" which define which parameter set
* is going to be variable when computing
* derivatives (leastSquareDerivatives) and minimizing
* (minimize()). whatToFit="Positions Distortions"
Expand Down Expand Up @@ -134,14 +134,8 @@ class AstrometryFit : public FitterBase {
void saveChi2RefContributions(std::string const &filename) const override;

private:
bool _fittingDistortions, _fittingPos, _fittingRefrac, _fittingPM;
bool _fittingDistortions, _fittingPos, _fittingPM;
std::shared_ptr<AstrometryModel> _astrometryModel;
double _referenceColor, _sigCol; // average and r.m.s color
double _refractionCoefficient; // fit parameter
Eigen::Index _refracPosInMatrix; // where it stands

// counts in parameter subsets.
std::size_t _nParRefrac;

double _epoch; // epoch to correct proper motion/parallax to (Julian Epoch year, e.g. J2000.0)
double _posError; // constant term on error on position (in pixel unit)
Expand All @@ -164,13 +158,11 @@ class AstrometryFit : public FitterBase {
*
* @param fittedStar The star to transform.
* @param sky2TP Transformation from sky coordinates to CcdImage tangent plane.
* @param refractionVector On-sky vector to refraction correct position.
* @param refractionCoeff Amount of refraction correction to apply.
* @param deltaYears Difference in years between FittedStar and MeasuredStar epochs.
* @return Corrected position of FittedStar.
*/
Point transformFittedStar(FittedStar const &fittedStar, AstrometryTransform const &sky2TP,
Point const &refractionVector, double refractionCoeff, double deltaYears) const;
double deltaYears) const;

/// Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) from one CcdImage.
void accumulateStatImage(CcdImage const &ccdImage, Chi2Accumulator &accum) const;
Expand Down
2 changes: 0 additions & 2 deletions include/lsst/jointcal/FittedStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,8 +100,6 @@ class FittedStar : public BaseStar {
//! Get the astrometric reference star associated with this star.
const RefStar* getRefStar() const { return _refStar; };

double color; // TODO: remove this holdover from PmBlock class

private:
Eigen::Index _indexInMatrix;
int _measurementCount;
Expand Down
43 changes: 0 additions & 43 deletions src/Associations.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -523,49 +523,6 @@ void Associations::collectMCStars(int realization) {
}
}

void Associations::setFittedStarColors(std::string dicStarListName, std::string color,
double matchCutArcSec) {
// decode color string in case it is x-y
size_t pos_minus = color.find('-');
bool compute_diff = (pos_minus != string::npos);
std::string c1, c2;
c1 = color.substr(0, pos_minus); // if pos_minus == npos, means "up to the end"
if (compute_diff) c2 = color.substr(pos_minus + 1, string::npos);
DicStarList cList(dicStarListName);
if (!cList.HasKey(c1))
throw(GastroException("Associations::SetFittedstarColors : " + dicStarListName +
" misses a key named \"" + c1 + "\""));
if (compute_diff && !cList.HasKey(c2))
throw(GastroException("Associations::SetFittedstarColors : " + dicStarListName +
" misses a key named \"" + c2 + "\""));
// we associate in some tangent plane. The reference catalog is expressed on the sky,
// but FittedStar's may be still in this tangent plane.
BaseStarList &l1 = (BaseStarList &)fittedStarList;
AstrometryTransformIdentity id;
TanRaDecToPixel proj(AstrometryTransformLinear(), getCommonTangentPoint());
// project or not ?
AstrometryTransform *id_or_proj = &proj;
if (fittedStarList.inTangentPlaneCoordinates) id_or_proj = &id;
// The color List is to be projected:
TStarList projected_cList((BaseStarList &)cList, proj);
// Associate
auto starMatchList = listMatchCollect(Fitted2Base(fittedStarList), (const BaseStarList &)projected_cList,
id_or_proj, matchCutArcSec / 3600);

LOGLS_INFO(_log, "Matched " << starMatchList->size() << '/' << fittedStarList.size()
<< " FittedStars to color catalog");
// Evaluate and assign colors.
for (auto i = starMatchList->begin(); i != starMatchList->end(); ++i) {
BaseStar *s1 = i->s1;
FittedStar *fs = dynamic_cast<FittedStar *>(s1);
BaseStar *s2 = i->s2;
const TStar *ts = dynamic_cast<const TStar *>(s2);
const DicStar *ds = dynamic_cast<const DicStar *>(ts->get_original());
fs->color = ds->getval(c1);
if (compute_diff) fs->color -= ds->getval(c2);
}
}

#endif /* TODO */
} // namespace jointcal
} // namespace lsst
84 changes: 10 additions & 74 deletions src/AstrometryFit.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,31 +64,15 @@ AstrometryFit::AstrometryFit(std::shared_ptr<Associations> associations,
std::shared_ptr<AstrometryModel> astrometryModel, double posError)
: FitterBase(associations),
_astrometryModel(astrometryModel),
_refractionCoefficient(0),
_nParRefrac(_associations->getNFilters()),
_epoch(_associations->getEpoch()),
_posError(posError) {
_log = LOG_GET("jointcal.AstrometryFit");
_referenceColor = 0;
_sigCol = 0;
std::size_t count = 0;
for (auto const &i : _associations->fittedStarList) {
_referenceColor += i->color;
_sigCol += std::pow(i->color, 2);
count++;
}
if (count) {
_referenceColor /= double(count);
if (_sigCol > 0) _sigCol = sqrt(_sigCol / count - std::pow(_referenceColor, 2));
}
LOGLS_INFO(_log, "Reference Color: " << _referenceColor << " sig " << _sigCol);
}

/* ! this routine is used in 3 instances: when computing
the derivatives, when computing the Chi2, when filling a tuple.
*/
Point AstrometryFit::transformFittedStar(FittedStar const &fittedStar, AstrometryTransform const &sky2TP,
Point const &refractionVector, double refractionCoeff,
double deltaYears) const {
Point fittedStarInTP;
if (fittedStar.getRefStar()) {
Expand All @@ -98,13 +82,6 @@ Point AstrometryFit::transformFittedStar(FittedStar const &fittedStar, Astrometr
} else {
fittedStarInTP = sky2TP.apply(fittedStar);
}

// account for atmospheric refraction: does nothing if color
// have not been assigned
// the color definition shouldbe the same when computing derivatives
double color = fittedStar.color - _referenceColor;
fittedStarInTP.x += refractionVector.x * color * refractionCoeff;
fittedStarInTP.y += refractionVector.y * color * refractionCoeff;
return fittedStarInTP;
}

Expand Down Expand Up @@ -143,9 +120,8 @@ void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage const &ccdImage,
// count parameters
std::size_t npar_mapping = (_fittingDistortions) ? mapping->getNpar() : 0;
std::size_t npar_pos = (_fittingPos) ? 2 : 0;
std::size_t npar_refrac = (_fittingRefrac) ? 1 : 0;
std::size_t npar_pm = (_fittingPM) ? 2 : 0;
std::size_t npar_tot = npar_mapping + npar_pos + npar_refrac + npar_pm;
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;
Expand All @@ -154,8 +130,6 @@ void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage const &ccdImage,

// FittedStar is "observed" epoch, MeasuredStar is "baseline"
double deltaYears = _epoch - ccdImage.getEpoch();
// refraction stuff
Point refractionVector = ccdImage.getRefractionVector();
// transformation from sky to TP
auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
// reserve matrices once for all measurements
Expand Down Expand Up @@ -206,8 +180,7 @@ void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage const &ccdImage,

std::shared_ptr<FittedStar const> const fs = ms.getFittedStar();

Point fittedStarInTP =
transformFittedStar(*fs, *sky2TP, refractionVector, _refractionCoefficient, deltaYears);
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
Expand All @@ -233,16 +206,6 @@ void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage const &ccdImage,
// indices[ipar + 1] = fs->getIndexInMatrix() + 3;
// ipar += npar_pm;
// }
if (_fittingRefrac) {
/* if the definition of color changes, it has to remain
consistent with transformFittedStar */
double color = fs->color - _referenceColor;
// sign checked
H(ipar, 0) = -refractionVector.x * color;
H(ipar, 1) = -refractionVector.y * color;
indices[ipar] = _refracPosInMatrix;
ipar += 1;
}

// We can now compute the residual
Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
Expand Down Expand Up @@ -367,8 +330,6 @@ void AstrometryFit::accumulateStatImage(CcdImage const &ccdImage, Chi2Accumulato
const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
// FittedStar is "observed" epoch, MeasuredStar is "baseline"
double deltaYears = _epoch - ccdImage.getEpoch();
// refraction stuff
Point refractionVector = ccdImage.getRefractionVector();
// transformation from sky to TP
auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
// reserve matrix once for all measurements
Expand All @@ -395,8 +356,7 @@ void AstrometryFit::accumulateStatImage(CcdImage const &ccdImage, Chi2Accumulato
transW(0, 1) = transW(1, 0) = -outPos.vxy / det;

std::shared_ptr<FittedStar const> const fs = ms->getFittedStar();
Point fittedStarInTP =
transformFittedStar(*fs, *sky2TP, refractionVector, _refractionCoefficient, deltaYears);
Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);

Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
double chi2Val = res.transpose() * transW * res;
Expand Down Expand Up @@ -461,12 +421,6 @@ void AstrometryFit::assignIndices(std::string const &whatToFit) {
LOGLS_INFO(_log, "assignIndices: Now fitting " << whatToFit);
_fittingDistortions = (_whatToFit.find("Distortions") != std::string::npos);
_fittingPos = (_whatToFit.find("Positions") != std::string::npos);
_fittingRefrac = (_whatToFit.find("Refrac") != std::string::npos);
if (_sigCol == 0 && _fittingRefrac) {
LOGLS_WARN(_log,
"Cannot fit refraction coefficients without a color lever arm. Ignoring refraction.");
_fittingRefrac = false;
}
_fittingPM = (_whatToFit.find("PM") != std::string::npos);
// When entering here, we assume that whatToFit has already been interpreted.

Expand All @@ -488,10 +442,6 @@ void AstrometryFit::assignIndices(std::string const &whatToFit) {
}
}
_nStarParams = ipar - _nModelParams;
if (_fittingRefrac) {
_refracPosInMatrix = ipar;
ipar += _nParRefrac;
}
_nTotal = ipar;
LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
<< " values: " << _nStarParams);
Expand Down Expand Up @@ -521,21 +471,9 @@ void AstrometryFit::offsetParams(Eigen::VectorXd const &delta) {
// }
}
}
if (_fittingRefrac) {
_refractionCoefficient += delta(_refracPosInMatrix);
}
}

void AstrometryFit::checkStuff() {
#if (0)
const char *what2fit[] = {"Positions",
"Distortions",
"Refrac",
"Positions Distortions",
"Positions Refrac",
"Distortions Refrac",
"Positions Distortions Refrac"};
#endif
const char *what2fit[] = {"Positions", "Distortions", "Positions Distortions"};
// DEBUG
for (unsigned k = 0; k < sizeof(what2fit) / sizeof(what2fit[0]); ++k) {
Expand All @@ -562,7 +500,7 @@ void AstrometryFit::saveChi2MeasContributions(std::string const &filename) const
ofile << "xErr" << separator << "yErr" << separator << "xyCov" << separator;
ofile << "xtpi" << separator << "ytpi" << separator;
ofile << "rxi" << separator << "ryi" << separator;
ofile << "color" << separator << "fsindex" << separator;
ofile << "fsindex" << separator;
ofile << "ra" << separator << "dec" << separator;
ofile << "chi2" << separator << "nm" << separator;
ofile << "chip" << separator << "visit" << std::endl;
Expand All @@ -574,7 +512,7 @@ void AstrometryFit::saveChi2MeasContributions(std::string const &filename) const
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 << "currently unused" << separator << "unique index of the fittedStar" << 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;
Expand All @@ -585,7 +523,6 @@ void AstrometryFit::saveChi2MeasContributions(std::string const &filename) const
const MeasuredStarList &cat = ccdImage->getCatalogForFit();
const AstrometryMapping *mapping = _astrometryModel->getMapping(*ccdImage);
const auto readTransform = ccdImage->getReadWcs();
const Point &refractionVector = ccdImage->getRefractionVector();
// FittedStar is "observed" epoch, MeasuredStar is "baseline"
double deltaYears = _epoch - ccdImage->getEpoch();
for (auto const &ms : cat) {
Expand All @@ -600,8 +537,7 @@ void AstrometryFit::saveChi2MeasContributions(std::string const &filename) const
FatPoint inputTpPos = readPixToTangentPlane->apply(inPos);
std::shared_ptr<FittedStar const> const fs = ms->getFittedStar();

Point fittedStarInTP =
transformFittedStar(*fs, *sky2TP, refractionVector, _refractionCoefficient, deltaYears);
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);
Expand All @@ -618,7 +554,7 @@ void AstrometryFit::saveChi2MeasContributions(std::string const &filename) const
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->color << separator << fs->getIndexInMatrix() << 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;
Expand All @@ -634,14 +570,14 @@ void AstrometryFit::saveChi2RefContributions(std::string const &filename) const
ofile << "rx" << separator << "ry" << separator;
ofile << "mag" << separator;
ofile << "xErr" << separator << "yErr" << separator << "xyCov" << separator;
ofile << "color" << separator << "fsindex" << 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 << "currently unused" << separator << "unique index of the fittedStar" << separator;
ofile << "unique index of the fittedStar" << separator;
ofile << "refStar contribution to Chi2 (2D dofs)" << separator
<< "number of measurements of this FittedStar" << std::endl;

Expand All @@ -663,7 +599,7 @@ void AstrometryFit::saveChi2RefContributions(std::string const &filename) const
ofile << rsProj.x << separator << rsProj.y << separator;
ofile << fs.getMag() << separator;
ofile << rsProj.vx << separator << rsProj.vy << separator << rsProj.vxy << separator;
ofile << fs.color << separator << fs.getIndexInMatrix() << separator;
ofile << fs.getIndexInMatrix() << separator;
ofile << chi2 << separator << fs.getMeasurementCount() << std::endl;
} // loop on FittedStars
}
Expand Down

0 comments on commit 102ebf1

Please sign in to comment.