/******************************************************************************/
/*                                                                            */
/*                      Fuzzy System                                             */
/*             2001.08.03 by CHEN Ching-Han                                   */
/*                                                                            */
/******************************************************************************/
#include <fstream.h>
#include <math.h>
#include <string>
#include <iomanip.h>
#include "array.h"
#define NB_RESOLUTION 20

class FuzzySet
{

public:
   FuzzySet (){ DomainResolution = 0; FunctionOrder=1;}
   void Initialize ( int ndx_resolu, float min, float max);
   void Initialize ( const FuzzySet& s ) ;
   void DefineMembershipFunction ( int npts , float *xpts , float *ypts ) ;
   void DefineMembershipFunction ( f1D &xpts , f1D &ypts ) ;
   void DefineMembershipFunction ( float x0,float y0,float x1,float y1 ,float x2,float y2) ;
   void DefineMembershipFunction ( float x0 , float y0 , float x1 , float y1 ,
              float x2 , float y2 , float x3 , float y3 ) ;
   ~FuzzySet () ;

 // conjunction and disjunction operation for fuzzy inferencing
   void AND(const FuzzySet& s);
   void OR(const FuzzySet& s);
   void Negate() ;
   void Scale(float fac ) ;
// defuzzification
   float Defuzzification();

// utilities for fuzzy set
   bool invalid(){return !DomainResolution ;}
   bool Zero();
   float membership ( float pt ) ;
	float FastMembership(float pt);
   void Print(ostream &out);
   void SetName(char *s){name = s;}
   void Labeling(int ndx){label = ndx;}
   int Label(){return label;}
   char* Name(){return name;}
	void SetFunctionOrder(int order){FunctionOrder=order;}

// operator overloading
   FuzzySet& operator= ( const FuzzySet& s ) ;

protected:
	float DomainMax;
	float DomainMin;
	int DomainResolution;
   f1D MembershipDegreeArray;
   f1D DomainValueArray;
   char *name;
   int label;
   int FunctionOrder;

private:
	float centroid (float *y, int n) ;
	float interpolation(float x0, float* xin, float* yin, int n, int order);
};


/******************************************************************************/
/*                                                                            */
/*               Linguistic Variable class                                    */
/*             1997.08.03 by CHEN Ching-Han                                   */
/*                                                                            */
/******************************************************************************/

class LinguisticVariable
{
public:
  FuzzySet* FzySet;
  int NbFzyTerm;
  float DomainMax;
  float DomainMin;
  f1D point;
  void Initialize(int nb_terms, float min, float max);
  void Print(ostream &os);
  void SetName(char name[]){Name=name;}
  int FindMaxDegreeFuzzySet(float t);
  float FindFuzzySetValue(int ndx_fzy_set, float t);
  FuzzySet GetFuzzySet(int i);
protected:
  void CreateFuzzySet();
  char *Name;
};

/******************************************************************************/
/*                                                                            */
/*                    Fuzzy Rule Primitive                                    */
/*                 1997.08.03 by CHEN Ching-Han                               */
/*                                                                            */
/******************************************************************************/

class FuzzyRule
{
public:
  FuzzySet* antecedent;
  FuzzySet conclusion;
  FuzzySet InferResult;
  int nb_antecedent;
  void Initialize(int nb);
  FuzzySet& Inferencing(f1D &para);
  float Effectiveness(f1D &para, float solution);
  void Print(ofstream &out);
};
/*-----------------------------------------------------------------------------*/




void FuzzySet::Initialize ( int ndx_resolu, float min , float max )
{
	if(ndx_resolu <1) DomainResolution = 1;else DomainResolution = ndx_resolu;
   DomainMax = max;
   DomainMin = min;
	MembershipDegreeArray.Initialize(DomainResolution);
   DomainValueArray.Initialize(DomainResolution);
   DomainValueArray[0] = DomainMin;
   float scale = (max-min)/DomainResolution;
   for(int i=1;i<DomainResolution;i++)
   	DomainValueArray[i] = DomainValueArray[i-1] + scale;
   for(int i=0;i<DomainResolution;i++)MembershipDegreeArray[i]=0.0;
}

void FuzzySet::DefineMembershipFunction ( int npts , float *xpts , float *ypts )
{
	for(int i=0;i<DomainResolution;i++)
   {
   	if(DomainValueArray[i]==xpts[npts-1])MembershipDegreeArray[i]=ypts[npts-1];
      else if (DomainValueArray[i]==xpts[0])MembershipDegreeArray[i]=ypts[0];
		else if((DomainValueArray[i]<xpts[npts-1])&&(DomainValueArray[i]>xpts[0]))
      {
			MembershipDegreeArray[i] = interpolation(DomainValueArray[i], xpts, ypts, npts, FunctionOrder);
      }
   }
}

void FuzzySet::DefineMembershipFunction (f1D &xpts , f1D &ypts )
{
	DefineMembershipFunction(xpts.nb, xpts.m,ypts.m);
}

void FuzzySet::DefineMembershipFunction (float x0 , float y0 , float x1 , float y1 ,
                     float x2 , float y2 , float x3 , float y3 )
{
   float xx[4], yy[4] ;
   xx[0]=x0;xx[1]=x1;xx[2]=x2;xx[3]=x3;
   yy[0]=y0;yy[1]=y1;yy[2]=y2;yy[3]=y3;
   DefineMembershipFunction(4, xx, yy);
}

void FuzzySet::DefineMembershipFunction(float x0, float y0, float x1, float y1,
               float x2, float y2)
{
   float xx[3], yy[3] ;
   xx[0]=x0;xx[1]=x1;xx[2]=x2;
   yy[0]=y0;yy[1]=y1;yy[2]=y2;
	DefineMembershipFunction(3, xx, yy);
}

void FuzzySet::Initialize ( const FuzzySet& s )
{
	DomainMax = s.DomainMax;
	DomainMin = s.DomainMin;
	DomainResolution = s.DomainResolution;
   MembershipDegreeArray = s.MembershipDegreeArray;
   DomainValueArray = s.DomainValueArray;
   FunctionOrder = s.FunctionOrder;
   Labeling(s.Label());
   SetName(s.Name());
}


FuzzySet::~FuzzySet ()
{
	MembershipDegreeArray.Finish();
	DomainValueArray.Finish();
}

FuzzySet& FuzzySet::operator= ( const FuzzySet& s )
{
	Initialize(s);
   return *this;
}

bool FuzzySet::Zero()
{
  bool veri=true;
  for(int i=0;i<DomainResolution;i++)if(MembershipDegreeArray[i]!=0)return false;
}

void FuzzySet::Print(ostream &out)
{
	out<<"Fuzzy Set "<<name<<endl;
	for(int i=0;i<DomainResolution;i++)
		out<<DomainValueArray[i]<<"\t"<<MembershipDegreeArray[i]<<endl;
}


float FuzzySet::membership ( float pt )
{
	return interpolation(pt, DomainValueArray.m, MembershipDegreeArray.m,
   	DomainValueArray.nb, FunctionOrder);
}

float FuzzySet::FastMembership ( float pt )
{
	for(int i=0;i<DomainResolution-1;i++)
   {
   	if((pt>=DomainValueArray[i])&&(pt<=DomainValueArray[i+1]))
		return (MembershipDegreeArray[i]*(pt-DomainValueArray[i+1])
      	- MembershipDegreeArray[i+1]*(pt-DomainValueArray[i]))
         /(DomainValueArray[i]-DomainValueArray[i+1]);
   }
}

void FuzzySet::Negate ()
{
	for(int i=0;i<DomainResolution;i++)
   	MembershipDegreeArray[i]=1.0-MembershipDegreeArray[i];
}

void FuzzySet::Scale ( float fac )
{
	for(int i=0;i<DomainResolution;i++)MembershipDegreeArray[i]*=fac;
}

void FuzzySet::AND(const FuzzySet& s)
{
	for(int i=0;i<DomainResolution;i++)if(MembershipDegreeArray[i]>s.MembershipDegreeArray[i])
     MembershipDegreeArray[i]=s.MembershipDegreeArray[i];
}

void FuzzySet::OR(const FuzzySet& s)
{
	for(int i=0;i<DomainResolution;i++)if(MembershipDegreeArray[i]<s.MembershipDegreeArray[i])
     MembershipDegreeArray[i]=s.MembershipDegreeArray[i];
}

float FuzzySet::Defuzzification()
{
	float pos = centroid(MembershipDegreeArray.m, MembershipDegreeArray.nb);
   return (DomainMin + pos*(DomainMax-DomainMin)/DomainResolution);
}

float FuzzySet::centroid (float *y, int n)
{
/*
   This routine centroids the input vector of data y[]. The assumption is
   that the true data represents a point process and that y[] is a quantized
   version of the process. It is further assumed that y[] has a single maximum
   and is monotonically decreasing below this value. The approach uses a 2-bin
   or a 3-bin centroid, depending on which is more appropriate for the specific
   shape of the data around the maximum.

   n ------ is the length of x[*]

   the returned value is the interpolated bin number of interest

   This routine is ideal for improving the frequency resolution of the power
   spectrum in single-tone frequency tests.
   Similarly, this routine is useful for improving the estimate of the time
   lag that best matches time series data in an autocorrelation procedure.
*/

	int ixp,ixn,ik,k,nav;
	unsigned imx;
	double s,d;
   float maxx=MembershipDegreeArray[0];
   for(int i=1;i<DomainResolution;i++)
   {
		if(MembershipDegreeArray[i] >= maxx)
		{
			maxx=MembershipDegreeArray[i];
			imx=i;
		}
	}
   ixp = imx+1;
   ixn = imx-1;
   if(ixn == -1) {
     ixp = 2;
     imx = 1;
     ixn = 0;
   }
   if(ixp == (int) n-1) {
     ixp = n-1;
     imx = n-2;
     ixn = n-3;
   }
   nav = 1;
   ik = ixn;
   if(y[ixp] < y[ixn])
   {
     if(y[ixp]*y[imx] < y[ixn]*y[ixn]) nav = 2;
   } else
   {
     if (y[ixp]*y[ixp] < y[ixn]*y[imx]) nav = 2;
     else ik = imx;
   }
   s = d = 0.0;
   for (k=ik; k<=nav+ik; k++)
   {
     s = s + k*y[k];
     d = d + y[k];
   }
   return s/d;
}

float FuzzySet::interpolation(float x0, float* xin, float* yin, int n, int order)
{
/*
Given the vectors of x-y points xin[*] and yin[*], this routine interpolates
between samples at x0 and returns the interpolated ordinate value y0.
It is assumed that the x values are ordered but not necessarily equally spaced.
That is, xin[0] < xin[1] < xin[2] < ... xin[n-1].
The method uses a Lagrange polynomial of degree "order" where
"order" = 1,2,3,...9 (i.e., any integer value less than 10).
The only other restriction on "order" is that "order" < "n", that is there
must be at least "order+1" data points.
(see Abramowitz and Stegun, pg. 878, eqns 25.2.1 and 25.2.2).
*/
	double p, h, y0;
	double tmp[10];
	int i, ihi, ilo, j, jj, npts, n2, nhi, nstart, nstop;

     if(x0 < xin[0] || x0 > xin[n-1]) return 0.0;
     if(order <= 0 || order >= 10) return  0.0;
     npts = order+1;
     if(order >= n) return 0.0;
     n2 = npts >> 1;
     nhi = n-n2-1;
     if(x0 <= xin[n2]) nstart = 0;
     else if (x0 >= xin[nhi]) { nstart = n-npts; }
     else {            /* find starting point through bisection               */
         ilo = n2;  ihi = nhi;
         while (ihi-ilo > 1) {
             i=(ihi+ilo) >> 1;  /* integer divide by two                      */
             if(xin[i] > x0) ihi = i;
             else ilo = i;
         }
         nstart = ihi - n2;
     }

     nstop = nstart + npts;

     for (j=0; j< npts; j++) { tmp[j] = x0-xin[j+nstart]; }

     y0 = 0.0;          /* start of Lagrange loop  */
     for(i=nstart; i<nstop; i++)
     {
         p = yin[i];
         jj = 0;
         for (j=nstart; j< nstop; j++)
         {
             if(j != i)
             {
                  h = xin[i] - xin[j];
                  if(h == 0.0)  return 0.0;
                  p *= tmp[jj]/h;
             }
             jj++;
         }
         y0 += p;
     }                /* end of Lagrange loop */
     return y0;
}

/*----------------------------------------------------------------------------*/

void LinguisticVariable::Initialize(int nb_terms, float min, float max)
{
  NbFzyTerm = nb_terms;
  DomainMax = max;
  DomainMin = min;
  point.Initialize(NbFzyTerm+2);
  FzySet = new FuzzySet[NbFzyTerm];
  float step = (max-min)/float(NbFzyTerm+1);
  point[0]= DomainMin;
  for(int i=1; i<NbFzyTerm+2; i++) point[i] = point[i-1] + step;
  for(int i=0;i<NbFzyTerm;i++)FzySet[i].Initialize(NB_RESOLUTION, DomainMin, DomainMax);
  CreateFuzzySet();
  Name=0;
}

void LinguisticVariable::CreateFuzzySet()
{
  float scale = (DomainMax-DomainMin)/(NB_RESOLUTION+1);
  for(int i=0;i<NbFzyTerm;i++)
  {
// create left shouder curve
    if(i == 0)FzySet[i].DefineMembershipFunction(point[i],1,point[i+1],1,point[i+2],0);
// create right shouder curve
    else if(i == NbFzyTerm-1)FzySet[i].DefineMembershipFunction(point[i],0,point[i+1],1,point[i+2],1);
// otherwise, create triangular curve
    else FzySet[i].DefineMembershipFunction(point[i],0,point[i+1]-scale,1,point[i+1]+scale,1,point[i+2],0);
  }
  for(int i=0;i<NbFzyTerm;i++)FzySet[i].Labeling(i);
}

void LinguisticVariable::Print(ostream &os)
{
  os<<"Linguistic Variable : "<<Name<<endl;
  os<<"There are "<<NbFzyTerm<<" terms(fuzzy set)"<<endl;
  for(int i=0;i<NbFzyTerm;i++)
  {
    FzySet[i].Print(os);
  }
}

int LinguisticVariable::FindMaxDegreeFuzzySet(float t)
{
  int i, ndx;
  float max=-1e20;
  float *d = new float[NbFzyTerm];
  for(int i=0;i<NbFzyTerm;i++)d[i] = FzySet[i].membership(t);
  for(int i=0;i<NbFzyTerm;i++)if(d[i]>max){max=d[i];ndx=i;}
  delete [] d;
  return ndx;
}

float LinguisticVariable::FindFuzzySetValue(int ndx_fzy_set, float t)
{
  float tt = FzySet[ndx_fzy_set].membership(t);
  return tt;
}

FuzzySet LinguisticVariable::GetFuzzySet(int ndx_fzy_set)
{
  FuzzySet fzy;
  fzy.Initialize(FzySet[ndx_fzy_set]);
  return fzy;
}

/*----------------------------------------------------------------------------*/

void FuzzyRule::Initialize(int nb)
{
  nb_antecedent = nb;
  antecedent = new FuzzySet[nb];
}

FuzzySet& FuzzyRule::Inferencing(f1D &para)
{
  int i;
  float *v = new float[nb_antecedent];
  for(i=0;i<nb_antecedent;i++)v[i] = antecedent[i].membership(para.m[i]);
  float min = 1.0;
  for(i=0;i<nb_antecedent;i++)if(v[i]<min)min=v[i];
  InferResult = conclusion;
  InferResult.Scale(min);
  return InferResult;
}

float FuzzyRule::Effectiveness(f1D &para, float solution)
{
  float eff = 1.0;
  for(int i=0;i<nb_antecedent;i++)eff *= antecedent[i].membership(para.m[i]);
  eff *= conclusion.membership(solution);
  return eff;
}

void FuzzyRule::Print(ofstream &out)
{
  out<<"IF ";
  for(int i=0;i<nb_antecedent-1;i++)out<<antecedent[i].Label()<<" AND ";
  out<<antecedent[nb_antecedent-1].Label()<<"  THEN ";
  out<<conclusion.Label()<<endl;
}

/*----------------------------------------------------------------------------*/