/****************************************************************/
/* Developpement : Stephane.Girard@imag.fr                      */
/* 2001                                                         */
/****************************************************************/

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


#define MAXN     1000
#define MAXNBZ   1000
#define MAXNBH   100

int	n,nbh,nbz,nbx,i,j,indh,indx,indz,indy,methode;
float	vectx[MAXN],vecty[MAXN],vecth[MAXNBH],vectz[MAXNBZ],critere[MAXNBH];
float	essx[MAXNBH],essy[MAXNBH];
float	qalpha[MAXN],q1malpha[MAXN],qalphalis[MAXN],q1malphalis[MAXN];
float 	deltax,hopt,xmin,xmax,hmin,hmax,ymin,ymax,ecart,fni,fni_j,criteremin;
float	ordre1,ordre2,h1opt,h2opt,hmean,htmp1,htmp2,temp,utmp,htmp,coeff,u01;
float	hfin1,hfin2,noyau;
int	alpha,nbessais,lissage;
double	cumul1,cumul2,cumuld1,cumuld2;
char	nomgrille[255],nomgrillebis[255];
char	nomfichier[255],nomfichier2[255],nomfichier3[255],nomfichier4[255];
char	nomfichiergen[255],nomparam[255],nomfichier5[255],nomfichierbis[255];
char	nomfichierter[255],repdon[255],repres[255],nomfichier3lis[255],nomfichier4lis[255];
float	fnitest[MAXNBZ],xtest,htest;


/**************************** Tri par insertion  *****************/
void tri_inser(U,V,k)
float U[],V[];
int k;

/* U tableau a trier ,et k = dimension du tableau */
{
        float min,aux;
        int i,j,indice;
        {
                for (j=1 ; j <= k; j = j+1)
                {
                        min = 1.0e8;
                        /* precondition = U[j] >= 0 */
                        for (i=j ; i <= k; i= i+1)
                        {
                                if (U[i] < min)
                                {
                                        indice = i;
                                        min = U[i];
                                }
                        }
                        aux = U[j];
                        U[j] = U[indice];
                        U[indice] = aux;
                        aux = V[j];
                        V[j] = V[indice];
                        V[indice] = aux;
                }
        }
}



/**************************** Lecture fichier parametres *****************/

void Lit_Fichier_Parametres(nomfich)

char    *nomfich;

{
        char    temp[255];
        FILE    *ff;

ff=fopen(nomfich,"r");
if (ff==0)
{
        printf("Nom fichier parametres incorrect.\n");
        exit(1);
}

fscanf(ff,"%s %d",temp,&methode);
fscanf(ff,"%s %d",temp,&lissage);
fscanf(ff,"%s %s",temp,repdon);
fscanf(ff,"%s %s",temp,repres);
fscanf(ff,"%s %s",temp,nomfichier);
fscanf(ff,"%s %d",temp,&alpha);
fscanf(ff,"%s %d",temp,&nbessais);
fscanf(ff,"%s %s",temp,nomgrille);
fscanf(ff,"%s %d",temp,&nbh);

if (nbh==0)
{ 
	if (methode==3)
		 fscanf(ff,"%s %f %f",temp,&h1opt,&h2opt);
	else 	 fscanf(ff,"%s %f",temp,&hopt);
}	
fclose(ff);



if (nbh>=MAXNBH)
{
        printf("Nb de H trop grand.\n");
        exit(1);
}

if (nbessais>=MAXNBH)
{
        printf("Nb essais trop grand.\n");
        exit(1);
}

if (methode!=1 && methode!=2 && methode!=3)
{
	printf("Methode %d inconnue.\n",methode);
	exit(1);
}

}

/********************* Ecriture fichier de 2 colonnes */

ecrit_fichier(nom,n,xtab,ytab)

	char *nom;
	float *xtab,*ytab;
	int n;

{
int i;
FILE *f,  *fopen();

f=fopen(nom,"w");
if (f==NULL)
{
      printf("probleme fichier sortie\n");
      exit(1);
}

for (i=1;i<=n;i++)
        fprintf(f,"%f %f\n",xtab[i],ytab[i]);

fclose(f);
}

/********************* Lecture fichier de 2 colonnes */

lit_fichier(nom,n,xtab,ytab)

	char *nom;
	int *n;
	float	*xtab,*ytab;

{
int i,nbpts;
float xtmp,ytmp;
FILE *f,*fopen() ;

f=fopen(nom,"r");
if (f==NULL)
{
        printf("probleme fichier entree %s\n",nom);
        exit(1);
}

nbpts=0;

while (fscanf(f,"%f %f",&xtmp,&ytmp)!=EOF)
        nbpts++;

fclose(f);

if (nbpts==0)
{
        printf("fichier vide!\n");
        exit(1);
}
if (nbpts>=MAXN)
{
        printf("nb maximum de pts autorise %d (ici %d)\n",MAXN,nbpts);
        exit(1);
}
f=fopen(nom,"r");
for (i=1;i<=nbpts;i++)
{
	fscanf(f,"%f %f",&xtmp,&ytmp);
	xtab[i]=xtmp;
	ytab[i]=ytmp;
}

fclose(f);
*n=nbpts;
}


/********************* Lecture fichier de 1 colonne */

lit_colonne(nom,n,xtab)

	char *nom;
	int *n;
	float	*xtab;

{
int i,nbpts;
float xtmp;
FILE *f,*fopen() ;

f=fopen(nom,"r");
if (f==NULL)
{
        printf("probleme fichier entree %s\n",nom);
        exit(1);
}

nbpts=0;

while (fscanf(f,"%f",&xtmp)!=EOF)
        nbpts++;

fclose(f);

if (nbpts==0)
{
        printf("fichier vide!\n");
        exit(1);
}
if (nbpts>=MAXNBZ)
{
        printf("nb maximum de pts autorise %d (ici %d)\n",MAXNBZ,nbpts);
        exit(1);
}
f=fopen(nom,"r");
for (i=1;i<=nbpts;i++)
{
	fscanf(f,"%f",&xtmp);
	xtab[i]=xtmp;
}

fclose(f);
*n=nbpts;
}


/******************** Fonction de repartition de N(0,1) ******************/

float Phi01(x)
float x ;
{
double p=0.2316419, b1=0.319381530,b2=-0.356563782,
b3=1.781477937, b4= -1.821255978, b5=1.330274429 ;
double t, aux, denaux ;
   if (x < 0.) {
               aux= Phi01(-x) ;
               return(1. - aux);

              }
   else
   {
      t= 1./(1.+p*x) ;

      denaux = (1./sqrt(6.28318))*exp(-pow(x,2.)/2 ) ;
      return(1. - denaux * (b1*t + b2*pow(t,2.) + b3*pow(t,3.) + b4*pow(t,4.) +
b5*pow(t,5.))) ;

   }
}

/******************** Fonction de repartition inverse de N(0,1) ******************/

float Phi01_1(x)
float x;
{
        float   milint,minint,maxint;
        int     it;

minint=-10;
maxint=10;
milint=(minint+maxint)/2;

for (it=1;it<=20;it++)
{
        if (Phi01(milint)<x)
        {
                minint=milint;
                milint=(minint+maxint)/2;
        }
        else
        {
                maxint=milint;
                milint=(minint+maxint)/2;
        }
}
return(milint);
}

/******************* Densite N(0,1) *********************************/

float	fi01(x)

	float x;

{
	float	tmp;

	tmp=exp(-(double)x*x/2)/sqrt(2*M_PI);
	return(tmp);
}


/********************** "Check function" **************/

float	rho(ordre,x)

	float	ordre,x;

{
	float 	res;

	if (x>=0) 
		res=ordre*x;
	else	res=(ordre-1)*x;

	return(res);
}

/******************* Calcul du noyau gaussien */

float	K(x)

	float x;

{
	float	tmp;

	tmp=fi01(x);
	return(tmp);
}

/******************* Calcul de la primitive du noyau uniforme */

float   Omega(x)

        float x;

{
        float   temp;

	if (x<=-1) temp=0;
	else	
	{
		if (x<=1)	temp=0.5*(x+1);
		else		temp=1;
	}
        return(temp);
}

/************** Fonction auxiliaire pour le calcul de l'estimateur 2.   */

float   Cout_estimateur_2(x,xdon,ydon,hn,nbpts,ordre,a)

        float   x,hn,a,ordre;
        float   *xdon,*ydon;
        int     nbpts;

{
        int     i;
        float   noyau,res,coeff;
        double  numer;

        numer=0;

	for (i=1;i<=nbpts;i++)
        {
                noyau=K((x-xdon[i])/hn);
		coeff=rho(ordre,ydon[i]-a);
                numer+=coeff*noyau;
        }
	res=(float)numer;
	
	return(res);
}

/************** Estimateur simple noyau fonction d'esperance conditionnelle    */
/************** (prive de l'observation numero "enleve")                       */

float   Estimateur_Regression(x,xdon,ydon,hn,nbpts,enleve)

        float   x,hn;
        float   *xdon,*ydon;
        int     nbpts,enleve;

{
        int     i;
        float   noyau,res;
        double  numer,denom;

        numer=0;
        denom=0;

        for (i=1;i<enleve;i++)
        {
                noyau=K((x-xdon[i])/hn);
                numer+=noyau*ydon[i];
                denom+=noyau;
        }
        for (i=enleve+1;i<=nbpts;i++)
        {
                noyau=K((x-xdon[i])/hn);
                numer+=noyau*ydon[i];
                denom+=noyau;
        }

        res=numer/denom;
        return(res);
}


/************** Estimateur simple noyau fonction de repartition conditionnelle    */
/************** (prive de l'observation numero "enleve")   		          */

float	EstimateurFn_1(y,x,xdon,ydon,hn,nbpts,enleve)

	float	x,y,hn;
	float 	*xdon,*ydon;
	int	nbpts,enleve;

{
	int	i;
	float	noyau,res;
	double	numer,denom;

	numer=0;
	denom=0;

	for (i=1;i<enleve;i++)
	{
		noyau=K((x-xdon[i])/hn);	
		if (ydon[i]<=y) numer+=noyau;
		denom+=noyau;
	}
	for (i=enleve+1;i<=nbpts;i++)
	{
		noyau=K((x-xdon[i])/hn);	
		if (ydon[i]<=y) numer+=noyau;
		denom+=noyau;
	}

	res=numer/denom;
	return(res);
}

/************** Estimateur double noyau fonction de repartition conditionnelle    */

float   EstimateurFn_2(y,x,xdon,ydon,h1n,h2n,nbpts)

        float   x,y,h1n,h2n;
        float   *xdon,*ydon;
        int     nbpts;

{
        int     i;
        float   noyau,res,coeff;
        double  numer,denom,temp;

        numer=0;
        denom=0;

        for (i=1;i<=nbpts;i++)
        {
                noyau=K((x-xdon[i])/h1n);
		temp=(y-ydon[i])/h2n;
		coeff=Omega(temp);
                numer+=coeff*noyau;
                denom+=noyau;
        }

        res=numer/denom;
        return(res);
}

/************* Calcul du quantile conditionnel, methode 1 ************************/

float calcul_quantile_1(ordre,h,n,xdon,ydontrie,xcalcul,nbit)

	float	ordre,h,xcalcul;
	int	n,nbit;
	float	*xdon,*ydontrie;


{
	float	milint,minint,maxint;
	int	it;

minint=ydontrie[1]-0.1;
maxint=ydontrie[n]+0.1;
milint=(minint+maxint)/2;

for (it=1;it<=nbit;it++)
{
	if (EstimateurFn_1(milint,xcalcul,xdon,ydontrie,h,n,0)<ordre)
	{
		minint=milint;
		milint=(minint+maxint)/2;
	}
	else
	{
		maxint=milint;
		milint=(minint+maxint)/2;
	}
}
return(milint);
		
}

/************* Calcul du quantile conditionnel, methode 2 ************************/

float calcul_quantile_2(ordre,h,n,xdon,ydontrie,xcalcul,nbit,essx,essy)

        float   ordre,h,xcalcul;
        int     n,nbit;
        float   *xdon,*ydontrie,*essx,*essy;


{
	int	i;
	float	crit,critmin,res;

	critmin=10000000;
	for (i=1;i<=nbit;i++)
	{
		crit=Cout_estimateur_2(xcalcul,xdon,ydontrie,h,n,ordre,essx[i]);
		essy[i]=crit;

		if (crit<=critmin)
		{
			critmin=crit;
			res=essx[i];
		}
	}
	return(res);
}

/************* Calcul du quantile conditionnel, methode 3 ************************/

float calcul_quantile_3(ordre,h1,h2,n,xdon,ydontrie,xcalcul,nbit)

        float   ordre,h1,h2,xcalcul;
        int     n,nbit;
        float   *xdon,*ydontrie;


{
        float   milint,minint,maxint;
        int     it;

minint=ydontrie[1]-0.1;
maxint=ydontrie[n]+0.1;
milint=(minint+maxint)/2;

for (it=1;it<=nbit;it++)
{
        if (EstimateurFn_2(milint,xcalcul,xdon,ydontrie,h1,h2,n)<ordre)
        {
                minint=milint;
                milint=(minint+maxint)/2;
        }
        else
        {
                maxint=milint;
                milint=(minint+maxint)/2;
        }
}
return(milint);
               
}

/************ Recherche H optimal, methode 1 *****************************/

float Recherche_hopt_1(nbh,n,vectx,vecty,vecth,critere)

	int	nbh,n;
	float	*vectx,*vecty,*critere,*vecth;

{
	int	i,indh,indx,indy;
	float	ecart,criteremin,hopt;
	

        criteremin=1e10;

	for (indh=1;indh<=nbh;indh++)
        {
                printf("Essai h No %d / %d = %f \n",indh,nbh,vecth[indh]);

                ecart=0;
                for (indx=1;indx<=n;indx++)
                {
                        printf("   Abscisse No %d / %d\r",indx,n);
			fflush(stdout);
                        for (indy=1;indy<=n;indy++)
                        {
                                fni=EstimateurFn_1(vecty[indy],vectx[indx],
                                vectx, vecty,vecth[indh],n,indx);

                                if (indx<indy)  fni_j=1;
                                else            fni_j=0;

                                ecart+=(fni-fni_j)*(fni-fni_j);
                        }
                }
                critere[indh]=ecart;

                printf("\n --> critere = %f\n",ecart);

                if (ecart<=criteremin)
                {
                        criteremin=ecart;
                        hopt=vecth[indh];
                }
        }
	return(hopt);
}

/************ Recherche H mean, methodes 2 et 3 *****************************/

float Recherche_hmean(nbh,n,vectx,vecty,vecth,critere)

        int     nbh,n;
        float   *vectx,*vecty,*critere,*vecth;

{
        int     i,indh,indx,indy;
        float   ecart,criteremin,hmean;


        criteremin=1e10;

        for (indh=1;indh<=nbh;indh++)
        {
                printf("Essai h No %d / %d = %f \n",indh,nbh,vecth[indh]);

                ecart=0;
                for (indx=1;indx<=n;indx++)
                {
                        printf("   Abscisse No %d / %d\r",indx,n);
			fflush(stdout);
                        fni=Estimateur_Regression(vectx[indx],
                                    vectx, vecty,vecth[indh],n,indx);

                        fni_j=vecty[indx];

                        ecart+=(fni-fni_j)*(fni-fni_j);
                }
                critere[indh]=ecart;

                printf("\n --> critere = %f \n",ecart);

                if (ecart<=criteremin)
                {
                        criteremin=ecart;
                        hmean=vecth[indh];
                }
        }
        return(hmean);
}


/************************************************************************/


main(argc,argv)

int     argc;
char    **argv;
{


        if (argc!=2)
        {
                printf("USAGE: nom_programme <nom_fichier_param>\n");
                exit(1);
        }

        strcpy(nomparam,argv[1]);
        printf("Utilisation de fichier de parametres %s\n",nomparam);
        Lit_Fichier_Parametres(nomparam);

	strcpy(nomgrillebis, repdon);
	strcat(nomgrillebis, nomgrille);
	lit_colonne(nomgrillebis,&nbz,vectz);
	if (nbz>=MAXNBZ)
	{
		printf("Nb de Z trop grand.\n");
		exit(1);
	}
	printf("Evaluation du quantile en %d points.\n",nbz);

	strcpy(nomfichierter, repdon);
	strcat(nomfichierter, nomfichier);
	strcat(nomfichierter,".dat");
	lit_fichier(nomfichierter,&n,vectx,vecty);
	printf("Fichier de %d donnees.\n",n);

	strcpy(nomfichierbis, repres);
	strcat(nomfichierbis, "est%d.");
	strcat(nomfichierbis,nomfichier);
	sprintf(nomfichier,nomfichierbis,methode);
	
	strcpy(nomfichier2, nomfichier);
	strcat(nomfichier2,".CV");

	strcpy(nomfichiergen, nomfichier);
	strcat(nomfichiergen,".Qn%d");
	sprintf(nomfichier3,nomfichiergen,alpha);
	sprintf(nomfichier4,nomfichiergen,100-alpha);

	strcpy(nomfichier5, nomfichier);
	strcat(nomfichier5,".rho");

	strcat(nomfichiergen,".lis");
	sprintf(nomfichier3lis,nomfichiergen,alpha);
	sprintf(nomfichier4lis,nomfichiergen,100-alpha);

	tri_inser(vecty,vectx,n);

	xmin=vectx[1];
	xmax=vectx[1];

	for (i=1;i<=n;i++)
	{
		if (vectx[i]<=xmin) xmin=vectx[i];
		if (vectx[i]>=xmax) xmax=vectx[i];
	}

	printf("Abscisses dans intervalle [%f,%f].\n",xmin,xmax);
	deltax=xmax-xmin;

	ymin=vecty[1];
	ymax=vecty[n];

	printf("Ordonnees dans intervalle [%f,%f].\n",ymin,ymax);

	ordre1=(float)(alpha)/100.;
	ordre2=1-ordre1;

	printf("======= Utilisation de la methode %d ===========\n",methode);

	if (nbh!=0)
	{
		hmin=deltax/n;
	        hmax=deltax/4;
		for (i=1;i<=nbh;i++)
			vecth[i]=hmin+(i-1)*(hmax-hmin)/(nbh-1);

		printf("ETAPE 1: Recherche du h optimal ...\n");

		if (methode==1)
			hopt= Recherche_hopt_1(nbh,n,vectx,vecty,vecth,critere);
		else 
		{
			hmean=Recherche_hmean(nbh,n,vectx,vecty,vecth,critere);
			u01=Phi01_1(ordre1);
			temp=fi01(u01);
			coeff=pow((double)(ordre1*ordre2/(temp*temp)),0.2);
			hopt=hmean*coeff;
			utmp=0;
			temp=fi01(utmp);
			coeff=pow((double)(0.25/(temp*temp)),0.2);
			htmp=hmean*coeff;
			h1opt=hopt;
			if (htmp>=1)
				h2opt=htmp*htmp*htmp*htmp/(hopt*hopt*hopt);
			else	
			{
				htmp1=htmp*htmp*htmp*htmp*htmp/(hopt*hopt*hopt);
				htmp2=hopt/10.;
				if (htmp1>htmp2) h2opt=htmp1;
				else		 h2opt=htmp2;
			}
		}

		ecrit_fichier(nomfichier2,nbh,vecth,critere);
	}
	else	printf("ETAPE 1: supprimee, utilisation valeurs utilisateurs ...\n");

	printf("ETAPE 2: Calcul des quantiles conditionnels ...\n");


	if (methode==1)
	{
		printf("Parametre de lissage : h = %f\n",hopt);

  		for (i=1;i<=nbz;i++)
		{
			qalpha[i]  =calcul_quantile_1(ordre1,hopt,n,vectx,vecty,
							vectz[i],nbessais);
			q1malpha[i]=calcul_quantile_1(ordre2,hopt,n,vectx,vecty,
							vectz[i],nbessais);
		}
	}

	if (methode==2)
	{
		printf("Parametre de lissage : h = %f\n",hopt);

		for (i=1;i<=nbessais;i++)
			essx[i]=ymin-0.1+(i-1)*(ymax-ymin+0.2)/(nbessais-1);

  		for (i=1;i<=nbz;i++)
		{
			qalpha[i]  =calcul_quantile_2(ordre1,hopt,n,vectx,vecty,
							vectz[i],nbessais,essx,essy);
			q1malpha[i]=calcul_quantile_2(ordre2,hopt,n,vectx,vecty,
							vectz[i],nbessais,essx,essy);
		}
		ecrit_fichier(nomfichier5,nbessais,essx,essy);
	}

	if (methode==3)
	{
		printf("Parametres de lissage : h1 = %f et h2= %f\n",h1opt,h2opt);

  		for (i=1;i<=nbz;i++)
		{
			qalpha[i]  =calcul_quantile_3(ordre1,h1opt,h2opt,n,vectx,
							vecty, vectz[i],nbessais);
			q1malpha[i]=calcul_quantile_3(ordre2,h1opt,h2opt,n,vectx,
							vecty, vectz[i],nbessais);
		}
	}

	ecrit_fichier(nomfichier3,nbz,vectz,qalpha);
	ecrit_fichier(nomfichier4,nbz,vectz,q1malpha);

	if (lissage==1)
	{
		printf("ETAPE 3: Lissage des quantiles conditionnels ...\n");

		hmin=deltax/nbz;
	        hmax=deltax/4;
		nbh=20;
		for (i=1;i<=nbh;i++)
			vecth[i]=hmin+(i-1)*(hmax-hmin)/(nbh-1);

		hfin1=Recherche_hmean(nbh,nbz,vectz,qalpha,vecth,critere);
		hfin2=Recherche_hmean(nbh,nbz,vectz,q1malpha,vecth,critere);

		printf("Parametres de lissage des resultats: h1 = %f et h2= %f\n",hfin1,hfin2);
	
		for (i=1;i<=nbz;i++)
		{
			qalphalis[i]=Estimateur_Regression(vectz[i],vectz,qalpha,hfin1,nbz,0);
			q1malphalis[i]=Estimateur_Regression(vectz[i],vectz,q1malpha,hfin2,nbz,0);
		}
		ecrit_fichier(nomfichier3lis,nbz,vectz,qalphalis);
		ecrit_fichier(nomfichier4lis,nbz,vectz,q1malphalis);
	}
}