// $Id$
//-----------------------------------------------------------------------
//       The GSI Online Offline Object Oriented (Go4) Project
//         Experiment Data Processing at EE department, GSI
//-----------------------------------------------------------------------
// Copyright (C) 2000- GSI Helmholtzzentrum fuer Schwerionenforschung GmbH
//                     Planckstr. 1, 64291 Darmstadt, Germany
// Contact:            http://go4.gsi.de
//-----------------------------------------------------------------------
// This software can be used under the license agreements as stated
// in Go4License.txt file which is part of the distribution.
//-----------------------------------------------------------------------

#include "TGo4FitModelGaussN.h"

#include <iostream>

#include "TMath.h"
#include "TGo4FitParameter.h"

TGo4FitModelGaussN::TGo4FitModelGaussN()
{
}

TGo4FitModelGaussN::TGo4FitModelGaussN(const char *iName, Int_t iNDimension)
   : TGo4FitModel(iName, "N-dimensional Gaussian", kTRUE), fxIndexes(iNDimension)
{
   for (Int_t n = 0; n < iNDimension; n++)
      fxIndexes[n] = n;

   for (Int_t n = 0; n < iNDimension; n++)
      NewParameter(GetPosParName(n), "position", 1.);

   for (Int_t n = 0; n < iNDimension; n++)
      NewParameter(GetWidthParName(n), "width", .5);

   for (Int_t n1 = 0; n1 < iNDimension; n1++)
      for (Int_t n2 = n1 + 1; n2 < iNDimension; n2++)
         NewParameter(GetCovarParName(n1, n2), "covariation", 0.);
}

TGo4FitModelGaussN::~TGo4FitModelGaussN()
{
   if (Vect_mu) {
      delete Vect_mu;
      Vect_mu = nullptr;
   }
   if (Matr_sig) {
      delete Matr_sig;
      Matr_sig = nullptr;
   }
   if (Vect_x) {
      delete Vect_x;
      Vect_x = nullptr;
   }
   if (Vect_dx) {
      delete Vect_dx;
      Vect_dx = nullptr;
   }
}

TString TGo4FitModelGaussN::GetPosParName(Int_t naxis)
{
   TString res;
   res.Form("Pos%d", naxis);
   return res;
}

TString TGo4FitModelGaussN::GetWidthParName(Int_t naxis)
{
   TString res;
   res.Form("Width%d", naxis);
   return res;
}

TString TGo4FitModelGaussN::GetCovarParName(Int_t naxis1, Int_t naxis2)
{
   TString res;
   res.Form("Cov%d_%d", naxis1, naxis2);
   return res;
}

Bool_t TGo4FitModelGaussN::SetAxisNumbers(Int_t naxis)
{
   Int_t oldnaxis = GetAxisNumbers();
   if ((naxis < 0) || (oldnaxis == naxis))
      return kFALSE;
   if (naxis < oldnaxis) {
      for (Int_t n = naxis; n < oldnaxis; n++) {
         RemovePar(GetPosParName(n));
         RemovePar(GetWidthParName(n));
      }
      for (Int_t n1 = 0; n1 < oldnaxis; n1++)
         for (Int_t n2 = n1 + 1; n2 < oldnaxis; n2++)
            if ((n1 >= naxis) || (n2 >= naxis))
               RemovePar(GetCovarParName(n1, n2));
      fxIndexes.Set(naxis);
   } else {
      Int_t lastindx = 0;
      for (Int_t n = oldnaxis; n < naxis; n++) {
         lastindx = GetParIndex(FindPar(GetPosParName(n - 1))) + 1;
         NewParameter(GetPosParName(n), "position", 1., kFALSE, lastindx);
         lastindx = GetParIndex(FindPar(GetWidthParName(n - 1))) + 1;
         NewParameter(GetWidthParName(n), "width", 5., kFALSE, lastindx);
      }
      lastindx++;
      for (Int_t n1 = 0; n1 < naxis; n1++)
         for (Int_t n2 = n1 + 1; n2 < naxis; n2++) {
            TGo4FitParameter *par = FindPar(GetCovarParName(n1, n2));
            if (par)
               lastindx = GetParIndex(par) + 1;
            else
               NewParameter(GetCovarParName(n1, n2), "covariation", 0., kFALSE, lastindx++);
         }

      fxIndexes.Set(naxis);
      for (Int_t n = oldnaxis; n < naxis; n++) {
         Int_t test = 0;
         Bool_t find;
         do {
            find = kFALSE;
            for (Int_t n1 = 0; n1 < n; n1++)
               if (test == fxIndexes[n1]) {
                  find = kTRUE;
                  test++;
                  break;
               }
         } while (find);
         fxIndexes[n] = test;
      }
   }

   return kTRUE;
}

Bool_t TGo4FitModelGaussN::ResortIndexes(Int_t leaveaxis)
{
   Bool_t res = kFALSE;
   for (Int_t n = fxIndexes.GetSize() - 1; n >= 0; n--) {
      if (n == leaveaxis)
         continue;
      Int_t test = fxIndexes[n];
      Bool_t find, first = kTRUE;
      do {
         find = kFALSE;
         for (Int_t n1 = 0; n1 < fxIndexes.GetSize(); n1++)
            if ((n != n1) && (test == fxIndexes[n1])) {
               if (first)
                  test = 0;
               else
                  test++;
               find = kTRUE;
               first = kFALSE;
               break;
            }
      } while (find);
      fxIndexes[n] = test;
      res = res || (!first);
   }
   return res;
}

Int_t TGo4FitModelGaussN::GetPosParIndex(Int_t naxis)
{
   for (Int_t n = 0; n < fxIndexes.GetSize(); n++)
      if (naxis == fxIndexes[n])
         return 1 + n;
   return -1;
}

Int_t TGo4FitModelGaussN::GetWidthParIndex(Int_t naxis)
{
   for (Int_t n = 0; n < fxIndexes.GetSize(); n++)
      if (naxis == fxIndexes[n])
         return 1 + fxIndexes.GetSize() + n;
   return -1;
}

void TGo4FitModelGaussN::FillMuVector(TVectorD &Mu)
{
   Int_t ndim = fxIndexes.GetSize();
   Mu.ResizeTo(ndim);
   for (Int_t n = 0; n < ndim; n++)
      Mu(n) = GetPar(1 + n)->GetValue();
}

void TGo4FitModelGaussN::FillSigmaMatrix(TMatrixD &Sigma)
{
   Int_t ndim = fxIndexes.GetSize();
   Sigma.ResizeTo(ndim, ndim);
   for (Int_t n = 0; n < ndim; n++)
      Sigma(n, n) = GetPar(1 + ndim + n)->GetValue() * GetPar(1 + ndim + n)->GetValue();
   Int_t indx = 1 + 2 * ndim;
   for (Int_t n1 = 0; n1 < ndim; n1++)
      for (Int_t n2 = n1 + 1; n2 < ndim; n2++) {
         Double_t zn = GetPar(indx++)->GetValue();
         Sigma(n1, n2) = zn;
         Sigma(n2, n1) = zn;
      }
}

Bool_t TGo4FitModelGaussN::BeforeEval(Int_t NDimension)
{
   Par_ndim = fxIndexes.GetSize();
   Par_indx = fxIndexes.GetArray();
   for (Int_t i = 0; i < Par_ndim; i++)
      if (Par_indx[i] >= NDimension) {
         std::cout << "TGo4FitModelGaussN::   invalid index " << std::endl;
         return kFALSE;
      }

   Vect_mu = new TVectorD(Par_ndim);
   Matr_sig = new TMatrixD(Par_ndim, Par_ndim);
   Vect_x = new TVectorD(Par_ndim);
   Vect_dx = new TVectorD(Par_ndim);

   FillMuVector(*Vect_mu);
   FillSigmaMatrix(*Matr_sig);
   Double_t determ;
   Matr_sig->Invert(&determ);
   if (determ == 0.) {
      std::cout << "TGo4FitModelGaussN::   Invalid sigma matrice " << std::endl;
      AfterEval();
      return kFALSE;
   }

   return kTRUE;
}

Double_t TGo4FitModelGaussN::EvalN(const Double_t *v)
{
   for (Int_t n = 0; n < Par_ndim; n++)
      (*Vect_x)(n) = v[Par_indx[n]];
   *Vect_x -= *Vect_mu;
   *Vect_dx = *Vect_x;
   *Vect_dx *= *Matr_sig;
   Double_t z = *Vect_dx * *Vect_x; // see 11.1
   return TMath::Exp(-0.5 * z);
}

void TGo4FitModelGaussN::AfterEval()
{
   if (Vect_mu) {
      delete Vect_mu;
      Vect_mu = nullptr;
   }
   if (Matr_sig) {
      delete Matr_sig;
      Matr_sig = nullptr;
   }
   if (Vect_x) {
      delete Vect_x;
      Vect_x = nullptr;
   }
   if (Vect_dx) {
      delete Vect_dx;
      Vect_dx = nullptr;
   }
}

void TGo4FitModelGaussN::Print(Option_t *option) const
{
   TGo4FitModel::Print(option);
   std::cout << "   N-dimensional Gauss" << std::endl;
   std::cout << "   axis indexes:";
   for (Int_t i = 0; i < fxIndexes.GetSize(); i++)
      std::cout << " " << fxIndexes[i];
   std::cout << std::endl;
}