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


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


#define MAXDIM   11
#define MAXN     1000
#define MAXNBZ   1000
#define MAXNBH   100

int	n,nbh,nbz,nbx,i,j,indh,indx,indz,indy,dim;
float	vectx[MAXN][MAXDIM],vecty[MAXN],vecth[MAXNBH],vectz[MAXNBZ][MAXDIM],critere[MAXNBH];
float	essx[MAXNBH],essy[MAXNBH];
float	qalpha[MAXN],q1malpha[MAXN],qalphalis[MAXN],q1malphalis[MAXN];
float 	deltax,hopt,xmin[MAXDIM],xmax[MAXDIM],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,deltatmp;
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],nomfichier4lis[255];
float	fnitest[MAXNBZ],xtest,htest;


/**************************** Tri par insertion  *****************/
void tri_inser_multi(k,dim)
int k;

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

				aux=vectx[j][jj];
				vectx[j][jj]=vectx[indice][jj];
				vectx[indice][jj]=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 %s",temp,repdon);
fscanf(ff,"%s %s",temp,repres);
fscanf(ff,"%s %s",temp,nomfichier);
fscanf(ff,"%s %d",temp,&dim);
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)
{ 
	 fscanf(ff,"%s %f",temp,&hopt);
}	
fclose(ff);


if (dim>=MAXDIM)
{
        printf("Dimension trop grande.\n");
        exit(1);
}


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

if (nbessais>=MAXNBH)
{
        printf("Nb essais trop grand.\n");
        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);
}

/********************* Ecriture fichier de (dim+1) colonnes */

ecrit_multifichier(nom,n,dim,quant)

	char *nom;
	float *quant;
	int n,dim;

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

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

for (i=1;i<=n;i++)
{
	for (j=1;j<=dim;j++)
	{
       		fprintf(f,"%f ",vectz[i][j]);
	}

    	fprintf(f,"%f\n",quant[i]);
}
fclose(f);
}

/********************* Lecture fichier de (dim+1) colonnes */

lit_multifichier(nom,dim,n)

	int  dim;
	char *nom;
	int *n;

{
int i,j,nbpts,ntmp;
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",&xtmp)!=EOF)
        nbpts++;

fclose(f);

if (nbpts==0)
{
        printf("fichier vide!\n");
        exit(1);
}
ntmp=nbpts/(dim+1);

if (ntmp>=MAXN)
{
        printf("nb maximum de pts autorise %d (ici %d)\n",MAXN,ntmp);
        exit(1);
}

f=fopen(nom,"r");
for (i=1;i<=ntmp;i++)
{
	for (j=1;j<=dim;j++)
	{
		fscanf(f,"%f",&xtmp);
		vectx[i][j]=xtmp;
	}
	fscanf(f,"%f",&ytmp);
	vecty[i]=ytmp;
}

fclose(f);

*n=ntmp;
}

/********************* Lecture fichier de dim colonnes */

lit_multifichier_bis(nom,dim,n)

	int  dim;
	char *nom;
	int *n;

{
int i,j,nbpts,ntmp;
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);
}
ntmp=nbpts/dim;

if (ntmp>=MAXNBZ)
{
        printf("nb maximum de pts autorise %d (ici %d)\n",MAXNBZ,ntmp);
        exit(1);
}

f=fopen(nom,"r");
for (i=1;i<=ntmp;i++)
{
	for (j=1;j<=dim;j++)
	{
		fscanf(f,"%f",&xtmp);
		vectz[i][j]=xtmp;
	}
}

fclose(f);
*n=ntmp;
}



/******************** 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);
}



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

float	K(x,dim)

	float *x;
	int dim;

{
	int	i;
	float	tmp;

	tmp=1;
	for (i=1;i<=dim;i++)
		tmp=tmp*fi01(x[i]);

	return(tmp);
}


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

float	EstimateurFn_1(y,x,hn,nbpts,dim,enleve)

	float	*x,y,hn;
	int	nbpts,enleve,dim;

{
	int	i,j;
	float	noyau,res;
	double	numer,denom;
	float 	xc[MAXDIM];

	numer=0;
	denom=0;

	for (i=1;i<enleve;i++)
	{
		for (j=1;j<=dim;j++)
		{
			xc[j]=(x[j]-vectx[i][j])/hn;
		}

		noyau=K(xc,dim);	
		if (vecty[i]<=y) numer+=noyau;
		denom+=noyau;
	}
	for (i=enleve+1;i<=nbpts;i++)
	{
		for (j=1;j<=dim;j++)
		{
			xc[j]=(x[j]-vectx[i][j])/hn;
		}

		noyau=K(xc,dim);	
		if (vecty[i]<=y) numer+=noyau;
		denom+=noyau;
	}

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

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

float calcul_quantile_1(ordre,h,n,dim,xcalcul,nbit)

	float	ordre,h,*xcalcul;
	int	n,nbit,dim;


{
	float	milint,minint,maxint;
	int	it;

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

for (it=1;it<=nbit;it++)
{
	if (EstimateurFn_1(milint,xcalcul,h,n,dim,0)<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,dim,vecth,critere)

	int	nbh,n,dim;
	float	*critere,*vecth;

{
	int	i,indh,indx,indy;
	float	ecart,criteremin,hopt;
	float	projx[MAXDIM];
	

        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("   Point No %d / %d\r",indx,n);
			fflush(stdout);
                        for (indy=1;indy<=n;indy++)
                        {
				for (j=1;j<=dim;j++)
				{
					projx[j]=vectx[indx][j];
				}
					
                                fni=EstimateurFn_1(vecty[indy],projx,vecth[indh],n,dim,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);
}


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


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 du fichier de parametres %s\n",nomparam);
        Lit_Fichier_Parametres(nomparam);
	printf("Covariable de dimension %d.\n",dim);

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

	strcpy(nomgrillebis, repdon);
	strcat(nomgrillebis, nomgrille);
	lit_multifichier_bis(nomgrillebis,dim,&nbz);
	printf("Evaluation du quantile en %d points.\n",nbz);

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

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

	tri_inser_multi(n,dim);

	deltax=0;
	for (j=1;j<=dim;j++)
	{
		xmin[j]=vectx[1][j];
		xmax[j]=vectx[1][j];

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

		printf("Covariable [%d] : valeurs dans intervalle [%f,%f].\n",j,xmin[j],xmax[j]);
		deltatmp=xmax[j]-xmin[j];
		if (deltatmp>deltax) deltax=deltatmp;
	}

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

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

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

	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");

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

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

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

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

	ecrit_multifichier(nomfichier3,nbz,dim,qalpha);
	ecrit_multifichier(nomfichier4,nbz,dim,q1malpha);

}