/*
 * Claude TOUZET - Mai 2001 - Perceptron multicouche et 
 * reconnaissance de caractères manuscrits.
 */

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

#define SEUIL 0.5
#define N_ENTREE 15
#define N_SORTIE 10

int   N_CACHEE;
int   NBIT;
float LAMDA;
int   nbproto;
float **poids12;
float **poids23;
float **poids13;
float *poids2seuil;
int   *chiffre;
float **x_proto;
float **y_des;
float **y_sortie;
int   phase;
FILE  *fp(), *res;
char  nomfic1[50];
unsigned short xsubi;

double erand48(unsigned short xsubi);
void codage();
void sortie_des(int rangpt,int rang);
void init_poids();
void sauve_poids();
void init_tab();
void lec_poids();
int classe(int rangpt);
void lec_fich(char nomfic[50]); 
void reseau1 (float **poids12,float **poids23,float **poids13,
	float *poids2seuil,float **x_proto,float **y_des,int *iter,float **y_sort );


double        erand48(unsigned short xsubi)
{
  unsigned short toto;
  toto=xsubi;
  return((double) rand()/RAND_MAX);
  toto=toto+1;
}

/********************** Procedure de codage des entrees *******************/

void codage()
{
  int i,j;

  for (i=0;i<nbproto;i++)
  {
    for (j=0;j<N_ENTREE;j++) x_proto[i][j]=(x_proto[i][j]-21.0)/21.0;
  }
}

/******************** Determination de la sortie desiree ****************/

void sortie_des(int rangpt,int rang)
{
  int i;

  for (i=0;i<N_SORTIE;i++) y_des[rangpt][i]=-1;
  y_des[rangpt][rang]=1;

}

/****** Procedure d'initialisation des poids pour l'apprentissage  ******/

void init_poids()
{
  int i,j;
  
  for (i=0;i<N_ENTREE;i++) {
    for (j=0;j<N_CACHEE;j++) poids12[i][j] = erand48(xsubi);
  }
  for (i=0;i<N_CACHEE;i++) {
    for (j=0;j<N_SORTIE;j++) poids23[i][j] =  erand48(xsubi);
  }
  for (i=0;i<N_ENTREE;i++) {
    for (j=0;j<N_SORTIE;j++) poids13[i][j] =  erand48(xsubi);
  }
  for (i=0;i<N_CACHEE;i++) poids2seuil[i] = erand48(xsubi);

}

/*************** Procedure de sauvegarde des poids *****************/

void sauve_poids()
{
  FILE *fp;
  int i,j,reponse;
  char nomfic[50];

  printf("\nVoulez-vous sauvegarder les poids, Non : 0 ; Oui : 1 ? ");
  scanf("%d", &reponse);
  if (reponse==1) {
    printf("\nNom du fichier des poids ? ");
    scanf("%s", nomfic);
    fp=fopen(nomfic,"w");
    fprintf(fp, "%d ", N_CACHEE);
    for (i=0;i<N_CACHEE;i++) fprintf(fp,"%14.12f ",poids2seuil[i]);
    for (i=0;i<N_ENTREE;i++) {
      for (j=0;j<N_CACHEE;j++) fprintf(fp,"%14.12f ",poids12[i][j]);
    }
    for (i=0;i<N_CACHEE;i++) {
      for (j=0;j<N_SORTIE;j++) fprintf(fp,"%14.12f ",poids23[i][j]);
    }
    for (i=0;i<N_ENTREE;i++) {
      for (j=0;j<N_SORTIE;j++) fprintf(fp,"%14.12f ",poids13[i][j]);
    }
    fclose(fp);
  }

}

/************************** Initialisation ******************************/

void init_tab()
{
  int i;

  poids12 = (float**) malloc( N_ENTREE * sizeof( float* ) );
  poids13 = (float**) malloc( N_ENTREE * sizeof( float* ) );
  for (i=0;i<N_ENTREE;i++) {
    *( poids12 + i ) = (float*) malloc( N_CACHEE * sizeof( float ) );
    *( poids13 + i ) = (float*) malloc( N_SORTIE * sizeof( float ) );
  }
  poids23 = (float**) malloc( N_CACHEE * sizeof( float* ) );
  for (i=0;i<N_CACHEE;i++) {
    *( poids23 + i ) = (float*) malloc( N_SORTIE * sizeof( float) );
  }
  poids2seuil = (float*) malloc( N_CACHEE * sizeof(float ) );

}

/********************** Procedure de lecture des poids  ****************/

void lec_poids()
{
  FILE *fp;
  int i,j;
  char nomfic[10];
  
  printf("\nNom du fichier des poids ? ");
  scanf("%s", nomfic);
  /*strcpy(nomfic,"poids");*/
  fp=fopen(nomfic,"r");
  fscanf(fp, "%d", &N_CACHEE);
  init_tab();
  for (i=0;i<N_CACHEE;i++) fscanf(fp,"%f ", &poids2seuil[i]);
  for (i=0;i<N_ENTREE;i++) {
    for (j=0;j<N_CACHEE;j++) fscanf(fp,"%f ", &poids12[i][j]);
  }
  for (i=0;i<N_CACHEE;i++) {
    for (j=0;j<N_SORTIE;j++) fscanf(fp,"%f ", &poids23[i][j]);
  }
  for (i=0;i<N_ENTREE;i++) {
    for (j=0;j<N_SORTIE;j++) fscanf(fp,"%f ", &poids13[i][j]);
  }
  fclose(fp);

}
        

/******************** Procedure de classement  **************/

int classe(int rangpt)
{
  int clas;
  int i,nb1;
  int sign[N_SORTIE];
  
  nb1=0;
  for (i=0;i<N_SORTIE;i++)  sign[i]=(y_sortie[rangpt][i]>=0.5)? 1 : -1;
  for (i=0;i<N_SORTIE;i++) {
    if (sign[i]==1)
    {
      clas=i;
      nb1++;
    }
  }
  if (nb1==1) return(clas);
  else return(-1);
}

/**********  Procedure de lecture des prototypes en fichier  *************/

void lec_fich(char nomfic[50])
{
  FILE  *fp;
  int i,j;

  x_proto = (float**) malloc( nbproto * sizeof( float* ) );
  y_des = (float**) malloc( nbproto * sizeof( float* ) );
  y_sortie = (float**) malloc( nbproto * sizeof( float* ) );
  for (i=0; i< nbproto; i++ )
  {
    *( x_proto + i ) = (float*) malloc( N_ENTREE * sizeof( float ) );
    *( y_des + i ) = (float*) malloc( N_SORTIE * sizeof( float ) );
    *( y_sortie + i ) = (float*) malloc( N_SORTIE * sizeof( float ) );
  }
  chiffre = (int*) malloc( nbproto * sizeof( int ) );
  if ((fp=fopen(nomfic,"r")) == NULL) printf(" Pas touvé le fichier des prototypes ! \n");
  for (i=0;i<nbproto;i++)
  {
    for (j=0;j<N_ENTREE;j++)
    {
    fscanf(fp,"%f", &x_proto[i][j]);
    }
    fscanf(fp,"%d", &chiffre[i]);
  }
  fclose(fp);

}


/********************** reseau *************************/

void reseau1 (float **poids12,float **poids23,float **poids13,
float *poids2seuil,
float **x_proto,
float **y_des,
int *iter,
float **y_sort )
{
  float *y_0;
  float *x_1, *y_1;
  float *x_2, *y_2;
  float *grad1, *grad2;
  int   i, j, nuproto;
  int   iteration;
  int   fin_app;
  float res_inter;
  float som_modif;
  float crit;

  double         erand48();

  printf("reseau \n");
  
  
/******************  Allocation des tableaux *****************/
  y_0   = (float *) malloc( N_ENTREE * sizeof( float ) );
  x_1   = (float *) malloc( N_CACHEE * sizeof( float ) );
  y_1   = (float *) malloc( N_CACHEE * sizeof( float ) );
  x_2   = (float *) malloc( N_SORTIE * sizeof( float ) );
  y_2   = (float *) malloc( N_SORTIE * sizeof( float ) );
  
  if (phase)
  {

/* apprentissage */

    grad1 = (float *) malloc( N_CACHEE * sizeof( float ) );
    grad2 = (float *) malloc( N_SORTIE * sizeof( float ) );

    iteration=0;
    fin_app=0;
    while ((fin_app==0) && (iteration++<NBIT))
    {
      fin_app=1;
      for (nuproto=0;nuproto<nbproto;nuproto++)
      {

/************** calcul de la sortie **********************/

	for (i=0;i<N_ENTREE;i++)
	{
	  res_inter=exp((double) x_proto[nuproto][i]);
	  y_0[i]=(res_inter-1)/(res_inter+1);
	}
	for (i=0;i<N_CACHEE;i++)
	{
	  x_1[i] = -poids2seuil[i];
	  for (j=0;j<N_ENTREE;j++)
	  {
	  x_1[i]+=poids12[j][i]*y_0[j];
	  }
	  res_inter=exp((double) x_1[i]);
	  y_1[i]=(res_inter-1)/(res_inter+1);
	}
	for (i=0;i<N_SORTIE;i++)
	{
	  x_2[i]=0.0;
	  for (j=0;j<N_CACHEE;j++)
	  {
	  x_2[i]+=poids23[j][i]*y_1[j];
	  }
	  for (j=0;j<N_ENTREE;j++)
	  {
	  x_2[i]+=poids13[j][i]*y_0[j];
	  }
	  res_inter=exp((double) x_2[i]);
	  y_2[i]=(res_inter-1)/(res_inter+1);
	}

/*************** modification des poids ****************/

	/* calcul du gradient sur la couche de sortie   */
	for (i=0;i<N_SORTIE;i++)
	{
	grad2[i]= 2.0*((1-y_2[i]*y_2[i])/2.0)*(y_des[nuproto][i]-y_2[i]);
	}
	/* calcul du gradient sur la couche cachee     */
	for (i=0;i<N_CACHEE;i++)
	{
	  som_modif=0.0;
	  for (j=0;j<N_SORTIE;j++)
	  {
	  som_modif+=poids23[i][j]*grad2[j];
	  }
	  grad1[i]=som_modif*(1-y_1[i]*y_1[i])/2.0;
	}

	/* modification des poids synaptiques                                                */
	for (i=0;i<N_ENTREE;i++)

	{
	  for (j=0;j<N_CACHEE;j++) poids12[i][j]+=LAMDA*grad1[j]*y_0[i];
	} 
	for (i=0;i<N_CACHEE;i++)
	{
	  for (j=0;j<N_SORTIE;j++)
	  {
	  poids23[i][j]+=LAMDA*grad2[j]*y_1[i];
	  }
	}
	for (i=0;i<N_ENTREE;i++)
	{
	  for (j=0;j<N_SORTIE;j++)
	  {
	  poids13[i][j]+=LAMDA*grad2[j]*y_0[i];
	  }
	}

	/* modification des poids du seuil*/
	for (j=0;j<N_CACHEE;j++)
	{
	poids2seuil[j]+=LAMDA*grad1[j];
	}
      }
      for (nuproto=0;nuproto<nbproto;nuproto++)
      {

	for (i=0;i<N_ENTREE;i++)
	{
	  res_inter=exp((double) x_proto[nuproto][i]);
	  y_0[i]=(res_inter-1)/(res_inter+1);
	}
	for (i=0;i<N_CACHEE;i++)
	{
	  x_1[i]=-poids2seuil[i];
	  for (j=0;j<N_ENTREE;j++)
	  {
	  x_1[i]+=poids12[j][i]*y_0[j];
	  }
	  res_inter=exp((double) x_1[i]);
	  y_1[i]=(res_inter-1)/(res_inter+1);
	}
	for (i=0;i<N_SORTIE;i++)
	{
	  x_2[i]=0.0;
	  for (j=0;j<N_CACHEE;j++)
	  {
	  x_2[i]+=poids23[j][i]*y_1[j];
	  }
	  for (j=0;j<N_ENTREE;j++)
	  {
	  x_2[i]+=poids13[j][i]*y_0[j];
	  }
	  res_inter=exp((double) x_2[i]);
	  y_2[i]=(res_inter-1)/(res_inter+1);
	}

/**************** calcul de l'erreur *****************/

	crit=0.0;
	for (i=0;i<N_SORTIE;i++)
	{
	crit+=(y_des[nuproto][i]-y_2[i])*(y_des[nuproto][i]-y_2[i]);
	}
	if (SEUIL<crit) fin_app=0;

      }
      iter[0]=iteration;
      printf("\n%d", iteration);
    }
    free(grad1);
    free(grad2);
  }

  else

  {
/* reconnaissance : calcul de la sortie*/

    for (nuproto=0;nuproto<nbproto;nuproto++)
    {
      fin_app=1;

      for (i=0;i<N_ENTREE;i++)
      {
	res_inter=exp((double) x_proto[nuproto][i]);
	y_0[i]=(res_inter-1)/(res_inter+1);
      }
      for (i=0;i<N_CACHEE;i++)
      {
	x_1[i]=-poids2seuil[i];
	for (j=0;j<N_ENTREE;j++)
	{
	x_1[i]+=poids12[j][i]*y_0[j];
	}
	res_inter=exp((double) x_1[i]);
	y_1[i]=(res_inter-1)/(res_inter+1);
      }
      for (i=0;i<N_SORTIE;i++)
      {
	x_2[i]=0.0;
	for (j=0;j<N_CACHEE;j++) x_2[i]+=poids23[j][i]*y_1[j];
	for (j=0;j<N_ENTREE;j++) x_2[i]+=poids13[j][i]*y_0[j];
	res_inter=exp((double) x_2[i]);
	y_2[i]=(res_inter-1)/(res_inter+1);
      }

      for (i=0;i<N_SORTIE;i++)
      {
      y_sort[nuproto][i]=y_2[i];
      }
    }
  }
  free(y_0);
  free(x_1);
  free(y_1);
  free(x_2);
  free(y_2);

}


/**************** main *****************/

int main(void)
{
  void  locate();
  int   nuproto;
  int   iteration;
  int   resultat;
  float pourcent;

  printf("\nPhase d'apprentissage (1) ou de reconnaissance (0) ou fin (2) : ");
  scanf("%d", &phase);
while (phase != 2)
{
  if (phase == 1) {
    printf("\nDimension de la couche cachee : ");
    scanf("%d", &N_CACHEE);
    init_tab();
    init_poids();
    printf("\nNom du fichier des prototypes : ");
    scanf("%s", nomfic1);
    printf("\nNombre de prototypes dans le fichier ( < 190 ) : ");
    scanf("%d", &nbproto);
    lec_fich(nomfic1);
    printf("\nEntrez la valeur de lamda : ");
    scanf("%f", &LAMDA);
    printf("\nNombre d'iterations maximum : ");
    scanf("%d", &NBIT);
    for (nuproto=0;nuproto<nbproto;nuproto++)
    {
    sortie_des(nuproto,chiffre[nuproto]);
    }
  }
  else {
    lec_poids();
    printf("\nNom du fichier des prototypes : ");
    scanf("%s", nomfic1);
    printf("\nNombre de prototypes dans le fichier ( < 190 ) : ");
    scanf("%d", &nbproto);
    lec_fich(nomfic1);
  }
  codage();
  reseau1( poids12, poids23, poids13,poids2seuil, x_proto, y_des, &iteration, y_sortie);

  if (phase == 1)
  {
    printf("\n\n%d iterations\n", iteration);
    sauve_poids();
  }
  else
  {
    resultat=0;
    for (nuproto=0;nuproto<nbproto;nuproto++)
    {
      printf("Le prototype est un %d, sa classe est %d\n", chiffre[nuproto], classe(nuproto));
      if (classe(nuproto)==chiffre[nuproto]) resultat++;
    }
    pourcent = 100.0* (float) resultat/nbproto;
    printf("\nTaux de reconnaissance : %f\n", pourcent);
  }
  printf("\n");
  printf("\nPhase d'apprentissage (1) ou de reconnaissance (0) ou fin (2) : ");
  scanf("%d", &phase);

}
	return 0;
}