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.

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

DM-31226: Add functionality needed for WCSFitTask in drp_tasks #9

Merged
merged 7 commits into from Dec 16, 2022

Conversation

cmsaunders
Copy link
Collaborator

@cmsaunders cmsaunders commented Oct 6, 2022

This PR adds the following functionality:

  • Conversion from SIP to TPV for the input WCS
  • More information about the pixel maps is exposed to the Python wrapper
  • Options to output catalogs with residuals, fit star positions, etc., in the form of dictionaries that can be trivially converted to Pandas Dataframes.
  • Improvements in interface for fitting proper motions and parallax.

Copy link

@natelust natelust left a comment

Choose a reason for hiding this comment

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

There is some work to be done here, and I would be happy to weight in on any responses, or look again if you wish, but I am also happy for you to take these as advisement and fix it up how you decide.

@@ -256,6 +257,56 @@ namespace astrometry {
virtual void write(YAML::Emitter& os) const;
#endif
virtual string getType() const {return type();}
Matrix33 getPixMatrix() {
cerr << "start get pix" << endl;
Copy link

Choose a reason for hiding this comment

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

Do we normally send things out to std err vs logging? Or is this just because our gbdes forking?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

There's a separate ticket to move to proper logging. However, in this case this should just be removed.

@@ -256,6 +257,56 @@ namespace astrometry {
virtual void write(YAML::Emitter& os) const;
#endif
virtual string getType() const {return type();}
Matrix33 getPixMatrix() {
Copy link

Choose a reason for hiding this comment

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

Why are you defining all of this in the header, instead of the signature here and the body in a src file?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I probably did this before I read that section in the dev guide...will fix.

auto gnom = dynamic_cast<const astrometry::Gnomonic*>(pix);
cerr << "did dyn cast" << endl;
Matrix33 tmpMatrix;
tmpMatrix.setToIdentity();
Copy link

Choose a reason for hiding this comment

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

would it make sense to do this with tmpMatrix only if it is actually going to be needed? (same below)

if (!icrs) {
Vector2 lonlat;
cerr << "not gnom or icrs" << endl;
return lonlat;
Copy link

Choose a reason for hiding this comment

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

I'm not familiar with this object, is lonlat a class level variable?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I just removed this function, because it was something I used for debugging, and I don't think I will need it again.

return icrs->getLonLat();
}
SphericalICRS pole = gnom->getOrient()->getPole();
Vector2 lonlat = pole.getLonLat();
Copy link

Choose a reason for hiding this comment

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

if lonlat is a class level variable name, maybe it would be a good idea to use a different variable name here for human readability

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

See above.

vv[astrometry::VX] *= units[i];
vv[astrometry::VY] *= units[i];
vv[astrometry::PAR] *= units[i];
vv[astrometry::VX] *= units[astrometry::VX];

Choose a reason for hiding this comment

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

I cant tell if this is a github rendering, but is this inside the for loop?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yes. I changed this because I saw that i was being reused as an index. I guess it wasn't affecting the i in the outer scope, because I think the vector was still being filled, but it seems like bad practice.

projection = catalogProjections[detptr->catalogNumber];
if (projection) break;
}
if (projection) {

Choose a reason for hiding this comment

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

couldn't you save this branch if you just did this work inside line 2772, and did break at the end?

// Fill ac and bc
for (int m = 0; m <= a_order; m++) {
for (int n = 0; n <= (a_order - m); n++) {
std::string keyString = "A_" + std::to_string(m) + "_" + std::to_string(n);

Choose a reason for hiding this comment

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

is addition the best way to do this? or is there a standard formatting?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This seems like the simplest way in c++17. I spent some time looking and didn't see better options, but please tell me if I'm missing something.

// Create a SubMap that owns a duplicate of these one or two maps (no sharing):
SubMap sm(pmlist, name, false);
// Delete the polymap if we made it:
delete pv;

Choose a reason for hiding this comment

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

I'm not sure I follow all the need for allocating a new object, and freeing like this. Couldn't this be refactored to not need the dynamic memory like this? perhaps by directly creating the SubMap?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I basically copied what is in readTPV. We can discuss how to refactor.

@@ -352,11 +365,15 @@ void WCSFit::setObjects(int i, const map<std::string, std::vector<double>> &tabl
const std::string &idKey, const std::string &pmCovKey, const std::string &magKey,
const int &magKeyElement, const std::string &magErrKey, const int &magErrKeyElement,
const std::string &pmRaKey, const std::string &pmDecKey,
const std::string &parallaxKey) {
const std::string &parallaxKey, const std::vector<std::vector<double>> &pmCov) {

Choose a reason for hiding this comment

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

Does this parameter have a default in a header somewhere?

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, it defaults to zero. Why do you ask, because it changes the API?

@cmsaunders cmsaunders merged commit ce05b32 into lsst-dev Dec 16, 2022
@cmsaunders cmsaunders deleted the tickets/DM-31226 branch December 16, 2022 20:00
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants