cpp nomad chk1 go

cmakelists.txt

@@ -8,6 +8,7 @@
 cmake_minimum_required(VERSION 3.12)
 
 project(starid)
+set(CMAKE_CXX_STANDARD 20)
 
 #add_subdirectory(libstarid)
 add_subdirectory(libstarid-debug)

libstarid/triangles/nomad_triangle.cpp

@@ -1,62 +1,56 @@
+#include <set>
 #include "nomad_triangle.h"
 
-starid::StartriangleNOMAD::StartriangleNOMAD(Eigen::MatrixXd &starvecs, starid::Starpairs &starpairs) :
-    side1{}, stara{0}, starb{0}, starc{0}, starvecs(starvecs), starpairs(starpairs) {
+StartriangleNOMAD::StartriangleNOMAD(MatrixXd &starvecs, Starpairs &starpairs) :
+    stara{0}, starb{0}, starc{0}, starvecs(starvecs), starpairs(starpairs) {
     sv1 << 0., 0., 0.;
     sv2 << 0., 0., 0.;
     sv3 << 0., 0., 0.;
 }
 
-void starid::StartriangleNOMAD::first() {
+void StartriangleNOMAD::first() {
     stara = 0; starb = 2; starc = 3;
     sv1 << 0., 0., 1.;
     sv2 = starvecs.row(starb);
     sv3 = starvecs.row(starc);
-    side1 = starid::Startriangleside2(sv1, sv2, starpairs);
-    stara = 0;
-//    side1 = starid::Startriangleside(starid::pairangle(starvecs, stara, starb), starpairs);
-//    side2 = starid::Startriangleside(starid::pairangle(starvecs, starb, starc), starpairs);
-//    side3 = starid::Startriangleside(starid::pairangle(starvecs, starc, stara), starpairs);
-//    this->constrain();
+    side1 = Startriangleside2(sv1, sv2, starpairs);
+    side2 = Startriangleside2(sv2, sv3, starpairs);
+    side3 = Startriangleside2(sv3, sv1, starpairs);
+    this->chk1();
 }
 
-void starid::StartriangleNOMAD::from_parent() {
+void StartriangleNOMAD::from_parent() {
 //    while (starc == stara || starc == starb) starc = rand() % starvecs.rows();
-//    side2 = starid::Startriangleside(starid::pairangle(starvecs, starb, starc), starpairs);
-//    side3 = starid::Startriangleside(starid::pairangle(starvecs, starc, stara), starpairs);
+//    side2 = Startriangleside(pairangle(starvecs, starb, starc), starpairs);
+//    side3 = Startriangleside(pairangle(starvecs, starc, stara), starpairs);
 //    this->constrain();
 }
 
-void starid::StartriangleNOMAD::constrain() {
-//    for (auto STAR1side1 = side1.stars.begin(), end = side1.stars.end(); STAR1side1 != end; ++STAR1side1) {
-//        auto STAR1side3 = side3.stars.find(STAR1side1->first);
-//        if (STAR1side3 == side3.stars.end()) continue;
-//        auto &pairs3 = STAR1side3->second;
-//        auto &pairs1 = STAR1side1->second;
-//        for (auto pairs1it = pairs1.begin(), end = pairs1.end(); pairs1it != end; ++pairs1it) {
-//            auto STAR2side2 = side2.stars.find(pairs1it->first);
-//            if (STAR2side2 == side2.stars.end()) continue;
-//            auto &pairs2 = STAR2side2->second;
-//            for (auto pairs2it = pairs2.begin(), end = pairs2.end(); pairs2it != end; ++pairs2it) {
-//                auto STAR3side3 = pairs3.find(pairs2it->first);
-//                if (STAR3side3 == pairs3.end()) continue;
-//                pairs1it->second = 1;
-//                pairs2it->second = 1;
-//                STAR3side3->second = 1;
-//            }
-//        }
-//    } // keep these pairs
-//    side1.drops();
-//    side2.drops();
-//    side3.drops();
+void StartriangleNOMAD::chk1() {
+    std::set<int> ok1, ok2, ok3;
+    for (auto & a1 : side1.stars) {
+        if (!side3.stars.contains(a1.first)) continue;
+        for (auto & b1 : a1.second) {
+            if (!side2.stars.contains(b1)) continue;
+            for (auto & c2 : side2.stars[b1]) {
+                if (!side3.stars[a1.first].contains(c2)) continue;
+                ok1.insert(a1.first); ok1.insert(b1);
+                ok2.insert(b1); ok2.insert(c2);
+                ok3.insert(c2); ok3.insert(a1.first);
+            }
+        }
+    }
+    side1.update_side(ok1);
+    side2.update_side(ok2);
+    side3.update_side(ok3);
 }
 
-void starid::StartriangleNOMAD::feedback(starid::Startriangleside &nextside1) {
+void StartriangleNOMAD::feedback(Startriangleside &nextside1) {
 //    side2.stars = nextside1.stars;
 //    this->constrain();
 }
 
-bool starid::StartriangleNOMAD::stop() {
+bool StartriangleNOMAD::stop() {
 //    if (s12.stars.size() <= 1) return true;
     return false;
 }

libstarid/triangles/nomad_triangle.h

@@ -4,28 +4,30 @@
 #include "Eigen/Core"
 #include "../sky/geometry.h"
 #include "star_triangle_side.h"
+using namespace Eigen;
+using namespace starid;
 
 namespace starid {
 
     class StartriangleNOMAD {
     public:
-        StartriangleNOMAD(Eigen::MatrixXd &starvecs, Starpairs &starpairs);
+        StartriangleNOMAD(MatrixXd &starvecs, Starpairs &starpairs);
 
         void first();
 
         void from_parent();
 
-        void feedback(starid::Startriangleside &next);
+        void feedback(Startriangleside &next);
 
-        void constrain();
+        void chk1();
 
         bool stop();
     private:
-        Eigen::MatrixXd starvecs;
+        MatrixXd starvecs;
         Starpairs starpairs;
-        Startriangleside2 side1;
+        Startriangleside2 side1, side2, side3;
         int stara, starb, starc;
-        Eigen::Vector3d sv1, sv2, sv3;
+        Vector3d sv1, sv2, sv3;
 
     };
 

libstarid/triangles/star_pairs.cpp

@@ -1,8 +1,8 @@
 #include "star_pairs.h"
 
-starid::Starpairs::Starpairs() {}
+Starpairs::Starpairs() {}
 
-void starid::Starpairs::generate(starid::Sky &sky) {
+void Starpairs::generate(Sky &sky) {
     int pairndx = 0;
     for (auto star: sky.stars) {
         std::vector<int> starndxs = sky.stars_near_point(star.x, star.y, star.z);
@@ -16,7 +16,7 @@ void starid::Starpairs::generate(starid::Sky &sky) {
                 double angle = acos((sky.stars[starndx1].x * sky.stars[starndx2].x) +
                                     (sky.stars[starndx1].y * sky.stars[starndx2].y) +
                                     (sky.stars[starndx1].z * sky.stars[starndx2].z));
-                if (std::abs(angle) > starid::star_pair_angle_limit) continue; // max pair angle
+                if (std::abs(angle) > star_pair_angle_limit) continue; // max pair angle
                 std::tuple<double, int, int> starpair{angle, starndx1, starndx2};
                 starpairs.push_back(starpair);
                 starpairsndxs.insert({key, pairndx}); // update map of unique pairs
@@ -28,17 +28,16 @@ void starid::Starpairs::generate(starid::Sky &sky) {
     angndxs.sort();
 }
 
-std::unordered_map<int, std::unordered_map<int, int>>
-starid::Starpairs::pairs_for_angle(double angle, double tol_radius) {
-    std::unordered_map<int, std::unordered_map<int, int>> stars;
+starsdict Starpairs::pairs_for_angle(double angle, double tol_radius) {
+    starsdict stars;
     double ang1 = angle - tol_radius;
     double ang2 = angle + tol_radius;
     double epsilon = 1.0;
     if (ang1 <= 0) ang1 = 0;
     if (ang2 <= epsilon * tol_radius) ang2 = epsilon * tol_radius;
-    if (ang1 >= starid::star_pair_angle_limit - epsilon * tol_radius)
-        ang1 = starid::star_pair_angle_limit - epsilon * tol_radius;
-    if (ang2 >= starid::star_pair_angle_limit) ang2 = starid::star_pair_angle_limit;
+    if (ang1 >= star_pair_angle_limit - epsilon * tol_radius)
+        ang1 = star_pair_angle_limit - epsilon * tol_radius;
+    if (ang2 >= star_pair_angle_limit) ang2 = star_pair_angle_limit;
     std::vector<int> pairndxs = angndxs.findndxs(ang1, ang2);
     for (auto ndx: pairndxs) {
         int star1 = std::get<1>(starpairs[ndx]);
@@ -65,7 +64,7 @@ starid::Starpairs::pairs_for_angle(double angle, double tol_radius) {
     return stars;
 }
 
-std::string starid::Starpairs::pair_labeler(int catndx1, int catndx2) {
+std::string Starpairs::pair_labeler(int catndx1, int catndx2) {
     if (catndx1 > catndx2) {
         int tmp = catndx1;
         catndx1 = catndx2;

libstarid/triangles/star_pairs.h

@@ -10,20 +10,22 @@
 #include <unordered_map>
 #include <map>
 #include "../sky/sky.h"
-
+using namespace Eigen;
+using namespace starid;
+using starsdict = std::unordered_map<int, std::unordered_map<int, int>>;
 namespace starid {
 
     class Starpairs {
 
     public:
         Starpairs();
 
-        void generate(starid::Sky &sky);
+        void generate(Sky &sky);
 
-        std::unordered_map<int, std::unordered_map<int, int>> pairs_for_angle(double angle, double tol_radius);
+        starsdict pairs_for_angle(double angle, double tol_radius);
 
     private:
-        starid::FloatsIndexer angndxs;
+        FloatsIndexer angndxs;
         std::vector<std::tuple<double, int, int>> starpairs;
         std::unordered_map<std::string, int> starpairsndxs;
 

libstarid/triangles/star_triangle_side.cpp

@@ -1,21 +1,22 @@
+#include <algorithm>
 #include "star_triangle_side.h"
 
-starid::Startriangleside::Startriangleside(double ang, starid::Starpairs &pairs) {
+Startriangleside::Startriangleside(double ang, Starpairs& pairs) {
     angtol = 2.0 * std::sqrt(500.0 * 500.0 + 500.00 * 500.0) * arcseconds_to_radians;
     stars = pairs.pairs_for_angle(ang, angtol);
     log_star_count.push_back(stars.size());
     log_pair_count.push_back(countpairs());
 }
 
-starid::Startriangleside::Startriangleside() {}
+Startriangleside::Startriangleside() {}
 
-void starid::Startriangleside::update(Startriangleside &side) {
+void Startriangleside::update(Startriangleside& side) {
     stars = side.stars;
     for (auto tmp: side.log_pair_count) log_pair_count.push_back(tmp);
     for (auto tmp: side.log_star_count) log_star_count.push_back(tmp);
 }
 
-void ::starid::Startriangleside::drops(bool dolog) {
+void Startriangleside::drops(bool dolog) {
     for (auto star1 = stars.begin(), end = stars.end(); star1 != end; ++star1) {
         auto &pairs = star1->second;
         for (auto star2 = pairs.begin(); star2 != pairs.end();) {
@@ -37,40 +38,48 @@ void ::starid::Startriangleside::drops(bool dolog) {
     }
 }
 
-int starid::Startriangleside::countpairs() {
+int Startriangleside::countpairs() {
     int result = 0;
     for (auto it1 = stars.begin(), end = stars.end(); it1 != end; ++it1) {
         result += it1->second.size();
     }
     return result;
 }
 
-starid::Startriangleside2::Startriangleside2(Eigen::Vector3d sv1_, Eigen::Vector3d sv2_, starid::Starpairs &pairs) {
-    angtol = .003; sv1 = sv1_; sv2 = sv2_; ang = std::acos(sv1.transpose() * sv2);
-    stars = pairs.pairs_for_angle(ang, angtol);
+Startriangleside2::Startriangleside2(Vector3d sv1, Vector3d sv2, Starpairs& pairs, float angtol) :
+    angtol{angtol}, sv1{sv1}, sv2{sv2}, ang{std::acos(sv1.transpose() * sv2)}, stars{} {
+    auto tmp = pairs.pairs_for_angle(ang, angtol);
+    for (auto & star1 : tmp) {
+        for (auto & star2 : star1.second) {
+            if (!stars.contains(star1.first)) stars.insert({star1.first, {star2.first}});
+            else stars[star1.first].insert(star2.first);
+        }
+    }
     log_star_count.push_back(stars.size());
     log_pair_count.push_back(countpairs());
 }
 
-starid::Startriangleside2::Startriangleside2() :
-    stars{{0, {{0, 0}}}}, ang{0}, angtol{0}, log_star_count{} {
+Startriangleside2::Startriangleside2() :
+    stars{}, ang{0}, angtol{0}, log_star_count{}, log_pair_count{} {}
+
+void Startriangleside2::update_side(std::set<int>& ok) {
+    std::set<int> drops;
+    for (auto & star : stars) {
+        if (!ok.contains(star.first)) drops.insert(star.first);
+        else {
+            std::set<int> tmp;
+            std::set_intersection(ok.begin(), ok.end(), star.second.begin(), star.second.end(), std::inserter(tmp, tmp.begin()));
+            star.second = tmp;
+            if (star.second.size() == 0) drops.insert(star.first);
+        }
+    }
+    for (auto & star : drops) stars.erase(star);
+    log_star_count.push_back(stars.size());
+    log_pair_count.push_back(countpairs());
 }
 
-int starid::Startriangleside2::countpairs() {
+int Startriangleside2::countpairs() {
     int result = 0;
-    for (auto it1 = stars.begin(), end = stars.end(); it1 != end; ++it1) {
-        result += it1->second.size();
-    }
+    for (auto & star : stars) result += star.second.size();
     return result;
 }
-
-//starid::Startriangleside2 &starid::Startriangleside2::operator=(starid::Startriangleside2 other) {
-//    log_pair_count = other.log_pair_count;
-//    log_star_count = other.log_star_count;
-//    stars = other.stars;
-//    ang = other.ang;
-//    angtol = other.angtol;
-//    sv1 = other.sv1;
-//    sv2 = other.sv2;
-//    return *this;
-//}

libstarid/triangles/star_triangle_side.h

@@ -4,14 +4,20 @@
 #include <unordered_map>
 #include <map>
 #include <vector>
+#include <set>
 #include "star_pairs.h"
 #include "../sky/geometry.h"
+using namespace Eigen;
+using namespace starid;
+using starsdict = std::unordered_map<int, std::unordered_map<int, int>>;
+using starsdict2 = std::unordered_map<int, std::set<int>>;
+using intvec = std::vector<int>;
 
 namespace starid {
 
     class Startriangleside {
     public:
-        Startriangleside(double ang, starid::Starpairs &pairs);
+        Startriangleside(double, Starpairs&);
 
         Startriangleside();
 
@@ -21,31 +27,29 @@ namespace starid {
 
         int countpairs();
 
-        std::unordered_map<int, std::unordered_map<int, int>> stars;
-        std::vector<int> log_star_count;
-        std::vector<int> log_pair_count;
+        starsdict stars;
+        intvec log_star_count;
+        intvec log_pair_count;
         double angtol;
     };
 
     class Startriangleside2 {
     public:
-        Startriangleside2(Eigen::Vector3d sv1_, Eigen::Vector3d sv2_, starid::Starpairs &pairs);
-
-//        Startriangleside2& operator=(Startriangleside2 other);
+        Startriangleside2(Vector3d, Vector3d, Starpairs&, float angtol = .003);
 
         Startriangleside2();
 
-//        void drops(bool dolog = true);
+        void update_side(std::set<int>&);
 
 //        void update(Startriangleside &side);
 
         int countpairs();
 
-        Eigen::Vector3d sv1, sv2;
+        Vector3d sv1, sv2;
         double ang;
-        std::unordered_map<int, std::unordered_map<int, int>> stars;
-        std::vector<int> log_star_count;
-        std::vector<int> log_pair_count;
+        starsdict2 stars;
+        intvec log_star_count;
+        intvec log_pair_count;
     private:
         double angtol;
     };

starid/triangles/nomad_triangle.py

@@ -29,6 +29,7 @@ def first(self):
         self.side2 = StarTriangleSide(self.sv2, self.sv3, self.starpairs, angtol=self.tol)
         self.side3 = StarTriangleSide(self.sv3, self.sv1, self.starpairs, angtol=self.tol)
         self.chk1()
+        return
 
     def from_parent(self, side2, starb, starc):
         """every triangle after the first is built on and extends its parent triangle. for a parent abca triangle,