Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-40652: Make SVD calculation optional #16

Merged
merged 2 commits into from
Dec 9, 2023
Merged

DM-40652: Make SVD calculation optional #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b137222
Make SVD calculation on failure optional
cmsaunders Dec 4, 2023
8963d62
Purge bad exposures after clipping
cmsaunders Dec 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion include/Match.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -325,8 +325,11 @@ class CoordAlign {
// Re-map Detections using current params - either all or just those being fit.
void remap(bool doAll = true);
// Fitting routine: returns chisq of current fit, updates params. and remaps.
// If Cholesky decomposition fails and doSVD=true, do singular value decomposition to report nature of
// degeneracies.
double fitOnce(bool reportToCerr = true,
bool inPlace = false); // Set inPlace to save space, but can't debug singularities
bool inPlace = false, // Set inPlace to save space, but can't debug singularities
bool doSVD = false);
// Conduct one round of sigma-clipping. If doReserved=true,
// then only clip reserved Matches. If =false, then
// only clip non-reserved Matches.
Expand Down
2 changes: 1 addition & 1 deletion include/WCSFitRoutine.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ class WCSFit {
void fit(double maxError = 100., int minFitExposures = 200, double reserveFraction = 0.2,
int randomNumberSeed = 1234, double minimumImprovement = 0.02, double clipThresh = 5.0,
double chisqTolerance = 0.001, bool clipEntireMatch = false, bool divideInPlace = false,
bool purgeOutput = false, double minColor = -10.0, double maxColor = 10.0);
bool purgeOutput = false, double minColor = -10.0, double maxColor = 10.0, bool calcSVD = false);

Astro::outputCatalog getOutputCatalog();

Expand Down
2 changes: 1 addition & 1 deletion pydir/wcsfit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@ PYBIND11_MODULE(wcsfit, m) {
py::arg("minimumImprovement") = 0.02, py::arg("clipThresh") = 5.0,
py::arg("chisqTolerance") = 0.001, py::arg("clipEntireMatch") = false,
py::arg("divideInPlace") = false, py::arg("purgeOutput") = false,
py::arg("minColor") = -10.0, py::arg("maxColor") = 10.0)
py::arg("minColor") = -10.0, py::arg("maxColor") = 10.0, py::arg("calcSVD") = false)
.def("saveResults", &WCSFit::saveResults)
.def("getOutputCatalog", [](WCSFit &self) {
Astro::outputCatalog outCat = self.getOutputCatalog();
Expand Down
128 changes: 65 additions & 63 deletions src/subs/Match.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1168,7 +1168,7 @@ void CoordAlign::operator()(const DVector &p, double &chisq, DVector &beta, DMat
} // End degenerate parameter check
}

double CoordAlign::fitOnce(bool reportToCerr, bool inPlace) {
double CoordAlign::fitOnce(bool reportToCerr, bool inPlace, bool calcSVD) {
DVector p = getParams();
// First will try doing Newton iterations, keeping a fixed Hessian.
// If it increases chisq or takes too long to converge, we will
Expand Down Expand Up @@ -1240,75 +1240,77 @@ double CoordAlign::fitOnce(bool reportToCerr, bool inPlace) {
// SVD, which is what we want anyway to find the non-pos-def values

if (choleskyFails) {
cerr << "Caught exception during Cholesky" << endl;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is dropping this line intentional?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, because you just get the same message again in the std::runtime_error(), so this std::cerr message isn't necessary.

if (inPlace) {
throw std::runtime_error("Cannot describe degeneracies while dividing in place");
}
int N = alpha.cols();
set<int> degen;
DMatrix U(N, N);
DVector S(N);

#ifdef USE_TMV
tmv::Eigen(symAlpha, U, S);
#elif defined USE_EIGEN
{
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(alpha);
U = eig.eigenvectors();
S = eig.eigenvalues();
}
#endif
// Both packages promise to return eigenvalues in increasing
// order, but let's not depend on that. Report largest/smallest
// abs values of eval's, and print them all
int imax = S.size() - 1; // index of largest, smallest eval
double smax = abs(S[imax]);
int imin = 0;
double smin = abs(S[imin]);
for (int i = 0; i < U.cols(); i++) {
cerr << i << " Eval: " << S[i] << endl;
double s = abs(S[i]);
if (s > smax) {
smax = s;
imax = i;
if (calcSVD) {
if (inPlace) {
throw std::runtime_error("Cannot describe degeneracies while dividing in place");
}
if (s < smin) {
smin = s;
imin = i;
int N = alpha.cols();
set<int> degen;
DMatrix U(N, N);
DVector S(N);

#ifdef USE_TMV
tmv::Eigen(symAlpha, U, S);
#elif defined USE_EIGEN
{
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(alpha);
U = eig.eigenvectors();
S = eig.eigenvalues();
}
if (S[i] < 1e-6) degen.insert(i);
}
cerr << "Largest abs(eval): " << smax << endl;
cerr << "Smallest abs(eval): " << smin << endl;
degen.insert(imin);
// Find biggest contributors to non-positive (or marginal) eigenvectors
const int ntop = MIN(N, 20);
for (int isv : degen) {
cerr << "--->Eigenvector " << isv << " eigenvalue " << S(isv) << endl;
// Find smallest abs coefficient
#endif
// Both packages promise to return eigenvalues in increasing
// order, but let's not depend on that. Report largest/smallest
// abs values of eval's, and print them all
int imax = S.size() - 1; // index of largest, smallest eval
double smax = abs(S[imax]);
int imin = 0;
for (int i = 0; i < U.rows(); i++)
if (abs(U(i, isv)) < abs(U(imin, isv))) imin = i;
vector<int> top(ntop, imin);
for (int i = 0; i < U.rows(); i++) {
for (int j = 0; j < ntop; j++) {
if (abs(U(i, isv)) >= abs(U(top[j], isv))) {
// Push smaller entries to right
for (int k = ntop - 1; k > j; k--) top[k] = top[k - 1];
top[j] = i;
break;
}
double smin = abs(S[imin]);
for (int i = 0; i < U.cols(); i++) {
cerr << i << " Eval: " << S[i] << endl;
double s = abs(S[i]);
if (s > smax) {
smax = s;
imax = i;
}
if (s < smin) {
smin = s;
imin = i;
}
if (S[i] < 1e-6) degen.insert(i);
}
for (int j : top) {
string badAtom = pmc.atomHavingParameter(j);
int startIndex, nParams;
pmc.parameterIndicesOf(badAtom, startIndex, nParams);
cerr << "Coefficient " << U(j, isv) << " at parameter " << j << " Map " << badAtom << " "
<< j - startIndex << " of " << nParams << endl;
cerr << "Largest abs(eval): " << smax << endl;
cerr << "Smallest abs(eval): " << smin << endl;
degen.insert(imin);
// Find biggest contributors to non-positive (or marginal) eigenvectors
const int ntop = MIN(N, 20);
for (int isv : degen) {
cerr << "--->Eigenvector " << isv << " eigenvalue " << S(isv) << endl;
// Find smallest abs coefficient
int imin = 0;
for (int i = 0; i < U.rows(); i++)
if (abs(U(i, isv)) < abs(U(imin, isv))) imin = i;
vector<int> top(ntop, imin);
for (int i = 0; i < U.rows(); i++) {
for (int j = 0; j < ntop; j++) {
if (abs(U(i, isv)) >= abs(U(top[j], isv))) {
// Push smaller entries to right
for (int k = ntop - 1; k > j; k--) top[k] = top[k - 1];
top[j] = i;
break;
}
}
}
for (int j : top) {
string badAtom = pmc.atomHavingParameter(j);
int startIndex, nParams;
pmc.parameterIndicesOf(badAtom, startIndex, nParams);
cerr << "Coefficient " << U(j, isv) << " at parameter " << j << " Map " << badAtom << " "
<< j - startIndex << " of " << nParams << endl;
}
}
}
throw std::runtime_error("Cholesky decomposition failed");
throw std::runtime_error("Cholesky decomposition failed, likely because data is not sufficient"
" to constrain the model");
}

// Now attempt Newton iterations to solution, with fixed alpha
Expand Down
19 changes: 17 additions & 2 deletions src/subs/WCSFitRoutine.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -381,7 +381,7 @@ void WCSFit::setObjects(int i, const map<std::string, std::vector<double>> &tabl

void WCSFit::fit(double maxError, int minFitExposures, double reserveFraction, int randomNumberSeed,
double minimumImprovement, double clipThresh, double chisqTolerance, bool clipEntireMatch,
bool divideInPlace, bool purgeOutput, double minColor, double maxColor) {
bool divideInPlace, bool purgeOutput, double minColor, double maxColor, bool calcSVD) {

PROGRESS(2, Purging defective detections and matches);

Expand Down Expand Up @@ -466,7 +466,7 @@ void WCSFit::fit(double maxError, int minFitExposures, double reserveFraction, i
}

// Do the fit here!!
double chisq = ca.fitOnce(verbose >= 1, divideInPlace); // save space if selected
double chisq = ca.fitOnce(verbose >= 1, divideInPlace, calcSVD); // save space if selected
// Note that fitOnce() remaps *all* the matches, including reserved ones.
double max;
int dof;
Expand Down Expand Up @@ -507,6 +507,21 @@ void WCSFit::fit(double maxError, int minFitExposures, double reserveFraction, i
}
if (verbose >= 0) std::cerr << "# Clipped " << nclip << " matches " << std::endl;

auto badExposures = findUnderpopulatedExposures<Astro>(minFitExposures, matches, exposures,
extensions, mapCollection);

PROGRESS(2, Purging bad exposure parameters and Detections after clipping);
// Freeze parameters of an exposure model and clip all
// Detections that were going to use it.
for (auto i : badExposures) {
cout << "WARNING: Shutting down exposure map " << i.first << " with only "
<< i.second << " fitted detections after clipping " << std::endl;
freezeMap<Astro>(i.first, matches, extensions, mapCollection);
}
if (purgeOutput) {
PROGRESS(2, Purging unfittable maps);
mapCollection.purgeInvalid();
}
} while (coarsePasses || nclip > 0);

// If there are reserved Matches, run sigma-clipping on them now.
Expand Down