/*******************************************************************/
/***************  NEURAL FUZZY NETWORK IDENTIFIER  *****************/
/*** GOAL:(1) on-line initial parameter choosing                ****/
/***      (2) parameter learning with  BP algorithm             ****/
/***      (3) add new rule ,if the current rules are lack       ****/
/*** FILE: FNNIP2.C (the dynamic system separate to an unknown  ****/
/**        unlinear part and a known linear part)               ****/
/*** This is example 2 -- a SISO system                         ****/
/********************************************  C. T. Chao 1993/3/8 */

#include   <stdio.h>
#include   <string.h>
#include   <stdlib.h>
#include   <alloc.h>
#include   <math.h>
#define TRUE 1
#define FALSE 0
#define maxseed 32767.0
#define pi 3.14159

FILE *fp1,*fp2,*fp3,*fp4,*fp5;
char initial[12],final[12],ferror[12],fout[12],foutp[12];
/**********************************************************************
final: SAVE (1). STRUCTRUE DEFINITION (2). MEAN & var (3). w34
FMFILE : SAVE MEANERRR
***********************************************************************/
/********** MEMBERSHIP FUNCTIONS **************/
/* gaussaian :  membership gaussian function  */
/**********************************************/
float gaussian(float,float,float);

main()
{

    float **mean,**var,*w34;
    float *inv,desireout;
    float output,uk;
    float **o2,*o3;
    float **dmean,**dvar,*dw34;
    float **delta2,*delta3,delta4;
    unsigned int sizeofin=2,rulenum;
    /****** sizeofin is the state number of non-linear function ********/

    float *meanerror;
    float errnow=0.0;
    long float den,index;

 /*******************************************************
     maxpass  : MAX ITERATION NUMBER
     pass     : PRESENT ITERATION NUMBER
     alphaw   : LEARNING RATE OF w34
     gammaw   : MOMENTUM CONSTANCE OF w34
     alpham   : LEARNING RATE OF MEAN AND var
     gammam   : MOMENTUM CONSTANCE OF MEAN AND var
 ********************************************************/

  int maxpass,pass=0,ptspan,trainning=TRUE;
  float alphaw,gammaw,alpham,gammam,et;

  /* temporary parameters */

  int i,j,k,L,l;
  float tempf,maxtempf,mintempf;
  float meanst,meannd,hold,overlapfactor;
  int maxrule;

  clrscr();

printf("\nINPUT THE OUTPUT FILE NAME :(.m)" );
scanf("%s",&final);
strcpy(initial,final);
strcpy(ferror,final);
strcpy(fout,final);
strcpy(foutp,final);
/*** add the extension for each file ***/
strcat(final,"f.m");
strcat(initial,"i.m");
strcat(ferror,"e.m");
strcat(fout,"o.m");
strcat(foutp,"p.m");

   printf("\n ################  THE STRUCTURE DEFINITION  #################");
/*   printf("\nHow many inputs are be fed into FUZZY IDENTIFIER? ");
     scanf("%d",&sizeofin);  */
   printf("\nHow many rules do you want for FUZZY IDENTIFIER? ");
   scanf("%d",&rulenum);
   printf("\nWhat is the overlap factor do you want? ");
   scanf("%f",&overlapfactor);

   printf("\n ###############  THE LEARNING ENIROMENT   ###################");
   printf("\nHow many iterations are needed in parameter learning: ");
   scanf("%d",&maxpass);
   printf("\nWhat is the error tolence? ");
   scanf("%f",&et);
   printf("\nWhat is the learning rate alphaw (0 < alpha < 1) : ");
   scanf("%f",&alphaw);
   printf("\nWhat is the momentum rate gammaw (0 < gamma < 1) : ");
   scanf("%f",&gammaw);
   printf("\nWhat is the learning rate alpham (0 < alpha < 1) : ");
   scanf("%f",&alpham);
   printf("\nWhat is the momentum rate gammam (0 < gamma < 1) : ");
   scanf("%f",&gammam);
   printf("\nDisplay data every ? iterations (say 50) : ");
   scanf("%d",&ptspan);

maxrule=40;

/*********************************************/
/* MEMORY ALLOCATION DECLARATION */
/*********************************************/
/*********************************************/
  /***** LAYER 1 *****/
  inv = (float *)calloc( sizeofin ,sizeof(float) );

  /***** LAYER 2 *****/
  o2 = (float **)calloc( maxrule , sizeof(float *) );
  for(i=0; i<maxrule ; i++)
     *(o2+i) = (float *)calloc( sizeofin , sizeof(float) );

  mean = (float **)calloc( maxrule , sizeof(float *) );
  dmean = (float **)calloc(maxrule , sizeof(float *) );
  delta2 = (float **)calloc(maxrule , sizeof(float *) );
  for(i=0; i<maxrule; i++)
  {
     *(mean+i)   = (float *)calloc( sizeofin , sizeof(float) );
     *(dmean+i)  = (float *)calloc( sizeofin , sizeof(float) );
     *(delta2+i) = (float *)calloc( sizeofin , sizeof(float) );
  }

  var  = (float **)calloc( maxrule , sizeof(float *) );
  dvar = (float **)calloc( maxrule , sizeof(float *) );
  for(i=0; i<maxrule; i++)
  {
     *(var+i)  = (float *)calloc( sizeofin , sizeof(float) );
     *(dvar+i) = (float *)calloc( sizeofin, sizeof(float) );
  }

  /***** LAYER 3 *****/
  o3 = (float *)calloc( maxrule , sizeof(float) );
  delta3 = (float *)calloc( maxrule , sizeof(float) );

  /***** Layer 4 *****/
  w34  = (float *)calloc( maxrule , sizeof(float) );
  dw34 = (float *)calloc( maxrule , sizeof(float) );


  meanerror = (float *)calloc( maxpass , sizeof(float) );
  /**********************************************************************/
  /*************              INITIAL  PHASE              ***************/
  /**********************************************************************/
  printf("\a%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
  printf("\n       THE INITIAL PARAMETER CHOOSING PHASE ....................");

  /************** DYNAMIC NONLINEAR PLANT *********************/
  /****   Y(k+1)=f( Y(k); Y(k-1) ) +U(k) ****/
  /****   Y(k+1)=f( x1,x2) +U(k)         ****/
  /*************** INITIAL CONTIDITION *************/
  /************ initialize input vector ************/

  for(i=0;i<sizeofin;i++)
  inv[i]=0.0;

  randomize();
  L=0;
  while( L <rulenum ) /** RULE SELLECTION **/
  {
      uk=4*(float)rand()/maxseed -2.0;
      if(uk>1.99)  uk=2.0;
      if(uk<-1.99) uk=-2.0;
      tempf=1.0 + pow(inv[0],2)+pow(inv[1],2) ;
      desireout=inv[0]*inv[1]*( inv[0] + 2.5 )/tempf +uk;
  /********* find the mean and variance for this data point **********/
      w34[L]=desireout;

      for(i=0;i<sizeofin;i++)
          mean[L][i]=inv[i];

      inv[1]=inv[0];
      inv[0]=desireout;

     L++;
  }  /****************** end initial rule sellection ***********************/
      printf("\nTHE RULE NUMBER IS %d.",L);

    /********* NEW METHOD FOR CALCULATE VARIANCE *************/
     for(k=0;k<sizeofin;k++)
     {
     for(i=0;i<rulenum;i++)
     {  meanst=10000;/** distance **/

        for(j=0;j<rulenum;j++)
        {
           if(mean[j][k]!=mean[i][k])
           {
             hold=fabs(mean[i][k]-mean[j][k]);
             if(hold<meanst)
             {  meanst=hold; l=j; }    }/*end if */    } /* end j-loop */


         if(mean[l][k]> mean[i][k])/*** find the second nearest meannd ***/
         {  meannd=10000;
            for(j=0;j<rulenum;j++)
            {
               if(mean[j][k]<mean[i][k])
               {  hold=fabs(mean[i][k]-mean[j][k]);
                  if(hold<meannd)
                    meannd=hold;   }         }/* end j-loop */

               if(meannd==10000)
                   meannd=meanst;            } /** end if */
          /*** find the second nearest meannd ***/
          else if(mean[l][k]< mean[i][k])
          {  meannd=10000;
             for(j=0;j<rulenum;j++)
             {
                if(mean[j][k]>mean[i][k])
                {  hold=fabs(mean[i][k]-mean[j][k]);
                   if(hold<meannd)
                     meannd=hold;  }          }  /* end j-loop */

                if(meannd==10000)
                   meannd=meanst;            } /** end if */

        var[i][k]=meannd/overlapfactor;
            }/**end i-loop FOR MEMBERSHIP */
              } /*** end k-loop ,FOR EACH INPUT*/

  /************************************************/
  /** WRITE THE INITIAL WEIGHTS ONTO ini.m FILE ***/
  /************************************************/

       if((fp1 = fopen(initial,"wt")) == NULL)
       {
        printf("\aERROR : Cannot open file %s",initial);
        exit(0);
       }

        fprintf(fp1,"\%%s : INITIAL PARAMETER",initial);
        fprintf(fp1,"\nrulenum=%d",rulenum);
        fprintf(fp1,"\nsizeofin=%d",sizeofin);
        fprintf(fp1,"\nRULE=[\n");
        for(i=0; i<rulenum; i++)
        {
            for(j=0; j<sizeofin; j++)
                fprintf(fp1,"%10.4f%10.4f",mean[i][j],var[i][j]);
            fprintf(fp1,"%10.4f",w34[i]);
            fprintf(fp1,"\n");       }

            fprintf(fp1,"];\n\n");

     fprintf(fp1,"\nalphaw34=%1.2f alpham=%1.2f",alphaw,alpham);
     fprintf(fp1,"\n gammaw34=%1.2f gammam=%1.2f",gammaw,gammam);

     fclose(fp1);


  /*** INITIAL DWEIGHT **/
  for(i=0;i<maxrule;i++)
     dw34[i]=0.0;


  for(i=0; i<maxrule ;i++)
  for(j=0; j<sizeofin;j++)
  {   dmean[i][j] = 0.0;
      dvar[i][j] = 0.0;   }

 /*****************************************************************/
 /*******************  PARAMETER LEARNING PHASE  ******************/
 /*****************************************************************/
L=0;
while( (pass < maxpass) && (trainning) )  /** PARAMETER LEARNING **/
{
  /*** INPUT OF DANAMIC SYSTEM **/
  uk=4*(float)rand()/maxseed -2.0;
  if(uk>1.99)
     uk=2.0;
  if(uk<-1.99)
     uk=-2.0;

  /*** THE OUTPUT OF DANAMIC SYSTEM **/
   tempf=1.0 + pow(inv[0],2)+pow(inv[1],2) ;
   desireout=inv[0]*inv[1]*( inv[0] + 2.5 )/tempf +uk;

  /******** FIND THE BOUND FOR MEMBERSHIP FUNCTION *********/

       for(k=0;k<sizeofin;k++)
       {
             maxtempf=-10000;
             mintempf=10000;

            for(i=0;i<rulenum;i++)
            {
                 if(mean[i][k]>maxtempf)
                    maxtempf=mean[i][k];
                 if(mean[i][k]<mintempf)
                    mintempf=mean[i][k];   }

            if(inv[k]<mintempf) inv[k]=mintempf;
            else if(inv[k]>maxtempf) inv[k]=maxtempf;    }

/**************************************************************/
/*** PROPAGATE FOWARD THROUGH SUCCEDING LAYERS            *****/
/***  AND CALCULATE OUTPUTS OF EACH LAYER      ****************/
/**************************************************************/

    /* LAYER 2 */

    for(i=0; i<rulenum; i++)
    for(j=0; j<sizeofin;j++)
        o2[i][j] = gaussian(mean[i][j],var[i][j],inv[j]);

  /* LAYER 3 */

  den=0.0;
  for(i=0; i<rulenum; i++)
  {   tempf = 1;
      for(j=0; j<sizeofin ; j++)
         tempf *= o2[i][j];
      o3[i] = tempf;
      den+=o3[i];               }

  /************************************************************************/
  /*************     Add new rule ,if sigma(o3)=0.0    ********************/
  /************************************************************************/
  index=pow(0.01,sizeofin);
  if( (den<index)&&(rulenum<maxrule) )
    {
      rulenum=rulenum+1;
      for(i=0;i<sizeofin;i++)
         mean[rulenum-1][i]=inv[i];
      w34[rulenum-1]=desireout;

     /********* NEW METHOD FOR CALCULATE VARIANCE *************/
     for(k=0;k<sizeofin;k++)
     {
       meanst=10000;/** distance **/

        for(j=0;j<rulenum;j++)
        {
           if(mean[j][k]!=mean[rulenum-1][k])
           {
             hold=fabs(mean[rulenum-1][k]-mean[j][k]);
             if(hold<meanst)
             {  meanst=hold; l=j; }    }/*end if */    } /* end j-loop */


         if(mean[l][k]> mean[rulenum-1][k])/*** find the second nearest meannd ***/
         {  meannd=10000;
            for(j=0;j<rulenum;j++)
            {
               if(mean[j][k]<mean[rulenum-1][k])
               {  hold=fabs(mean[rulenum-1][k]-mean[j][k]);
                  if(hold<meannd)
                    meannd=hold;   }         }/* end j-loop */

               if(meannd==10000)
                   meannd=meanst;            } /** end if */
          /*** find the second nearest meannd ***/
          else if(mean[l][k]< mean[rulenum-1][k])
          {  meannd=10000;
             for(j=0;j<rulenum;j++)
             {
                if(mean[j][k]>mean[rulenum-1][k])
                {  hold=fabs(mean[rulenum-1][k]-mean[j][k]);
                   if(hold<meannd)
                     meannd=hold;  }          }  /* end j-loop */

                if(meannd==10000)
                   meannd=meanst;            } /** end if */

        var[rulenum-1][k]=meannd/overlapfactor;
              } /*** end k-loop ,FOR EACH INPUT*/


  /**** calculate layer 2 ****/
   for(i=0;i<rulenum;i++)
   for(j=0;j<sizeofin;j++)
       o2[i][j]=gaussian(mean[i][j],var[i][j],inv[j]);

  /**** calculate layer 3 ****/
  den=0.0;
  for(i=0;i<rulenum;i++)
  {  tempf=1.0;
     for(j=0; j<sizeofin ; j++)
        tempf *= o2[i][j];
     o3[i]=tempf;
   den+=o3[i];    }

   printf("\a Add one rule ,the current rule number is %d",rulenum);
   printf("\n Hit any key to continue.");getche();
   } /*** END IF(add rule) ***/

/************************ continue BP algorithm *****************************/
  /* LAYER 4 */

   tempf = 0.0;
   for(i=0; i<rulenum; i++)
      tempf += o3[i]*w34[i];

   /*** THE OUTPUT OF FNN IDENTIFIER **/
   /*** Linear part + Non-linear part ***/
   output = tempf/den +uk;

   delta4 = (desireout-output)/den;

  /* COMPUTE THE DELTAS OF THE OUTPUT LAYER */

  for(i=0; i<rulenum; i++)
     delta3[i]=delta4*( w34[i]-output);


  for(i=0; i<rulenum; i++)
  for(j=0; j<sizeofin; j++)
       delta2[i][j]=delta3[i]*o3[i];

  /* UPDATE WEIGHTS */
  for(i=0; i<rulenum; i++)
    {  dw34[i] = alphaw * delta4 * o3[i] + gammaw * dw34[i];
       w34[i] += dw34[i];  }


    for(i=0; i<rulenum; i++)
    for(j=0; j<sizeofin; j++)
     {
      tempf = 2 * (inv[j]-mean[i][j]) / pow(var[i][j],2.0);
      dmean[i][j] = alpham * delta2[i][j] * tempf+gammam * dmean[i][j];
      mean[i][j] += dmean[i][j];
      tempf = 2*pow(inv[j]-mean[i][j],2.0)/pow(var[i][j],3.0);
      dvar[i][j] = alpham * delta2[i][j] * tempf + gammam * dvar[i][j];
      var[i][j] += dvar[i][j];
     }

  /***** SQUARE ERROR MEASURE *****/

     errnow +=( desireout-output )*( desireout-output );
  /** Display results during iteration */
  if ((pass % ptspan)==0)
  {
   for(i=0; i<sizeofin; i++)
       printf("\nin[%d] = %f ",i+1,inv[i]);
   printf("\n");
   printf("out = %1.5f\tdesire = %1.5f\n",output,desireout);

   }

  /* ITERATION PASS UPDATE */

  if(L==200)
  {
      /* MEAN-SQUARE CRITERION (CRITERION = 2) */
   meanerror[pass] = errnow/200;

   if((meanerror[pass] < et) & ((meanerror[pass] * et) >= 0))
      trainning = FALSE;
   printf("\nIteration = %d\tError = %f\n",pass+1,meanerror[pass]);
   pass++;
   L=0;
   errnow = 0.0;
  }
  L++;

     /**** INPUT VECTOR SHIFP ****/
      inv[1]=inv[0];
      inv[0]=desireout;

   } /* end while ( LEARNING PHASE )*/

/******************************************************************/
  /***** WRITE RESULTS INTO DATA FILE ************/

    if((fp2= fopen(final,"wt")) == NULL)
    {
     printf("\aERROR : Cannot open file %s",final);
     exit(0);
     }

     fprintf(fp2,"%d\n",sizeofin);
     fprintf(fp2,"%d\n",rulenum);

     for(i=0;i<rulenum;i++)
     {  for(j=0;j<sizeofin;j++)
           fprintf(fp2,"%10.4f %10.4f",mean[i][j],var[i][j]);
        fprintf(fp2,"%10.4f",w34[i]);
        fprintf(fp2,"\n");                                   }

fclose(fp2);
       /*************************************************/
       /** WRITE THE LEARNING ERROR  ONTO MATLAB FILE ***/
       /*************************************************/

    if((fp3= fopen(ferror,"wt")) == NULL)
    {
     printf("\aERROR : Cannot open file %s",ferror);
     exit(0);
     }
     fprintf(fp3,"\nmeanerr=[\n");
     for(i=0; i<pass; i++)
        fprintf(fp3,"%f\n",meanerror[i]);
     fprintf(fp3,"];");

fclose(fp3);
/*****************************************************************/
/***********  TEST THE MODELED IDENTIFIER  RESPONSE **************/
/***********                   use paralle-series model  *********/
/*****************************************************************/

for(i=0;i<sizeofin;i++)
inv[i]=0.0;
desireout=0.0;

if( (fp4= fopen(fout,"wt") ) == NULL)
{
 printf("\aERROR : Cannot open file %s",fout);
 exit(0);
 }

pass=0;
while( (pass < 100) )
{
    inv[1]=inv[0];
    inv[0]=desireout;
   /***** the input of plant *****/
    uk= sin(2*pi*pass/25);
   /*** find the plant output ***/
   tempf=1.0 + pow(inv[0],2)+pow(inv[1],2) ;
   desireout=inv[0]*inv[1]*( inv[0] + 2.5 )/tempf +uk;

   /******** FIND THE BOUND FOR MEMBERSHIP FUNCTION *********

        for(k=0;k<sizeofin;k++)
        {
              maxtempf=-10000;
              mintempf=10000;

             for(i=0;i<rulenum;i++)
             {
                  if(mean[i][k]>maxtempf)
                     maxtempf=mean[i][k];
                  if(mean[i][k]<mintempf)
                     mintempf=mean[i][k];   }

             if(inv[k]<mintempf) inv[k]=mintempf;
             else if(inv[k]>maxtempf) inv[k]=maxtempf;    }

/**************************************************************/
/***      PROPAGATE FOWARD THROUGH SUCCEDING LAYERS         ***/
/***      AND CALCULATE OUTPUTS OF EACH LAYER               ***/
/**************************************************************/

  /* LAYER 2 */

  for(i=0; i<rulenum; i++)
  for(j=0; j<sizeofin;j++)
      o2[i][j] = gaussian(mean[i][j],var[i][j],inv[j]);

  /* LAYER 3 */

  den=0.0;
  for(i=0; i<rulenum; i++)
  {   tempf = 1;
      for(j=0; j<sizeofin ; j++)
         tempf *= o2[i][j];
      o3[i] = tempf;
      den+=o3[i];               }

  /* LAYER 4 */

   tempf = 0.0;

   for(i=0; i<rulenum; i++)
      tempf += o3[i]*w34[i];

   output = tempf/den +uk;

   printf("\noutput(%d)=%f",pass+1,output);

fprintf(fp4,"out(%3d)=%10.4f;desire(%3d)=%10.4f;\n",pass+1,output,pass+1,desireout);

   pass++;
                } /* end while */

 fclose(fp4);

/*****************************************************************/
/***********  TEST THE MODELED IDENTIFIER  RESPONSE **************/
/***********                   use paralle model    **************/
/*****************************************************************/

for(i=0;i<sizeofin;i++)
inv[i]=0.0;
output=0.0;

if((fp5= fopen(foutp,"wt")) == NULL)
{
 printf("\aERROR : Cannot open file %s",foutp);
 exit(0);
 }

pass=0;
while( (pass < 100) )
{
   inv[1]=inv[0];
   inv[0]=output;

   /***** the input of plant *****/

    uk= sin(2*pi*pass/25);
   /******** FIND THE BOUND FOR MEMBERSHIP FUNCTION *********

        for(k=0;k<sizeofin;k++)
        {
              maxtempf=-10000;
              mintempf=10000;

             for(i=0;i<rulenum;i++)
             {
                  if(mean[i][k]>maxtempf)
                     maxtempf=mean[i][k];
                  if(mean[i][k]<mintempf)
                     mintempf=mean[i][k];   }

             if(inv[k]<mintempf) inv[k]=mintempf;
             else if(inv[k]>maxtempf) inv[k]=maxtempf;    }

/**************************************************************/
/***      PROPAGATE FOWARD THROUGH SUCCEDING LAYERS         ***/
/***      AND CALCULATE OUTPUTS OF EACH LAYER               ***/
/**************************************************************/

     /* LAYER 2 */

  for(i=0; i<rulenum; i++)
  for(j=0; j<sizeofin;j++)
      o2[i][j] = gaussian(mean[i][j],var[i][j],inv[j]);

  /* LAYER 3 */

  den=0.0;
  for(i=0; i<rulenum; i++)
  {   tempf = 1;
      for(j=0; j<sizeofin ; j++)
         tempf *= o2[i][j];
      o3[i] = tempf;
      den+=o3[i];               }

  /* LAYER 4 */

   tempf = 0.0;

   for(i=0; i<rulenum; i++)
      tempf += o3[i]*w34[i];

   output = tempf/den +uk;

   printf("\noutput(%d)=%f",pass+1,output);
   fprintf(fp5,"outp(%3d)=%10.4f;\n",pass+1,output);
   pass++;
          } /* end while */

   fclose(fp5);
/* Outfile message */
printf("\n%s , %s are generated. Check it !!",final,ferror);
printf("\n%s , %s are generated. Check it !!",fout,foutp);

getche();
} /* END MAIN */

float gaussian(float m, float v, float x)
{
return(exp(-pow(x-m,2.0)/pow(v,2.0)));
}