#ifndef REGIME_DETECTOR_MQH
#define REGIME_DETECTOR_MQH
#include "AgentBase.mqh"
#include "../Core/PeriodCalculator.mqh"

class RegimeDetector : public IAgent {
private:
   int    hurstPeriod;
   int    userPeriod;
   int    minPeriod, maxPeriod;
   double prevZ;
   int    warmup;
   int    targetWindows;
   
    int LogReturns(const double &close[], int len, double &ret[]) const {
       int n = len - 1;
       ArrayResize(ret, n);
       for(int i=0; i<n; i++) {
          double r = close[i] / close[i+1];
          if(r <= 0) { ret[i] = 0; continue; }
          ret[i] = MathLog(r);
       }
       return n;
    }
   
     double ComputeDFA(const double &close[], int len) {
        targetWindows = MathMax(4, MathMin(14, len / 50));
       int minN = MathMax(3, targetWindows);
        if(len < targetWindows * minN * 2) return 0.5;
        
        int maxN = MathMax(minN * 2, len / (targetWindows / 2));
        if(maxN < minN * 2) return 0.5;
        
        double returns[];
        int nRet = LogReturns(close, len, returns);
        if(nRet < maxN) return 0.5;
        
         // Integra: profilo (somma cumulativa dei rendimenti)
         double profile[];
         ArrayResize(profile, nRet);
         profile[0] = returns[0];
         for(int i=1; i<nRet; i++)
            profile[i] = profile[i-1] + returns[i];
         
         // Varianza dei rendimenti per soglia data-scaled DFA
         double retVar = 0;
         for(int i=0; i<nRet; i++) retVar += returns[i] * returns[i];
          retVar = MathMax(1e-15, retVar / nRet);
         
         // Step derivato dal numero di window target
         double step = MathPow((double)maxN / minN, 1.0 / (targetWindows - 1));
         step = MathMax(1.3, MathMin(2.0, step));
        
        double logF[], logN[];
        ArrayResize(logF, targetWindows);
        ArrayResize(logN, targetWindows);
        int pts = 0;
        
        for(int n = minN; n <= maxN; n = (int)(n * step) + 1) {
           int m = nRet / n;
           if(m < 3) continue;
           if(pts >= targetWindows) break;
          
          double sumF2 = 0;
          int validWin = 0;
          
          for(int j=0; j<m; j++) {
             int base = j * n;
             
             // OLS detrend lineare della finestra
             double sx=0, sy=0, sxx=0, sxy=0;
             for(int k=0; k<n; k++) {
                double x = k;
                double y = profile[base + k];
                sx  += x;   sy  += y;
                sxx += x*x; sxy += x*y;
             }
              double slope = (n * sxy - sx * sy) / (n * sxx - sx * sx + 1e-15);
              double intercept = (sy - slope * sx) / n;
              
              // Varianza del residuo (dopo detrend)
              double var = 0;
              for(int k=0; k<n; k++) {
                 double fit = intercept + slope * k;
                 double res = profile[base + k] - fit;
                 var += res * res;
              }
              var /= n;
               // Soglia: varianza attesa per unbiased RW = retVar * n
               // DATA_EPS: soglia numerica scalata con la varianza attesa
               double epsVar = DATA_EPS(retVar * n);
              if(var < epsVar) continue;
              sumF2 += var;
              validWin++;
          }
          
          if(validWin < 2) continue;
          double F = MathSqrt(sumF2 / validWin);
          logF[pts] = MathLog(F);
          logN[pts] = MathLog(n);
          pts++;
       }
       
       if(pts < 3) return 0.5;
       
       double sumX=0, sumY=0, sumXY=0, sumX2=0;
       for(int i=0; i<pts; i++) {
          sumX  += logN[i];
          sumY  += logF[i];
          sumXY += logN[i] * logF[i];
          sumX2 += logN[i] * logN[i];
       }
       double H = (pts * sumXY - sumX * sumY) / (pts * sumX2 - sumX * sumX);
       
       H = MathMax(0.01, MathMin(1.50, H));
       return H;
    }
    
    // Fallback a R/S se DFA non converge
    double ComputeHurst(const double &close[], int len) {
      double H = ComputeDFA(close, len);
      double hSe = MathSqrt(12.0 / len);  // SE approssimato di Hurst per unbiased RW
      if(H < hSe || H > 1.0 - hSe || MathAbs(H - 0.5) < hSe)
         H = ComputeRS(close, len);
      return MathMax(hSe, MathMin(1.0 - hSe, H));
    }
    
     double ComputeRS(const double &close[], int len) {
        targetWindows = MathMax(3, MathMin(10, len / 60));
       int minN = MathMax(3, targetWindows);
        if(len < targetWindows * minN * 2) return 0.5;
        
        int maxN = MathMax(minN * 2, len / (targetWindows / 2));
        if(maxN < minN * 2) return 0.5;
        
         double returns[];
         int nRet = LogReturns(close, len, returns);
         if(nRet < maxN) return 0.5;
         
         // Varianza di riferimento per soglia data-scaled
          double retVarRef = 0;
          for(int i=0; i<nRet; i++) retVarRef += returns[i] * returns[i];
           retVarRef = DATA_EPS(retVarRef / nRet);
          double epsVarRS = DATA_EPS(retVarRef);
         
         double step = MathPow((double)maxN / minN, 1.0 / (targetWindows - 1));
         step = MathMax(1.3, MathMin(2.5, step));
        
        double logRS[], logN[];
        ArrayResize(logRS, targetWindows);
        ArrayResize(logN, targetWindows);
        int pts = 0;
        
        for(int n = minN; n <= maxN; n = (int)(n * step) + 1) {
           int m = nRet / n;
           if(m < 2) continue;
           if(pts >= targetWindows) break;
          
          double sumRS = 0;
          int validSub = 0;
          
          for(int j=0; j<m; j++) {
             int base = j * n;
             
             double sum = 0, sumSq = 0;
             for(int k=0; k<n; k++) {
                double r = returns[base + k];
                sum += r;
                sumSq += r * r;
             }
              double mean = sum / n;
              double var = sumSq / n - mean * mean;
              double std = (var > epsVarRS) ? MathSqrt(var) : 0;
              if(std < MathSqrt(epsVarRS)) continue;
             
             double cumDev[];
             ArrayResize(cumDev, n);
             cumDev[0] = returns[base] - mean;
             for(int k=1; k<n; k++)
                cumDev[k] = cumDev[k-1] + returns[base + k] - mean;
             
             int maxIdx = 0, minIdx = 0;
             for(int k=1; k<n; k++) {
                if(cumDev[k] > cumDev[maxIdx]) maxIdx = k;
                if(cumDev[k] < cumDev[minIdx]) minIdx = k;
             }
             double R = cumDev[maxIdx] - cumDev[minIdx];
             sumRS += R / std;
             validSub++;
          }
          
          if(validSub < 1) continue;
          double avgRS = sumRS / validSub;
          logRS[pts] = MathLog(avgRS);
          logN[pts]  = MathLog(n);
          pts++;
       }
       
       if(pts < 3) return 0.5;
       
       double sumX=0, sumY=0, sumXY=0, sumX2=0;
       for(int i=0; i<pts; i++) {
          sumX  += logN[i];
          sumY  += logRS[i];
          sumXY += logN[i] * logRS[i];
          sumX2 += logN[i] * logN[i];
       }
       double H = (pts * sumXY - sumX * sumY) / (pts * sumX2 - sumX * sumX);
       return MathMax(0.01, MathMin(0.99, H));
    }
   
public:
   RegimeDetector(string n="Hurst", double w=1.0, int hp=0)
      : IAgent(n, w), userPeriod(hp), hurstPeriod(0), minPeriod(40), maxPeriod(200), prevZ(0), warmup(0) { signalStats.SetR(0.05); }
   
    double Analyze(const MarketData &data) override {
       warmup++;
       
       // Periodo: fisso se utente lo specifica, altrimenti data-driven + EWMA
       if(userPeriod > 0) {
          hurstPeriod = userPeriod;
       } else {
           int cycle = PeriodCalculator::DominantCycle(data.close, data.count, 20, 100);
           // Periodo: max(2x ciclo, minWindows * targetWindows)
           int minForWindows = targetWindows * MathMax(3, targetWindows);
           int newP = MathMax(cycle * 2, minForWindows);
           newP = MathMax(20, MathMin(200, newP));
           if(hurstPeriod <= 0) hurstPeriod = newP;
           else {
               double alpha = 1.0 / (1.0 + warmup * 0.1);
               double minAlpha = 1.0 / MathMax(2.0, (double)MathMax(1, hurstPeriod));
               alpha = MathMax(minAlpha, alpha);  // solo floor, niente max clamp
              hurstPeriod = (int)MathRound(alpha * newP + (1.0 - alpha) * hurstPeriod);
           }
          if(hurstPeriod < minPeriod) hurstPeriod = minPeriod;
       }
       
       double H = ComputeHurst(data.close, MathMin(hurstPeriod, data.count));
       
       signalStats.Update(H);
       double zRaw = signalStats.ZScore(H);
       
        // EWMA con alpha che scala con il numero di osservazioni
        double alpha = 1.0 / (1.0 + signalStats.Count() * 0.1);
        double minAlpha = 1.0 / MathMax(2.0, (double)MathMax(1, hurstPeriod));
        alpha = MathMax(minAlpha, alpha);  // solo floor, niente max clamp
       prevZ = (1.0 - alpha) * prevZ + alpha * zRaw;
       double calibrated = CalibrateZ(prevZ);
       lastZScore = MathTanh(calibrated);
       lastRawSignal = H;
      
      SHARED_regimeH = H;
      
      return lastZScore;
   }
   
   void Interact(IAgent *&allAgents[], int count) override {}
   
   void Learn(double predictedZ, double actualReturnZ) override {}
   
   void Save(int fh) const override {
      IAgent::Save(fh);
      FileWriteDouble(fh, prevZ);
   }
    void Load(int fh) override {
       IAgent::Load(fh);
       prevZ = FileReadDouble(fh);
       warmup = signalStats.Count();  // ripristina warmup dal conteggio statistiche
    }
   
   void Reset() override {
      IAgent::Reset();
      prevZ = 0;
      warmup = 0;
   }
   
    string SignalInfo() const override {
       return name + " z=" + StringFormat("%+.3f", lastZScore)
            + " H=" + StringFormat("%.3f", SHARED_regimeH)
            + " p=" + (string)hurstPeriod
            + " " + signalStats.ToString();
     }
};
#endif