/* ts_tools defines basis functions used to compute jacobians, outputs */
#include <math.h>
#include "ts_tools.h"


void process(	double	*f,
   				double	*xMinC,
			 	double	*memb,
			 	double	*y,
			 	double 	*sumMemb,
			 	double	*x,
			 	double	*cp,
			 	double	*vp,
			 	double	*lp,
			 	int 	nIn,
				int		nOut,
				int		nRules,
				int		norm,
				int		memb_type,
				double		m)
{
	int i,j,r,s;
	int iOffIn,iOffOut,iOff2,rOff;
	int nInSq,nInOut;
	
						
	/* Compute the differences between x and the centers */
	for (i=0,iOffIn=0;i<nRules;i++,iOffIn+=nIn)
 	{
 		for(j=0;j<nIn;j++)
 		{
 			xMinC[iOffIn+j]=x[j]-cp[iOffIn+j];
 		}
 	}
 	
 	/* Compute the membership functions*/
 	nInSq=nIn*nIn;
	for (i=0,iOffIn=0,iOff2=0;i<nRules;i++,iOffIn+=nIn,iOff2+=nInSq)
 	{
 		memb[i] = 0;
 		/* Non diagonal elements */ 		
 		for(r=0,rOff=0;r<nIn;r++,rOff+=nIn)
 		{
 			for(s=r+1;s<nIn;s++)
 			{
 				memb[i] += xMinC[iOffIn+r]*vp[iOff2+rOff+s]*xMinC[iOffIn+s];
 			}
 		}
 		
 		memb[i] *= 2;
 		
 		/* Diagonal elements */		
		for(r=0,rOff=0;r<nIn;r++,rOff+=nIn+1)
 		{
 			memb[i] += xMinC[iOffIn+r]*vp[iOff2+rOff]*xMinC[iOffIn+r];
 		}
 		
 		switch (memb_type)
 		{
 			case 1:
 				memb[i]=exp(-memb[i]);
 			break;
 			
 			case 2:
 				memb[i]= pow((memb[i] + 1e-100),-1/(m-1));
 		}
 				
 				
 	}
	
	/* Compute the sum of the memberships */
	*sumMemb = 2.2251e-308;
	for (i=0;i<nRules;i++)
 	{	
 		*sumMemb += memb[i];
 	}
 	
 	/* Compute the local outputs */
 	
 	nInOut=(nIn+1)*nOut;
 	
 	for (i=0,iOffOut=0,iOff2=0;i<nRules;i++,iOff2+=nInOut,iOffOut+=nOut)
 	{	
	 	for (j=0;j<nOut;j++)
	 	{
	 		y[iOffOut+j] = 0;
	 		
	 		for(r=0,rOff=0;r<nIn;r++,rOff+=nOut)
	 		{
	 			y[iOffOut+j] += x[r]*lp[iOff2+rOff+j];
	 		}
	 		
	 		y[iOffOut+j] += lp[iOff2+nOut*nIn+j];
	 	}
 	}
 	
 	/* Compute the global outputs */
	for (j=0; j<nOut; j++)
	{
		f[j] = 0;
		for (i=0,iOffOut=0;i<nRules;i++,iOffOut+=nOut)
	 	{	
	 		f[j] += memb[i] * y[iOffOut+j];
	 	}
	 	if (norm)
	 	{
	 		f[j] /= *sumMemb;
	 	}
	 }
}
