//+------------------------------------------------------------------+
//| HSMMFilter.mqh                                                   |
//| Copyright 2026, Gbadebo Adewumi                                  |
//| https://www.mql5.com                                             |
//+------------------------------------------------------------------+
#property copyright "Gbadebo Adewumi"
#property link      "https://www.mql5.com"

//+------------------------------------------------------------------+
//| Compile-time model constants (three-way contract)                |
//+------------------------------------------------------------------+
//--- these four constants define the fixed shape of the HSMM model.
//--- they must match the corresponding fields inside hsmm_manifest.json
//--- EXACTLY, or the model trained offline in Python will not line up
//--- with the arrays allocated here. ValidateContract() enforces this
//--- at OnInit() time and refuses to run on any mismatch, so a stale
//--- manifest can never silently drive the EA with garbage indices.
#define NUM_STATES        4   // TrendUp, TrendDown, Range, HighVolChop
#define NUM_FEATURES       3  // realized_vol_dev, efficiency_ratio, return_skew
#define MAX_DURATION       60 // longest modelled regime sojourn, in bars
#define FEATURE_LOOKBACK   20 // rolling window (bars) used to build each feature

//+------------------------------------------------------------------+
//| Human-readable state labels (for logs and chart text)            |
//+------------------------------------------------------------------+
string HSMM_DefaultStateNames[NUM_STATES] =
  {
   "TrendUp",    // state 0 -- sustained directional up-move
   "TrendDown",  // state 1 -- sustained directional down-move
   "Range",      // state 2 -- low-efficiency mean-reverting chop
   "HighVolChop" // state 3 -- volatility expansion without direction
  };

//+------------------------------------------------------------------+
//| Class CJSONArrayUtils                                            |
//| Minimal bracket-aware JSON array parsing helpers                 |
//+------------------------------------------------------------------+
class CJSONArrayUtils
  {
public:
   static void       SplitTopLevel(string s, string &parts_out[]);
   static void       ParseDoubleArray(string raw, double &out[]);
   static string     ExtractValueRaw(string json, string key);
   static string     StripBrackets(string s);
  };

//+------------------------------------------------------------------+
//| CJSONArrayUtils::StripBrackets                                   |
//+------------------------------------------------------------------+
string CJSONArrayUtils::StripBrackets(string s)
  {
//--- trims whitespace and, if present, one outer pair of [ ] or { }
//--- this is called before splitting so callers never have to worry
//--- about leading/trailing bracket characters throwing off the
//--- depth counter used by SplitTopLevel().
   StringTrimLeft(s);
   StringTrimRight(s);
   int len = StringLen(s);
   if(len >= 2)
     {
      ushort first = StringGetCharacter(s, 0);
      ushort last  = StringGetCharacter(s, len - 1);
      if((first == '[' && last == ']') || (first == '{' && last == '}'))
         s = StringSubstr(s, 1, len - 2);
     }
   return(s);
  }

//+------------------------------------------------------------------+
//| CJSONArrayUtils::SplitTopLevel                                   |
//+------------------------------------------------------------------+
void CJSONArrayUtils::SplitTopLevel(string s, string &parts_out[])
  {
//--- splits a comma-separated JSON list at commas that sit at
//--- bracket-depth zero only. A manifest field such as
//--- "[[1,2],[3,4]]" must NOT be split at the commas that separate
//--- 1 from 2, only at the comma that separates the two inner
//--- arrays -- so we track '[' / '{' depth as we scan and only
//--- treat a comma as a delimiter when depth == 0.
   ArrayResize(parts_out, 0);
   int depth = 0;
   int start = 0;
   int len = StringLen(s);
   for(int i = 0; i < len; i++)
     {
      ushort c = StringGetCharacter(s, i);
      if(c == '[' || c == '{')
         depth++;
      else if(c == ']' || c == '}')
         depth--;
      else if(c == ',' && depth == 0)
        {
         string token = StringSubstr(s, start, i - start);
         StringTrimLeft(token);
         StringTrimRight(token);
         int n = ArraySize(parts_out);
         ArrayResize(parts_out, n + 1);
         parts_out[n] = token;
         start = i + 1;
        }
     }
//--- flush the final token -- there is no trailing comma after the
//--- last element, so the loop above never emits it on its own.
   string last_token = StringSubstr(s, start, len - start);
   StringTrimLeft(last_token);
   StringTrimRight(last_token);
   if(StringLen(last_token) > 0)
     {
      int n2 = ArraySize(parts_out);
      ArrayResize(parts_out, n2 + 1);
      parts_out[n2] = last_token;
     }
  }

//+------------------------------------------------------------------+
//| CJSONArrayUtils::ParseDoubleArray                                |
//+------------------------------------------------------------------+
void CJSONArrayUtils::ParseDoubleArray(string raw, double &out[])
  {
//--- flattens a JSON numeric array of ANY nesting depth (1-D, 2-D or
//--- 3-D) into a single flat double array, row-major. Flattening
//--- instead of using MQL5 dynamic multi-dimensional arrays lets the
//--- caller index with simple compile-time strides (state*COLS+col),
//--- which is both faster and easier to reason about than nested
//--- dynamic arrays.
   string body = CJSONArrayUtils::StripBrackets(raw);
   string tokens[];
   CJSONArrayUtils::SplitTopLevel(body, tokens);
   ArrayResize(out, 0);
   for(int i = 0; i < ArraySize(tokens); i++)
     {
      string t = tokens[i];
      StringTrimLeft(t);
      StringTrimRight(t);
      if(StringLen(t) == 0)
         continue;
      ushort c0 = StringGetCharacter(t, 0);
      if(c0 == '[')
        {
         //--- nested array -- recurse and append the child's flat
         //--- output onto our own flat output buffer.
         double child[];
         CJSONArrayUtils::ParseDoubleArray(t, child);
         int base = ArraySize(out);
         ArrayResize(out, base + ArraySize(child));
         for(int k = 0; k < ArraySize(child); k++)
            out[base + k] = child[k];
        }
      else
        {
         //--- leaf value -- a plain JSON number.
         int base2 = ArraySize(out);
         ArrayResize(out, base2 + 1);
         out[base2] = StringToDouble(t);
        }
     }
  }

//+------------------------------------------------------------------+
//| CJSONArrayUtils::ExtractValueRaw                                 |
//+------------------------------------------------------------------+
string CJSONArrayUtils::ExtractValueRaw(string json, string key)
  {
//--- returns the raw (still JSON-encoded) value text that follows
//--- "\"key\":" inside json, respecting bracket depth so a value that
//--- is itself an object/array is captured whole rather than cut off
//--- at its first internal comma. This is a targeted extractor (not
//--- a general parser) because the manifest has a small, fixed key
//--- set -- a full recursive-descent grammar would be more code for
//--- no practical benefit here.
   string needle = "\"" + key + "\":";
   int pos = StringFind(json, needle);
   if(pos < 0)
      return("");
   int cursor = pos + StringLen(needle);
   int len = StringLen(json);
//--- skip whitespace between the colon and the value proper.
   while(cursor < len && StringGetCharacter(json, cursor) == ' ')
      cursor++;
   ushort first = StringGetCharacter(json, cursor);
   int start = cursor;
   if(first == '[' || first == '{')
     {
      //--- bracketed value: walk forward counting depth until the
      //--- matching close bracket is found.
      ushort open_ch  = first;
      ushort close_ch = (first == '[') ? ']' : '}';
      int depth = 0;
      for(int i = cursor; i < len; i++)
        {
         ushort c = StringGetCharacter(json, i);
         if(c == open_ch)
            depth++;
         else if(c == close_ch)
           {
            depth--;
            if(depth == 0)
               return(StringSubstr(json, start, i - start + 1));
           }
        }
      return(StringSubstr(json, start, len - start));
     }
   else if(first == '"')
     {
      //--- quoted string value: find the closing quote.
      int i = cursor + 1;
      while(i < len && StringGetCharacter(json, i) != '"')
         i++;
      return(StringSubstr(json, start, i - start + 1));
     }
   else
     {
      //--- bare number/literal: ends at the next comma or closing
      //--- brace at the current (top) level.
      int i = cursor;
      while(i < len)
        {
         ushort c = StringGetCharacter(json, i);
         if(c == ',' || c == '}' || c == '\n' || c == '\r')
            break;
         i++;
        }
      return(StringSubstr(json, start, i - start));
     }
  }

//+------------------------------------------------------------------+
//| Class CHSMMFilter                                                |
//| Duration-aware (explicit-duration / residual-time) HSMM filter   |
//+------------------------------------------------------------------+
class CHSMMFilter
  {
private:
   //--- fitted model parameters, loaded from hsmm_manifest.json
   double            m_emission_mean[NUM_STATES * NUM_FEATURES];   // per-state feature means (z-scored space)
   double            m_emission_var[NUM_STATES * NUM_FEATURES];    // per-state feature variances (diagonal cov)
   double            m_transition[NUM_STATES * NUM_STATES];        // off-diagonal transition-on-exit matrix
   double            m_duration_pmf[NUM_STATES * MAX_DURATION];    // P(total sojourn = k) at index [state*MAX_DURATION+(k-1)]
   double            m_feat_norm_mean[NUM_FEATURES];               // training-set feature means, for z-scoring live features
   double            m_feat_norm_std[NUM_FEATURES];                // training-set feature std devs, for z-scoring live features
   string            m_state_names[NUM_STATES];                    // human-readable labels, read from manifest

   //--- live belief state -- this is what actually moves every bar
   double            m_alpha[NUM_STATES * MAX_DURATION];           // alpha[state*MAX_DURATION+d] = P(state=j, remaining duration=d | O_1..t)

   //--- values parsed straight from the manifest, used only to run
   //--- ValidateContract() against the compile-time constants above
   int               m_manifest_num_states;
   int               m_manifest_num_features;
   int               m_manifest_max_duration;

   bool              m_loaded;                                     // true once LoadManifest() has succeeded

   void              NormalizeFeatures(const double &raw_features[], double &norm_features[]);
   double            GaussianLogLik(int state, const double &norm_features[]);
   bool              ReadWholeFile(string filename, string &out_text);

public:
                     CHSMMFilter(void);
   bool              LoadManifest(string filename);
   bool              ValidateContract(void);
   void              ResetBelief(void);
   void              PredictStep(void);
   void              UpdateStep(const double &raw_features[]);
   void              StepFilter(const double &raw_features[]);
   int               GetFilteredState(double &confidence);
   double            GetExpectedRemainingDuration(int state);
   string            StateName(int state);
   bool              ComputeFeatures(const double &close[], int shift, double atr_value, double &features_out[]);
  };

//+------------------------------------------------------------------+
//| CHSMMFilter::CHSMMFilter (constructor)                           |
//+------------------------------------------------------------------+
CHSMMFilter::CHSMMFilter(void)
  {
   m_loaded = false;
   m_manifest_num_states   = 0;
   m_manifest_num_features = 0;
   m_manifest_max_duration = 0;
   ArrayInitialize(m_alpha, 0.0);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::ReadWholeFile                                       |
//+------------------------------------------------------------------+
bool CHSMMFilter::ReadWholeFile(string filename, string &out_text)
  {
//--- FileReadString() on a FILE_TXT handle silently returns only the
//--- FIRST line of the file (a known MQL5 gotcha), so a naive single
//--- call would truncate the manifest to its opening brace. Instead
//--- we open in FILE_ANSI mode explicitly and loop line-by-line with
//--- FileIsEnding(), concatenating everything before handing the
//--- full text off to the JSON key extractor.
   int handle = FileOpen(filename, FILE_READ | FILE_TXT | FILE_ANSI | FILE_COMMON);
   if(handle == INVALID_HANDLE)
     {
      Print("CHSMMFilter: failed to open manifest file: ", filename, " err=", GetLastError());
      return(false);
     }
   out_text = "";
   while(!FileIsEnding(handle))
     {
      string line = FileReadString(handle);
      out_text += line;
     }
   FileClose(handle);
   return(true);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::LoadManifest                                        |
//+------------------------------------------------------------------+
bool CHSMMFilter::LoadManifest(string filename)
  {
   string json;
   if(!ReadWholeFile(filename, json))
      return(false);

//--- scalar contract fields -- read first so ValidateContract() can
//--- reject a mismatched manifest before we waste time parsing the
//--- (potentially large) numeric arrays below.
   m_manifest_num_states   = (int)StringToInteger(CJSONArrayUtils::ExtractValueRaw(json, "num_states"));
   m_manifest_max_duration = (int)StringToInteger(CJSONArrayUtils::ExtractValueRaw(json, "max_duration"));

//--- feature_names is an array; its element COUNT (not its content)
//--- is what the three-way contract actually checks against
//--- NUM_FEATURES -- there is no separate scalar "feature count"
//--- field in the manifest, since that would just be a second source
//--- of truth that could itself drift out of sync with the array.
   string feat_names_raw = CJSONArrayUtils::ExtractValueRaw(json, "feature_names");
   string feat_name_tokens[];
   CJSONArrayUtils::SplitTopLevel(CJSONArrayUtils::StripBrackets(feat_names_raw), feat_name_tokens);
   m_manifest_num_features = ArraySize(feat_name_tokens);

   if(!ValidateContract())
      return(false);

//--- numeric parameter blocks -- each parsed into a flat double
//--- array via CJSONArrayUtils::ParseDoubleArray, then copied into
//--- the fixed-size member arrays index-for-index (row-major, so
//--- emission_mean[state][feature] lives at state*NUM_FEATURES+feature).
   double tmp[];

   CJSONArrayUtils::ParseDoubleArray(CJSONArrayUtils::ExtractValueRaw(json, "emission_mean"), tmp);
   for(int i = 0; i < ArraySize(tmp) && i < NUM_STATES * NUM_FEATURES; i++)
      m_emission_mean[i] = tmp[i];

   CJSONArrayUtils::ParseDoubleArray(CJSONArrayUtils::ExtractValueRaw(json, "emission_var"), tmp);
   for(int i = 0; i < ArraySize(tmp) && i < NUM_STATES * NUM_FEATURES; i++)
      m_emission_var[i] = tmp[i];

   CJSONArrayUtils::ParseDoubleArray(CJSONArrayUtils::ExtractValueRaw(json, "transition_matrix"), tmp);
   for(int i = 0; i < ArraySize(tmp) && i < NUM_STATES * NUM_STATES; i++)
      m_transition[i] = tmp[i];

   CJSONArrayUtils::ParseDoubleArray(CJSONArrayUtils::ExtractValueRaw(json, "duration_pmf"), tmp);
   for(int i = 0; i < ArraySize(tmp) && i < NUM_STATES * MAX_DURATION; i++)
      m_duration_pmf[i] = tmp[i];

   CJSONArrayUtils::ParseDoubleArray(CJSONArrayUtils::ExtractValueRaw(json, "feature_norm_mean"), tmp);
   for(int i = 0; i < ArraySize(tmp) && i < NUM_FEATURES; i++)
      m_feat_norm_mean[i] = tmp[i];

   CJSONArrayUtils::ParseDoubleArray(CJSONArrayUtils::ExtractValueRaw(json, "feature_norm_std"), tmp);
   for(int i = 0; i < ArraySize(tmp) && i < NUM_FEATURES; i++)
      m_feat_norm_std[i] = tmp[i];

//--- state_names is a string array, not numeric, so it gets its own
//--- lightweight split-and-strip-quotes pass rather than reusing
//--- ParseDoubleArray.
   string names_raw = CJSONArrayUtils::ExtractValueRaw(json, "state_names");
   string name_tokens[];
   CJSONArrayUtils::SplitTopLevel(CJSONArrayUtils::StripBrackets(names_raw), name_tokens);
   for(int i = 0; i < ArraySize(name_tokens) && i < NUM_STATES; i++)
     {
      string t = name_tokens[i];
      StringTrimLeft(t);
      StringTrimRight(t);
      StringReplace(t, "\"", "");
      m_state_names[i] = t;
     }

   m_loaded = true;
   ResetBelief();
   Print("CHSMMFilter: manifest loaded OK from ", filename);
   return(true);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::ValidateContract                                    |
//+------------------------------------------------------------------+
bool CHSMMFilter::ValidateContract(void)
  {
//--- the three-way contract: compile-time constants, manifest
//--- fields, and (implicitly, via array sizes above) the actual
//--- parsed array lengths must all agree. We fail loudly here rather
//--- than silently truncating or padding a mismatched array, because
//--- a silent truncation would let the EA run with a corrupted model
//--- and produce plausible-looking but meaningless regime calls.
   bool ok = true;
   if(m_manifest_num_states != NUM_STATES)
     {
      Print("CHSMMFilter: CONTRACT MISMATCH -- manifest num_states=", m_manifest_num_states,
            " but compiled NUM_STATES=", NUM_STATES);
      ok = false;
     }
   if(m_manifest_num_features != NUM_FEATURES)
     {
      Print("CHSMMFilter: CONTRACT MISMATCH -- manifest feature_names count=", m_manifest_num_features,
            " but compiled NUM_FEATURES=", NUM_FEATURES);
      ok = false;
     }
   if(m_manifest_max_duration != MAX_DURATION)
     {
      Print("CHSMMFilter: CONTRACT MISMATCH -- manifest max_duration=", m_manifest_max_duration,
            " but compiled MAX_DURATION=", MAX_DURATION);
      ok = false;
     }
   return(ok);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::ResetBelief                                         |
//+------------------------------------------------------------------+
void CHSMMFilter::ResetBelief(void)
  {
//--- cold-start prior: spread belief mass uniformly across all
//--- (state, duration) cells. This is deliberately uninformative --
//--- after a few dozen bars the recursive update in StepFilter()
//--- will have concentrated the mass onto whichever state/duration
//--- combination the observed features actually support, so a crude
//--- starting prior does not meaningfully bias the converged filter.
   double uniform = 1.0 / (double)(NUM_STATES * MAX_DURATION);
   for(int i = 0; i < NUM_STATES * MAX_DURATION; i++)
      m_alpha[i] = uniform;
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::PredictStep                                         |
//+------------------------------------------------------------------+
void CHSMMFilter::PredictStep(void)
  {
//--- explicit-duration (residual-time) HSMM time update. Each belief
//--- cell alpha[j,d] holds P(state=j, d bars remaining in this state).
//--- Two things can happen to the mass in a cell between bar t-1 and
//--- bar t:
//---   1) d > 0: the state simply continues -- remaining duration
//---      counts down by exactly one, deterministically. No branch,
//---      no probability multiply, because the duration was already
//---      drawn when the state was entered.
//---   2) d == 0: this bar was the LAST bar of the state's sojourn,
//---      so the process transitions OUT to some new state j' (drawn
//---      from the off-diagonal transition matrix) and simultaneously
//---      draws a brand-new total duration k for j' from that state's
//---      duration pmf, landing in cell (j', k-1).
   double pred[NUM_STATES * MAX_DURATION];
   ArrayInitialize(pred, 0.0);

   double finishing[NUM_STATES];
   for(int j = 0; j < NUM_STATES; j++)
      finishing[j] = m_alpha[j * MAX_DURATION + 0];

//--- case 1: deterministic countdown for every cell with d > 0.
   for(int j = 0; j < NUM_STATES; j++)
     {
      for(int d = 1; d < MAX_DURATION; d++)
        {
         double mass = m_alpha[j * MAX_DURATION + d];
         if(mass <= 0.0)
            continue;
         pred[j * MAX_DURATION + (d - 1)] += mass;
        }
     }

//--- case 2: mass finishing state j is redistributed into every
//--- other state j', weighted by the transition probability and by
//--- that state's own duration pmf (which sets its NEW remaining
//--- duration on entry). Truncating k at MAX_DURATION means a tiny
//--- amount of probability mass for exceptionally long regimes is
//--- dropped rather than wrapped or clipped into the last bucket --
//--- an acceptable approximation as long as the duration pmf was fit
//--- so its tail is already small beyond MAX_DURATION bars.
   for(int j = 0; j < NUM_STATES; j++)
     {
      if(finishing[j] <= 0.0)
         continue;
      for(int jp = 0; jp < NUM_STATES; jp++)
        {
         if(jp == j)
            continue;
         double outflow = finishing[j] * m_transition[j * NUM_STATES + jp];
         if(outflow <= 0.0)
            continue;
         for(int k = 1; k <= MAX_DURATION; k++)
           {
            double p_dur = m_duration_pmf[jp * MAX_DURATION + (k - 1)];
            if(p_dur <= 0.0)
               continue;
            pred[jp * MAX_DURATION + (k - 1)] += outflow * p_dur;
           }
        }
     }

   ArrayCopy(m_alpha, pred);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::NormalizeFeatures                                   |
//+------------------------------------------------------------------+
void CHSMMFilter::NormalizeFeatures(const double &raw_features[], double &norm_features[])
  {
//--- the emission Gaussians were fit in Python on z-scored features
//--- (see train_hsmm.py), so live features must be put through the
//--- SAME training-set mean/std before scoring -- feeding raw,
//--- un-normalized features into GaussianLogLik would silently shift
//--- every state's likelihood by an arbitrary, feature-dependent
//--- amount and break the whole filter.
   ArrayResize(norm_features, NUM_FEATURES);
   for(int f = 0; f < NUM_FEATURES; f++)
     {
      double std = m_feat_norm_std[f];
      if(std < 1e-8)
         std = 1e-8;   // guard against a degenerate (near-zero variance) training feature
      norm_features[f] = (raw_features[f] - m_feat_norm_mean[f]) / std;
     }
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::GaussianLogLik                                      |
//+------------------------------------------------------------------+
double CHSMMFilter::GaussianLogLik(int state, const double &norm_features[])
  {
//--- diagonal-covariance multivariate Gaussian log-density, summed
//--- feature-by-feature. Working in log space (rather than
//--- multiplying raw densities together) avoids underflow to exactly
//--- zero when several features are simultaneously several standard
//--- deviations from a state's mean -- a real risk once NUM_FEATURES
//--- and the number of states both grow.
   double log_lik = 0.0;
   for(int f = 0; f < NUM_FEATURES; f++)
     {
      double var = m_emission_var[state * NUM_FEATURES + f];
      if(var < 1e-8)
         var = 1e-8;    // guard: a state that saw almost no feature variance during training
      double mean = m_emission_mean[state * NUM_FEATURES + f];
      double diff = norm_features[f] - mean;
      log_lik += -0.5 * MathLog(2.0 * M_PI * var) - 0.5 * (diff * diff) / var;
     }
   return(log_lik);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::UpdateStep                                          |
//+------------------------------------------------------------------+
void CHSMMFilter::UpdateStep(const double &raw_features[])
  {
//--- Bayesian measurement update: multiply every (state,duration)
//--- cell by that state's emission likelihood for the bar just
//--- closed, then renormalize. Duration does not affect the
//--- emission -- by construction the observation model only reads
//--- the discrete regime label, not how long it has been active --
//--- so every cell belonging to state j gets the SAME multiplier.
   double norm_features[];
   NormalizeFeatures(raw_features, norm_features);

   double log_lik[NUM_STATES];
   double max_log_lik = -DBL_MAX;
   for(int j = 0; j < NUM_STATES; j++)
     {
      log_lik[j] = GaussianLogLik(j, norm_features);
      if(log_lik[j] > max_log_lik)
         max_log_lik = log_lik[j];
     }

//--- numerically stable exponentiation: subtract the max log-lik
//--- before calling MathExp() so the largest term becomes exp(0)=1
//--- and every other term is a well-behaved value in (0,1], instead
//--- of risking exp() underflow/overflow on raw log-likelihoods that
//--- can run to several hundred in magnitude.
   double lik[NUM_STATES];
   for(int j = 0; j < NUM_STATES; j++)
      lik[j] = MathExp(log_lik[j] - max_log_lik);

   double total = 0.0;
   for(int j = 0; j < NUM_STATES; j++)
     {
      for(int d = 0; d < MAX_DURATION; d++)
        {
         m_alpha[j * MAX_DURATION + d] *= lik[j];
         total += m_alpha[j * MAX_DURATION + d];
        }
     }

//--- renormalize so the belief stays a proper probability
//--- distribution. If total collapses to (near) zero -- e.g. after
//--- many bars of extremely unlikely observations under every state
//--- -- fall back to a fresh uniform prior rather than dividing by
//--- (near) zero and injecting NaNs into the belief state.
   if(total < 1e-300)
     {
      ResetBelief();
      return;
     }
   for(int i = 0; i < NUM_STATES * MAX_DURATION; i++)
      m_alpha[i] /= total;
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::StepFilter                                          |
//+------------------------------------------------------------------+
void CHSMMFilter::StepFilter(const double &raw_features[])
  {
//--- convenience wrapper: advance the belief by exactly one closed
//--- bar. Predict must run BEFORE Update -- Predict propagates last
//--- bar's posterior forward through the countdown/transition
//--- dynamics to form this bar's prior, and only then does Update
//--- fold in this bar's actual observation. Calling these in the
//--- opposite order would condition on an observation before the
//--- state had even been allowed to evolve, which is not the filter
//--- this class is documented to implement.
   PredictStep();
   UpdateStep(raw_features);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::GetFilteredState                                    |
//+------------------------------------------------------------------+
int CHSMMFilter::GetFilteredState(double &confidence)
  {
//--- marginalizes the belief over duration to get P(state=j) for
//--- each state, then returns the argmax and its probability mass as
//--- a confidence score the caller can threshold against before
//--- acting on a regime call.
   double state_prob[NUM_STATES];
   for(int j = 0; j < NUM_STATES; j++)
     {
      double s = 0.0;
      for(int d = 0; d < MAX_DURATION; d++)
         s += m_alpha[j * MAX_DURATION + d];
      state_prob[j] = s;
     }
   int best = 0;
   for(int j = 1; j < NUM_STATES; j++)
      if(state_prob[j] > state_prob[best])
         best = j;
   confidence = state_prob[best];
   return(best);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::GetExpectedRemainingDuration                        |
//+------------------------------------------------------------------+
double CHSMMFilter::GetExpectedRemainingDuration(int state)
  {
//--- E[remaining bars | state=j] = sum_d d * alpha[j,d] / sum_d alpha[j,d].
//--- This is the number that makes the model "duration-aware" in a way
//--- a plain HMM's transition-matrix diagonal never gives you directly:
//--- a caller can gate entries on the regime STILL having runway left,
//--- not just on which regime is currently most likely.
   double num = 0.0, den = 0.0;
   for(int d = 0; d < MAX_DURATION; d++)
     {
      double mass = m_alpha[state * MAX_DURATION + d];
      num += d * mass;
      den += mass;
     }
   if(den < 1e-12)
      return(0.0);
   return(num / den);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::StateName                                           |
//+------------------------------------------------------------------+
string CHSMMFilter::StateName(int state)
  {
   if(state < 0 || state >= NUM_STATES)
      return("Unknown");
   return(m_state_names[state]);
  }

//+------------------------------------------------------------------+
//| CHSMMFilter::ComputeFeatures                                     |
//+------------------------------------------------------------------+
bool CHSMMFilter::ComputeFeatures(const double &close[], int shift, double atr_value, double &features_out[])
  {
//--- close[] must be a series-indexed array (ArraySetAsSeries(...,true))
//--- so that close[shift] is the bar we are scoring and close[shift+1]
//--- is one bar older. FEATURE_LOOKBACK+1 closes are required to form
//--- FEATURE_LOOKBACK log-returns, so callers must guarantee that many
//--- bars exist before calling this -- callers should check
//--- Bars(...) >= FEATURE_LOOKBACK+2 before entering their bar loop.
   ArrayResize(features_out, NUM_FEATURES);

   double returns[FEATURE_LOOKBACK];
   for(int i = 0; i < FEATURE_LOOKBACK; i++)
     {
      double c_new = close[shift + i];
      double c_old = close[shift + i + 1];
      if(c_old <= 0.0 || c_new <= 0.0)
         return(false);   // guard: corrupt/zero price would poison MathLog below
      returns[i] = MathLog(c_new / c_old);
     }

//--- feature 0: realized_vol_dev -- how far the SHORT-WINDOW realized
//--- volatility of returns sits from the volatility implied by the
//--- broker's own ATR reading. Zero means "about what ATR would
//--- predict"; a large positive value flags a volatility expansion
//--- that ATR (a slower, smoothed measure) has not yet caught up to
//--- -- exactly the kind of regime break this model is meant to
//--- detect early.
   double mean_r = 0.0;
   for(int i = 0; i < FEATURE_LOOKBACK; i++)
      mean_r += returns[i];
   mean_r /= FEATURE_LOOKBACK;

   double var_r = 0.0;
   for(int i = 0; i < FEATURE_LOOKBACK; i++)
     {
      double d = returns[i] - mean_r;
      var_r += d * d;
     }
   var_r /= FEATURE_LOOKBACK;
   double realized_sigma = MathSqrt(var_r);

   double atr_frac = atr_value / close[shift];
   if(atr_frac < 1e-8)
      atr_frac = 1e-8;   // guard against a zero/near-zero ATR reading on a very quiet symbol
   features_out[0] = (realized_sigma / atr_frac) - 1.0;

//--- feature 1: Kaufman efficiency ratio -- net directional travel
//--- divided by total path length walked to get there, over the same
//--- lookback window. Close to 1.0 means the price moved in a
//--- straight line (a clean trend); close to 0.0 means it wandered
//--- back and forth and ended up near where it started (chop/range).
   double net_change = MathAbs(close[shift] - close[shift + FEATURE_LOOKBACK]);
   double path_sum = 0.0;
   for(int i = 0; i < FEATURE_LOOKBACK; i++)
      path_sum += MathAbs(close[shift + i] - close[shift + i + 1]);
   if(path_sum < 1e-12)
      path_sum = 1e-12;   // guard: a perfectly flat window would otherwise divide by zero
   features_out[1] = net_change / path_sum;

//--- feature 2: return skew -- third standardized moment of the
//--- return window. Distinguishes a slow grinding trend (mildly
//--- skewed) from a panic move (heavily skewed in one direction),
//--- which is useful because two regimes can share similar realized
//--- volatility and efficiency ratio yet have very different skew
//--- signatures.
   double m3 = 0.0;
   for(int i = 0; i < FEATURE_LOOKBACK; i++)
     {
      double d = returns[i] - mean_r;
      m3 += d * d * d;
     }
   m3 /= FEATURE_LOOKBACK;
   double sigma3 = MathPow(realized_sigma, 3.0);
   if(sigma3 < 1e-12)
      sigma3 = 1e-12;   // guard: near-zero variance window would otherwise blow up the ratio
   features_out[2] = m3 / sigma3;

   return(true);
  }