//+------------------------------------------------------------------+
//|                  Robust_Location_Estimators.mqh                  |
//|                                                                  |
//|  ROBUST LOCATION (center) ESTIMATOR library.                     |
//|  Taken from Robust_ATR_Pivot (pure array functions, NO state     |
//|  and NO globals). Use it for the STARC MIDLINE and the Z-Score   |
//|  CENTER: pass a window of N prices and get the robust center     |
//|  of that window.                                                 |
//|                                                                  |
//|  WHERE IT PLUGS IN:                                              |
//|   - STARC: add new cases to the "Starc average type" selector.   |
//|       For each bar j, build the window w[] with the N prices     |
//|       ending at j and call the estimator:                        |
//|         midline[j] = BiweightOf(w, N);   // Tukey biweight       |
//|   - Z-Score: same in the center selector (CenterMethod).         |
//|                                                                  |
//|  Element order in w[] does not matter (these are window          |
//|  statistics). No external deps: only Math* and ArrayResize.      |
//|                                                                  |
//|  RECOMMENDED for counter-trend / crypto scale trading (a stable, |
//|  wick-robust band): MedianOf, HodgesLehmannOf, BiweightOf        |
//|  (Tukey), HuberOf, MidhingeOf, TrimeanOf, and the high-breakdown |
//|  MMLocOf / LmsLocOf / TauLocOf / Hampel3LocOf.                   |
//|  Avoid low-lag (EMA/Hull/ZLEMA - which you already have) on the  |
//|  midline: they cling to price and worsen the band.              |
//|                                                                  |
//|  NOTE: taken from a file that compiles; still COMPILE it in your |
//|  project. If any name collides with yours, rename it.            |
//+------------------------------------------------------------------+
#ifndef __ROBUST_LOCATION_ESTIMATORS_MQH__
#define __ROBUST_LOCATION_ESTIMATORS_MQH__

void InsSort(double &a[], int n)
{
   for(int x = 1; x < n; x++) { double v = a[x]; int y = x - 1; while(y >= 0 && a[y] > v) { a[y+1] = a[y]; y--; } a[y+1] = v; }
}
double MeanOf(const double &a[], int n)
{
   double s = 0.0; for(int k = 0; k < n; k++) s += a[k]; return (n > 0) ? s / n : 0.0;
}
double MedianOf(const double &a[], int n)
{
   double t[]; ArrayResize(t, n); for(int k = 0; k < n; k++) t[k] = a[k];
   InsSort(t, n);
   if(n % 2 == 1) return t[n/2];
   return 0.5 * (t[n/2 - 1] + t[n/2]);
}
double QsortedAt(const double &a[], int n, double q)
{
   if(n <= 1) return a[0];
   double pos = q * (n - 1); int lo = (int)MathFloor(pos); double fr = pos - lo;
   if(lo >= n - 1) return a[n-1]; if(lo < 0) return a[0];
   return a[lo] + fr * (a[lo+1] - a[lo]);
}
// Linear-interpolated quantile of an UNSORTED array (copies, sorts, interpolates).
double QuantileOf(const double &a[], int n, double q)
{
   if(n <= 0) return 0.0;
   if(n == 1) return a[0];
   double t[]; ArrayResize(t, n); for(int k = 0; k < n; k++) t[k] = a[k];
   InsSort(t, n);
   if(q <= 0.0) return t[0];
   if(q >= 1.0) return t[n-1];
   double pos = q * (n - 1); int lo = (int)MathFloor(pos); double fr = pos - lo;
   if(lo >= n - 1) return t[n-1];
   return t[lo] + fr * (t[lo+1] - t[lo]);
}
// Scaled MAD: 1.4826 * median(|x - median(x)|) -> robust dispersion estimate.
double MadScaleOf(const double &a[], int n)
{
   double med = MedianOf(a, n);
   double d[]; ArrayResize(d, n);
   for(int k = 0; k < n; k++) d[k] = MathAbs(a[k] - med);
   return MedianOf(d, n) * 1.4826;
}
// Population standard deviation.
double StdevPopOf(const double &a[], int n)
{
   double m = MeanOf(a, n), v = 0.0;
   for(int k = 0; k < n; k++) { double e = a[k] - m; v += e * e; }
   return MathSqrt(v / n);
}
// Tukey biweight location (robust mean; redescending, c = 4.685*MAD).
double BiweightOf(const double &a[], int n)
{
   double m = MedianOf(a, n), d[]; ArrayResize(d, n);
   for(int it = 0; it < 8; it++)
   {
      for(int j = 0; j < n; j++) d[j] = MathAbs(a[j] - m);
      double mad = MedianOf(d, n) * 1.4826; if(mad <= 0.0) break;
      double c = 4.685 * mad, num = 0.0, den = 0.0;
      for(int j = 0; j < n; j++)
      {
         double u = (a[j] - m) / c;
         if(MathAbs(u) < 1.0) { double ww = (1.0 - u*u); ww = ww*ww; num += ww * a[j]; den += ww; }
      }
      if(den <= 0.0) break;
      double mn = num / den; if(MathAbs(mn - m) < 1e-9) { m = mn; break; } m = mn;
   }
   return m;
}
// Huber M-estimator location (robust mean; k = 1.345*MAD).
double HuberOf(const double &a[], int n)
{
   double m = MedianOf(a, n), d[]; ArrayResize(d, n);
   for(int it = 0; it < 10; it++)
   {
      for(int j = 0; j < n; j++) d[j] = MathAbs(a[j] - m);
      double mad = MedianOf(d, n) * 1.4826; if(mad <= 0.0) break;
      double k = 1.345 * mad, num = 0.0, den = 0.0;
      for(int j = 0; j < n; j++)
      {
         double r = a[j] - m, ar = MathAbs(r);
         double ww = (ar <= k) ? 1.0 : (k / ar);
         num += ww * a[j]; den += ww;
      }
      if(den <= 0.0) break;
      double mn = num / den; if(MathAbs(mn - m) < 1e-9) { m = mn; break; } m = mn;
   }
   return m;
}
// 10% trimmed mean.
double TrimMeanOf(const double &a[], int n)
{
   double w[]; ArrayResize(w, n); for(int k = 0; k < n; k++) w[k] = a[k];
   InsSort(w, n);
   int k = (int)MathFloor(0.1 * n); double sum = 0.0; int cnt = 0;
   for(int j = k; j < n - k; j++) { sum += w[j]; cnt++; }
   return (cnt > 0) ? sum / cnt : MedianOf(a, n);
}
// Winsorized mean (clamp the g/1-g tails, then average).
double WinsMeanOf(const double &a[], int n, double g)
{
   double w[]; ArrayResize(w, n); for(int k = 0; k < n; k++) w[k] = a[k];
   InsSort(w, n);
   double gg = (g > 0.0 && g < 0.5) ? g : 0.1;
   double lo = QsortedAt(w, n, gg), hi = QsortedAt(w, n, 1.0 - gg);
   double sum = 0.0;
   for(int t = 0; t < n; t++) { double v = w[t]; if(v < lo) v = lo; if(v > hi) v = hi; sum += v; }
   return sum / n;
}
// Least Median of Squares location = midpoint of the shortest half (50% breakdown).
double LmsLocOf(const double &a[], int n)
{
   if(n <= 1) return (n == 1) ? a[0] : 0.0;
   double t[]; ArrayResize(t, n); for(int k = 0; k < n; k++) t[k] = a[k];
   InsSort(t, n);
   int h = n/2 + 1;
   if(h >= n) return MedianOf(a, n);
   int bi = 0; double bg = t[h-1] - t[0];
   for(int i = 1; i + h - 1 < n; i++)
   {
      double g = t[i + h - 1] - t[i];
      if(g < bg) { bg = g; bi = i; }
   }
   return 0.5 * (t[bi] + t[bi + h - 1]);
}
// Tukey biweight M-scale (~50% breakdown; c=1.547, b=0.5). Used by MM and Tau.
double SScaleOf(const double &a[], int n)
{
   double m = MedianOf(a, n);
   double s = MadScaleOf(a, n);
   if(s <= 0.0) { s = StdevPopOf(a, n); if(s <= 0.0) s = 1.0; }
   double c = 1.547, b = 0.5;
   for(int it = 0; it < 40; it++)
   {
      double sumRho = 0.0;
      for(int j = 0; j < n; j++)
      {
         double u = (a[j] - m) / (c * s);
         double rho = (MathAbs(u) < 1.0) ? (1.0 - MathPow(1.0 - u*u, 3.0)) : 1.0;
         sumRho += rho;
      }
      double sn = s * MathSqrt(MathMax(sumRho / n, 1e-9) / b);
      if(MathAbs(sn - s) < 1e-8 * MathMax(1.0, s)) { s = sn; break; }
      s = sn;
   }
   return MathMax(s, 1e-9);
}
// MM-estimator location: high-breakdown S-scale (fixed) + 95%-efficient redescending biweight M.
double MMLocOf(const double &a[], int n)
{
   double s = SScaleOf(a, n);
   double m = MedianOf(a, n);
   double c = 4.685;
   for(int it = 0; it < 25; it++)
   {
      double num = 0.0, den = 0.0;
      for(int j = 0; j < n; j++)
      {
         double u = (a[j] - m) / (c * s);
         if(MathAbs(u) < 1.0) { double ww = (1.0 - u*u); ww = ww*ww; num += ww * a[j]; den += ww; }
      }
      if(den <= 0.0) break;
      double mn = num / den;
      if(MathAbs(mn - m) < 1e-9 * MathMax(1.0, MathAbs(m))) { m = mn; break; }
      m = mn;
   }
   return m;
}
// Hampel three-part redescending M-estimator location (a,b,c = 1.7,3.4,8.5 * MAD).
double Hampel3LocOf(const double &a[], int n)
{
   double m = MedianOf(a, n);
   double A = 1.7, B = 3.4, C = 8.5;
   for(int it = 0; it < 25; it++)
   {
      double s = MadScaleOf(a, n);
      if(s <= 0.0) break;
      double num = 0.0, den = 0.0;
      for(int j = 0; j < n; j++)
      {
         double r  = (a[j] - m) / s;
         double ar = MathAbs(r);
         double sg = (r >= 0.0) ? 1.0 : -1.0;
         double psi;
         if(ar <= A)      psi = r;
         else if(ar <= B) psi = A * sg;
         else if(ar <= C) psi = A * sg * (C - ar) / (C - B);
         else             psi = 0.0;
         double ww = (ar > 1e-12) ? (psi / r) : 1.0;
         num += ww * a[j]; den += ww;
      }
      if(den <= 0.0) break;
      double mn = num / den;
      if(MathAbs(mn - m) < 1e-9 * MathMax(1.0, MathAbs(m))) { m = mn; break; }
      m = mn;
   }
   return m;
}
// Data-driven adaptive trimmed mean: pick the trim fraction that minimizes the
// estimated variance of the trimmed mean over a small grid (Jaeckel-style).
double AdTrimMeanOf(const double &a[], int n)
{
   double t[]; ArrayResize(t, n); for(int k = 0; k < n; k++) t[k] = a[k];
   InsSort(t, n);
   double grid[6]; grid[0]=0.0; grid[1]=0.05; grid[2]=0.10; grid[3]=0.15; grid[4]=0.20; grid[5]=0.25;
   double best = MedianOf(a, n), bestV = 1e300;
   double wv[]; ArrayResize(wv, n);
   for(int gi = 0; gi < 6; gi++)
   {
      double al = grid[gi];
      int g = (int)MathFloor(al * n);
      if(n - 2*g < 1) continue;
      double mu = 0.0; int cnt = 0;
      for(int j = g; j < n - g; j++) { mu += t[j]; cnt++; }
      mu /= cnt;
      double wmean = 0.0;
      for(int j = 0; j < n; j++)
      {
         double v = t[j];
         if(g > 0) { if(j < g) v = t[g]; else if(j >= n - g) v = t[n - g - 1]; }
         wv[j] = v; wmean += v;
      }
      wmean /= n;
      double sw2 = 0.0;
      for(int j = 0; j < n; j++) { double e = wv[j] - wmean; sw2 += e*e; }
      sw2 /= (n - 1);
      double dd = (1.0 - 2.0*al); dd = dd*dd;
      double vv = sw2 / dd / n;
      if(vv < bestV) { bestV = vv; best = mu; }
   }
   return best;
}
// Tau-estimator location: high-breakdown S-scale + one redescending (biweight) reweight.
double TauLocOf(const double &a[], int n)
{
   double m = MedianOf(a, n);
   double s = SScaleOf(a, n);
   if(s <= 0.0) return m;
   double c1 = 4.5, num = 0.0, den = 0.0;
   for(int j = 0; j < n; j++)
   {
      double u = (a[j] - m) / (c1 * s);
      if(MathAbs(u) < 1.0) { double ww = (1.0 - u*u); ww = ww*ww; num += ww * a[j]; den += ww; }
   }
   return (den > 0.0) ? (num / den) : m;
}
// Hodges-Lehmann location = median of all pairwise (Walsh) averages. Smoother than
// the plain median (~96% efficiency) yet robust (~29% breakdown): a steady pivot
// for placing scale-in ladders. Used by the CENTER_HL center option.
double HodgesLehmannOf(const double &a[], int n)
{
   if(n <= 1) return (n == 1) ? a[0] : 0.0;
   int m = n * (n + 1) / 2;
   double w[]; ArrayResize(w, m);
   int idx = 0;
   for(int i = 0; i < n; i++)
      for(int j = i; j < n; j++)
         w[idx++] = 0.5 * (a[i] + a[j]);
   return MedianOf(w, m);
}
// True if a same-direction arrow exists within the last InpArrowMargin bars
// (series indexing: higher index = older). Suppresses repeat arrows inside the

//==================== EXTRAS (not in the original library) ====================
// Midhinge = (Q1 + Q3)/2. Robust and cheap.
double MidhingeOf(const double &a[], int n)
{
   if(n <= 0) return 0.0;
   return 0.5 * (QuantileOf(a, n, 0.25) + QuantileOf(a, n, 0.75));
}
// Tukey trimean = (Q1 + 2*Q2 + Q3)/4. Robust.
double TrimeanOf(const double &a[], int n)
{
   if(n <= 0) return 0.0;
   return 0.25 * (QuantileOf(a, n, 0.25) + 2.0 * QuantileOf(a, n, 0.50) + QuantileOf(a, n, 0.75));
}

#endif // __ROBUST_LOCATION_ESTIMATORS_MQH__