#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#ifdef __MWERKS__
#include "unixmac.h"
#endif


#include "fmath.h"

#include "bruit.h"
#include "nbrandom.h"
// #include "rancern.h"
// #include "hbook.h"
#ifdef OS_MACOSX
#include <limits.h>
#endif

//  Le code des classes   NoiseGenerator  RWalkNoise

//++ 
// Class	NoiseGenerator
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links	Childs
//
//    RWalkNoise OOFNoise EXPNoise MemNoise  SumNoise
//     
//-- 
//++
// Titre	Constructors
//--

// Definition de fonctions pour float - a changer en utilisant cmath
//                              Reza Juillet 2002
inline float _myfabsf(float x) { return ((float)(fabs((double)x)) ); }
inline float _mysqrtf(float x) { return ((float)(sqrt((double)x)) ); }
 
/* --Methode-- */
//++
NoiseGenerator::NoiseGenerator(float sigma)
//
//--
{
if (sigma < 0.)  sigma = 1.;
mNCoups = 0;
mSigma = sigma;
//printf("-- NoiseGenerator::NoiseGenerator(%g) (Constructeur) ---\n", sigma);
}

//++
// Titre	Destructor
//--
/* --Methode-- */
//++
NoiseGenerator::~NoiseGenerator()
//
//--
{
    //printf("-- NoiseGenerator::~NoiseGenerator() (Destructeur) --- \n");
}


//++
//
// inline  unsigned long int  NoiseGenerator::NCoups() 
//--
//++
// Titre	Public Methods
//--

/* --Methode-- */
//++
float NoiseGenerator::Noise()
//
//--
{
mNCoups++;
return(NorRand()*mSigma);
}
//++
// Class	RWalkNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links	Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre	Constructor
//--

/* --Methode-- */
//++
RWalkNoise::RWalkNoise(float sigma)
    : NoiseGenerator(sigma)
//
//--
{
mState = 0.;
//printf("--  RWalkNoise::RWalkNoise(%g) (Constructeur) ---\n", sigma);
}


//++
// Titre	Destructor
//--
/* --Methode-- */
//++
RWalkNoise::~RWalkNoise()
//
//--
{
//printf("--  RWalkNoise::~RWalkNoise (Destructeur) ---\n");

}

//++
// Titre	Public Methods
//--
/* --Methode-- */
//++
float RWalkNoise::Noise()
//
//--
{
mState += NoiseGenerator::Noise();
return(mState);
}

//++
// Class	OOFNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links	Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre	Constructor
//--

/* --Methode-- */
//++
OOFNoise::OOFNoise(float sigma, int typacf, int mem, float tau)
    : NoiseGenerator(sigma)
//
//--
{
if (typacf != ACF_Exp)   typacf = ACF_Exp;
// Fonction d'autocorrelation donnant un bruit en 1/f
//                  
mTypACF = typacf;
if (mem < 1)  mem = 1;
mMemL = mem;
if (tau < 1.e-6)  tau = 1.e-6;
mTau = tau;
mState = new double[mem];
mDyn = new double[mem];
int i;
mDyn[0]=1.;
mState[0]=0.;
for(i=1; i<mem; i++) {
  mState[i] = 0.;   mDyn[i] = 1./sqrt((double)i);
}
//printf("--  OOFNoise::OOFNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
//       sigma, typacf, mem, tau);
}


//++
// Titre	Destructor
//--
/* --Methode-- */
//++
OOFNoise::~OOFNoise()
//
//--
{
delete[] mState;
delete[] mDyn;
// printf("--  OOFNoise::~OOFNoise (Destructeur) ---\n");
}


//++
// Titre	Public Methods
//--
/* --Methode-- */
//++
float OOFNoise::Noise()
//
//--
{
int i;
double rn;

for(i=(mMemL-1); i>0; i--)   mState[i] = mState[i-1];
mState[0] = NoiseGenerator::Noise();
rn = 0.;
for(i=0; i< mMemL; i++)  rn += mDyn[i]*mState[i];
return(rn);
}

/* --Methode-- */
//++
void OOFNoise::Print()
//
//--
{
int i,j;
printf("OOFNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
printf("OOFNoise::Print() Vecteur de Dynamique / Etat : \n");
for(i=0; i<mMemL-1; i+=2) 
  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
return;
}
//++
// Class	EXPNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links	Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre	Constructor
//--

/* --Methode-- */
//++
EXPNoise::EXPNoise(float sigma, int typacf, int mem, float tau)
    : NoiseGenerator(sigma)
//
//--
{
if (typacf != ACF_Exp)   typacf = ACF_Exp;
// Fonction d'autocorrelation exponentiel disponible uniquement
//                  Reza , Juillet 97
mTypACF = typacf;
if (mem < 1)  mem = 1;
mMemL = mem;
if (tau < 1.e-6)  tau = 1.e-6;
mTau = tau;
mState = new double[mem];
mDyn = new double[mem];
int i;
for(i=0; i<mem; i++) {
  mState[i] = 0.;   mDyn[i] = exp(-(double)i/tau);
}
//printf("--  EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
//       sigma, typacf, mem, tau);
}


//++
// Titre	Destructor
//--
/* --Methode-- */
//++
EXPNoise::~EXPNoise()
//
//--
{
delete[] mState;
delete[] mDyn;
//printf("--  EXPNoise::~EXPNoise (Destructeur) ---\n");
}


//++
// Titre	Public Methods
//--
/* --Methode-- */
//++
float EXPNoise::Noise()
//
//--
{
int i;
double rn;

for(i=(mMemL-1); i>0; i--)   mState[i] = mState[i-1];
mState[0] = NoiseGenerator::Noise();
rn = 0.;
for(i=0; i< mMemL; i++)  rn += mDyn[i]*mState[i];
return(rn);
}

/* --Methode-- */
//++
void EXPNoise::Print()
//
//--
{
int i,j;
printf("EXPNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
printf("EXPNoise::Print() Vecteur de Dynamique / Etat : \n");
for(i=0; i<mMemL-1; i+=2) 
  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   

return;
}

//++
// Class	MemNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links	Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre	Constructor
//--


// ----------------      $CHECK$   Reza 1/12/99   ------------
// ----- MAJ MemNoise et SumNoise / version F. Couchot (~couchot/CoSa/Samba/bruit.cc  )
// -----------------------------------------------------------

/* --Methode-- */
//++
MemNoise::MemNoise(float sigma, int mem, float tau, int ava)
    : NoiseGenerator(_mysqrtf(_myfabsf(tau))*sigma)
//
//--
{
/* on tire les instants des impulsions successives selon une
loi de distance exponentielle de duree tau
*/
const int place = 1000;
mTypACF = ACF_Exp;
if (mem < 1)  mem = 1;
mMemL = place;
mCoupHaut=mem;
if (tau < 1.)  tau = 1.;
mTau = tau;
mStPos = new float[place];
mStNeg = new float[place];
mTePos = new float[place];
mTeNeg = new float[place];
// printf("apres reservation\n");
int i;
for(i=0; i<place; i++)
{
mStNeg[i] = 0.; mStPos[i] = 0.;
}
// printf("apres initialisation tableaux \n");

mNappel=0;
mNtirage=1;
mNappLast=0;
mMemPos=1;
mMemNeg=1;
mduree=0.;
mTePos[0]=0.;
mTeNeg[0]=0.;
mTdernier=0.;

Avance(ava);
// printf("fin du createur MemNoise \n");
//printf("--  EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
//       sigma, typacf, mem, tau);
}


//++
// Titre	Destructor
//--
/* --Methode-- */
//++
MemNoise::~MemNoise()
//
//--
{
delete[] mStPos;
delete[] mStNeg;
delete[] mTePos;
delete[] mTeNeg;
//printf("--  EXPNoise::~EXPNoise (Destructeur) ---\n");
}

//++
// Titre	Public Methods
//--
/* --Methode-- */
//++
float MemNoise::Noise()
//
//--
{
return(Avance(0));
}

/* --Methode-- */
//++
float MemNoise::Avance(long Asauter)
//
//--
{
int i,j;
float boum;
float rn = 0;
float Somme,SbaryT,Poids,BaryT,Approx,Erreur;
int nReduit;
const float epsilon=.05;
// printf("entree dans methode Noise de MemNoise \n");
long mNsaut;
for (mNsaut = Asauter ; mNsaut >= 0 ; mNsaut--)
    {  // debut mNsaut
	// printf("mNsaut = %ld\n",mNsaut);
mNappel++;
mTdernier += 1.;
for(i=0; i<mMemPos; i++) mTePos[i] += 1.;
for(i=0; i<mMemNeg; i++) mTeNeg[i] += 1.;

if (mduree < mTdernier)
    {
	/* la duree depuis le dernier tirage est maintenant ecoulee
	   on commence par remettre a jour le temps ecoule entre la derniere
	   impulsion (celle qu'on va tirer maintenant) et l'instant present  */
	   
	   mTdernier = mTdernier - mduree;
        
	/* puis on tire une nouvelle impulsion de bruit */

	   boum =  NoiseGenerator::Noise();

        /* on decale la memoire d'une case vers le passe 
	   puis on met a jour la case [0] */

if(boum>0)
    {
	for(i=mMemPos; i>0; i--)
	    {
		mStPos[i] = mStPos[i-1];
		mTePos[i] = mTePos[i-1];
	    }
	mTePos[0] = mTdernier;
	mStPos[0] = boum;
	if (mMemPos<mMemL) mMemPos++;
/* a chaque tirage, on regroupe les impulsions precedentes */
for(i=mMemPos-2 ; i>0 ; i--)
    {
	if(mStPos[i]==0.)continue;
	//	printf("i,temps,bruit %d,%g,%g \n",i,mTePos[i],mStPos[i]);
	Somme = mStPos[i]/_mysqrtf(mTePos[i]);
	    SbaryT = mStPos[i] * mTePos[i];
	    Poids = mStPos[i];
	    for(j=i-1 ; j>=0 ; j--)
		{
		    //		printf("  j,temps,bruit %d,%g,%g \n",j,mTePos[j],mStPos[j]);
		Somme += mStPos[j]/_mysqrtf(mTePos[j]);
		    Poids += mStPos[j];
		    SbaryT += mStPos[j] * mTePos[j];
		    BaryT = SbaryT / Poids;
		    Approx = Poids / _mysqrtf(BaryT);
		    Erreur = (Approx - Somme) / Somme;
		    //		    		    printf("i,j = %d %d Approx = %g \n",i,j,_myfabsf(Erreur));
		    //		    printf("Temps,Bruit = %g %g \n",BaryT,Poids);
		    if(_myfabsf(Erreur)<epsilon)
			{
			    mStPos[i]=Poids;
			    mTePos[i]=BaryT;
			    mStPos[j]=0.;
			}
		    else break;
		}
    }
nReduit=0;
for(i=0; i<mMemPos-1 ; i++)
    {
	if(mStPos[i]!=0)
	    {
		mStPos[nReduit]=mStPos[i];
		mTePos[nReduit]=mTePos[i];
		nReduit++;
	    }
    }
mMemPos=nReduit+1;

    }
else
    {
	for(i=mMemNeg; i>0; i--)
	    {
		mStNeg[i] = mStNeg[i-1];
		mTeNeg[i] = mTeNeg[i-1];
	    }
	mTeNeg[0] = mTdernier;
	mStNeg[0] = boum;
	if (mMemNeg<mMemL) mMemNeg++;
    

/* idem pour negatifs */
for(i=mMemNeg-2 ; i>0 ; i--)
    {
	if(mStNeg[i]==0.)continue;
	//	printf("i,temps,bruit %d,%g,%g \n",i,mTeNeg[i],mStNeg[i]);
	Somme = mStNeg[i]/_mysqrtf(mTeNeg[i]);
	    SbaryT = mStNeg[i] * mTeNeg[i];
	    Poids = mStNeg[i];
	    for(j=i-1 ; j>=0 ; j--)
		{
		    //		printf("  j,temps,bruit %d,%g,%g \n",j,mTeNeg[j],mStNeg[j]);
		Somme += mStNeg[j]/_mysqrtf(mTeNeg[j]);
		    Poids += mStNeg[j];
		    SbaryT += mStNeg[j] * mTeNeg[j];
		    BaryT = SbaryT / Poids;
		    Approx = Poids / _mysqrtf(BaryT);
		    Erreur = (Approx - Somme) / Somme;
		    //		    printf("i,j = %d %d Approx = %g \n",i,j,Erreur);
		    //		    printf("Temps,Bruit = %g %g \n",BaryT,Poids);
		    if(_myfabsf(Erreur)<epsilon)
			{
			    mStNeg[i]=Poids;
			    mTeNeg[i]=BaryT;
			    mStNeg[j]=0.;
			}
		    else break;
		}
    }
nReduit=0;
for(i=0; i<mMemNeg-1 ; i++)
    {
	if(mStNeg[i]!=0)
	    {
		mStNeg[nReduit]=mStNeg[i];
		mTeNeg[nReduit]=mTeNeg[i];
		nReduit++;
	    }
    }
mMemNeg=nReduit+1;
    }
/* on tire la duree a attendre avant la prochaine impulsion de bruit */
// $CHECK$ avec Francois - rndm() remplace par Reza 10/01/99 
//  passage en drand le 1/12/99
mduree = -mTau*log(drand01());  

// if(mNtirage<mMemL) mNtirage++;
mNtirage++;
/*
int idn=1;
float pds=1.;
float nul=0.;
hfill_(&idn,&mduree,&nul,&pds);
*/
// printf("Ntirage %d,mMemPos,Neg %d,%d, Suppress = %d \n",mNtirage,mMemPos,mMemNeg, Suppress);
    }
/* On calcule le bruit total */

//printf("mMemPos,Neg,%d %d\n",mMemPos,mMemNeg);
//for(i=0;i<mMemPos;i++) printf("mTePos %g \n",mTePos);
//for(i=0;i<mMemNeg;i++) printf("mTeNeg %g \n",mTeNeg);
if(mNsaut == 0)
    {
	// printf("coucou\n");
rn = 0.;
int indlim=0;
for(i=0; i < mMemPos-1; i++) 
{
if(mTePos[i] > 1.) break; 		// Modif F.C et D.Y. 03/2000
rn += mStPos[i];
indlim++;
}
// float ajout=0.;
for(i=indlim; i < mMemPos-1; i++)  rn += mStPos[i]/_mysqrtf(mTePos[i]);
indlim=0;
for(i=0; i < mMemNeg-1; i++) 
{
if(mTeNeg[i] > 1.) break;		// Modif F.C et D.Y. 03/2000
rn += mStNeg[i];
indlim++;
}
for(i=indlim; i < mMemNeg-1; i++) rn += mStNeg[i]/_mysqrtf(mTeNeg[i]);
}
}
return(rn);
}
/* --Methode-- */
//++
int MemNoise::Print()
//
//--
{
int i,j,rc=0;
printf("MemNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
printf("mNtirage= %d , mMemPos= %d , mMemNeg= %d \n",mNtirage,mMemPos,mMemNeg);
// printf("MemNoise::Print() Vecteur de Dynamique / Etat : \n");
/* for(i=0; i<mMemL-1; i+=2) 
  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
	 */
return(rc);
}

//++
// Class	SumNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//    Generateur de bruit blanc + 1/f avec fknee

//--
//++
//
// Links	Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre	Constructor
//--
// Generateur de bruit blanc + 1/f avec fknee


/* --Methode-- */
//++
SumNoise::SumNoise(float fknee, float sig, int mem, float tau)
    : NoiseGenerator(sig)
//
//--
{
// Reza 27/01/98 :
// Je calcule les sigma du MemNoise d'apres la formule de Francois 
// mn = new MemNoise(sig*fknee);  
// mn = new MemNoise(1.4142*sig*_mysqrtf(fknee),1,10.,1000000); // modif appel FC 13 mai 98.  
mn = new MemNoise(1.4142*sig*sqrt((double)fknee),mem,tau,12000); // FC 13 mai 98. -> Reza 1/12/99 
}


//++
// Titre	Destructor
//--
/* --Methode-- */
//++
SumNoise::~SumNoise()
//
//--
{
delete mn;
}

//++
// Titre	Public Methods
//--
/* --Methode-- */
//++
float SumNoise::Noise()
//
//--
{
return(NoiseGenerator::Noise() + mn->Noise());
}