#property copyright "Copyleft © 2007, GammaRat Forex"
#property link      "http://www.gammarat.com/Forex/"

//Hoe I use this:
// the indicators are currently set up for a 1 hour type chart, to produce
// the trend, 1 deviation on either side, 2 deviations on either side, and 2 deviations 
// on either side calculated off the highs and lows.
// The "Samples" parameter is used to set the 
// die off of the exponentially averaged deviations. Bigger numbers give
// smoother edges.
// By setting
// the deviation levels to zero, they don't get plotted.
// don't bother with the linear gain, leave it at sero.
//
// For the acceleration de-gain (AccDegain) decreasing it increases the 
// oscillations, increasing it decreases the oscillations. Confusing, no?
// Old habits die hard.
//
//I've found it works best set up as follows -
// One filter with the de-gain set at 0, and all the deviations 
// turned on
// Another with the deviation levels set to 0, and the degain set to 1.
// this gives you slower oscillations.
//A third with the deviations set to 0, and the degain set to -1. This gives yu faster oscillations.
// then use it as you would running averages.
//
//I've set up a primitive matrix liberary and pasted it below. It's very weak, I hope to incorporate
// DLL's along the line of LINPACK later, so I haven't put much time into it.
//
//It is currently set up to allow it to be consistent across time frames
// normalized to 1 hour views, for frames of 4 hours or less. For days and longer, or for 
// independent looks at different time frames.
// turn off the normalize_to_1hr variable.
//
//the first several hundred points are garbage, that's when the filter is ramping up.
//
//Known bugs - I don't have the end point set correctly, so the filter continues running on
// ticks. For now, just go back and forth between (say) your primary chart and another one
//without changing anything. That clears the buffers.
//
//and sorry about the matrices, I want to leave the code open for
//higher order filtering.
//
// - oddity - I use geometric averaging on the point (MathPow(HLCC,0.25)) rather than the
// usual sums. See the funtion "get_avg" below, and change it
// if you like. Again, just a personal pref....

//comments & questions - brobeck@ns.sympatico.ca
//
//cheers;  




//#include <MatrixMath.mqh>
#property indicator_chart_window
#property indicator_buffers 7
#property indicator_color1 Yellow
#property indicator_color2 Yellow
#property indicator_color3 Yellow
#property indicator_color4 Red
#property indicator_color5 Red
#property indicator_color6 Red
#property indicator_color7 Red
#property indicator_style1 0
#property indicator_style2 2
#property indicator_style3 2
#property indicator_style4 2
#property indicator_style5 2
#property indicator_style6 0
#property indicator_style7 0


//---- input parameters
extern int  Samples=120;
extern double DevLevel1 = 1;
extern double DevLevel2 = 2;
extern double ZGain = 0;
extern double AccDeGain = 0;
extern bool Normalize_to_1hr = true;

//---- buffers
double KalmanBuffer[];
double KalmanBufferPlus1[];
double KalmanBufferNeg1[];
double KalmanBufferPlus2[];
double KalmanBufferNeg2[];
double KalmanBufferPlus3[];
double KalmanBufferNeg3[];
double xkk[2][1],xkk1[2][1],xk1k1[2][1],yk[1][1],zk[1][1];
double Pkk[2][2],Pkk1[2][2],Pk1k1[2][2];
double Qkk[2][2],Hk[1][2],Hk_t[2][1],Sk[1][1],Sk_inv[1][1],Rk[1][1],Kk[2][1];
double Fk[2][2],Fk_t[2][2],GGT[2][2];
double Eye[2][2];
double temp22[2][2],temp21[2][1],temp12[1][2],temp11[1][1];
double LookAhead=0;

//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
   if(LookAhead <0)LookAhead=0;
   SetIndexStyle(0,DRAW_LINE);
   SetIndexBuffer(0,KalmanBuffer);
   SetIndexShift(0,LookAhead);
   SetIndexDrawBegin(0,LookAhead);
   SetIndexLabel(0,"Kalman Trend");
   if(MathAbs(DevLevel1) > 0) {
      SetIndexStyle(1,DRAW_LINE);
      SetIndexBuffer(1,KalmanBufferPlus1);
      SetIndexShift(1,LookAhead);
      SetIndexDrawBegin(1,LookAhead);
      SetIndexLabel(1,"Kalman +" + DoubleToStr(DevLevel1,1)  + " STD");
      SetIndexStyle(2,DRAW_LINE);
      SetIndexBuffer(2,KalmanBufferNeg1);
      SetIndexShift(2,LookAhead);
      SetIndexDrawBegin(2,LookAhead+1);
      SetIndexLabel(2,"Kalman -" + DoubleToStr(DevLevel1,1) + " STD");
    }
   if(MathAbs(DevLevel2) > 0){
      SetIndexStyle(3,DRAW_LINE);
      SetIndexBuffer(3,KalmanBufferPlus2);
      SetIndexShift(3,LookAhead);
      SetIndexDrawBegin(3,LookAhead);
      SetIndexLabel(3,"Kalman +" + DoubleToStr(DevLevel2,1) + " STD");
      SetIndexStyle(4,DRAW_LINE);
      SetIndexBuffer(4,KalmanBufferNeg2);
      SetIndexShift(4,LookAhead);
      SetIndexDrawBegin(4,LookAhead);
      SetIndexLabel(4,"Kalman -" + DoubleToStr(DevLevel2,1) + " STD");
      
      SetIndexStyle(5,DRAW_LINE);
      SetIndexBuffer(5,KalmanBufferPlus3);
      SetIndexShift(5,LookAhead);
      SetIndexDrawBegin(5,LookAhead);
      SetIndexLabel(5,"Kalman High +" + DoubleToStr(DevLevel2,1) + " STD");
      SetIndexStyle(6,DRAW_LINE);
      SetIndexBuffer(6,KalmanBufferNeg3);
      SetIndexShift(6,LookAhead);
      SetIndexDrawBegin(6,LookAhead);
      SetIndexLabel(6,"Kalman Low -" + DoubleToStr(DevLevel2,1) + " STD");
   }
//----
   return(0);
  }
//| Point and figure                                |
//+------------------------------------------------------------------+
int start()
{
   compute();
   return(0);
}
int compute(){
   int    i,j,counted_bars=IndicatorCounted();
   static bool initiated=false;
   double c0;
   static double c0_above=0, c0_below = 0;
   static double c1_above=0, c1_below = 0;
   static double acc_sigma2;
   static double z_sigma2;
   static double normalization_time;
   static double working_samples;
   double z,diff; 
   double ts[3][1];
   static double dev_level1,dev_level2; 
  
   //----
   if(Bars<5) return(0);
   if(counted_bars >= Bars-1) return(0);
   if( !initiated){
      if(Period() < 241 && Normalize_to_1hr){
         normalization_time = Period()/60.;
      }else{
         normalization_time=1;
      }      
      dev_level1 = DevLevel1*MathSqrt(2);
      dev_level2 = DevLevel2*MathSqrt(2);
      xkk[0][0] = get_avg(Bars-1);
      xkk[1][0] = get_avg(Bars-1)-get_avg(Bars-15);
      //MatrixPrint(xkk);
      MatrixEye(Pkk);
      Fk[0][0] = 1;
      Fk[0][1]=1;
      Fk[1][0]=0;
      Fk[1][1] =1;
      MatrixTranspose(Fk,Fk_t);
      Qkk[0][0] = 1;
      Qkk[0][1] = 0;
      Qkk[1][0]=0;
      Qkk[1][1] = 1;
      GGT[0][0]=.25;
      GGT[1][0]=.5;
      GGT[0][1]=.5;
      GGT[1][1]=1;
      Hk[0][0] = 1;
      Hk[0][1] = 0;
      MatrixTranspose(Hk,Hk_t);
      
      Rk[0][0] = .1;
      MatrixEye(Eye);
      z_sigma2 = MathPow(10,-ZGain);
      acc_sigma2 = Point*MathPow(10,-AccDeGain)*normalization_time;
      working_samples = Samples/normalization_time;
      acc_sigma2 *= acc_sigma2;
      initiated = true;
      MatrixScalarMul(acc_sigma2,Pkk);
      //Print(1./MathSqrt(acc_sigma2));
      //Print(" ");
      
   }

   
   for(i=Bars-counted_bars; i>0;i--){
  
      
      //make copies of the last iteration;
      ArrayCopy(Pk1k1,Pkk);
      ArrayCopy(xk1k1,xkk);
      //predict the state
      //Print("Doing Predict");
      MatrixMul(Fk,xk1k1,xkk1);
      //MatrixPrint(xkk1);
      MatrixMul(Fk,Pk1k1,temp22);
      MatrixMul(temp22,Fk_t,Pk1k1);
      MatrixAdd(Pk1k1,Qkk,Pkk1,1);
      //Print("Predict Ended");
      KalmanBuffer[i-1] = (xkk1[0][0]*Point);
      //update cycle
      //update
      //innovation
     
      MatrixZero(yk);
      zk[0][0] = get_avg(i);
      
      diff = (zk[0][0]*Point-KalmanBuffer[i]);
      if(diff>=0) {
         diff = diff*diff;
         c0_above = (c0_above*(working_samples-1)+MathPow(get_avg(i)*Point-KalmanBuffer[i],2))/working_samples;
         c0_below = c0_below*(working_samples-1)/working_samples;
         c1_above = (c1_above*(working_samples-1)+MathPow(High[i]-KalmanBuffer[i],2))/working_samples;
         c1_below = c1_below*(working_samples-1)/working_samples;
      }else{
         diff = diff*diff;
         c0_below = (c0_below*(working_samples-1)+MathPow(get_avg(i)*Point-KalmanBuffer[i],2))/working_samples;
         c0_above = c0_above*(working_samples-1)/working_samples;
         c1_below = (c1_below*(working_samples-1)+MathPow(Low[i]-KalmanBuffer[i],2))/working_samples;
         c1_above = c1_above*(working_samples-1)/working_samples;
      }  
      if(MathAbs(DevLevel1)>0){
         KalmanBufferPlus1[i-1] = KalmanBuffer[i]+dev_level1*MathSqrt(c0_above);
         KalmanBufferNeg1[i-1] = KalmanBuffer[i]-dev_level1*MathSqrt(c0_below);
      }
      if(MathAbs(DevLevel2)>0){   
         KalmanBufferPlus2[i-1] = KalmanBuffer[i]+dev_level2*MathSqrt(c0_above);
         KalmanBufferNeg2[i-1] = KalmanBuffer[i]-dev_level2*MathSqrt(c0_below);
         KalmanBufferPlus3[i-1] = KalmanBuffer[i]+dev_level2*MathSqrt(c1_above);
         KalmanBufferNeg3[i-1] = KalmanBuffer[i]-dev_level2*MathSqrt(c1_below);
      } 
      zk[1][0]= get_avg(i)-get_avg(i+1);
      Rk[0][0] = z_sigma2;
      
      //acc += zk[1][0];
      //acc=MathAbs(acc);
      //acc_avg= (acc+(Samples-1)*acc_avg)/Samples;
      //      acc_sigma2 = (MathPow(acc-acc_avg,2)/Samples+(Samples-1)*acc_sigma2)/Samples;
      //Print(acc_sigma2);
      ArrayCopy(Qkk,GGT);
      MatrixScalarMul(acc_sigma2,Qkk);
      MatrixMul(Hk,xkk1,temp11);
      //MatrixPrint(xkk1);
      MatrixAdd(zk,temp11,yk, -1); 
      MatrixMul(Hk,Pkk1,temp12);
      MatrixMul(temp12,Hk_t,temp11);
      MatrixAdd(temp11,Rk,Sk,1);
      MatrixInvert(Sk,Sk_inv);
      MatrixMul(Pkk1,Hk_t,temp21);
      MatrixMul(temp21,Sk_inv,Kk);
      MatrixMul(Kk,yk,temp21);
      MatrixAdd(temp21,xkk1,xkk,1);
      
      MatrixMul(Kk,Hk,temp22);
      MatrixAdd(Eye,temp22,temp22,-1);
      MatrixMul(temp22,Pkk1,Pkk);
      
    
   }
}
 
 double get_avg(int k){
   return(MathPow((High[k]*Low[k]*Close[k]*Close[k]),1/4.)/Point);
}        
double determinant_fixed(double a[][], int k){
   double z;
    switch(k){
      case 1: return(a[0][0]);
         break;
      case 2: return( a[0][0]*a[1][1]-a[1][0]*a[0][1]);
         break;
      case 3:
         z = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + 
             a[0][2]*a[1][0]*a[2][1] -
            (a[2][0]*a[1][1]*a[0][2] + a[1][0]*a[0][1]*a[2][2] + 
             a[0][0]*a[2][1]*a[1][2]);
         return(z);
         break;
      default:      
          Print("array_range too large for determinant calculation");
         return(0);
   }  
}
double determinant(double a[][]){
   double z;
   if(ArrayRange(a,0) != ArrayRange(a,1)){
      Print("Non-square array entered into determinant");
      return(0);
    }
    switch(ArrayRange(a,0)){
      case 1: return(a[0][0]);
         break;
      case 2: return( a[0][0]*a[1][1]-a[1][0]*a[0][1]);
         break;
      case 3:
         z = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + 
             a[0][2]*a[1][0]*a[2][1] -
            (a[2][0]*a[1][1]*a[0][2] + a[1][0]*a[0][1]*a[2][2] + 
             a[0][0]*a[2][1]*a[1][2]);
         return(z);
         break;
      default:      
          Print("array_range too large for determinant calculation");
         return(0);
   }  
}
  
int MatrixAdd(double a[][], double b[][],double & c[][],double w=1){
   int i,j;
   for(i=0;i<ArrayRange(c,0);i++){
      for(j=0;j<ArrayRange(c,1);j++){
         c[i][j] = a[i][j]+w*b[i][j];
      }
   }
    return(0);   
}
int MatrixMul(double a[][], double b[][],double & c[][]){
   int i,j, k;
   if(ArrayRange(a,1) != ArrayRange(b,0)){
      Print("array Range Mismatch in Matrix Mul");
      return(-1);
   }
   for(i=0;i<ArrayRange(c,0);i++){
      for(j=0;j<ArrayRange(c,1);j++){
         c[i][j] = 0;
         for(k=0;k<ArrayRange(a,1);k++){
           c[i][j] += a[i][k]*b[k][j];
         }  
      }
   }
  return(0);   
}
int MatrixScalarMul(double b,double & a[][]){
   int i,j;
   for(i =0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         a[i][j] *= b;
      }
   }
}         
int MatrixInvert(double a[][],double & b[][]){
   double d;
   double temp1[4][4],temp2[4][4];
   int i,j,i1,j1,k;
   
   if(ArrayRange(a,0) != ArrayRange(a,1)){
      Print("Matrix Inverse attempted on non square matrix. Not implemented yet");
      return(-1);
   }
   if(ArrayRange(a,1) > 3){
      Print("Bugger off, not implemented yet");
      return(-1);
   }
   d = determinant(a);
   if(MathAbs(d)<0.000001){
    Print("unstable matrix");
    d = 1;
   }
   k = ArrayRange(a,0);
   d = determinant(a);
   for(i1=0;i1<k;i1++){
      for(j1=0;j1<k;j1++){
         temp1[i1][j1] = a[i1][j1];
      }
   }   
   for(i=0;i<k;i++){
      for(j=0;j<k;j++){
         //copy the arry. There must be a better way!!
         ArrayCopy(temp2,temp1);
         for(j1=0;j1<k;j1++){
            if(j1==j){
               temp2[i][j1] = 1;
            }else{
               temp2[i][j1] = 0;
            }
         }   
         b[i][j]=determinant_fixed(temp2,k)/d;
      }
   }
}   
int MatrixTranspose(double a[][],double & b[][]){
   int i,j;
   for(i=0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         b[j][i] = a[i][j];
      }
   }
   return(0);
 }
int MatrixZero(double & a[][]){
   int i,j; 
   for(i= 0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         a[i][j]=0;
      }
   }
   return(0);
} 
int MatrixEye(double & a[][]){
   int i,j; 
   for(i= 0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         a[i][j]=0;
         if(i==j)a[i][j]=1;
      }
   }
   return(0);
} 
int MatrixPrint(double a[][]){
   int i,j; 
   for(i= 0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
      Print(i," ",j," ",a[i][j]);
      }
   }
   return(0);
}           
//+------------------------------------------------------------------+