Merge pull request #10 from lsst/tickets/DM-37737

DM-37737: Fix non-deterministic behavior

astrometry/include/PixelMapCollection.h

@@ -204,6 +204,8 @@ namespace astrometry {
     // Set/get the master parameter vector for all PixelMaps
     void setParams(const DVector& p);
     DVector getParams() const;
+    // Get the parameter vector for one map, or return an empty DVector if there are no parameters
+    DVector getParams(string mapName) const;
     std::map<std::string, astrometry::DVector> getParamDict() const;
     int nParams() const {return parameterCount;}
 

astrometry/src/PixelMapCollection.cpp

@@ -135,6 +135,25 @@ PixelMapCollection::getParams() const {
   return p;
 }
 
+DVector
+PixelMapCollection::getParams(string mapName) const {
+  auto melpair = mapElements.find(mapName);
+  const MapElement& mel = melpair->second;
+  // Return an empty vector if the map is fixed
+  if (mel.isFixed) {
+    return DVector();
+  }
+  if (mel.atom) {
+    int nSub = mel.nParams;
+    DVector mparams = mel.atom->getParams().subVector(0, nSub);
+    return mparams;
+  }
+  // Map is not atomic and doesn't have its own parameters; return an empty vector
+  else {
+    return DVector();
+  }
+}
+
 std::map<std::string, astrometry::DVector>
 PixelMapCollection::getParamDict() const {
   map<string, linalg::Vector<double>> outMap;

include/FoF.h

@@ -10,6 +10,7 @@
 #include <set>
 #include <vector>
 #include <algorithm>
+#include <numeric>
 
 #include <iostream>
 
@@ -82,6 +83,24 @@ class Match : public list<const P *> {
         this->push_back(&point);
         cells.insert(cellsContainingPoint.begin(), cellsContainingPoint.end());
     }
+    const std::array<double, DIM> calculateAverage() {
+        // Calculate the average of the x values for the points in Match
+
+        // Return zero if there are no points in the match
+        if (this->size() == 0) {
+            std::array<double, DIM> average = std::array<double, DIM>{0};
+            return std::array<double, DIM>{0};
+        } 
+        // Calculate the average along each dimension
+        std::array<double, DIM> average;
+        for (int d=0; d < DIM; d++) {
+            double d_average = std::accumulate(this->begin(), this->end(), 0.0,
+                                               [d](double a, P const *p) { return a + p->getX()[d]; });
+            d_average /= this->size();
+            average[d] = d_average;
+        }
+        return average;
+    }
 
 private:
     // Upper and lower coords of bounding rectangle of points

src/subs/WCSFoFRoutine.cpp

@@ -346,8 +346,19 @@ void FoFClass::sortMatches(int fieldNumber, int minMatches, bool allowSelfMatche
         if (pcat.empty()) continue;
 
         long matches = 0;
+
+        // Copy Matches into a vector to facilitate sorting
+        std::vector<fof::Match<Point, 2> *> matchVector(pcat.size());
+        matchVector.assign(pcat.begin(), pcat.end());
+
+        // Sort the matches by their ra position to ensure deterministic behavior downstream.
+        const auto matchSort = [](fof::Match<Point, 2> *a, fof::Match<Point, 2> *b) {
+            return a->calculateAverage()[0] < b->calculateAverage()[0];
+        };
+        std::sort(matchVector.begin(), matchVector.end(), matchSort);
+
         // Now loop through matches in this catalog
-        for (PointCat::const_iterator j = pcat.begin(); j != pcat.end(); ++j) {
+        for (auto j = matchVector.cbegin(); j != matchVector.cend(); ++j) {
             // Skip any Match that is below minimum match size
 
             if ((*j)->size() < minMatches) continue;

src/subs/WcsSubs.cpp

@@ -9,6 +9,8 @@
 #include <random>
 using namespace astrometry;
 
+std::mt19937 g(12345);
+
 void fitDefaulted(PixelMapCollection &pmc, set<Extension *> extensions,
                   const vector<unique_ptr<Instrument>> &instruments,
                   const vector<unique_ptr<Exposure>> &exposures, bool logging) {
@@ -85,8 +87,6 @@ void fitDefaulted(PixelMapCollection &pmc, set<Extension *> extensions,
         vector<int> vx(nGridPoints);
         for (int i = 0; i < vx.size(); i++) vx[i] = i;
         vector<int> vy = vx;
-        std::random_device rd;
-        std::mt19937 g(rd());
         std::shuffle(vy.begin(), vy.end(), g);
         double xstep = (b.getXMax() - b.getXMin()) / nGridPoints;
         double ystep = (b.getYMax() - b.getYMin()) / nGridPoints;