//+------------------------------------------------------------------+
//|                                                         frel.mqh |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<Math\Stat\Uniform.mqh>
#include<Math\Stat\Math.mqh>
#include<Powells.mqh>
#include<np.mqh>


//+--------------------------------------------------------------------------------+
//|Class Implementing Feature weighting as regularized energy based learning : FREL|
//+--------------------------------------------------------------------------------+
class FREL:public PowellsMethod
  {
private:
   int               m_num_boot;        //number of bootstrap operations
   int               m_bootsize;        //sample size of each bootstrap operation
   int               m_indices[];       //indices of all observations in m_data (rows)
   matrix            m_data;            //matrix containing candidate predictors with target in the last column
   int               m_target_bin[];    //array containing index of bin each target value belongs to
   double            m_trial_weights[]; //intermediary collection of preliminary weights
   double            m_work_weights[];  //suplementary collenction intermediary weights
   double            m_reg_factor;      //regularization factor controlling degree of regularization
   double            m_loss;            //loss value
   bool              m_memory_allocated;//flag indicating successful memory allocation
   //+------------------------------------------------------------------+
   //| divide data in array into categories of roughly equal size       |
   //+------------------------------------------------------------------+
   bool              bin_array(double &data[],int &num_partitions,double &upperbound_thresholds[],int &categories[])
     {
      int parts;
      bool is_equal ;
      int idata[], indices[], cat_stop[];
      int ibegin, iend, nforward, nbackward, noptimal;
      int best_thresh, best_split ;
      best_thresh = best_split = -1;
      double xdata[] ;
      int n = int(data.Size());

      if(num_partitions > n)
         num_partitions = n ;

      parts = num_partitions ;

      if(ArrayResize(xdata,n)!=n || ArrayResize(idata,n)!=n ||
         ArrayResize(indices,n)!=n || ArrayResize(cat_stop,parts)!=parts)
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }

      for(int i=0 ; i<n ; i++)
        {
         xdata[i] = data[i] ;
         indices[i] = i ;
        }

      MathQuickSortAscending(xdata,indices,0, n-1);
      int k = 0;
      idata[0] = k;
      for(int i=1 ; i<n ; i++)
        {
         if(xdata[i] - xdata[i-1] >= 1.e-12 * (1.0 + fabs(xdata[i]) + fabs(xdata[i-1])))
            ++k ;
         idata[i] = k ;
        }

      k = 0 ;
      int j;
      for(int i=0 ; i<parts ; i++)
        {
         j = (n - k) / (parts - i) ;
         k += j ;
         cat_stop[i] = k-1 ;
        }


      while(true)
        {
         is_equal = false;

         for(int ithresh=0 ; ithresh<parts-1 ; ithresh++)
           {
            if(idata[cat_stop[ithresh]] == idata[cat_stop[ithresh]+1])
              {
               for(int i=ithresh+1 ; i<parts ; i++)
                  cat_stop[i-1] = cat_stop[i] ;
               --parts ;
               is_equal = true ;
               break ;
              }
           }

         if(! is_equal)
            break ;

         ibegin = 0 ;
         noptimal = -1 ;
         for(int ithresh=0 ; ithresh<parts ; ithresh++)
           {
            iend = cat_stop[ithresh] ;
            for(int i=ibegin ; i<iend ; i++)
              {
               if(idata[i] == idata[i+1])
                  continue ;
               nforward = i - ibegin + 1 ;
               nbackward = iend - i ;
               if(nforward < nbackward)
                 {
                  if(nforward > noptimal)
                    {
                     noptimal = nforward ;
                     best_thresh = ithresh ;
                     best_split = i ;
                    }
                 }
               else
                 {
                  if(nbackward > noptimal)
                    {
                     noptimal = nbackward ;
                     best_thresh = ithresh ;
                     best_split = i ;
                    }
                 }
              }
            ibegin = iend + 1 ;
           }

         if(noptimal < 0)
            continue ;

         for(int ithresh=parts-1 ; ithresh>=best_thresh ; ithresh--)
            cat_stop[ithresh+1] = cat_stop[ithresh] ;
         cat_stop[best_thresh] = best_split ;
         ++parts ;
        }

      num_partitions = parts ;

      if(upperbound_thresholds.Size()<uint(parts) && ArrayResize(upperbound_thresholds,parts)!=parts)
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }

      for(int ithresh=0 ; ithresh<parts ; ithresh++)
         upperbound_thresholds[ithresh] = xdata[cat_stop[ithresh]] ;

      ibegin = 0 ;
      for(int ithresh=0 ; ithresh<parts ; ithresh++)
        {
         iend = cat_stop[ithresh] ;
         for(int i=ibegin ; i<=iend ; i++)
            categories[indices[i]] =  ithresh ;
         ibegin = iend + 1 ;
        }

      return true;
     }
   //+------------------------------------------------------------------+
   //|  loss over all data                                              |
   //+------------------------------------------------------------------+

   double            total_loss(double &w[])
     {
      int category,first, other ;
      double  distance, top, bottom, loss ;

      loss = 0.0 ;
      for(int i=0; i<m_bootsize; i++)
        {
         other = m_indices[i] ;
         category = m_target_bin[other] ;
         top = bottom = DBL_MAX ;

         for(int iother=0 ; iother<m_bootsize; iother++)
           {
            first = m_indices[iother] ;
            if(first == other)
               continue ;

            distance = 0.0 ;
            for(ulong v=0 ; v<m_data.Cols()-1; v++)
              {
               distance += w[v] * fabs(m_data[other][v] - m_data[first][v]) ;
              }

            if(m_target_bin[first] == category)
              {
               if(distance < top)
                  top = distance ;
              }
            else
              {
               if(distance < bottom)
                  bottom = distance ;
              }
           }

         distance = top - bottom ;
         if(distance > 30.0)
            loss += distance ;
         else
            loss += log(1.0 + exp(distance));
        }

      return loss ;
     }
   //+------------------------------------------------------------------+
   //| calculates the loss function                                     |
   //+------------------------------------------------------------------+
   double            loss(double &w[])
     {
      double totaloss = total_loss(w);

      totaloss/=double(m_data.Rows());

      if(m_reg_factor>0.0)
        {
         for(ulong i=0; i<m_data.Cols()-1;i++)
            totaloss+=m_reg_factor*pow(w[i],2.0);
        }

      return totaloss;
     }
   //+------------------------------------------------------------------+
   //| function minimized by Powells optimization method                |
   //+------------------------------------------------------------------+
   virtual double          func(const double& p[])
     {
      double pen = 0.0 ;

      for(ulong i=0 ; i<m_data.Cols()-1 ; i++)
        {
         if(p[i] > 4.0)
           {
            m_work_weights[i] = exp(4.0) + p[i] - 4.0 ;
            pen += (p[i] - 4.0) * (p[i] - 4.0) ;
           }
         else
            if(p[i] < -3.0)
              {
               m_work_weights[i] = exp(-3.0) + p[i] + 3.0 ;
               pen += (p[i] + 3.0) * (p[i] + 3.0) ;
              }
            else
               m_work_weights[i] = exp(p[i]) ;
        }

      return (loss(m_work_weights) + pen) ;
     }
   //+------------------------------------------------------------------+
   //| calculates the optimal weights of candidate variables            |
   //+------------------------------------------------------------------+

   int               calc_wt(int num_bins_target,double &loss_value, double &w[])
     {
      int ret,rand_error, class_count[] ;

      ret = 0;

      if(ArrayResize(class_count,num_bins_target)!=num_bins_target || (w.Size()!=uint(m_data.Cols()-1) && ArrayResize(w,int(m_data.Cols()-1))!=int(m_data.Cols()-1)))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return -1;
        }

      ArrayInitialize(w,0.0);
      loss_value = 0.0 ;

      for(ulong i=0 ; i<m_data.Rows() ; i++)
         m_indices[i] = int(i) ;
     
      for(int ibootstrap=0 ; ibootstrap<m_num_boot; ibootstrap++)
        { 
         Comment(" Bootstrap iteration ", ibootstrap+1);
           
         ArrayInitialize(class_count,0);

         int ii, j, k, m;

         ii = int (m_data.Rows()) ;
         while(ii > 1)
           {
            m = int (m_data.Rows()) - ii ;
            if(m >= m_bootsize)
               break ;
            j = (int)(MathRandomUniform(0.0,1.0,rand_error) * ii) ;
            if(j >= ii)
               j = ii - 1 ;
            k = m_indices[m] ;
            m_indices[m] = m_indices[m+j] ;
            m_indices[m+j] = k ;
            --ii ;
            ++class_count[m_target_bin[m_indices[m]]] ;
           }
           

         for(int i=0 ; i<num_bins_target ; i++)
           {
            if(class_count[i] < 2)
               Print(__FUNCTION__, "  class at ", i, " has less than 2 members. Consider adjusting Frel parameters. (number of partitions or bootstrap sample size)");
           }
            
         ArrayInitialize(m_trial_weights,0.0);

         ret += Optimize(m_trial_weights);
         loss_value += PowellsMethod::GetFret() ;

         for(ulong i=0 ; i<m_data.Cols()-1 ; i++)
            w[i] += m_trial_weights[i] ; 

        }

      for(ulong i=0 ; i<m_data.Cols()-1; i++)
        w[i] /= double(m_num_boot) ;

      return ret ;

     }

public:
   //+------------------------------------------------------------------+
   //| constructor                                                      |
   //+------------------------------------------------------------------+

                     FREL(matrix &in_data,int numboot=1, int bootsize=0)
     {
      m_data = in_data;
      m_num_boot=(numboot>0)?numboot:1;
      m_bootsize=(bootsize>2 && bootsize<=int(m_data.Rows()) && m_num_boot>1)?bootsize:int(m_data.Rows());
      
      
      if(ArrayResize(m_indices, int(m_data.Rows()))!=int(m_data.Rows()) ||
         ArrayResize(m_target_bin, int(m_data.Rows()))!=int(m_data.Rows()) ||
         ArrayResize(m_trial_weights, int(m_data.Cols()-1))!=int(m_data.Cols()-1) ||
         ArrayResize(m_work_weights, int(m_data.Cols()-1))!=int(m_data.Cols()-1)
        )
        {
         Print(__FUNCTION__, " error ", GetLastError());
         m_memory_allocated = false;
        }
      else
         m_memory_allocated = true;
     }
   //+------------------------------------------------------------------+
   //| destructor                                                       |
   //+------------------------------------------------------------------+

                    ~FREL(void)
     {
      Comment("");
     }
   //+-----------------------------------------------------------------------+
   //| Find the most relevant variables from a dataset of candidate variables|
   //+-----------------------------------------------------------------------+

   bool               WeighVars(int num_bins_target, double reg_factor,int &index[],double &weights[])
     {
      
      if(!m_memory_allocated)
        {
         Print(" INTERNAL ERROR ");
         return false;
        }
        
        
      if(num_bins_target<=1 || num_bins_target>int(m_data.Rows()))
       {
        Print(__FUNCTION__, " invalid function parameter: num_bins_target. Parameter should be >=2 ");
        return false;
       }  

      int ret=0;
      double target[], target_thresholds[] ;
      double sum ;

      int n_cases = int(m_data.Rows());
      m_reg_factor = MathAbs(reg_factor);
      m_loss = 0.0;


      if(ArrayResize(index,int(m_data.Cols()-1))!=int(m_data.Cols()-1) ||
         !np::vecAsArray(m_data.Col(m_data.Cols()-1),target)
        )
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }

      int k = num_bins_target ;
      if(!bin_array(target, k, target_thresholds, m_target_bin))
         return false;

      if(k<num_bins_target)
        {
         Print("error bins of target vector ", num_bins_target," : ", k);
         return false;
        }

      for(int i=0 ; i<n_cases ; i++)
        {
         if(m_target_bin[i] >= num_bins_target)
           {
            Print("error m_target_bin array at index ", i, " is ",m_target_bin[i], " should be less than ", num_bins_target);
            return false;
           }
        }

      ret = calc_wt(num_bins_target,m_loss,weights);

      if(ret<0)
         return false;

      sum = 0.0 ;
      for(ulong var=0 ; var<m_data.Cols()-1 ; var++)
        {
         weights[var] = m_data.Col(var).Std() * exp(weights[var]);
         sum += weights[var] ;
        }

      for(ulong var=0 ; var<m_data.Cols()-1 ; var++)
        {
         weights[var] *= 100.0 / sum ;
         index[var] = int(var) ;
        }

      MathQuickSortDescending(weights,index,0,int(weights.Size()-1)) ;

      return true;
     }
  };
//+------------------------------------------------------------------+