#ifndef STATISTICS_MQH
#define STATISTICS_MQH

// Precisione macchina double
#define DBL_EPS 2.2204460492503131e-016
// Soglia numerica data-scaled: |x| * DBL_EPS * 1e9 = |x| * 2.22e-7 ≈ |x| * 1e-6
// Epsilon data-driven: sqrt(machine_epsilon) × |x|, mai sotto sqrt(machine_epsilon)
// sqrt(DBL_EPS) è la soglia standard per confronti floating-point (Num. Recipes)
#define DATA_EPS(x) MathMax(MathSqrt(DBL_EPS), MathAbs(x) * MathSqrt(DBL_EPS))

class EWMA {
   double value;
   double alpha;
   bool   init;
public:
   EWMA(double a=0.1) : alpha(a), init(false), value(0) {}
   
   void SetAlpha(double a) { alpha = a; }
   double Alpha() const { return alpha; }
   
   double Update(double x) {
      if(!init) { value = x; init = true; }
      else value = alpha * x + (1.0 - alpha) * value;
      return value;
   }
   
   double Value() const { return init ? value : 0; }
   bool   IsInit() const { return init; }
   void   Reset() { init = false; value = 0; }
   
   void Save(int fh) const {
      FileWriteDouble(fh, value);
      FileWriteInteger(fh, init ? 1 : 0);
   }
   void Load(int fh) {
      value = FileReadDouble(fh);
      init  = FileReadInteger(fh) == 1;
   }
};

class RunningStats {
   EWMA  mean;
   EWMA  var;
   int   minSamples;
   int   count;
   int   freezeAfter;
   
   double AdaptiveEpsilon() const {
      double m = MathAbs(mean.Value());
      return (m > 0) ? DATA_EPS(m) : 1e-15;
   }
public:
   RunningStats(double a=0.05, int minS=20, int freeze=0)
      : mean(a), var(a), minSamples(minS), count(0), freezeAfter(freeze) {}
   
   void SetFreezeAfter(int n) { freezeAfter = n; }
   bool IsFrozen() const { return freezeAfter > 0 && count >= freezeAfter; }
   
   void Update(double x) {
      if(IsFrozen()) return;
      if(count == 0) {
         mean.Update(x);
         var.Update(0);
         count = 1;
         return;
      }
      double prevMean = mean.Value();
      mean.Update(x);
      double diff = x - prevMean;
      var.Update(diff * diff);
      count++;
   }
   
   double ZScore(double x) {
      double m = mean.Value();
      double s = MathSqrt(var.Value());
      if(s < AdaptiveEpsilon() || count < minSamples) return 0;
      return (x - m) / s;
   }
   
   double RawZScore(double x) {
      double m = mean.Value();
      double s = MathSqrt(var.Value());
      if(s < AdaptiveEpsilon()) return 0;
      return (x - m) / s;
   }
   
   double Mean() const { return mean.Value(); }
   double Std()  const { return MathSqrt(var.Value()); }
   bool   Ready() const { return count >= minSamples; }
   void   Reset() { mean.Reset(); var.Reset(); count = 0; }
   int    Count() const { return count; }
   
   void Save(int fh) const {
      mean.Save(fh);
      var.Save(fh);
      FileWriteInteger(fh, count);
      FileWriteInteger(fh, freezeAfter);
   }
   void Load(int fh) {
      mean.Load(fh);
      var.Load(fh);
      count = FileReadInteger(fh);
      freezeAfter = FileReadInteger(fh);
   }
   
   string ToString() const {
      string s = "mean=" + StringFormat("%.5f", Mean())
               + " std=" + StringFormat("%.5f", Std())
               + " n=" + (string)count + "/" + (string)minSamples;
      if(IsFrozen()) s += " FROZEN";
      return s;
   }
};

class RunningCorrelation {
   double alpha;
   double meanX, meanY;
   double cov, varX, varY;
   int    count;
   int    minSamples;
   
   double AdaptiveEpsilon() const {
      double mx = MathAbs(meanX), my = MathAbs(meanY);
      double ref = (mx + my) * 0.5;
      return (ref > 0) ? DATA_EPS(ref) : 1e-15;
   }
public:
   RunningCorrelation(double a=0.05, int minS=10) 
      : alpha(a), minSamples(minS), count(0),
        meanX(0), meanY(0), cov(0), varX(0), varY(0) {}
   
   void Update(double x, double y) {
      count++;
      if(count == 1) {
         meanX = x; meanY = y;
         return;
      }
      
      double dx = x - meanX;
      double dy = y - meanY;
      
      meanX += alpha * dx;
      meanY += alpha * dy;
      
      double dxNew = x - meanX;
      double dyNew = y - meanY;
      
      cov  = (1.0 - alpha) * cov  + alpha * dx * dyNew;
      varX = (1.0 - alpha) * varX + alpha * dx * dxNew;
      varY = (1.0 - alpha) * varY + alpha * dy * dyNew;
   }
   
   double Correlation() {
      double denom = MathSqrt(varX * varY);
      if(denom < AdaptiveEpsilon() || count < minSamples) return 0;
      double r = cov / denom;
      double maxObserved = 1.0;
      return MathMax(-maxObserved, MathMin(maxObserved, r));
   }
   
    bool   Ready() const { return count >= minSamples; }
    void   Reset() { count = 0; meanX = meanY = cov = varX = varY = 0; }
    int    Count() const { return count; }

    void Save(int fh) const {
       FileWriteDouble(fh, meanX);
       FileWriteDouble(fh, meanY);
       FileWriteDouble(fh, cov);
       FileWriteDouble(fh, varX);
       FileWriteDouble(fh, varY);
       FileWriteInteger(fh, count);
    }
    void Load(int fh) {
       meanX = FileReadDouble(fh);
       meanY = FileReadDouble(fh);
       cov   = FileReadDouble(fh);
       varX  = FileReadDouble(fh);
       varY  = FileReadDouble(fh);
       count = FileReadInteger(fh);
    }
};

// --- Kalman Filter Normalizer con Q adattivo ---
class KalmanNormalizer {
   double x;
   double P;
   double Q;
   double R;
   int    count;
   int    minSamples;
   
   double AdaptiveEpsilon() const {
      double ax = MathAbs(x);
      return (ax > 0) ? DATA_EPS(ax) : 1e-15;
   }
public:
   KalmanNormalizer(double q=0.001, double r=0.1, int minS=20)
      : x(0), P(1.0), Q(q), R(r), count(0), minSamples(minS) {}
   
   void SetQ(double q) { Q = q; }
   void SetR(double r) { R = r; }
   
   void Update(double obs) {
      if(count == 0) {
         x = obs;
         P = R;
         count = 1;
         return;
      }
      double innov = obs - x;
      
      P += Q;
      
      double K = P / (P + R);
      x += K * innov;
      P = (1.0 - K) * P;
      
      double innovVar = innov * innov;
      // AdaptRate: innovVar/(R+innovVar) → [0, 1), nessun clamp
      double adaptRate = innovVar / (R + innovVar);
      Q = (1.0 - adaptRate) * Q + adaptRate * innovVar;
      double qMin = DATA_EPS(R);
      double qMax = R - DATA_EPS(R);  // Q < R garantito → K < 0.5
      Q = MathMax(qMin, MathMin(qMax, Q));
      
      count++;
   }
   
   double Mean() const { return x; }
   double Std()  const { return MathSqrt(P); }
   bool   Ready() const { return count >= minSamples; }
   int    Count() const { return count; }
   void   Reset() { x = 0; P = 1.0; count = 0; }
   
   double ZScore(double obs) {
      double s = Std();
      if(s < AdaptiveEpsilon() || count < minSamples) return 0;
      return (obs - x) / s;
   }
   
   double RawZScore(double obs) {
      double s = Std();
      if(s < AdaptiveEpsilon()) return 0;
      return (obs - x) / s;
   }
   
   void Save(int fh) const {
      FileWriteDouble(fh, x);
      FileWriteDouble(fh, P);
      FileWriteInteger(fh, count);
   }
   void Load(int fh) {
      x = FileReadDouble(fh);
      P = FileReadDouble(fh);
      count = FileReadInteger(fh);
   }
   
   string ToString() const {
      return "μ=" + StringFormat("%.5f", x)
           + " σ=" + StringFormat("%.5f", Std())
           + " n=" + (string)count + "/" + (string)minSamples;
   }
};
#endif