Find a stable orbit with cgd

_code/notes/diffsim/main.cpp

@@ -26,6 +26,32 @@ struct LaunchObjective {
         initial_conditions_.segment<dim_>(3 * dim_) = Vector2<Scalar>(0, 0);  // velocity of earth
     }
 
+    NBody& nbody() { return nbody_; }
+
+    void prep_initial_conditions(VectorX<Scalar> const& velocity) {
+        initial_conditions_.segment<dim_>(3 * dim_) = velocity;
+        nbody_.initialize(dim_, masses_, initial_conditions_, delta_t_);
+    }
+
+    void operator()(VectorX<Scalar> const& velocity, Scalar& loss) {
+        // prepare initial conditions
+        initial_conditions_.segment<dim_>(3 * dim_) = velocity;
+
+        // run the simulation
+        nbody_.initialize(dim_, masses_, initial_conditions_, delta_t_);
+        for (uint32_t i = 0; i < n_steps_; ++i) {
+            nbody_.step();
+        }
+
+        // compute loss
+        Vector2<Scalar> const goal = Vector2<Scalar>(0, 1);
+        Vector2<Scalar> const difference = nbody_.state().segment<dim_>(dim_) - goal;
+        loss = difference.squaredNorm();
+
+        // print stuff
+        std::cout << "loss = " << loss << "\n";
+    }
+
     void operator()(VectorX<Scalar> const& velocity, Scalar& loss, VectorX<Scalar>& gradient) {
         // prepare initial conditions
         initial_conditions_.segment<dim_>(3 * dim_) = velocity;
@@ -50,8 +76,8 @@ struct LaunchObjective {
         std::cout << "loss = " << loss << "\n";
         std::cout << "loss gradient\n";
         std::cout << gradient << '\n';
-        std::cout << "state jacobian\n";
-        std::cout << nbody_.state_jacobian() << '\n';
+        //std::cout << "state jacobian\n";
+        //std::cout << nbody_.state_jacobian() << '\n';
     }
 
     constexpr static uint32_t dim_ = 2;
@@ -68,67 +94,66 @@ struct LaunchObjective {
 
 //--------------------------------------------------------------------------------------------------
 int main() {
-    optimization::Rosenbrock<-1> rosenbrock;
-    optimization::ConjugateGradientDescent<-1> cgd(
-        1e-3,   // gradient threshold
-        -1,     // abs threshold
-        100000, // max evaluations
-        10000   // max iterations
-    );
-    optimization::NelderMead<-1> nelder_mead(
-        10000, // max evaluations
-        1e-6   // relative y tolerance
-    );
-
-    VectorX<Scalar> initial_point;
-    uint32_t const dim = 10;
-    initial_point.resize(dim);
-    initial_point[0] = -1;
-    initial_point[1] = 2;
-    for (uint32_t i = 2; i < dim; ++i) {
-        initial_point[i] = -1;
-    }
-
-    cgd.optimize(rosenbrock, initial_point);
-    std::cout << "n evaluations: " << cgd.n_evaluations() << '\n';
-    std::cout << "n iterations: " << cgd.n_iterations() << '\n';
-    std::cout << "final point: " << cgd.point() << '\n';
-    std::cout << '\n';
-
-    nelder_mead.optimize(rosenbrock, initial_point);
-    std::cout << "n evaluations: " << nelder_mead.n_evaluations() << '\n';
-    std::cout << "final point: " << nelder_mead.point() << '\n';
-    std::cout << '\n';
-
-    return 0;
-
-
     LaunchObjective objective;
 
     VectorX<Scalar> satellite_velocity;
     satellite_velocity.resize(2);
     satellite_velocity[0] = 0;
-    satellite_velocity[1] = 1;
-
-    Scalar loss;
-    VectorX<Scalar> gradient;
-
-    objective(satellite_velocity, loss, gradient);
-
-    /*
-    grapher.initialize(
-        &nbody,
-        BoundingBox<Vector2<Scalar>>(
-            {-2, -2},
-            {2, 2}
-        )
-    );
-    grapher.init_glut();
-    glutDisplayFunc(display);
-    glutMouseFunc(mouse);
-    glutIdleFunc(idle);
-    grapher.run();
-    */
+    satellite_velocity[1] = 0.4;
+
+    //Scalar loss;
+    //VectorX<Scalar> gradient;
+    //objective(satellite_velocity, loss, gradient);
+    //std::cout << '\n';
+
+    bool do_nm = false;
+    bool do_cgd = true;
+    bool do_vis = true;
+
+    if (do_nm) {
+        optimization::NelderMead<-1> nelder_mead(
+            10000, // max evaluations
+            1e-6   // relative y tolerance
+        );
+        nelder_mead.optimize(objective, satellite_velocity, 0.5);
+        satellite_velocity = nelder_mead.point();
+        std::cout << "n evaluations: " << nelder_mead.n_evaluations() << '\n';
+        std::cout << "final point: " << nelder_mead.point() << '\n';
+        std::cout << "final loss: " << nelder_mead.value() << '\n';
+        std::cout << '\n';
+    }
+
+    if (do_cgd) {
+        optimization::ConjugateGradientDescent<-1> cgd(
+            1e-3,   // gradient threshold
+            -1,     // abs threshold
+            100000, // max evaluations
+            10000   // max iterations
+        );
+        cgd.optimize(objective, satellite_velocity);
+        satellite_velocity = cgd.point();
+        std::cout << "n evaluations: " << cgd.n_evaluations() << '\n';
+        std::cout << "n iterations: " << cgd.n_iterations() << '\n';
+        std::cout << "final point: " << cgd.point() << '\n';
+        std::cout << '\n';
+    }
+
+    if (do_vis) {
+        objective.prep_initial_conditions(satellite_velocity);
+
+        grapher.initialize(
+            &objective.nbody(),
+            BoundingBox<Vector2<Scalar>>(
+                {-2, -2},
+                {2, 2}
+            )
+        );
+        grapher.init_glut();
+        glutDisplayFunc(display);
+        glutMouseFunc(mouse);
+        glutIdleFunc(idle);
+        grapher.run();
+    }
 
     return 0;
 }

_code/notes/diffsim/n_body.cpp

@@ -110,8 +110,8 @@ void NBody::initialize(
     points_.resize(law_.n() * law_.d);
     points_ = state_.segment(0, law_.n() * law_.d);
 
-    std::cout << "state.size() = " << state_.size() << '\n';
-    std::cout << "points.size() = " << points_.size() << '\n';
+    //std::cout << "state.size() = " << state_.size() << '\n';
+    //std::cout << "points.size() = " << points_.size() << '\n';
 }
 
 //..................................................................................................