#property strict
#property indicator_separate_window
#property indicator_buffers 1
#property indicator_color1 clrAqua
#property indicator_width1 2
#property indicator_level1 2.0
#property indicator_level2 -2.0
#property indicator_level3 0.0

enum SDIPreSmoothMethod
{
   SDI_SMA  = 0,
   SDI_EMA  = 1,
   SDI_SMMA = 2,
   SDI_LWMA = 3,
   SDI_VWMA = 4
};

input int SDIPeriod = 19;
input int SDIPreSmoothPeriod = 3;
input SDIPreSmoothMethod SDIPreSmoothType = SDI_SMA;
input int ZScorePeriod = 36;
input int MaxBarsToProcess = 3000;

double ZScoreBuffer[];
double VolumeEMA[];
double RangeEMA[];
double BuyRaw[];
double SellRaw[];
double BuyEMA[];
double SellEMA[];
double SDIRaw[];
double SDISmoothed[];

int OnInit()
{
   if(SDIPeriod < 2 || SDIPreSmoothPeriod < 1 ||
      ZScorePeriod < 2 || MaxBarsToProcess < 10)
      return(INIT_PARAMETERS_INCORRECT);

   SetIndexBuffer(0,ZScoreBuffer);
   ArraySetAsSeries(ZScoreBuffer,true);
   SetIndexStyle(0,DRAW_LINE,STYLE_SOLID,2,clrAqua);
   SetIndexLabel(0,"SDI Z-Score");
   SetIndexEmptyValue(0,EMPTY_VALUE);

   IndicatorDigits(2);
   IndicatorShortName("Sibbet SDI Adaptive Z-Score PreSmooth");
   SetLevelStyle(STYLE_DOT,1,clrSilver);

   return(INIT_SUCCEEDED);
}

void PrepareArray(double &array[],const int size)
{
   if(ArraySize(array) != size)
      ArrayResize(array,size);
   ArraySetAsSeries(array,true);
}

double SafeExp(double value)
{
   if(value > 50.0)
      value = 50.0;
   if(value < -50.0)
      value = -50.0;
   return(MathExp(value));
}

double GetWeightedPrice(const double highValue,
                        const double lowValue,
                        const double closeValue)
{
   return(highValue + lowValue + 2.0 * closeValue);
}

double GetTwoBarRange(const double &high[],
                      const double &low[],
                      const int shift,
                      const int rates_total)
{
   double rangeHigh = high[shift];
   double rangeLow = low[shift];

   if(shift + 1 < rates_total)
   {
      if(high[shift + 1] > rangeHigh)
         rangeHigh = high[shift + 1];
      if(low[shift + 1] < rangeLow)
         rangeLow = low[shift + 1];
   }

   return(rangeHigh - rangeLow);
}

double CalculatePreSmooth(const int shift,
                          const int rates_total,
                          const int calcStart,
                          const long &tick_volume[])
{
   if(SDIPreSmoothPeriod <= 1)
      return(SDIRaw[shift]);

   if(SDIPreSmoothType == SDI_EMA)
   {
      double alpha = 2.0 / (SDIPreSmoothPeriod + 1.0);
      if(shift == calcStart)
         return(SDIRaw[shift]);

      return(alpha * SDIRaw[shift]
             + (1.0 - alpha) * SDISmoothed[shift + 1]);
   }

   if(SDIPreSmoothType == SDI_SMMA)
   {
      if(shift == calcStart)
         return(SDIRaw[shift]);

      return((SDISmoothed[shift + 1] * (SDIPreSmoothPeriod - 1.0)
              + SDIRaw[shift]) / SDIPreSmoothPeriod);
   }

   double weightedSum = 0.0;
   double weightSum = 0.0;

   for(int offset = 0; offset < SDIPreSmoothPeriod; offset++)
   {
      int index = shift + offset;
      if(index >= rates_total || index > calcStart)
         return(EMPTY_VALUE);

      double weight = 1.0;
      if(SDIPreSmoothType == SDI_LWMA)
         weight = SDIPreSmoothPeriod - offset;
      else if(SDIPreSmoothType == SDI_VWMA)
         weight = (double)tick_volume[index];

      weightedSum += SDIRaw[index] * weight;
      weightSum += weight;
   }

   if(weightSum <= 0.0)
   {
      if(SDIPreSmoothType == SDI_VWMA)
      {
         double simpleSum = 0.0;
         for(int fallbackOffset = 0; fallbackOffset < SDIPreSmoothPeriod; fallbackOffset++)
            simpleSum += SDIRaw[shift + fallbackOffset];
         return(simpleSum / SDIPreSmoothPeriod);
      }
      return(SDIRaw[shift]);
   }

   return(weightedSum / weightSum);
}

int OnCalculate(const int rates_total,
                const int prev_calculated,
                const datetime &time[],
                const double &open[],
                const double &high[],
                const double &low[],
                const double &close[],
                const long &tick_volume[],
                const long &volume[],
                const int &spread[])
{
   ArraySetAsSeries(high,true);
   ArraySetAsSeries(low,true);
   ArraySetAsSeries(close,true);
   ArraySetAsSeries(tick_volume,true);

   ArrayInitialize(ZScoreBuffer,EMPTY_VALUE);

   int minimumBars = SDIPeriod + SDIPreSmoothPeriod + ZScorePeriod + 2;
   if(rates_total < minimumBars)
      return(rates_total);

   PrepareArray(VolumeEMA,rates_total);
   PrepareArray(RangeEMA,rates_total);
   PrepareArray(BuyRaw,rates_total);
   PrepareArray(SellRaw,rates_total);
   PrepareArray(BuyEMA,rates_total);
   PrepareArray(SellEMA,rates_total);
   PrepareArray(SDIRaw,rates_total);
   PrepareArray(SDISmoothed,rates_total);

   ArrayInitialize(VolumeEMA,0.0);
   ArrayInitialize(RangeEMA,0.0);
   ArrayInitialize(BuyRaw,0.0);
   ArrayInitialize(SellRaw,0.0);
   ArrayInitialize(BuyEMA,0.0);
   ArrayInitialize(SellEMA,0.0);
   ArrayInitialize(SDIRaw,0.0);
   ArrayInitialize(SDISmoothed,0.0);

   int outputBars = MathMin(MaxBarsToProcess,rates_total);
   int warmupBars = SDIPeriod * 10 + SDIPreSmoothPeriod + ZScorePeriod + 10;
   int calcStart = outputBars - 1 + warmupBars;

   if(calcStart > rates_total - 2)
      calcStart = rates_total - 2;

   double alpha = 2.0 / (SDIPeriod + 1.0);
   double epsilon = 1.0e-12;

   for(int bar = calcStart; bar >= 0; bar--)
   {
      double weightedPrice = GetWeightedPrice(high[bar],low[bar],close[bar]);
      double twoBarRange = GetTwoBarRange(high,low,bar,rates_total);
      double currentVolume = (double)tick_volume[bar];

      if(bar == calcStart)
      {
         VolumeEMA[bar] = currentVolume;
         RangeEMA[bar] = twoBarRange;
      }
      else
      {
         VolumeEMA[bar] = alpha * currentVolume + (1.0 - alpha) * VolumeEMA[bar + 1];
         RangeEMA[bar] = alpha * twoBarRange + (1.0 - alpha) * RangeEMA[bar + 1];
      }

      double previousWeightedPrice = GetWeightedPrice(high[bar + 1],low[bar + 1],close[bar + 1]);
      double volumeNormalized = 0.0;
      if(VolumeEMA[bar] > epsilon)
         volumeNormalized = currentVolume / VolumeEMA[bar];

      double buyPower = volumeNormalized;
      double sellPower = volumeNormalized;

      if(weightedPrice < previousWeightedPrice)
      {
         if(RangeEMA[bar] > epsilon && MathAbs(weightedPrice) > epsilon)
         {
            double exponent = 0.375
                           * (weightedPrice + previousWeightedPrice)
                           / RangeEMA[bar]
                           * (previousWeightedPrice - weightedPrice)
                           / weightedPrice;
            buyPower = volumeNormalized / SafeExp(exponent);
         }
      }
      else if(weightedPrice > previousWeightedPrice)
      {
         if(RangeEMA[bar] > epsilon && MathAbs(previousWeightedPrice) > epsilon)
         {
            double exponent = 0.375
                           * (weightedPrice + previousWeightedPrice)
                           / RangeEMA[bar]
                           * (weightedPrice - previousWeightedPrice)
                           / previousWeightedPrice;
            sellPower = volumeNormalized / SafeExp(exponent);
         }
      }

      BuyRaw[bar] = buyPower;
      SellRaw[bar] = sellPower;

      if(bar == calcStart)
      {
         BuyEMA[bar] = BuyRaw[bar];
         SellEMA[bar] = SellRaw[bar];
      }
      else
      {
         BuyEMA[bar] = alpha * BuyRaw[bar] + (1.0 - alpha) * BuyEMA[bar + 1];
         SellEMA[bar] = alpha * SellRaw[bar] + (1.0 - alpha) * SellEMA[bar + 1];
      }

      double divisor = MathMax(BuyEMA[bar],SellEMA[bar]);
      double dividend = MathMin(BuyEMA[bar],SellEMA[bar]);
      double normalizedDifference = 0.0;
      if(divisor > epsilon)
         normalizedDifference = 1.0 - (dividend / divisor);

      if(SellEMA[bar] > BuyEMA[bar])
         SDIRaw[bar] = -100.0 * normalizedDifference;
      else
         SDIRaw[bar] = 100.0 * normalizedDifference;
   }

   for(int smoothBar = calcStart; smoothBar >= 0; smoothBar--)
   {
      double smoothValue = CalculatePreSmooth(smoothBar,rates_total,calcStart,tick_volume);
      if(smoothValue == EMPTY_VALUE)
         SDISmoothed[smoothBar] = SDIRaw[smoothBar];
      else
         SDISmoothed[smoothBar] = smoothValue;
   }

   for(int outBar = outputBars - 1; outBar >= 0; outBar--)
   {
      double sum = 0.0;
      double sumSquares = 0.0;
      bool valid = true;

      for(int window = 0; window < ZScorePeriod; window++)
      {
         int index = outBar + window;
         if(index >= rates_total || index > calcStart)
         {
            valid = false;
            break;
         }

         double value = SDISmoothed[index];
         sum += value;
         sumSquares += value * value;
      }

      if(!valid)
      {
         ZScoreBuffer[outBar] = EMPTY_VALUE;
         continue;
      }

      double mean = sum / ZScorePeriod;
      double variance = (sumSquares / ZScorePeriod) - (mean * mean);
      if(variance < 0.0)
         variance = 0.0;

      double standardDeviation = MathSqrt(variance);
      if(standardDeviation > epsilon)
         ZScoreBuffer[outBar] = (SDISmoothed[outBar] - mean) / standardDeviation;
      else
         ZScoreBuffer[outBar] = 0.0;
   }

   return(rates_total);
}