//+------------------------------------------------------------------+
//|                                           PolyFitVelocity_v1.mq4 |
//|                                Copyright © 2007, TrendLaboratory |
//|            http://finance.groups.yahoo.com/group/TrendLaboratory |
//|                                   E-mail: igorad2003@yahoo.co.uk |
//+------------------------------------------------------------------+
#property copyright "Copyright © 2007, TrendLaboratory"
#property link      "http://finance.groups.yahoo.com/group/TrendLaboratory"

#property indicator_separate_window
#property indicator_buffers 2
#property indicator_color1  LightBlue
#property indicator_width1  2
#property indicator_color2  OrangeRed
#property indicator_width2  1
#property indicator_style2  2
#property indicator_level1  0
//---- external variables
extern int Price        =   0;   //Apply to Price(0-Close;1-Open;2-High;3-Low;4-Median price;5-Typical price;6-Weighted Close)
extern int Length       =  14;   //Number of bars for an evaluation
extern int Degree       =   2;   //Degree of a Polynomial(no more 12)
extern int MA_Length    =   1;   //Period of Preliminary Smoothing
extern int MA_Mode      =   0;   //Mode of MA:0-SMA,1-EMA,2-Wilders(SMMA),3-LWMA
extern int Signal       =   3;   //Period of Signal MA
extern int Sig_Mode     =   0;   //Mode of Signal MA:0-SMA,1-EMA,2-Wilders(SMMA),3-LWMA
extern int FitMode      =   1;   //Fitting Mode: 0-Fitting,1-Moving
extern int CountBars    =1000;   //Maximum Number of Bars(0-all Bars)
//---- indicator buffers
double PolyVel[];
double SigVel[];

//---- global variables
double PriceArray[];
int cBars;
//+---------------------------------------------------------------------------+
//| Custom indicator initialization function                                  |
//+---------------------------------------------------------------------------+
int init()
{
   IndicatorDigits(MarketInfo(Symbol(), MODE_DIGITS )+2);
   //IndicatorBuffers(3);
   SetIndexStyle(0, DRAW_LINE);
   SetIndexStyle(1, DRAW_LINE);
   SetIndexBuffer(0, PolyVel);
   SetIndexBuffer(1, SigVel);
  
   if(CountBars == 0) cBars = Bars;
   else cBars = CountBars;
   
   SetIndexDrawBegin( 0,Bars-cBars+(Length+MA_Length));
   SetIndexDrawBegin( 1,Bars-cBars+(Length+MA_Length+Signal));  
   
   string short_name="PolyFitVelocity("+Price+","+Length+","+Degree+")";
   IndicatorShortName(short_name);
   SetIndexLabel(0,"PolyFitVelocity");
   SetIndexLabel(1,"Signal");
   
   ArrayResize(PriceArray, Length);
   
   return(0);
}

double PolyFit(double price[],int deg,int len,int bar)
{
double   result=0,Sum=0;
double   AX[13,13],BX[13],CX[13],ZX[13],Pow[26];
int      Row[13];  
// FILL MATRIX FIRST
   if (len<=1)	Sum=price[0];
   else
   {	
   
   for (int j=1;j<=deg+1;j++) BX[j]=0;
		
		for (int k=1;k<=len;k++) 
  	   {
  	   double YK=price[len-k];
  	   double XK=k;
  	   double Prod=1;
  	      for (j=1;j<=deg+1;j++)
    	   {
    	   BX[j]+= YK*Prod;
    	   Prod*=XK;
  		   }
      }
      
   for (j=0;j<=2*deg;j++) Pow[j]=0;
	
	Pow[0]=len;
	
      for (k=1;k<=len;k++)
	   {
  	   XK=k;
  	   Prod=k;
  	      for (j=1;j<=2*deg;j++)
  	      {
    	   Pow[j]+= Prod;    	   
    	   Prod*= XK;
  		   }
      }	
	
	   for (j=1;j<=deg+1;j++)
	   {	
         for (int l=1;l<=deg+1;l++) 
  	      AX[j,l]=Pow[j+l-2];
      }	
//NOW SOLVE FOR COEFFICIENTS

	for (j=1;j<=deg+1;j++) Row[j]=j;
  	
      for (int i=1;i<=deg;i++)
  	   {
  	      for (k=i+1;k<=deg+1;k++)
         {
            if(MathAbs(AX[Row[k],i]) > MathAbs(AX[Row[i],i]))
      	   {	
      	   int Temp=Row[i];
      	   Row[i]=Row[k];
      	   Row[k]=Temp;
            }
         }

         for (k=i+1;k<=deg+1;k++)
    	   {
    	   if(AX[Row[i],i]!=0) AX[Row[k],i]=AX[Row[k],i]/AX[Row[i],i];
    		   for (l=i+1;l<=deg+1;l++)
      	   AX[Row[k],l]=AX[Row[k],l]-AX[Row[k],i]*AX[Row[i],l];
    		}
		}

   ZX[1]=BX[Row[1]];
      
      for (k=2;k<=deg+1;k++)
  	   {
  	   Sum=0;
      for (l=1;l<=k-1;l++) Sum+=AX[Row[k],l]*ZX[l];
  	   
  	   ZX[k]=BX[Row[k]]-Sum;
	   }
     
	if(AX[Row[deg+1],deg+1] != 0) CX[deg+1]=ZX[deg+1]/AX[Row[deg+1],deg+1];
      	
      for (k=deg;k>=1;k--)
  	   {
  	   Sum=0;
  	   //l=k+1;  		
  	   for (l=k+1;l<=deg+1;l++) Sum += AX[Row[k],l]*CX[l]; //AX[Row[KX],Col]*CX[Col];
  	   CX[k]=(ZX[k]-Sum)/AX[Row[k],k];
	   }
    
//NOW COMPUTE NEXT POINT IN SERIES AND RETURN}

	Sum=deg*CX[deg+1];
	
	if(deg > 1) 
	  for (k=deg;k>=2;k--) {Sum=(k-1)*CX[k]+Sum*(len+bar);}
 
   }
   return (Sum);
}	
//+------------------------------------------------------------------+
//| PolyFitVelocity_v1                                               |
//+------------------------------------------------------------------+
int start()
{
   int    i, j, k, m, shift, counted_bars=IndicatorCounted(), limit;
        
   
   if ( counted_bars < 0 )  return(0);
   if ( counted_bars ==0 )  limit=Bars-1; 
         
   if (FitMode==1 && counted_bars > 0) {limit = Bars - counted_bars; int len = Length+MA_Length;}
   else
   if (FitMode==0) {limit = Length; len = cBars - Length+1;} 
   
   if ( counted_bars < 1 || FitMode==0)
   { 
      for(i=1;i<len;i++)
      { 
      PolyVel[cBars-i]= EMPTY_VALUE;    
      SigVel[cBars-i] = EMPTY_VALUE;
      }
   }   

   for(shift = 0;shift<limit;shift++)
   {
      for(j = 0; j < Length; j++) 
      PriceArray[j] = iMA(NULL,0,MA_Length,0,MA_Mode,Price,shift*FitMode + j);
      
   PolyVel[shift] = PolyFit(PriceArray,Degree,Length,shift*(FitMode-1));
   }
   if ( FitMode == 1)   
      for(shift = limit;shift>=0;shift--)
      SigVel[shift] = iMAOnArray(PolyVel,0,Signal,0,Sig_Mode,shift); 
        
   return(0);
}
//+---------------------------------------------------------------------------+