#include "sopnamsp.h"
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#ifdef __MWERKS__
#include "mwerksmath.h" // Portage mac D. Y.
#include "unixmac.h"
#endif
#include <string.h>
#include <string>

#include "pexceptions.h"
#include "generalfunc.h"

//================================================================
// GeneralFunction
//================================================================

/*!
  \class SOPHYA::GeneralFunction
  \ingroup NTools
   Classe de fonctions parametrees a plusieurs variables:
  \f$ F[x1,x2,x3,...:a1,a2,a3,...] \f$
*/

//////////////////////////////////////////////////////////////////////
/*!
 Creation d'une fonction de `nVar' variables et `nPar' parametres:
 \f$ F[x(1),x(2),x(3),...x(nVar) : a(1),a(2),a(3),...,a(nPar)] \f$
*/
GeneralFunction::GeneralFunction(unsigned int nVar, unsigned int nPar)
  : mNVar(nVar), mNPar(nPar)
{
 if(nVar==0 || nPar==0)
   throw SzMismatchError("GeneralFunction::GeneralFunction_Error: nVar==0 || nPar==0 !");
 deltaParm = new double[nPar];
 tmpParm   = new double[nPar];
}

GeneralFunction::~GeneralFunction()
{
 if(deltaParm != NULL) delete[] deltaParm;
 if(tmpParm != NULL) delete[] tmpParm;
}

//////////////////////////////////////////////////////////////////////
/*!
  Valeur et Derivees de la fonction (fct virtuelle par defaut).
*/
double GeneralFunction::Val_Der(double const xp[], double const* parm
                               , double *DgDpar)
{
 for(int i=0;i<mNPar;i++) tmpParm[i] = parm[i];
 {for(int i=0;i<mNPar;i++) {
   double d = deltaParm[i];
   if(d==0.) { DgDpar[i] = 0.; continue;}
   tmpParm[i] -= d/2.;
   double vg = Value(xp,tmpParm);
   tmpParm[i] += d;
   double vd = Value(xp,tmpParm);
   DgDpar[i] = (vd - vg)/d;
   tmpParm[i] = parm[i];
 }}
 return Value(xp, parm);
}

//////////////////////////////////////////////////////////////////////
/*!
  Definition de la variation du parametre numPar
  pour calculer la derivee automatiquement.
*/
void GeneralFunction::SetDeltaParm(int numPar, double d)
{
 if(numPar<0 || numPar >= mNPar)
   throw SzMismatchError("GeneralFunction::SetDeltaParm_Error: numPar<0 || numPar>=mNPar !");
 deltaParm[numPar] = d;
}


/*!
  Idem precedente fonction mais pour tous les parametres
*/
void GeneralFunction::SetDeltaParm(double const* dparam)
{
 for(int i=0;i<mNPar;i++) deltaParm[i] = dparam[i];
}

//================================================================
// GeneralFunc
//================================================================

/*!
  \class SOPHYA::GeneralFunc
  \ingroup NTools
  Classe de fonctions parametrees a plusieurs variables
  derivant de ``GeneralFunction''. Permet de definir
  une fonction a fiter sans passer par une classe derivee
  en utilisant l'ecriture courante du C. La fonction
  retournant les derivees par rapport aux parametres du fit
  peut etre egalement fournie (optionnel).
*/

/////////////////////////////////////////////////////////////////
/*!
  Createur, on passe le nom ``fun'' de la fonction a la mode C.
  On peut optionellement egalement passer le nom de la fonction
  ``funder'' qui retourne les valeurs des derivees par rapport
  aux parametres du fit.
  \verbatim
    ----------------------
    Exemple d'utilisation:
    ----------------------
    include "generalfit.h"
    ...
    double   gaussc(double const* x,double const* p);
    double d_gaussc(double const* x,double const* p,double* dp);
    ...
    main {
     ...
     // Fit SANS calcul automatique des derivees
     GeneralFunc      myfunc(2,7,gaussc);
     GeneralFit       myfit(&myfunc);
     ...
     myfit.Fit();
     ...
     // Fit AVEC calcul automatique des derivees
     GeneralFunc      myfunc(2,7,gaussc,d_gaussc);
     GeneralFit       myfit(&myfunc);
     ...
     myfit.Fit();
    }
    // Definition de la fonction a fitter a la mode C
    double gaussc(double const* x,double const* p)
    // Fonction: X=(x[0]-p[1])/p[3], Y=(x[1]-p[2])/p[4],
    //  f = p[0]*exp{-0.5*[X^2+Y^2-2*p[5]*X*Y]} + p[6]
    {
     double X = (x[0]-p[1])/p[3];
     double Y = (x[1]-p[2])/p[4];
     return p[0]*exp(-(X*X+Y*Y-2*p[5]*X*Y)/2)+p[6];
    }
    // Definition de la fonction des derivees / parametres
    // Cette fonction retourne aussi la valeur de la fonction a fitter.
    double d_gaussc(double const* x,double const* p,double* dp)
    {
     dp[0] = derivee de gaussc par rapport au parametre p[0]
     ...
     dp[6] = derivee de gaussc par rapport au parametre p[6]
     return gaussc(x,p);
    }
  \endverbatim
*/
GeneralFunc::GeneralFunc(unsigned int nvar, unsigned int npar
             , double (*fun) (double const*, double const*)
             , double (*funder) (double const*, double const*, double*) )
: GeneralFunction(nvar,npar), tmpFun(fun), tmpFunDer(funder)
{
}

GeneralFunc::~GeneralFunc()
{
}

double GeneralFunc::Value(double const xp[], double const* Par)
{
return tmpFun(xp,Par);
}

double GeneralFunc::Val_Der(double const xp[],double const* parm, double* DgDpar)
{
if(tmpFunDer) return tmpFunDer(xp,parm,DgDpar);
  else        return GeneralFunction::Val_Der(xp,parm,DgDpar);
}