#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

input int SDIPeriod = 19;
input int ZScorePeriod = 36;
input int MaxBarsToProcess = 3000;

double ZScoreBuffer[];
double VolumeEMA[];
double RangeEMA[];
double BuyRaw[];
double SellRaw[];
double BuyEMA[];
double SellEMA[];
double SDI[];

int OnInit()
{
   if(SDIPeriod < 2 || 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");
   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 currentHigh = high[shift];
   double currentLow = low[shift];
   double rangeHigh = currentHigh;
   double rangeLow = currentLow;

   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);
}

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);

   if(rates_total < ZScorePeriod + SDIPeriod + 2)
      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(SDI,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(SDI,0.0);

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

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

   const double alpha = 2.0 / (SDIPeriod + 1.0);
   const 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])
         SDI[bar] = -100.0 * normalizedDifference;
      else
         SDI[bar] = 100.0 * normalizedDifference;
   }

   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 = SDI[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] = (SDI[outBar] - mean) / standardDeviation;
      else
         ZScoreBuffer[outBar] = 0.0;
   }

   return(rates_total);
}