connor_solver.h

Go to the documentation of this file.

#ifndef CONNOR_SOLVER_H
#define CONNOR_SOLVER_H
#include "cost_function.h"
#include "convergence_check.h"
#include "observer.h"
#include <blitz/array.h>
#include <boost/shared_ptr.hpp>
#include <boost/utility.hpp>

namespace FullPhysics {

// Note that we have the various equations here in Latex. This is
// fairly unreadable as text comments, but if you view the doxgyen
// documentation created by this it gives a nicely formatted output.


/****************************************************************/
class ConnorSolverState;

class ConnorSolver : public Observable<ConnorSolver>, 
                     public Printable<ConnorSolver>,
                     boost::noncopyable {
public:
  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  ConnorSolver(const boost::shared_ptr<CostFunction>& Cf,
               const boost::shared_ptr<ConvergenceCheck>& Conv,
               double Gamma_initial = 0.0,
               const std::string& Save_test_data = "")
    : save_test_data_(Save_test_data), 
      cost_function_(Cf), convergence_check_(Conv), gamma_initial(Gamma_initial)
  { }
  virtual ~ConnorSolver() {}

  virtual void add_observer(Observer<ConnorSolver>& Obs)
  { add_observer_do(Obs, *this);}
  virtual void remove_observer(Observer<ConnorSolver>& Obs)
  { remove_observer_do(Obs, *this);}

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  void save_test_data(const std::string& Fname) { save_test_data_ = Fname;}

  void to_stream(std::ostream& os) const;
  void from_stream(std::istream& is);

  boost::shared_ptr<ConnorSolverState> state() const;
  void state(const ConnorSolverState& S);

  void test_do_inversion(const std::string& Fname, 
                         blitz::Array<double, 1>& Dx, 
                         blitz::Array<double, 2>& Kt_se_m1_k);
  virtual bool solve(const blitz::Array<double, 1>& Initial_guess, 
                     const blitz::Array<double, 1>& Apriori, 
                     const blitz::Array<double, 2>& Apriori_cov);
  virtual blitz::Array<double, 2> aposteriori_covariance_scaled() const;
  virtual blitz::Array<double, 2> aposteriori_covariance() const;
  blitz::Array<double, 1> x_solution_uncertainty() const;

  virtual blitz::Array<double, 2> averaging_kernel() const;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double gamma_last_step() const {return gamma_last_step_;}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  boost::shared_ptr<CostFunction> cost_function() const
  { return cost_function_; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  boost::shared_ptr<ConvergenceCheck> convergence_check() const 
  { return convergence_check_; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  int number_iteration() const {return fit_statistic().number_iteration; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  int number_divergent() const {return fit_statistic().number_divergent; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  int outcome_flag() const { return (int) fit_statistic().outcome; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  blitz::Array<double, 1> x_apriori() const { return x_a;}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  blitz::Array<double, 1> x_apriori_uncertainty() const 
  { return sigma_ap; }
  virtual blitz::Array<double, 1> x_solution() const { return x_i; }
  int x_solution_size() const { return x_solution().rows();}
  blitz::Array<double, 1> x_solution_zero_unused() const;
  blitz::Array<double, 2> apriori_covariance() const;
  blitz::Array<double, 1> x_update() const { return dx; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 2> jacobian() const {return k;}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual FitStatistic fit_statistic() const {return fstat; }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 1> residual() const {return residual_;}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 1> residual_covariance_diagonal() const 
  {return se;}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 2> apriori_covariance_inv_norm() const 
  {return sa_m1_scaled;}

  void print(std::ostream& Os) const { Os << "ConnorSolver";}

  // These members are "protected" just so we get doxygen
  // documentation on them (private members don't appear in the
  // documentation in general).  
protected:
  blitz::Array<double, 1> zero_unused_parm() const;
  void do_inversion();

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  std::string save_test_data_;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  blitz::Array<double, 1> se;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  blitz::Array<double, 2> k;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  blitz::Array<double, 2> kt_se_m1_k;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  double gamma;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  double gamma_last_step_;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  FitStatistic fstat;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  blitz::Array<double, 2> apriori_cov_scaled;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  blitz::Array<double, 2> sa_m1_scaled;
  
  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  blitz::Array<double, 1> sigma_ap;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  blitz::Array<double, 1> x_a;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  blitz::Array<double, 1> x_i;

  //-----------------------------------------------------------------------
  // This gets updated each iteration in "solve".
  //-----------------------------------------------------------------------

  blitz::Array<double, 1> residual_;

  //-----------------------------------------------------------------------
  // This gets updated each iteration in "do_inversion".
  //-----------------------------------------------------------------------

  blitz::Array<double, 1> dx;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  boost::shared_ptr<CostFunction> cost_function_;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------

  boost::shared_ptr<ConvergenceCheck> convergence_check_;

  //-----------------------------------------------------------------------
  //-----------------------------------------------------------------------
  double gamma_initial;

  static double rcond;
};

/****************************************************************/
class ConnorSolverState: public Printable<ConnorSolverState>,
                         boost::noncopyable {
public:
  ConnorSolverState(const blitz::Array<double, 1>& X_i, 
                    const blitz::Array<double, 1>& X_a, 
                    const blitz::Array<double, 2>& Apriori_cov_scaled, 
                    const blitz::Array<double, 2>& Sa_m1_scaled,
                    const blitz::Array<double, 1>& Sigma_ap, 
                    double Gamma, double Gamma_last_step, double Gamma_intial,
                    const blitz::Array<double, 1>& Residual, 
                    const blitz::Array<double, 1>& Se, 
                    const blitz::Array<double, 2>& K, 
                    const blitz::Array<double, 2>& Kt_se_m1_k, 
                    const blitz::Array<double, 1>& Dx, 
                    const FitStatistic& Fstat)
    : x_i_(X_i.copy()), x_a_(X_a.copy()),
      apriori_cov_scaled_(Apriori_cov_scaled.copy()), 
      sa_m1_scaled_(Sa_m1_scaled.copy()),
      sigma_ap_(Sigma_ap.copy()), gamma_(Gamma), 
      gamma_last_step_(Gamma_last_step), gamma_initial_(Gamma_intial),
      residual_(Residual.copy()), se_(Se.copy()), k_(K.copy()), 
      kt_se_m1_k_(Kt_se_m1_k.copy()), dx_(Dx.copy()), fstat_(Fstat)
  { }
  virtual ~ConnorSolverState() {}
  void print(std::ostream& Os) const
  { Os << "ConnorSolverState"; }
  const blitz::Array<double, 1>& x_i() const { return x_i_;}
  const blitz::Array<double, 1>& x_a() const { return x_a_;}
  const blitz::Array<double, 2>& apriori_cov_scaled() const 
  { return apriori_cov_scaled_;}
  const blitz::Array<double, 2>& sa_m1_scaled() const { return sa_m1_scaled_;}
  const blitz::Array<double, 1>& sigma_ap() const { return sigma_ap_;}
  double gamma() const { return gamma_;}
  double gamma_last_step() const { return gamma_last_step_;}
  double gamma_initial() const { return gamma_initial_;}
  const blitz::Array<double, 1>& residual() const { return residual_;}
  const blitz::Array<double, 1>& se() const { return se_;}
  const blitz::Array<double, 2>& k() const { return k_;}
  const blitz::Array<double, 2>& kt_se_m1_k() const { return kt_se_m1_k_;}
  const blitz::Array<double, 1>& dx() const { return dx_;}
  const FitStatistic& fstat() const { return fstat_;}
private:
  blitz::Array<double, 1> x_i_;
  blitz::Array<double, 1> x_a_;
  blitz::Array<double, 2> apriori_cov_scaled_;
  blitz::Array<double, 2>sa_m1_scaled_;
  blitz::Array<double, 1> sigma_ap_;
  double gamma_;
  double gamma_last_step_;
  double gamma_initial_;
  blitz::Array<double, 1> residual_;
  blitz::Array<double, 1> se_;
  blitz::Array<double, 2> k_;
  blitz::Array<double, 2> kt_se_m1_k_;
  blitz::Array<double, 1> dx_;
  FitStatistic fstat_;
};

inline std::istream& operator>>(std::istream& Is, ConnorSolver& Solve)
{
  Solve.from_stream(Is);
  return Is;
}

}
#endif