#ifndef ROOT_TModel
#define ROOT_TModel
 
#include <TBuffer.h>
#include <TError.h>
#include <Rtypes.h>
#include <Math/GSLRndmEngines.h>
 
#include <iostream>
#include <math.h>

class TModel
{

private:
   static ROOT::Math::GSLRngRanLux48* fRng;

protected:
   UInt_t 	fDeg;		// degree of Legendre expansion
   Double_t	*fLeg;		// array of Legendre coefficients
 
public:
   TModel(UInt_t deg=5);
   TModel(const TModel& model);
   virtual ~TModel();

   Double_t* const operator[](UInt_t i) const;
   Double_t operator()(Double_t x, Double_t y) const;
   TModel& operator=(const TModel &source);
   TModel& operator+=(const TModel &source);
   TModel& operator-=(const TModel &source);
   TModel& operator*=(Double_t factor);
   TModel& operator/=(Double_t factor);
   TModel operator-() const;
   Bool_t operator==(const TModel &other) const;
   Bool_t operator!=(const TModel &other) const;
   UInt_t GetDegree() const;
   void SetDegree(const UInt_t deg);
   void Randomize(Double_t scale=1);
   void Mutate(Double_t scale=0.1, UInt_t order=UINT_MAX);
   void Mutate(const TModel& parent, Double_t scale=0.1, UInt_t order=UINT_MAX);

   friend TModel operator+(const TModel &v1, const TModel &v2);
   friend TModel operator-(const TModel &v1, const TModel &v2);
   friend TModel operator*(const TModel &v1, Double_t factor);
   friend TModel operator*(Double_t factor, const TModel &v1);
   friend TModel operator/(const TModel &v1, Double_t factor);
   friend TBuffer& operator>>(TBuffer &buf, TModel* &obj);
   friend TBuffer& operator<<(TBuffer &buf, const TModel* obj);
   friend std::istream& operator>>(std::istream &in, TModel* &obj);
   friend std::ostream& operator<<(std::ostream &out, const TModel* obj);
   void Print(Option_t *option="");
 
   ClassDef(TModel,1)  // Real 2D function on ([-1,1],[-1,1]) expressed by its Legendre expansion
};


//----- inlines ----------------------------------------------------------------
 
inline void TModel::SetDegree(const UInt_t deg)
{
   Double_t *leg = fLeg;
   fLeg = new Double_t[deg*deg];
   for (UInt_t i=0; i<deg; i++) {
      for (UInt_t j=0,k=i*fDeg; j<deg; j++) {
         fLeg[i*deg+j] = (i < fDeg && j < fDeg)? leg[k++] : 0;
      }
   }
   delete [] leg;
   fDeg = deg;
}

inline Double_t* const TModel::operator[](UInt_t i) const
{
   return &fLeg[i*fDeg];
}

inline TModel& TModel::operator=(const TModel &source)
{
   delete [] fLeg;
   fDeg = source.fDeg;
   UInt_t len=fDeg*fDeg;
   fLeg = new Double_t[len];
   for (UInt_t k=0; k<len; k++)
      fLeg[k] = source.fLeg[k];
   return *this;
}

inline TModel& TModel::operator+=(const TModel &source)
{
   UInt_t deg = fDeg;
   if (source.fDeg > fDeg) {
      deg = source.fDeg;
      SetDegree(deg);
   }
   for (UInt_t i=0,k=0; i<source.fDeg; i++) {
      for (UInt_t j=0; j<source.fDeg; j++) {
         fLeg[i*deg+j] += source.fLeg[k++];
      }
   }
   return *this;
}

inline TModel& TModel::operator-=(const TModel &source)
{
   UInt_t deg = fDeg;
   if (source.fDeg > fDeg) {
      deg = source.fDeg;
      SetDegree(deg);
   }
   for (UInt_t i=0,k=0; i<source.fDeg; i++) {
      for (UInt_t j=0; j<source.fDeg; j++) {
         fLeg[i*deg+j] -= source.fLeg[k++];
      }
   }
   return *this;
}

inline TModel& TModel::operator*=(Double_t factor)
{
   UInt_t len = fDeg*fDeg;
   for (UInt_t k=0; k<len; k++)
      fLeg[k] *= factor;
   return *this;
}

inline TModel& TModel::operator/=(Double_t factor)
{
   UInt_t len = fDeg*fDeg;
   for (UInt_t k=0; k<len; k++)
      fLeg[k] /= factor;
   return *this;
}

inline TModel TModel::operator-() const
{
   TModel model(*this);
   UInt_t len=fDeg*fDeg;
   for (UInt_t k=0; k<len; k++)
      model.fLeg[k] *= -1;
   return model;
}

inline Bool_t TModel::operator==(const TModel &other) const
{
   UInt_t len=fDeg*fDeg;
   if (fDeg != other.fDeg)
      return kFALSE;
   else for (UInt_t k=0; k<len; k++)
      if (fLeg[k] != other.fLeg[k])
         return kFALSE;
   return kTRUE;
}

inline Bool_t TModel::operator!=(const TModel &other) const
{
   return !(*this == other);
}

inline UInt_t TModel::GetDegree() const
{
   return fDeg;
}

inline TModel operator+(const TModel &v1, const TModel &v2)
{
   TModel res(v1);
   return res += v2;
}

inline TModel operator-(const TModel &v1, const TModel &v2)
{
   TModel res(v1);
   return res -= v2;
}

inline TModel operator*(const TModel &v1, Double_t factor)
{
   TModel res(v1);
   return res *= factor;
}

inline TModel operator*(Double_t factor, const TModel &v1)
{
   return v1*factor;
}

inline TModel operator/(const TModel &v1, Double_t factor)
{
   TModel res(v1);
   return res /= factor;
}


inline TBuffer& operator>>(TBuffer &buf, TModel* &obj)
{
   UInt_t deg;
   buf >> deg;
   if (deg) {
      obj = new TModel(deg);
      buf.ReadStaticArray((Double_t *)&obj->fLeg);
   }
   return buf;
}

inline TBuffer& operator<<(TBuffer &buf, const TModel* obj)
{
   buf << obj->fDeg;
   buf.WriteArray((Double_t *)&obj->fLeg, obj->fDeg*obj->fDeg);
   return buf;
}

inline std::istream& operator>>(std::istream &in, TModel* &obj)
{
   UInt_t deg;
   in >> deg;
   if (deg) {
      obj = new TModel(deg);
      UInt_t k=0;
      for (UInt_t i=0; i<deg; i++) {
         for (UInt_t j=0; j<deg; j++,k++) {
            in >> obj->fLeg[k];
         }
      }
   }
   return in;
}

inline std::ostream& operator<<(std::ostream &out, const TModel* obj)
{
   out << " " << obj->fDeg;
   UInt_t k=0;
   for (UInt_t i=0; i<obj->fDeg; i++) {
      for (UInt_t j=0; j<obj->fDeg; j++,k++) {
         out << " " << obj->fLeg[k];
      }
   }
   return out;
}

#endif