//---Author:  Richard T. Jones  07/12/2007
//---Version:  1.0

/*************************************************************************
 * Copyright(c) 1998, University of Connecticut, All rights reserved.    *
 * Author: Richard T. Jones, Asst. Prof. of Physics                      *
 *                                                                       *
 * Permission to use, copy, modify and distribute this software and its  *
 * documentation for non-commercial purposes is hereby granted without   *
 * fee, provided that the above copyright notice appears in all copies   *
 * and that both the copyright notice and this permission notice appear  *
 * in the supporting documentation. The author makes no claims about the *
 * suitability of this software for any purpose.                         *
 * It is provided "as is" without express or implied warranty.           *
 *************************************************************************/
//////////////////////////////////////////////////////////////////////////
//									//
// Michelson Interferometric Topography Analysis Package		//
//									//
// The present package implements tools for modeling the shape of a     //
// general 2D surface in terms of its Legendre expansion.		//
//									//
//////////////////////////////////////////////////////////////////////////
 
#include <iostream>
#include <Math/SpecFuncMathMore.h>

using namespace std;

#include <TModel.h>

ClassImp(TModel)

ROOT::Math::GSLRngRanLux48* TModel::fRng = 0;


TModel::TModel(UInt_t deg)
{
   fDeg = deg;
   fLeg = new Double_t[deg*deg];
   if (fRng == 0) {
      fRng = new ROOT::Math::GSLRngRanLux48;
      fRng->Initialize();
      fRng->SetSeed(0);
   }
}

TModel::TModel(const TModel& model)
{
   fDeg = model.fDeg;
   fLeg = new Double_t[fDeg*fDeg];
   UInt_t k=0;
   for (UInt_t i=0; i<fDeg; i++) {
      for (UInt_t j=0; j<fDeg; j++,k++) {
         fLeg[k] = model.fLeg[k];
      }
   }
   if (fRng == 0) {
      fRng = new ROOT::Math::GSLRngRanLux48;
      fRng->Initialize();
      fRng->SetSeed(0);
   }
}

TModel::~TModel()
{
   delete [] fLeg;
}

Double_t TModel::operator()(Double_t x, Double_t y) const
{
   Double_t *xfact = new Double_t[fDeg];
   Double_t *yfact = new Double_t[fDeg];
   for (UInt_t i=0; i<fDeg; i++) {
      xfact[i] = ROOT::Math::legendre(i,x);
      yfact[i] = ROOT::Math::legendre(i,y);
   }
   Double_t result = 0;
   UInt_t k=0;
   for (UInt_t i=0; i<fDeg; i++) {
      for (UInt_t j=0; j<fDeg; j++,k++) {
         result += fLeg[k]*xfact[i]*yfact[j];
      }
   }
   delete [] xfact;
   delete [] yfact;
   return result;
}

void TModel::Randomize(Double_t scale)
{
   UInt_t k=0;
   for (UInt_t i=0; i<fDeg; i++) {
      for (UInt_t j=0; j<fDeg; j++,k++) {
         fLeg[k] = fRng->Gaussian(scale);
      }
   }
}

void TModel::Mutate(Double_t scale, UInt_t order)
{
   UInt_t k=0;
   for (UInt_t i=0; i<fDeg; i++) {
      for (UInt_t j=0; j<fDeg; j++,k++) {
         if (i < order && j < order) {
            fLeg[k] += fRng->Gaussian(scale);
         }
      }
   }
}

void TModel::Mutate(const TModel& parent, Double_t scale, UInt_t order)
{
   UInt_t k=0;
   if (fDeg == parent.fDeg) {
      for (UInt_t i=0; i<fDeg; i++) {
         for (UInt_t j=0; j<fDeg; j++,k++) {
            if (i < order && j < order) {
               fLeg[k] = parent.fLeg[k] + fRng->Gaussian(scale);
            }
         }
      }
   }
   else {
      cerr << "TModel::Mutate error: parent and daughter models "
           << "have different dimensions!\n";
   }
}

void TModel::Print(Option_t *option)
{
   // PrUInt_t contents of the Legendre coefficient array as a square matrix

   cout << "TModel(";
   UInt_t k=0;
   for (UInt_t i=0; i<fDeg; i++) {
      UInt_t j;
      for (j=0; j<fDeg-1; j++,k++) {
         cout << fLeg[k] << ", ";
      }
      cout << fLeg[k] << endl;
   }
   cout << "); dimension " << fDeg << endl;
}

void TModel::Streamer(TBuffer &buf)
{
   // Put/get fDeg*fDeg Double_t values to/from stream buffer buf.

   if (buf.IsReading()) {
      buf >> fDeg;
      buf.ReadStaticArray(fLeg);
   } else {
      buf << fDeg;
      buf.WriteArray(fLeg, fDeg*fDeg);
   }
}