//+------------------------------------------------------------------+
//|                                            BOSSvsDTWBenchmark.mq5|
//|                                      MMQ — Muhammad Minhas Qamar |
//|                                   www.mql5.com/en/articles/23491 |
//+------------------------------------------------------------------+
#property copyright   "MMQ — Muhammad Minhas Qamar"
#property link        "https://www.mql5.com/en/articles/23491"
#property description "Head-to-head on real bars: dictionary-based BOSS versus elastic DTW, both as 1-nearest-neighbour classifiers of the SAME regime labels on the SAME data."
#property version     "1.00"
#property strict
#property script_show_inputs

#include <BOSS\BOSS.mqh>
#include <BOSS\RegimeLabeler.mqh>
#include <dtw.mqh>

input string Symbol_    = "BTCUSD";   // symbol to pull (falls back to chart symbol)
input ENUM_TIMEFRAMES TF = PERIOD_H1; // timeframe
input int    Bars_      = 50000;      // history bars to load
input int    Segment    = 120;        // bars per classified segment
input int    Window     = 24;         // SFA sliding-window length (inside a segment)
input int    WordLen    = 2;          // SFA coefficients kept (robust low-pass point)
input int    Alphabet   = 4;          // SFA letters per coefficient slot
input double NoiseSigma = 0.4;        // stdev of injected noise (in window-stdev units)
//--- labeler thresholds (tune these to balance the regime classes):
input double TrendACF   = 0.05;       // lag-1 autocorrelation above this -> TREND
input double VolHighMult= 1.3;        // volatility above mult*reference -> VOLATILE

//--- forward declarations (functions call each other out of order).
void   RunBenchmark(const CBossClassifier &boss,
                    CBossEnsemble &ens,bool ensOK,
                    const double &segs[],const int &labels[],
                    int offset,int ntest,int seg,
                    const double &protos[],const int &plab[],int nproto,
                    double noise);
double Recall(int hit,int tot);

//+------------------------------------------------------------------+
//| Load 'count' closes for the chosen symbol/timeframe into 'out',  |
//| oldest first. Returns the number actually loaded.                |
//+------------------------------------------------------------------+
int LoadCloses(double &out[])
  {
   string sym=Symbol_;
   if(sym=="" || !SymbolSelect(sym,true))
      sym=_Symbol;
   double c[];
   int got=CopyClose(sym,TF,1,Bars_,c);
   if(got<=0)
     {
      PrintFormat("CopyClose failed for %s: %d",sym,GetLastError());
      return 0;
     }
   ArraySetAsSeries(c,false);           // oldest first
   ArrayResize(out,got);
   for(int i=0;i<got;i++)
      out[i]=c[i];
   PrintFormat("loaded %d closes for %s %s",got,sym,EnumToString(TF));
   return got;
  }

//+------------------------------------------------------------------+
//| Copy a length-'len' slice out of 'series' starting at 'start'.   |
//+------------------------------------------------------------------+
void Slice(const double &series[],int start,int len,double &out[])
  {
   ArrayResize(out,len);
   for(int i=0;i<len;i++)
      out[i]=series[start+i];
  }

//+------------------------------------------------------------------+
//| Add Gaussian noise scaled to the segment's own return spread, so |
//| 'sigma' is a relative noise level comparable across regimes.     |
//| Uses a 12-uniform sum as a cheap standard-normal draw.           |
//+------------------------------------------------------------------+
void AddNoise(double &seg[],int len,double sigma)
  {
   double mean=0.0;
   for(int i=0;i<len;i++)
      mean+=seg[i];
   mean/=len;
   double var=0.0;
   for(int i=0;i<len;i++)
     {
      double d=seg[i]-mean;
      var+=d*d;
     }
   double sd=MathSqrt(var/len);
   double amp=sigma*sd;

   for(int i=0;i<len;i++)
     {
      double g=0.0;
      for(int k=0;k<12;k++)
         g+=(MathRand()/32767.0);
      g-=6.0;                            // ~N(0,1)
      seg[i]+=amp*g;
     }
  }

//+------------------------------------------------------------------+
//| Z-normalize src[0..len-1] into a vector, for DTW shape matching. |
//| Without this DTW would match on absolute price level, which is   |
//| meaningless on an instrument whose price ranges over tens of     |
//| thousands - so this keeps BOSS and DTW on equal (shape) footing. |
//+------------------------------------------------------------------+
void ZNormVec(const double &src[],int len,vector &dst)
  {
   double mean=0.0;
   for(int i=0;i<len;i++)
      mean+=src[i];
   mean/=len;
   double var=0.0;
   for(int i=0;i<len;i++)
     {
      double d=src[i]-mean;
      var+=d*d;
     }
   double sd=MathSqrt(var/len);
   if(sd<1.0e-12)
      sd=1.0;
   dst.Resize(len);
   for(int i=0;i<len;i++)
      dst[i]=(src[i]-mean)/sd;
  }

//+------------------------------------------------------------------+
//| DTW 1-NN over prototype segments.                                |
//|  Classifies query segment 'q' (length seg) by the label of the   |
//|  nearest prototype segment under DTW distance. Both are          |
//|  z-normalized first. Prototypes are the flat array 'protos'      |
//|  (nproto segments of 'seg'), labels in 'plab'.                   |
//+------------------------------------------------------------------+
int DTWClassify(const double &q[],int seg,const double &protos[],
                const int &plab[],int nproto)
  {
   vector vq;
   ZNormVec(q,seg,vq);

   Cdtw    engine;
   double  best=DBL_MAX;
   int     bestlab=REGIME_UNDEFINED;
   double  pbuf[];
   ArrayResize(pbuf,seg);
   vector  vp;

   for(int p=0;p<nproto;p++)
     {
      for(int i=0;i<seg;i++)
         pbuf[i]=protos[p*seg+i];
      ZNormVec(pbuf,seg,vp);
      if(!engine.dtw(vq,vp,DIST_EUCLIDEAN,STEP_SYMM2))
         continue;
      double d=engine.distance();
      if(d<best)
        {
         best=d;
         bestlab=plab[p];
        }
     }
   return bestlab;
  }

//+------------------------------------------------------------------+
//| Script entry: load bars, cut into segments, label, split,        |
//| train, benchmark.                                                |
//+------------------------------------------------------------------+
void OnStart()
  {
   MathSrand(20260712);

   double closes[];
   int n=LoadCloses(closes);
   if(n<4*Segment)
     { Print("not enough bars for the chosen segment size"); return; }
   if(Segment<=Window+4)
     { Print("Segment must be larger than Window"); return; }

//--- STEP 1: cut non-overlapping segments and label each by the
//--- regime of the whole segment.
   int nseg=n/Segment;
   double  segbuf[];
   double  segs[];
   ArrayResize(segs,nseg*Segment);
   int     labels[];
   ArrayResize(labels,nseg);

//--- reference volatility for the labeler: mean segment volatility.
   CRegimeLabeler labeler;
   labeler.SetTrendACF(TrendACF);
   labeler.SetVolHighMult(VolHighMult);
   double volsum=0.0;
   for(int s=0;s<nseg;s++)
     {
      Slice(closes,s*Segment,Segment,segbuf);
      volsum+=labeler.WindowVolatility(segbuf,Segment);
     }
   labeler.SetReferenceVol((nseg>0)?volsum/nseg:0.0);

   int kept=0, cnt[4]= {0,0,0,0};
   for(int s=0;s<nseg;s++)
     {
      Slice(closes,s*Segment,Segment,segbuf);
      ENUM_REGIME r=labeler.Label(segbuf,Segment);
      if(r==REGIME_UNDEFINED)
         continue;
      for(int i=0;i<Segment;i++)
         segs[kept*Segment+i]=segbuf[i];
      labels[kept]=(int)r;
      cnt[(int)r]++;
      kept++;
     }
   ArrayResize(segs,kept*Segment);
   ArrayResize(labels,kept);
   PrintFormat("cut %d segments of %d bars -> trend=%d range=%d volatile=%d",
               kept,Segment,cnt[REGIME_TREND],cnt[REGIME_RANGE],cnt[REGIME_VOLATILE]);
   if(kept<8)
     { Print("too few labelled segments; load more bars or shrink Segment"); return; }

//--- STEP 2: split into train (first half) and test (second half).
   int split=kept/2;
   int ntrain=split, ntest=kept-split;

//--- prototype segments for DTW (the train half), flat + labels.
   double protos[];
   int plab[];
   ArrayResize(protos,ntrain*Segment);
   ArrayResize(plab,ntrain);
   for(int i=0;i<ntrain;i++)
     {
      for(int j=0;j<Segment;j++)
         protos[i*Segment+j]=segs[i*Segment+j];
      plab[i]=labels[i];
     }

//--- STEP 3: harvest the MCB training pool. Every sliding window
//--- inside every training segment feeds the breakpoint estimate,
//--- so the bins see the true spread of coefficients BOSS will meet.
   int perSeg=Segment-Window+1;
   int poolN=ntrain*perSeg;
   double pool[];
   ArrayResize(pool,poolN*Window);
   int pk=0;
   for(int i=0;i<ntrain;i++)
      for(int start=0;start<=Segment-Window;start++)
        {
         for(int j=0;j<Window;j++)
            pool[pk*Window+j]=protos[i*Segment+start+j];
         pk++;
        }

//--- STEP 4: train BOSS. MCB fit on the window pool; each training
//--- SEGMENT added as one labelled bag (many words, the real thing).
   CBossClassifier boss;
   if(!boss.Build(WordLen,Alphabet,pool,poolN,Window))
     { Print("BOSS build failed"); return; }
   for(int i=0;i<ntrain;i++)
     {
      double ex[];
      Slice(protos,i*Segment,Segment,ex);
      boss.AddExample(ex,Segment,CRegimeLabeler::RegimeName((ENUM_REGIME)plab[i]));
     }
   PrintFormat("trained single-BOSS: %d segments, MCB pool=%d windows, word=%d letters",
               ntrain,poolN,WordLen*2);

//--- STEP 4b: train the BOSS ENSEMBLE over several window sizes.
//--- The ensemble is handed the raw training segments and builds
//--- each member's own pool/bags internally, keeps members within
//--- 92% of the best LOO accuracy, and majority-votes at classify.
   int candWins[]= {12,18,24,32,48,64};
   string classNames[3];
   classNames[REGIME_TREND]   =CRegimeLabeler::RegimeName(REGIME_TREND);
   classNames[REGIME_RANGE]   =CRegimeLabeler::RegimeName(REGIME_RANGE);
   classNames[REGIME_VOLATILE]=CRegimeLabeler::RegimeName(REGIME_VOLATILE);

   CBossEnsemble ens;
   bool ensOK=ens.Train(protos,plab,ntrain,Segment,WordLen,Alphabet,
                        candWins,ArraySize(candWins),classNames);
   if(ensOK)
     {
      string kept="";
      for(int i=0;i<ens.MemberCount();i++)
         kept+=StringFormat("w%d(%.0f%%) ",ens.MemberWindow(i),100.0*ens.MemberAcc(i));
      PrintFormat("trained ensemble: %d members kept -> %s",ens.MemberCount(),kept);
     }
   else
      Print("ensemble training failed");

//--- STEP 5: benchmark on the test half, clean then noisy.
   Print("=== clean queries ===");
   RunBenchmark(boss,ens,ensOK,segs,labels,split,ntest,Segment,protos,plab,ntrain,0.0);

   PrintFormat("=== noisy queries (sigma=%.2f) ===",NoiseSigma);
   RunBenchmark(boss,ens,ensOK,segs,labels,split,ntest,Segment,protos,plab,ntrain,NoiseSigma);

   Print("=== benchmark complete ===");
  }

//+------------------------------------------------------------------+
//| Run BOSS and DTW over the test segments and print accuracy,      |
//| wall-clock, a majority-class baseline and per-class recall.      |
//| 'offset' is where the test half begins inside segs[]. 'noise'>0  |
//| injects noise into each query segment before classifying.        |
//+------------------------------------------------------------------+
void RunBenchmark(const CBossClassifier &boss,
                  CBossEnsemble &ens,bool ensOK,
                  const double &segs[],const int &labels[],
                  int offset,int ntest,int seg,
                  const double &protos[],const int &plab[],int nproto,
                  double noise)
  {
   int    bossCorrect=0, ensCorrect=0, dtwCorrect=0;
   ulong  bossUs=0, ensUs=0, dtwUs=0;
   double q[];

   int perTot[4]= {0,0,0,0};
   int bossHit[4]= {0,0,0,0};
   int ensHit[4]= {0,0,0,0};
   int dtwHit[4]= {0,0,0,0};

   for(int t=0;t<ntest;t++)
     {
      int idx=offset+t;
      Slice(segs,idx*seg,seg,q);
      if(noise>0.0)
         AddNoise(q,seg,noise);
      int truth=labels[idx];
      string truthName=CRegimeLabeler::RegimeName((ENUM_REGIME)truth);
      perTot[truth]++;

      //--- single BOSS: 1-NN over word-frequency bags
      ulong t0=GetMicrosecondCount();
      double bd;
      string blab=boss.Classify(q,seg,bd);
      bossUs+=GetMicrosecondCount()-t0;
      if(blab==truthName)
        { bossCorrect++; bossHit[truth]++; }

      //--- BOSS ensemble: majority vote over kept window sizes
      if(ensOK)
        {
         ulong te=GetMicrosecondCount();
         string elab=ens.Classify(q,seg);
         ensUs+=GetMicrosecondCount()-te;
         if(elab==truthName)
           { ensCorrect++; ensHit[truth]++; }
        }

      //--- DTW 1-NN over the same prototype segments
      ulong t1=GetMicrosecondCount();
      int dlab=DTWClassify(q,seg,protos,plab,nproto);
      dtwUs+=GetMicrosecondCount()-t1;
      if(dlab==truth)
        { dtwCorrect++; dtwHit[truth]++; }
     }

//--- majority-class baseline: always guess the most common test label
   int major=0;
   for(int c=1;c<3;c++)
      if(perTot[c]>perTot[major])
         major=c;
   double baseline=(ntest>0)?100.0*perTot[major]/ntest:0.0;

//--- balanced (macro) accuracy: mean of per-class recalls
   double bMacro=0.0, eMacro=0.0, dMacro=0.0;
   int seen=0;
   for(int c=0;c<3;c++)
      if(perTot[c]>0)
        {
         bMacro+=100.0*bossHit[c]/perTot[c];
         eMacro+=100.0*ensHit[c]/perTot[c];
         dMacro+=100.0*dtwHit[c]/perTot[c];
         seen++;
        }
   if(seen>0)
     {
      bMacro/=seen;
      eMacro/=seen;
      dMacro/=seen;
     }

   PrintFormat("  baseline (always '%s'): %.1f%%",
               CRegimeLabeler::RegimeName((ENUM_REGIME)major),baseline);
   PrintFormat("  single BOSS : %d/%d = %.1f%%   macro=%.1f%%   (%.1f us/q)",
               bossCorrect,ntest,100.0*bossCorrect/ntest,bMacro,(double)bossUs/ntest);
   if(ensOK)
      PrintFormat("  BOSS ENSEM  : %d/%d = %.1f%%   macro=%.1f%%   (%.1f us/q)",
                  ensCorrect,ntest,100.0*ensCorrect/ntest,eMacro,(double)ensUs/ntest);
   PrintFormat("  DTW         : %d/%d = %.1f%%   macro=%.1f%%   (%.1f us/q)",
               dtwCorrect,ntest,100.0*dtwCorrect/ntest,dMacro,(double)dtwUs/ntest);
   PrintFormat("  recall  BOSS[tr=%.0f rg=%.0f vo=%.0f] ENS[tr=%.0f rg=%.0f vo=%.0f] DTW[tr=%.0f rg=%.0f vo=%.0f]",
               Recall(bossHit[0],perTot[0]),Recall(bossHit[1],perTot[1]),Recall(bossHit[2],perTot[2]),
               Recall(ensHit[0],perTot[0]), Recall(ensHit[1],perTot[1]), Recall(ensHit[2],perTot[2]),
               Recall(dtwHit[0],perTot[0]), Recall(dtwHit[1],perTot[1]), Recall(dtwHit[2],perTot[2]));
  }

//+------------------------------------------------------------------+
//| Per-class recall as a percentage, 0 if the class is absent.      |
//+------------------------------------------------------------------+
double Recall(int hit,int tot)
  {
   return (tot>0)?100.0*hit/tot:0.0;
  }
//+------------------------------------------------------------------+