//+------------------------------------------------------------------+
//|                                                     stepwise.mqh |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<np.mqh>
#define IA 16807
#define IM 2147483647
#define IQ 127773
#define IR 2836
//+------------------------------------------------------------------+
//| standardize a matrix                                             |
//+------------------------------------------------------------------+
matrix            stdmat(matrix &in)
  {
   vector mean = in.Mean(0);
   vector std = in.Std(0);
   std+=1e-10;
   matrix out = in;

   for(ulong row =0; row<out.Rows(); row++)
      if(!out.Row((in.Row(row)-mean)/std,row))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return matrix::Zeros(in.Rows(), in.Cols());
        }

   return out;
  }
//+------------------------------------------------------------------+
//| stepwise feature selection base class                            |
//+------------------------------------------------------------------+
class CStepwisebase
  {
protected:
   matrix            m_data;
   matrix            m_decision_matrix;
   vector            m_target;
   ulong             m_vars;
   ulong             m_samples;
   ulong             m_keptvars;
   ulong             m_folds;
   ulong             m_min_num_preds;
   ulong             m_max_num_preds;
   ulong             m_reps;
   ulong             m_var_indices[];

   ulong             m_trial_length;
   int               m_trial_vars[];
   int               m_already_tried[];
   int               m_all_best_trials[];
   int               m_prior_best_trials[];

   matrix            m_all_best_crits;
   vector            m_prior_best_crits;
   bool              m_verbose;

   //+------------------------------------------------------------------+
   //|   reset the size of all data structures                          |
   //+------------------------------------------------------------------+

   void              reset(void)
     {
      if(m_all_best_trials.Size())
        {
         ArrayFree(m_all_best_trials);
         ArrayFree(m_already_tried);
         ArrayFree(m_trial_vars);
         ArrayFree(m_prior_best_trials);
         ArrayFree(m_var_indices);
        }

      if(m_prior_best_crits.Size())
        {
         m_prior_best_crits = vector::Zeros(0);
         m_all_best_crits = m_decision_matrix = matrix::Zeros(0,0);
        }

      return;
     }
   //+------------------------------------------------------------------+
   //|   Generates  uniform value in [0, 1).                            |
   //+------------------------------------------------------------------+
   double            unif(int &iparam)
     {
      long k ;

      k = iparam / IQ ;
      iparam = IA * (iparam - int(k) * IQ) - IR * int(k) ;
      if(iparam < 0)
         iparam += IM ;
      return double(iparam) / (double) IM ;
     }
   //+------------------------------------------------------------------+
   //|  coordinates main stepwise search operation                      |
   //+------------------------------------------------------------------+

   int               stepwiseSearch(void)
     {
      int done = 1;
      int btrials;
      m_trial_length = m_max_num_preds*m_vars*m_keptvars;

      reset();

      if(!m_prior_best_crits.Resize(m_keptvars) || !m_all_best_crits.Resize(m_reps,m_keptvars))
        {
         Print(__FUNCTION__," data structure resizing error ", GetLastError());
         return done;
        }

      if(ArrayResize(m_all_best_trials,int(m_reps*m_max_num_preds*m_keptvars))<0  ||
         ArrayResize(m_already_tried,int(m_trial_length))<0 ||
         ArrayResize(m_trial_vars,int(m_max_num_preds))<0 ||
         ArrayResize(m_prior_best_trials,int(m_max_num_preds*m_keptvars))<0)
        {
         Print(__FUNCTION__,"  error resizing array ", GetLastError());
         return done;
        }


      ArrayInitialize(m_all_best_trials,-1);
      ArrayInitialize(m_already_tried,-1);
      ArrayInitialize(m_trial_vars,-1);
      ArrayInitialize(m_prior_best_trials,-1);

      vector target;
      int n_so_far =0;
      int rval, rem, j,n_this_rep = 0;
      int mcpt_mod_count, mcpt_change_count;
      mcpt_mod_count = mcpt_mod_count = mcpt_change_count = -1;
      double temp,prior_crit;
      double original_crit, original_change, new_crit;
      original_crit = original_change = new_crit = -DBL_MIN;

      if(m_verbose)
        {
         if(m_reps>1)
            Print("Criterion  || New pval  || Diff pval  || Column Indices");
         else
            Print("Criterion  || Column Indices");
        }
      while(true)
        {

         target = m_target;

         for(ulong rep=0; rep<m_reps; rep++)
           {


            if(rep)
              {
               rval = 17*int(rep) + 11;
               unif(rval);
               unif(rval);
               rem = int(m_samples);
               while(rem>1)
                 {
                  j = (int)(unif(rval))*rem;
                  if(j>=rem)
                     j = rem-1;
                  temp = target[--rem];
                  target[rem] = target[j];
                  target[j] = temp;
                 }
              }

            btrials = int(rep*m_max_num_preds*m_keptvars);
            if(n_so_far == 0)
               prior_crit = -1.e60;
            else
               prior_crit = m_all_best_crits[rep][0];
            n_this_rep = n_so_far;

            done = addFeature(target,n_this_rep,btrials,rep);

            if(done || n_this_rep!=(n_so_far+1))
              {
               Print(__FUNCTION__, " internal error ");
               return 1;
              }

            if(m_all_best_crits[rep][0]<=prior_crit && n_so_far>=int(m_min_num_preds) && !rep)
              {
               for(ulong i = 0 ; i<m_keptvars; i++)
                 {
                  m_all_best_crits[rep][i] = m_prior_best_crits[i];
                  if(ArrayCopy(m_all_best_trials,m_prior_best_trials,btrials+int(i*n_so_far),int(i*n_so_far),n_so_far)<0)
                    {
                     Print(__FUNCTION__, " ArrayCopy error ", GetLastError());
                     return 1;
                    }
                 }
               if(m_verbose)
                   Print(__FUNCTION__, " procedure terminated early because adding a new variable caused performance degradation");
               return 0;
              }

            if(prior_crit < 0.0)
               prior_crit = 0.0;

            new_crit = m_all_best_crits[rep][0];

            if(new_crit<0.0)
               new_crit = 0.0;

            if(rep == 0)
              {
               original_crit = new_crit;
               original_change = new_crit - prior_crit;
               mcpt_mod_count = mcpt_change_count = 1;
              }
            else
              {
               if(new_crit >= original_crit)
                  ++mcpt_mod_count;
               if(new_crit - prior_crit >= original_change)
                  ++mcpt_change_count;
              }

           }

         if(n_so_far == 0)
            mcpt_change_count = mcpt_mod_count;

         if(!m_decision_matrix.Resize(n_so_far+1,m_reps>1?3:1,100) || ArrayResize(m_var_indices,int(m_var_indices.Size())+n_this_rep,100)<0)
           {
            Print(__FUNCTION__, " container resize error ", GetLastError());
            return 1;
           }

         string msg;
         if(m_reps>1)
           {
            double mod_pval = (double) mcpt_mod_count / (double) m_reps;
            double change_pval = (double) mcpt_change_count / (double) m_reps;

            if(m_verbose)
               msg = StringFormat("%10.8lf %10.8lf %10.8lf  ", original_crit,mod_pval,change_pval);

            m_decision_matrix[n_so_far][0] = original_crit;
            m_decision_matrix[n_so_far][1] = mod_pval;
            m_decision_matrix[n_so_far][2] = change_pval;
           }
         else
           {
            msg = m_verbose?StringFormat("%10.8lf", original_crit):NULL;
            m_decision_matrix[n_so_far][0] = original_crit;
           }
         if(m_verbose)
           {
            for(int i = 0; i<n_this_rep; i++)
              {
               msg += StringFormat(" %s",IntegerToString(m_all_best_trials[i]));
              }
           }
         if(ArrayCopy(m_var_indices,m_all_best_trials,(n_so_far*(n_so_far+1))/2,0,n_this_rep)<0)
           {
            Print(__FUNCTION__, " array copy error ", GetLastError());
            return 1;
           }

         if(m_verbose)
            Print(msg);

         ++n_so_far;
         if(n_so_far == int(m_max_num_preds))
           {
            if(m_verbose)
              Print(" Stepwise selection successfully completed ");
            break;
           }
        }

      return 0;

     }
   //+------------------------------------------------------------------+
   //|   adds predictor to collection                                   |
   //+------------------------------------------------------------------+
   int               addFeature(vector &targs,int &index,int btindex, ulong rep_index)
     {
      int done = 1;

      if(index == 0)
        {
         index = 1;
         done = addFirstFeature(targs,btindex,rep_index);
        }
      else
        {
         done = addNewFeature(targs,index,btindex,rep_index);
        }

      return done;
     }
   //+------------------------------------------------------------------+
   //|     finds the first predictor in the set                         |
   //+------------------------------------------------------------------+
   int               addFirstFeature(vector &targ, int bindex,ulong rep_index)
     {
      double crit;

      for(ulong i = 0; i<m_keptvars; i++)
        {
         m_all_best_crits[rep_index][i] = -1.e60;
         m_all_best_trials[bindex+i] = -1;
        }

      for(int i = 0; i<int(m_vars); i++)
        {
         m_trial_vars[0] = i;
         if(crossEval(1,targ,crit))
           {
            Print(__FUNCTION__, " xval error ");
            return 1;
           }
         ulong j;
         for(j = 0; j<m_keptvars; j++)
           {
            if(crit > m_all_best_crits[rep_index][j])
               break;
           }
         if(j<m_keptvars)
           {
            for(long k = long(m_keptvars-2); k>=long(j); k--)
              {
               m_all_best_trials[bindex+k+1] = m_all_best_trials[bindex+k];
               m_all_best_crits[rep_index][k+1] = m_all_best_crits[rep_index][k];
              }
            m_all_best_trials[bindex+j] = i;
            m_all_best_crits[rep_index][j] = crit;
           }
        }

      return 0;
     }
   //+------------------------------------------------------------------+
   //|  looks for next predictor to include                             |
   //+------------------------------------------------------------------+
   int               addNewFeature(vector& target,int &ind, int bindex,ulong rep_index)
     {
      int i, j, k, ivar, ir ;
      int npred, n_already_tried,nbest, rootvars;
      nbest  = -1;
      double crit;

      ind = npred = ind+1;

      for(i = 0; i<int(m_keptvars); i++)
        {
         m_prior_best_crits[i] = m_all_best_crits[rep_index][i];
         if(ArrayCopy(m_prior_best_trials,m_all_best_trials,int(i*(npred-1)),bindex+int(i*(npred-1)),npred-1)<0)
           {
            Print(__FUNCTION__, " ArrayCopy error ", GetLastError());
            return 1;
           }
         m_all_best_crits[rep_index][i] = -1.e60;
        }
      n_already_tried = 0;
      for(ir = 0; ir<int(m_keptvars); ir++)
        {
         if(m_prior_best_crits[ir] < -1.e59)
            break;
        }
      nbest = ir;

      for(ir = 0; ir<nbest; ir++)
        {
         rootvars = ir*(npred-1);
         for(ivar=0; ivar<int(m_vars); ivar++)
           {
            for(i = 0; i<npred-1; i++)
              {
               if(m_prior_best_trials[rootvars+i] == ivar)
                  break;
              }
            if(i < npred-1)
               continue;
            k = 0;

            for(i = 0; i<int(m_vars); i++)
              {
               for(j=0; j<npred-1; j++)
                 {
                  if(m_prior_best_trials[rootvars+j] == i)
                    {
                     m_trial_vars[k++] = i;
                     break;
                    }
                 }
               if(ivar == i)
                  m_trial_vars[k++] = i;
              }
            for(i = 0; i<n_already_tried; i++)
              {
               for(j = 0; j<npred; j++)
                 {
                  if(m_trial_vars[j] != m_already_tried[i*npred+j])
                     break;
                 }
               if(j == npred)
                  break;
              }
            if(i < n_already_tried)
               continue;

            for(i =0; i<npred; i++)
               m_already_tried[n_already_tried*npred+i]=m_trial_vars[i];
            ++n_already_tried;

            if(crossEval(ulong(npred),target,crit))
              {
               Print(__FUNCTION__, " xval error ");
               return 1;
              }

            for(i = 0; i< int(m_keptvars); i++)
              {
               if(crit > m_all_best_crits[rep_index][i])
                  break;
              }

            if(i<int(m_keptvars))
              {
               for(j = int(m_keptvars-2); j>=i; j--)
                 {
                  m_all_best_crits[rep_index][j+1] = m_all_best_crits[rep_index][j];
                  if(ArrayCopy(m_all_best_trials,m_all_best_trials,bindex+int((j+1)*npred),bindex+int(j*npred),npred)<0)
                    {
                     Print(__FUNCTION__, " array copy error ", GetLastError());
                     return 1;
                    }
                 }
               m_all_best_crits[rep_index][i] = crit;
               if(ArrayCopy(m_all_best_trials,m_trial_vars,bindex+int(i*npred),0,npred)<0)
                 {
                  Print(__FUNCTION__, " array copy error ", GetLastError());
                  return 1;
                 }
              }
           }
        }
      return 0;
     }
   //+------------------------------------------------------------------+
   //|  evaluates a predictor set using cross validation                |
   //+------------------------------------------------------------------+
   int               crossEval(ulong num_vars,vector &target,double &criterion)
     {
      ulong ifold, n_remaining, test_start, test_stop ;
      double error, evalresult ;

      n_remaining = m_samples;
      test_start = 0;

      error = 0.0;

      for(ifold = 0; ifold<m_folds; ifold++)
        {
         test_stop = test_start + (n_remaining/(m_folds-ifold));
         if(!fitModel(num_vars,test_start,test_stop,target))
           {
            criterion = -DBL_MIN;
            Print(__FUNCTION__, " fit() failed ");
            return 1;
           }
         evalresult = evalModel(num_vars,test_start,test_stop,target);
         if(evalresult==EMPTY_VALUE)
           {
            criterion = -DBL_MIN;
            Print(__FUNCTION__, " evalModel() failed ");
            return 1;
           }
         error+=evalresult;
         n_remaining -= test_stop - test_start;
         test_start = test_stop;
        }

      criterion = 1.0 - error/double(m_samples);
      return 0;
     }
   //+------------------------------------------------------------------+
   //| fit a model                                                      |
   //+------------------------------------------------------------------+
   virtual bool      fitModel(ulong num_preds,ulong row_start,ulong row_stop, vector& targets) { return false;}
   //+------------------------------------------------------------------+
   //| evaluate a model                                                 |
   //+------------------------------------------------------------------+
   virtual double    evalModel(ulong num_preds,ulong row_start,ulong row_stop, vector& targets) { return EMPTY_VALUE;}
   //+------------------------------------------------------------------+
   //| constructor                                                      |
   //+------------------------------------------------------------------+

public:
                     CStepwisebase(void)
     {
      m_keptvars = 1;
      m_folds = 5;
      m_min_num_preds = m_max_num_preds = ULONG_MAX;
      m_reps = 1;
     }
   //+------------------------------------------------------------------+
   //| destructor                                                       |
   //+------------------------------------------------------------------+

                    ~CStepwisebase(void)
     {
      ArrayFree(m_all_best_trials);
      ArrayFree(m_already_tried);
      ArrayFree(m_trial_vars);
      ArrayFree(m_prior_best_trials);
     }

   //+------------------------------------------------------------------+
   //|   set the stepwise selection parameters                          |
   //+------------------------------------------------------------------+
   void              setParams(ulong num_kept, ulong num_folds, ulong min_num_predictors, ulong max_num_predictors,long num_replications,bool verbose_output)
     {
      m_keptvars = num_kept>0?num_kept:1;
      m_folds = num_folds>0?num_folds:ULONG_MAX;
      m_min_num_preds = min_num_predictors;
      m_max_num_preds = max_num_predictors;
      m_reps = num_replications>0?num_replications:1;
      m_verbose = verbose_output;
      return;
     }

   //+------------------------------------------------------------------+
   //|  main method to execute stepwise feature selection               |
   //+------------------------------------------------------------------+
   bool              stepwiseSelection(matrix& predictors, vector& targets)
     {

      if(targets.Size()!=predictors.Rows() || !targets.Size())
        {
         Print(__FUNCTION__, " invalid inputs ");
         return false;
        }

      m_data = stdmat(predictors);

      m_target = (targets-targets.Mean())/(targets.Std()+1e-10);

      m_vars = predictors.Cols();
      m_samples = predictors.Rows();
      
      if(m_keptvars>m_vars || m_keptvars<1)
         m_keptvars = m_vars;
      
      if(m_folds<1)
         m_folds = 1;
         
      if(m_folds>m_samples)
         m_folds = m_samples;

      if(m_min_num_preds==0 || m_min_num_preds>m_vars || m_min_num_preds>m_max_num_preds)
         m_min_num_preds = m_vars;

      if(m_max_num_preds==0 || m_max_num_preds>m_vars || m_max_num_preds<m_min_num_preds)
         m_max_num_preds = m_vars;
         
      if(m_reps<1)
         m_reps = 1;
         
      if(m_verbose)
        {
         Print(" STEPWISE SELECTION PARAMETERS ");
         Print(" Number of retained candidate variables for each iteration ", m_keptvars);
         Print(" Number of folds in cross validation procedure ", m_folds);
         Print(" Final minimum number of selected predictors ", m_min_num_preds);
         Print(" Final maximum number of selected predictors ", m_max_num_preds);
         Print(" Number of replications for MCP test ", m_reps);
        }
        
      return (stepwiseSearch()==0);
     }

   //+------------------------------------------------------------------+
   //|  get the selected features at a specific iteration of the process|
   //|  iterations denoted as zero based indices                        |
   //|  eg, first selected feature located at iteration index 0         |
   //+------------------------------------------------------------------+

   bool              getFeatureSetAt(ulong iteration_index,ulong &selectedvars[])
     {
      if(ArrayCopy(selectedvars,m_var_indices,0,int((iteration_index*(iteration_index+1))/2),int(iteration_index+1))>0)
         return true;
      else
        {
         Print(__FUNCTION__, " ArrayCopy error ", GetLastError());
         return false;
        }
     }
   //+---------------------------------------------------------------------+
   //| get maximum number of iterations cycled through in selection process|
   //+---------------------------------------------------------------------+
   ulong             getNumFeatureSets(void)
     {
      return m_decision_matrix.Rows();
     }
   //+------------------------------------------------------------------+
   //| get all crit values and corresponding p-values                   |
   //+------------------------------------------------------------------+
   matrix            getCriticalVals(void)
     {
      return m_decision_matrix;
     }

  };
//+------------------------------------------------------------------+