#include<iostream>
using namespace std;

#include "mat.h"
using namespace ariadne;

class state_diff: public Optimisable_function
{
   public:
   eos *an_eos;
   double preq, treq; // requested state
   double pwt, twt; // weights
   double rho, e; // primitive state
   int verbose;
   static const int nparams = 2;

   Array<double> cmin, cmax; // parameter limits
   Array_fixed_dimension<double> dc;

   void Read_param_limits(istream& stream)
   {
      cmin.Recreate(nparams); cmax.Recreate(nparams);
      for (int i=0;i<nparams;i++) stream >> cmin[i] >> cmax[i];
   }
   void Read(istream& stream)
   {
      stream >> an_eos;
      stream >> preq >> treq;
      stream >> pwt >> twt;
      stream >> rho >> e;
      dc.Recreate(nparams); stream >> dc;
      stream >> verbose;
   }
   void Write(ostream& stream) const
   {
      stream << "EOS:\n" << an_eos << endl;
      stream << "preq = " << preq << ", Treq = " << treq << endl;
      stream << "pwt = " << pwt << ", Twt = " << twt << endl;
      stream << "rho_est = " << rho << ", e_est = " << e << endl;
      stream << "parameter scales = " << dc << endl;
      stream << "verbose flag = " << verbose << endl;
   }

   Vector params() const
   {
      Vector param_val(nparams);
      param_val[0] = rho;
      param_val[1] = e;
      return param_val;
   }

   void delta(const Vector& c, Vector& delta_c) const { delta_c = dc; }

   double evaluate(const Vector& c) const
   {
      eos::state st;
      st.rho = c[0]; st.e = c[1];
      if (verbose) cout << "trying state " << st << ": ";
      double p = an_eos->p(st);
      double t = an_eos->t(st);
      double diff = pwt * pow(p - preq,2) + twt * pow(t - treq,2);
      if (verbose) cout << "p = " << p << ", T = " << t << ", diff = " << diff
         << endl;
      return diff;
   }
};

int main()
{
   cout << "state function:\n"; state_diff s; cin >> s; cout << s;
   cout << "optimiser:\n"; Optimiser *opt; cin >> opt; cout << opt;
/*
   if (!strcmp(opt->the_class_name(),"monte-carlo"))
      s.Read_param_limits(cin);
*/

   Vector c = s.params();
   if (!opt->minimise(s,c)) cout << "Optimiser failed!\n";
   cout.precision(10);
   cout << "final parameters: "; c.Write_elements(cout); cout << endl;
}