DiffitModel.hpp

#ifndef __MODELPUB_H__
#define __MODELPUB_H__

//#include "nr.hpp"
//#include "def.hpp"

#include <CapeOpen.hpp>
#include <Basic/DianaBasic.hpp>
#include <Diana/DianaDAESO.hpp>
#include <Basic/DianaArray.hpp>
#include <Diana.hpp>
#include <Basic/CapeModel.hpp>
//#include "DianaSensDAESO.hpp"
#include <string>
#include <math.h>
#include <vector>



using namespace Numeric::Solvers::Model;

class Glob;

using namespace std;

class ExperimentSystem
{
public:

    ExperimentSystem(const long _nvar, const long _npar, const long _nobs);
    ExperimentSystem(const ExperimentSystem& ex);
    virtual ~ExperimentSystem();

    void updateResidualsFrom(const  ExperimentSystem& _sys);
    int * y0fix;
    int * doP;

  long nobs;
  long nvar;    
  long npar;    

  // the actual data
  long nMeasure;
  double *xMeasure;
  double **sigma;
  double fitstart, fitend; 
  long firstMeasure;
  long lastMeasure;
  // multiple shooting data
  long nPoints;
  double *mesh;   
  long me;
  long mg;
  
  // data obtained by integration
  double objF;


  double **h;

  double **residues;

  double *r2;

  double *r3;

  double ***dyds;

  double ***dydp;
  double ***dmds;
  double ***dmdp;
  double ***dR2ds;
  double **dR2dp;
  double ***dR3ds;
  double **dR3dp;
  double *ua,**Ea,**Pa;
  double *ue,**Ee,**Pe;
  double *ug,**Eg,**Pg;
    double *par;
  double **dS;
  double *dP;

  double *y0;

  double **yComp;

  double **yPred;

    double **yTry;//,**yTrySave; 

  double **yMeasure;


  char *fileName;
  string *splineFile;
  double **splineNodes;
  double **splineY;
  double **splineGam;
  unsigned long *nNodes;
  int splinesDefined;
    long NSPLINES;


};



typedef ExperimentSystem GlobExp;
typedef std::vector<ExperimentSystem> ExperimentVector;


#ifndef TRUE
#define TRUE 1
#define FALSE 0
#endif

class Glob{
public:
  //flags
    long nvar;
  int noGnu;
  int noMeasurements;
  int wait;
  int usesig;
  int odessa;
  int stiff;
  int nowait;
  int reg;
  int simInit;
  int nodamp;
  //all the rest
  long npar;
  long maxit;
  long nIter;
  int nrExp;
  int *doP;
  double **covar;
  double chisq;
  double eps;
  long rkqs_ign;
  int maxstp;
  double minimp;
  double elastic;
  double epsilon;
  double lambda;
  double dt;
  double sig;
  double pert;
  int *y0fix;
  //damping section
  double wquer;

  //GNUplot
  FILE **gnuFp;
  long ngnu;
  long gnuindex;

    long maxPoints;

    long maxVar;
    Glob();
    virtual ~Glob();

};



namespace Diana {
  class DiffitESO : public DianaDAESO {
  protected:
      CapeBoolean           bInitialized;
      CapeBoolean           bSymbolicJacobian;
      virtual DianaDAESO*   createClone(const DianaDAESO& org);
  public:

      //DiffitESO();

    DiffitESO(const Common::Types::CapeString& _name,
              const Common::Types::CapeString& _description,
              Common::Types::CapeLong          _nVars,
              Common::Types::CapeLong          _nEqns,
              Common::Types::CapeLong          _nRealParams
              );
    DiffitESO(const DiffitESO& org);

    virtual ~DiffitESO();

    void Initialize()
            throw (Common::Error::ECapeUnknown, Common::Error::ECapeFailedInitialisation,
             Common::Error::ECapeOutOfResources, Common::Error::ECapeLicenceError,
             Common::Error::ECapeBadInvOrder);
      const Diana::DianaSparseArray& GetAllDiffJacobianValues()
          throw (Common::Error::ECapeUnknown, Common::Error::ECapeNoImpl);
 

    
protected:

    virtual void calculateResiduals();

  };

}

class SimpleModel{

protected:
    virtual void ode (double t, double *y, double *f, double *p)=0;    

public:

    virtual void ode_replace(double t, double *y, double *f, double *p)=0;
};
//using namespace Common::Parameter;




class DiffitModel : public Diana::DiffitESO //,public SimpleModel
{
protected:
const int NEQNS;       // number of ODEs    
const int NPARAMS;     // number of parameters
const char *DefModelDescription;
const double *DefParameters;
const double * DefYValues;
string * ParameterNames;
string * VariableNames;

public:
    //DiffitModel();
/* CTOR */
    DiffitModel(const char* description,
                const int NEQNS,
                //const int NOBS,
                const int NPARAMS,
                //const int NCONDS,
                const double * DefParameters,
                const double * MinParameters,
                const double * MaxParameters,
                const double * DefYValues,
                const double * MinYValues,
                const double * MaxYValues,
                const char ** _ParameterNames,
                const char ** _VariableNames
        );
    DiffitModel(const DiffitModel& org);


    virtual ~DiffitModel();
    virtual void Initialize()throw (Common::Error::ECapeBadInvOrder, Common::Error::ECapeLicenceError, Common::Error::ECapeOutOfResources, Common::Error::ECapeFailedInitialisation, Common::Error::ECapeUnknown);
    virtual void calculateResiduals();


    const Diana::DianaSparseArray& GetAllJacobianValues()
            throw (Common::Error::ECapeUnknown, Common::Error::ECapeNoImpl);

    const Diana::DianaSparseArray& GetAllParJacobianValues(const Common::Types::CapeArrayLong * parIndices=NULL)
        throw (Common::Error::ECapeUnknown, Common::Error::ECapeNoImpl);

virtual int observation (Glob *globs,ExperimentSystem *ex,double t, double *y,double *gy, double *p, double **dgdy,double **dgdp)=0;


//virtual void initInt(Glob *globs,ExperimentSystem *ex)=0;

/*
protected:
    long ode_fkmax,ode_fkount,ode_kmax,ode_kount,rkqs_ign;
    double ode_dxsav,*ode_fxp,*ode_xp,**ode_yp;
*/
protected:

    Diana::DianaDAESO* createClone();

virtual void jacobi (double t, double *y, double **jac, double *p)=0;
virtual void inhomo (double t, double *y, double **inh, double *p)=0;
virtual void inhomo1 (double t, double *y, double *inh1, double *p, int &jp)=0;



/* constraints */

//virtual void setNrConstraints(ExperimentSystem *data, int nP, double *parameters)=0;

    virtual void R2(Glob *globs,ExperimentSystem *ex,int computeDerivs)=0;

    virtual void R3(Glob *globs,ExperimentSystem *ex,int computeDerivs)=0;

    virtual void ode (double t, double *y, double *f, double *p)=0;        


    const DianaSparseArray& GuessJacobianStructure();
   
};






#endif