writing back model parameters

src/root_model.hpp

@@ -32,6 +32,7 @@ class Root_model {
    virtual void clear_problems() {}
    virtual std::size_t get_number_of_parameters() const = 0;
    virtual Eigen::ArrayXd get_parameters() const = 0;
+   virtual void set_parameters(const Eigen::ArrayXd&) = 0;
    virtual std::string name() const { return "unnamed"; }
    virtual void run_to(double, double eps = -1.0) = 0;
    virtual void print(std::ostream& out) const { out << "Model: " << name(); }

src/root_solver.cpp

@@ -159,6 +159,8 @@ int RGFlow<Root>::find_root()
    printf("\tRoot_finder status = %s\n", gsl_strerror(status));
 #endif
 
+   set_model_parameters(to_eigen_array(solver->x));
+
    gsl_vector_free(x_init);
    gsl_multiroot_fsolver_free(solver);
 
@@ -300,6 +302,19 @@ Eigen::ArrayXd RGFlow<Root>::calculate_function(const Eigen::ArrayXd& x)
    return f;
 }
 
+void RGFlow<Root>::set_model_parameters(const Eigen::ArrayXd& p)
+{
+   std::size_t offset = 0;
+
+   for (std::size_t m = 0; m < models.size(); m++) {
+      const std::size_t dim = models[m]->model->get_number_of_parameters();
+      models[m]->model->set_parameters(p.block(offset, 0, dim, 1));
+      offset += dim;
+   }
+
+   assert(offset == get_dimension());
+}
+
 Eigen::ArrayXd RGFlow<Root>::to_eigen_array(const gsl_vector* v)
 {
    const std::size_t dim = v->size;

src/root_solver.hpp

@@ -109,6 +109,7 @@ class RGFlow<Root> {
    std::size_t get_dimension() const;  ///< calculate dimension
    unsigned int get_max_iterations() const; ///< returns max. number of iterations
    void initial_guess();               ///< initial guess
+   void set_model_parameters(const Eigen::ArrayXd&); ///< set model parameters
    void starting_point(gsl_vector*);   ///< find starting point
 };
 

test/test_root_solver.cpp

@@ -20,7 +20,7 @@ class Static_model: public Root_model {
    virtual void calculate_spectrum() {}
    virtual std::string name() const { return "Static_model"; }
    virtual void run_to(double, double) {}
-   virtual void set_parameters(const Eigen::Array2d& v) { parameters = v; }
+   virtual void set_parameters(const Eigen::ArrayXd& v) { parameters = v; }
    virtual Eigen::ArrayXd get_parameters() const { return parameters; }
    virtual std::size_t get_number_of_parameters() const { return 2; }
    virtual void set_precision(double) {}