//+------------------------------------------------------------------+
//|                                                 SAXTransform.mqh |
//|                                      MMQ — Muhammad Minhas Qamar |
//|                                   www.mql5.com/en/articles/23484 |
//+------------------------------------------------------------------+
#property copyright "MMQ — Muhammad Minhas Qamar"
#property link      "https://www.mql5.com/en/articles/23484"
#property version   "1.00"
#property strict

//+------------------------------------------------------------------+
//| Symbolic Aggregate approXimation (SAX) — core transform.         |
//|                                                                  |
//|  A raw window of prices is turned into a short string of letters |
//|  in three steps:                                                 |
//|    1. z-normalize the window          (shape, not price level)   |
//|    2. PAA: average it into 'w' cells  (dimensionality reduction) |
//|    3. map each cell to a letter using equal-probability          |
//|       breakpoints of the standard normal   (discretization)      |
//|                                                                  |
//|  Letters are stored as int codes 0..a-1 (0 = lowest band). The   |
//|  MINDIST routine returns a distance between two words that       |
//|  provably lower-bounds the true Euclidean distance of the        |
//|  z-normalized series — this is what makes SAX search sound.      |
//+------------------------------------------------------------------+

//--- status of a single encode attempt
enum ENUM_SAX_STATUS
  {
   SAX_OK,               // word produced and usable
   SAX_BAD_PARAMS,       // w/a out of range, or len < w
   SAX_FLAT_WINDOW       // stdev ~ 0: window is flat, no meaningful shape
  };

//--- limits kept deliberately small; SAX words are meant to be short
#define SAX_MIN_ALPHABET 2
#define SAX_MAX_ALPHABET 10
#define SAX_FLAT_EPS     1.0e-12   // stdev below this -> treat window as flat

//+------------------------------------------------------------------+
//| CSAXTransform - encodes windows and compares words.              |
//|                                                                  |
//|  One instance is configured once (word length w, alphabet size   |
//|  a) and then reused to encode many windows. The Gaussian         |
//|  breakpoints for the chosen alphabet are built once at Configure |
//|  time and cached.                                                |
//+------------------------------------------------------------------+
class CSAXTransform
  {
private:
//--- configuration
   int               m_w;              // word length (number of PAA cells / letters)
   int               m_a;              // alphabet size (number of symbols)
   bool              m_ready;          // Configure() succeeded

//--- cached breakpoints: a-1 cut points that split N(0,1) into
//--- 'a' equal-probability bands. m_beta[i] is the i-th cut.
   double            m_beta[];

//--- helpers
   void              BuildBreakpoints(void);
   double            NormInv(double p) const;

public:
                     CSAXTransform(void);

//--- set word length and alphabet; rebuilds breakpoints. Returns
//--- false if parameters are out of the supported range.
   bool              Configure(int word_len,int alphabet);

   bool              IsReady(void)   const { return m_ready; }
   int               WordLength(void)const { return m_w; }
   int               Alphabet(void)  const { return m_a; }

//--- z-normalize src[0..len-1] into dst[]. Returns false on a flat
//--- window (stdev ~ 0), leaving dst untouched.
   bool              ZNormalize(const double &src[],int len,double &dst[]) const;

//--- PAA: average 'src' (length len) into 'w' cells -> paa[0..w-1].
//--- Handles len not divisible by w via fractional cell boundaries.
   bool              PAA(const double &src[],int len,double &paa[]) const;

//--- full encode: raw window -> integer symbol codes word[0..w-1].
//--- 'raw' holds len prices in chronological order (oldest first).
   ENUM_SAX_STATUS   Encode(const double &raw[],int len,int &word[]) const;

//--- map a single z-normalized value to its symbol code 0..a-1.
   int               SymbolOf(double zval) const;

//--- MINDIST between two words of equal length. Lower-bounds the
//--- Euclidean distance of the underlying z-normalized series when
//--- scaled by sqrt(n/w) (see Lin & Keogh, 2003). 'orig_len' is n,
//--- the length of the original window before PAA.
   double            MinDist(const int &wordA[],const int &wordB[],int orig_len) const;

//--- distance between two symbol codes under the current alphabet:
//--- 0 for adjacent/equal bands, otherwise the gap between the
//--- outer breakpoints they straddle. This is the cell table used
//--- by MINDIST; exposed for inspection/teaching.
   double            CellDist(int c1,int c2) const;

//--- render a word as a lowercase letter string ("cbaabdcc") for
//--- logging and on-chart display.
   string            WordToString(const int &word[]) const;
  };

//+------------------------------------------------------------------+
//| Construct an unconfigured transform.                             |
//+------------------------------------------------------------------+
CSAXTransform::CSAXTransform(void)
  {
   m_w=0;
   m_a=0;
   m_ready=false;
  }

//+------------------------------------------------------------------+
//| Configure word length and alphabet and cache the breakpoints.    |
//+------------------------------------------------------------------+
bool CSAXTransform::Configure(int word_len,int alphabet)
  {
   m_ready=false;
   if(word_len<1)
      return false;
   if(alphabet<SAX_MIN_ALPHABET || alphabet>SAX_MAX_ALPHABET)
      return false;

   m_w=word_len;
   m_a=alphabet;
   BuildBreakpoints();
   m_ready=true;
   return true;
  }

//+------------------------------------------------------------------+
//| Inverse standard-normal CDF (probit) via Acklam's rational       |
//| approximation. Accurate to ~1e-9 over p in (0,1), which is far   |
//| more than the breakpoints need. Returns the z with P(Z<z)=p.     |
//+------------------------------------------------------------------+
double CSAXTransform::NormInv(double p) const
  {
//--- coefficients for Acklam's algorithm
   static const double a1=-3.969683028665376e+01;
   static const double a2= 2.209460984245205e+02;
   static const double a3=-2.759285104469687e+02;
   static const double a4= 1.383577518672690e+02;
   static const double a5=-3.066479806614716e+01;
   static const double a6= 2.506628277459239e+00;

   static const double b1=-5.447609879822406e+01;
   static const double b2= 1.615858368580409e+02;
   static const double b3=-1.556989798598866e+02;
   static const double b4= 6.680131188771972e+01;
   static const double b5=-1.328068155288572e+01;

   static const double c1=-7.784894002430293e-03;
   static const double c2=-3.223964580411365e-01;
   static const double c3=-2.400758277161838e+00;
   static const double c4=-2.549732539343734e+00;
   static const double c5= 4.374664141464968e+00;
   static const double c6= 2.938163982698783e+00;

   static const double d1= 7.784695709041462e-03;
   static const double d2= 3.224671290700398e-01;
   static const double d3= 2.445134137142996e+00;
   static const double d4= 3.754408661907416e+00;

   const double plow =0.02425;
   const double phigh=1.0-plow;

   if(p<=0.0)
      return -DBL_MAX;
   if(p>=1.0)
      return  DBL_MAX;

   double q,r;
   if(p<plow)
     {
      //--- lower tail
      q=MathSqrt(-2.0*MathLog(p));
      return (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) /
             ((((d1*q+d2)*q+d3)*q+d4)*q+1.0);
     }
   if(p>phigh)
     {
      //--- upper tail
      q=MathSqrt(-2.0*MathLog(1.0-p));
      return -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) /
             ((((d1*q+d2)*q+d3)*q+d4)*q+1.0);
     }
//--- central region
   q=p-0.5;
   r=q*q;
   return (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q /
          (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1.0);
  }

//+------------------------------------------------------------------+
//| Build the a-1 breakpoints that split N(0,1) into 'a' bands of    |
//| equal probability 1/a. m_beta[i] = probit((i+1)/a).              |
//+------------------------------------------------------------------+
void CSAXTransform::BuildBreakpoints(void)
  {
   ArrayResize(m_beta,m_a-1);
   for(int i=0;i<m_a-1;i++)
     {
      double p=(double)(i+1)/(double)m_a;
      m_beta[i]=NormInv(p);
     }
  }

//+------------------------------------------------------------------+
//| z-normalize: dst = (src - mean) / stdev. Population stdev used.  |
//| (divide by N), matching the SAX literature. Flat window->false.  |
//+------------------------------------------------------------------+
bool CSAXTransform::ZNormalize(const double &src[],int len,double &dst[]) const
  {
   if(len<1)
      return false;

   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;
     }
   var/=len;                       // population variance
   double sd=MathSqrt(var);
   if(sd<SAX_FLAT_EPS)
      return false;                // flat window: no shape to encode

   ArrayResize(dst,len);
   for(int i=0;i<len;i++)
      dst[i]=(src[i]-mean)/sd;
   return true;
  }

//+------------------------------------------------------------------+
//| PAA: reduce 'src' (length len) to 'w' cell averages.             |
//|  When len is a multiple of w each cell is a clean block. When it |
//|  is not, we accumulate with fractional weights so every cell     |
//|  covers exactly len/w points — the standard general-case PAA.    |
//+------------------------------------------------------------------+
bool CSAXTransform::PAA(const double &src[],int len,double &paa[]) const
  {
   if(len<m_w || m_w<1)
      return false;
   ArrayResize(paa,m_w);

   if(len%m_w==0)
     {
      //--- clean blocks: plain averaging
      int blk=len/m_w;
      for(int j=0;j<m_w;j++)
        {
         double s=0.0;
         int base=j*blk;
         for(int k=0;k<blk;k++)
            s+=src[base+k];
         paa[j]=s/blk;
        }
      return true;
     }

//--- general case: each output cell spans len/w input points, but the
//--- span boundaries fall between samples, so points are split by the
//--- fraction of the cell they belong to.
   double cell=(double)len/(double)m_w;
   for(int j=0;j<m_w;j++)
     {
      double lo=j*cell;            // cell start in fractional index units
      double hi=(j+1)*cell;        // cell end
      double s=0.0;
      int    ilo=(int)MathFloor(lo);
      int    ihi=(int)MathCeil(hi)-1;
      for(int i=ilo;i<=ihi && i<len;i++)
        {
         double left =MathMax((double)i,lo);
         double right=MathMin((double)(i+1),hi);
         double wgt  =right-left;  // overlap of sample i with this cell
         if(wgt>0.0)
            s+=src[i]*wgt;
        }
      paa[j]=s/cell;
     }
   return true;
  }

//+------------------------------------------------------------------+
//| Map a z-normalized value to its symbol code in 0..a-1.           |
//|  Band 0 is (-inf, beta0), band a-1 is [beta_{a-2}, +inf).        |
//+------------------------------------------------------------------+
int CSAXTransform::SymbolOf(double zval) const
  {
   int c=0;
   while(c<m_a-1 && zval>=m_beta[c])
      c++;
   return c;
  }

//+------------------------------------------------------------------+
//| Full encode: raw window -> integer symbol codes.                 |
//|  raw[0..len-1] is oldest-first. The window is z-normalized, PAA  |
//|  reduced, then each cell mapped to a symbol.                     |
//+------------------------------------------------------------------+
ENUM_SAX_STATUS CSAXTransform::Encode(const double &raw[],int len,int &word[]) const
  {
   if(!m_ready || len<m_w)
      return SAX_BAD_PARAMS;

   double zn[];
   if(!ZNormalize(raw,len,zn))
      return SAX_FLAT_WINDOW;

   double paa[];
   if(!PAA(zn,len,paa))
      return SAX_BAD_PARAMS;

   ArrayResize(word,m_w);
   for(int j=0;j<m_w;j++)
      word[j]=SymbolOf(paa[j]);
   return SAX_OK;
  }

//+------------------------------------------------------------------+
//| Distance between two symbol codes (the MINDIST cell table).      |
//|  Adjacent or equal cells contribute 0; otherwise the distance is |
//|  the gap between the breakpoints that separate them:             |
//|     cell(r,c) = 0                       if |r-c| <= 1            |
//|              = beta[max-1] - beta[min]   otherwise               |
//+------------------------------------------------------------------+
double CSAXTransform::CellDist(int c1,int c2) const
  {
   int hi=(c1>c2)?c1:c2;
   int lo=(c1<c2)?c1:c2;
   if(hi-lo<=1)
      return 0.0;
   return m_beta[hi-1]-m_beta[lo];
  }

//+------------------------------------------------------------------+
//| MINDIST between two SAX words of equal length.                   |
//|  MINDIST = sqrt(n/w) * sqrt( sum_j cell(wordA[j],wordB[j])^2 )   |
//|  where n = orig_len. This provably lower-bounds the Euclidean    |
//|  distance of the two z-normalized series, so a MINDIST above a   |
//|  cutoff guarantees the true distance is above it too — the basis |
//|  for sound pruning in analog search. Returns -1 on mismatch.     |
//+------------------------------------------------------------------+
double CSAXTransform::MinDist(const int &wordA[],const int &wordB[],int orig_len) const
  {
   int n1=ArraySize(wordA);
   int n2=ArraySize(wordB);
   if(n1!=n2 || n1!=m_w || orig_len<m_w)
      return -1.0;

   double sum=0.0;
   for(int j=0;j<m_w;j++)
     {
      double d=CellDist(wordA[j],wordB[j]);
      sum+=d*d;
     }
   return MathSqrt((double)orig_len/(double)m_w)*MathSqrt(sum);
  }

//+------------------------------------------------------------------+
//| Render a word as lowercase letters: code 0 -> 'a', 1 -> 'b', ... |
//+------------------------------------------------------------------+
string CSAXTransform::WordToString(const int &word[]) const
  {
   int n=ArraySize(word);
   string s="";
   for(int j=0;j<n;j++)
      s+=CharToString((uchar)('a'+word[j]));
   return s;
  }
//+------------------------------------------------------------------+