Changeset 646

Timestamp:

10/01/09 17:55:03 (15 years ago)

Author:

mido

Message:

arx_straux od sarky

Location:

library/bdm/estim

Files:

• arx_straux.cpp (modified) (4 diffs)
• arx_straux.h (modified) (1 diff)

library/bdm/estim/arx_straux.cpp

@@ +1 @@
+
+
+#include <limits>
 #include "arx.h"
 
+
+#include <fstream>
+#include<iostream>
+#include<iomanip>
+
 namespace bdm {
 
-  struct str_aux {
+  
+
+        
+/*      
+struct str_aux {
         vec d0;
         double nu0;
@@ -15 +27 @@
         int posit1;                // regressand position
         int nbits;                                 // number of bits available in double
-        int bitstr; 
+        bvec bitstr; 
         double loglik;          // loglikelihood
   };
   
-ivec straux1(ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, ivec &rgrsout){
+ */
+
+  
+ /* bvec str_bitset( bvec in,ivec ns,int nbits){
+      //int index, bitindex,n;
+          bvec out = in;
+          int n;
+          
+         for (int i = 0; i < ns.length(); i++){
+          n = ns(i);
+          out(n-2) = 1;
+          cout << out; 
+         
+         }
+        
+        return out; 
+  
+  }*/
+ 
+  
+  void str_bitset( bvec &out,ivec ns,int nbits){
+     for (int i = 0; i < ns.length(); i++){
+          out(ns(i)-2) = 1;
+         }
+}
+
+
+
+
+
+
+double seloglik1(str_aux in){
+// This is the loglikelihood (non-constant part) - this should be used in
+// frequent computation
+   int len = length(in.d);
+   int p1  = in.posit1 - 1;
+   
+   double i1 = -0.5*in.nu *log(in.d(p1)) -0.5*sum(log(in.d.right(len - p1 -1)));
+   double i0 = -0.5*in.nu0*log(in.d0(p1)) -0.5*sum(log(in.d0.right(len - p1 -1)));   
+   return i1-i0;
+ //DEBUGGing print:
+ //fprintf('SELOGLIK1: str=%s loglik=%g\n', strPrintstr(in), l);*/
+ 
+  }
+
+
+void sedydr(mat &r,mat &f,double &Dr,double &Df,int R/*,int jl,int jh ,mat &rout, mat &fout, double &Drout, double &Dfoutint &kr*/){
+/*SEDYDR dyadic reduction, performs transformation of sum of 2 dyads
+%
+% [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh);
+% [rout, fout, Drout, Dfout] = sedydr(r,f,Dr,Df,R);
+%
+% Description: dyadic reduction, performs transformation of sum of 
+%  2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed by R is zeroed
+%
+%     r    : column vector of reduced dyad
+%     f    : column vector of reducing dyad
+%     Dr   : scalar with weight of reduced dyad
+%     Df   : scalar with weight of reducing dyad
+%     R    : scalar number giving 1 based index to the element of r,
+%            which is to be reduced to 
+%            zero; the corresponding element of f is assumed to be 1.
+%     jl   : lower index of the range within which the dyads are 
+%            modified (can be omitted, then everything is updated)
+%     jh   : upper index of the range within which the dyads are 
+%            modified (can be omitted then everything is updated)
+%     rout,fout,Drout,dfout : resulting two dyads
+%     kr   : coefficient used in the transformation of r
+%            rnew = r + kr*f
+%
+% Description: dyadic reduction, performs transformation of sum of 
+%            2 dyads r*Dr*r'+ f*Df*f' so that the element of r indexed by R is zeroed
+% Remark1: Constant mzero means machine zero and should be modified
+%           according to the precision of particular machine
+% Remark2: jl and jh are, in fact, obsolete. It takes longer time to
+%           compute them compared to plain version. The reason is that we
+%           are doing vector operations in m-file. Other reason is that
+%           we need to copy whole vector anyway. It can save half of time for
+%           c-file, if you use it correctly. (please do tests)
+%
+% Note: naming:
+%       se = structure estimation
+%       dydr = dyadic reduction
+%
+% Original Fortran design: V. Peterka 17-7-89
+% Modified for c-language: probably R. Kulhavy
+% Modified for m-language: L. Tesar 2/2003
+% Updated: Feb 2003
+% Project: post-ProDaCTool
+% Reference: none*/
+   
+/*if nargin<6;
+   update_whole=1;
+else
+   update_whole=0;
+end;*/
+ 
+double mzero = 1e-32;
+ 
+if (Dr<mzero){
+   Dr=0;
+}
+ 
+double r0   = r(R,0);
+double kD   = Df;
+double kr   = r0 * Dr;
+
+
+Df   = kD + r0 * kr;
+ 
+if (Df > mzero){
+    kD = kD / Df;
+    kr = kr / Df;
+}else{
+    kD = 1;
+    kr = 0;
+}
+ 
+Dr = Dr * kD;
+ 
+// Try to uncomment marked stuff (*) if in numerical problems, but I don't
+// think it can make any difference for normal healthy floating-point unit
+//if update_whole;
+   r = r - r0*f;
+//   rout(R) = 0;   // * could be needed for some nonsense cases(or numeric reasons?), normally not
+   f = f + kr*r;
+//   fout(R) = 1;   // * could be needed for some nonsense cases(or numeric reasons?), normally not
+/*else;  
+   rout = r;
+   fout = f;
+   rout(jl:jh) = r(jl:jh) - r0 * f(jl:jh);
+   rout(R) = 0;   
+   fout(jl:jh) = f(jl:jh) + kr * rout(jl:jh);
+end;*/
+  }
+
+
+  
+ /*mat*/ void seswapudl(mat &L, vec &d , int i/*, vec &dout*/){
+/*%SESWAPUDL swaps information matrix in decomposition V=L^T diag(d) L
+%
+%   [Lout, dout] = seswapudl(L,d,i);
+%
+% L : lower triangular matrix with 1's on diagonal of the decomposistion
+% d : diagonal vector of diagonal matrix of the decomposition
+% i : index of line to be swapped with the next one 
+% Lout : output lower triangular matrix
+% dout : output diagional vector of diagonal matrix D
+%
+% Description:
+%  Lout' * diag(dout) * Lout = P(i,i+1) * L' * diag(d) * L * P(i,i+1);
+% 
+%  Where permutation matrix P(i,j) permutates columns if applied from the
+%  right and line if applied from the left.
+%   
+% Note: naming:
+%       se = structure estimation
+%       lite = light, simple
+%       udl = U*D*L, or more precisely, L'*D*L, also called as ld
+%   
+% Design  : L. Tesar
+% Updated : Feb 2003
+% Project : post-ProDaCTool
+% Reference: sedydr*/
+ 
+int j = i+1;
+ 
+double pomd = d(i);
+d(i) = d(j);
+d(j) = pomd;
+ 
+/*vec pomL   = L.get_row(i);
+L.set_row(i, L.get_row(j));   
+L.set_row(j,pomL);*/
+ 
+L.swap_rows(i,j);
+L.swap_cols(i,j);
+
+/*pomL   = L.get_col(i);
+L.set_col(i, L.get_col(j));
+L.set_col(j,pomL);*/
+ 
+//% We must be working with LINES of matrix L !
+ 
+
+
+mat r = L.get_row(i);
+r = r.transpose();
+r = r.transpose();
+//????????????????
+mat f = L.get_row(j);
+f = f.transpose();
+f = f.transpose();
+
+
+
+
+double Dr = d(i);
+double Df = d(j);
+ 
+sedydr(r, f, Dr, Df, j);
+ 
+
+
+double r0 = r(i,0);
+Dr = Dr*r0*r0;
+r  = r/r0;
+
+
+
+mat pom_mat = r.transpose();
+L.set_row(i, pom_mat.get_row(0)); 
+pom_mat = f.transpose();
+L.set_row(j, pom_mat.get_row(0)); 
+
+d(i)   = Dr;
+d(j)   = Df;
+ 
+L(i,i) = 1;
+L(j,j) = 1;
+ 
+  
+  }
+  
+
+ void str_bitres(bvec &out,ivec ns,int nbits){
+    
+        
+        for (int i = 0; i < ns.length(); i++){
+          out(ns(i)-2) = 0;
+   }
+     
+  
+  }
+  
+str_aux sestrremove(str_aux in,ivec removed_elements){
+//% Removes elements from regressor
+   int n_strL = length(in.strL);
+    str_aux out = in;
+   for (int i = 0; i < removed_elements.length();i++){
+           
+          int f = removed_elements(i);
+          int posit1 = (find(out.strL==1))(0); 
+      int positf = (find(out.strL==f))(0);
+          int pom_strL;
+          for (int g = positf-1; g >posit1 -1; g--) {   
+             //% BEGIN: We are swapping g and g+1 NOW!!!!
+         seswapudl(out.L, out.d, g);
+         seswapudl(out.L0, out.d0, g);
+       
+                 pom_strL = out.strL(g);
+                 out.strL(g)= out.strL(g+1);
+                 out.strL(g+1) = pom_strL;
+                         
+                 //% END
+              }
+   }
+   out.posit1 = (find(out.strL==1))(0)+1;
+   out.strRgr = out.strL.right(n_strL - out.posit1);
+   out.strMis = out.strL.left(out.posit1-1);
+   str_bitres(out.bitstr,removed_elements,out.nbits);
+   out.loglik = seloglik1(out);
+  
+  return out;
+  }
+  
+
+  ivec setdiff(ivec a, ivec b){
+    ivec pos;
+          
+        for (int i = 0; i < b.length(); i++){
+        pos = find(a==b(i));
+                for (int j = 0; j < pos.length(); j++){
+        a.del(pos(j)-j);
+                }
+        }
+  return a;
+  }
+
+
+  
+/*
+  
+  Array<str_aux> add_new(Array<str_aux> global_best,str_aux newone,int nbest){
+// Eventually add to global best, but do not go over nbest values
+// Also avoids repeating things, which makes this function awfully slow
+          
+          Array<str_aux> global_best_out;
+          if (global_best.length() >= nbest){
+      //logliks = [global_best.loglik];
+      
+          vec logliks(1); 
+          logliks(0) =  global_best(0).loglik;
+          for (int j = 1; j < global_best.length(); j++)
+          logliks = concat(logliks, global_best(j).loglik);
+          
+          int i, addit;
+          double loglik = min(logliks, i);
+      global_best_out = global_best;
+          if (loglik < newone.loglik){
+         //         if ~any(logliks == new.loglik);
+          addit=1;
+                 
+                 
+         
+                  if (newone.bitstr.length() == 1) {
+          for (int j = 0; j < global_best.length(); j++){
+                  for(int i = 0; i < global_best(j).bitstr.length(); i++){
+
+                          if (newone.bitstr(0) == global_best(j).bitstr(i)){
+               addit = 0;
+               break;
+                          }
+                          }   
+          }   
+                  }
+                 if (addit){
+            global_best_out(i) = newone;
+            // DEBUGging print:
+            // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik);
+             }        
+         }  
+}
+   else
+       global_best_out = concat(global_best, newone);
+
+  return global_best_out;
+  
+  }
+  
+*/
+   
+void add_new(Array<str_aux> &global_best,str_aux newone,int nbest){
+// Eventually add to global best, but do not go over nbest values
+// Also avoids repeating things, which makes this function awfully slow
+          
+          int  addit, i = 0;
+if (global_best.length() >= nbest){
+      //logliks = [global_best.loglik];
+      
+          
+          for (int j = 1; j < global_best.length(); j++){
+                          if (global_best(j).loglik < global_best(i).loglik) {
+                                   i = j;
+                          }
+                  }
+         
+          if (global_best(i).loglik < newone.loglik){
+         //         if ~any(logliks == new.loglik);
+          addit=1;
+                 
+                 
+         //????????????????????????????????????????????
+                  // V MATLABU SE MISTO DVOU A VICE DOUBLU
+                  // POROVNAVA KAZDY ZVLAST, COZ DISKRIMINUJE
+                  // VICEROZMERNE (>52) MATICE.. KOD V MATLABU
+                  // JE ZREJME SPATNE .. TODO
+        if (newone.bitstr.length() == 1) {
+          for (int j = 0; j < global_best.length(); j++){
+                  for(int i = 0; i < global_best(j).bitstr.length(); i++){
+
+                          if (newone.bitstr(0) == global_best(j).bitstr(i)){
+               addit = 0;
+               break;
+                          }
+                }   
+          }   
+        }
+                 
+                  
+                //?????????????????????????????????????????????????  
+                  
+                  if (addit){
+            global_best(i) = newone;
+            // DEBUGging print:
+            // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik);
+                  }        
+         }  
+}
+   else
+       global_best = concat(global_best, newone);
+
+}
+   
+  
+  
+  
+  
+ 
+  str_aux sestrinsert(str_aux in,ivec inserted_elements){
+// Moves elements into regressor
+   int n_strL = in.strL.length();
+   str_aux out = in;
+   for (int j = 0;j < inserted_elements.length(); j++){
+      int f = inserted_elements(j);
+      int posit1 = (find(out.strL==1))(0);
+      int positf = (find(out.strL==f))(0);
+          for (int g = positf; g <= posit1-1; g++ ){
+                  
+          // BEGIN: We are swapping g and g+1 NOW!!!!
+          seswapudl(out.L,  out.d,  g);
+          seswapudl(out.L0, out.d0, g);
+          
+          
+                  int pom_strL = out.strL(g);
+                 out.strL(g)= out.strL(g+1);
+                 out.strL(g+1) = pom_strL;
+                  
+                  // END
+          }
+   }         
+
+    out.posit1 = (find(out.strL==1))(0)+1;
+   out.strRgr = out.strL.right(n_strL - out.posit1);
+   out.strMis = out.strL.left(out.posit1-1);
+   str_bitset(out.bitstr,inserted_elements,out.nbits);
+   
+   out.loglik = seloglik1(out);
+   
+   
+   
+   return out;
+  
+  }
+  
+  double seloglik2(str_aux in){
+// This is the loglikelihood (constant part) - this should be added to
+// everything at the end. It needs some computation, so it is useless to
+// make it for all the stuff
+   double logpi = log(pi);
+ 
+   double i1 = lgamma(in.nu /2) - 0.5*in.nu *logpi;
+   double i0 = lgamma(in.nu0/2) - 0.5*in.nu0*logpi;
+   return i1-i0;
+  }
+
+ 
+
+  
+  ivec straux1(ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, Array<str_aux> &rgrsout/*, stat &statistics*/){
 // see utia_legacy/ticket_12/ implementation and str_test.m
 
@@ -30 +481 @@
 int n_data = d.length();
 
-ivec belief_out = find(belief==4)+1; // we are avoiding to put this into regressor
-ivec belief_in  = find(belief==1)+1; // we are instantly keeping this in regressor
+ivec belief_out = find(belief==4)+2; // we are avoiding to put this into regressor
+ivec belief_in  = find(belief==1)+2; // we are instantly keeping this in regressor
+
 
 str_aux full;
@@ -40 +492 @@
 full.L   = L;
 full.d   = d;
+
+
+
+/*full.d0  = "0.012360650875200       0.975779169502626        1.209840558439000";
+         
+full.L0  = "1.0000         0         0;"
+           "0.999427690134298    1.0000         0;"
+           "0.546994424043659   0.534335486953833    1.0000";
+
+
+full.L   = "1.0000         0         0;"
+           "0.364376353850780    1.0000         0;"
+           "1.222141096674815   1.286534510946323    1.0000";
+
+full.d   =  "0.001610356837691          3.497566869589465       3.236640487818002";
+*/
+
+fstream F;
+
+
+
+
+F.open("soubor3.txt", ios::in);
+F << setiosflags(ios::scientific);
+F << setprecision(16);
+F.flags(0x1);
+
+
+for (int i = 0; i < n_data ; i++){
+F >> full.d0(i);
+}
+
+
+
+for (int i =0; i < n_data; i++){
+        for (int j = 0 ; j < n_data  ; j++){
+F >> full.L0(j,i);
+}
+}
+
+for (int i = 0; i < n_data ; i++){
+F >> full.d(i);
+}
+
+
+for (int i =0; i < n_data; i++){
+        for (int j = 0 ; j < n_data  ; j++){
+F >> full.L(j,i);
+}
+}
+
+
+full.nu0 = nu0;
 full.nu  = nu;
 full.strL = linspace(1,n_data);
 full.strRgr = linspace(2,n_data);
 full.strMis = ivec(0);                     
-full.posit1 = 1;   
-//full.nbits  = floor(log2(bitmax))-1;  //!!!!!!!
-
-return ivec(0); // 
-}
-
-
-
-}
+full.posit1 = 1;
+full.bitstr.set_size(n_data-1);
+full.bitstr.clear();
+str_bitset(full.bitstr,full.strRgr,full.nbits);
+//full.nbits  = std::numeric_lim its<double>::digits-1;  // number of bits available in double
+/*bvec in(n_data-1);
+in.clear();
+full.bitstr = str_bitset(in,full.strRgr,full.nbits);*/
+
+
+
+
+full.loglik = seloglik1(full);       // % loglikelihood
+   
+
+
+
+
+//% construct full and empty structure
+full = sestrremove(full,belief_out);
+str_aux empty = sestrremove(full,setdiff(full.strRgr,belief_in));
+ 
+//% stopping rule calculation:
+
+
+
+
+bmat local_max(0,0);
+int to, muto = 0;
+ 
+//% statistics:
+//double cputime0 = cputime; 
+//if nargout>=3;
+   
+ CPU_Timer timer;
+ timer.start(); 
+ 
+   ivec mutos(max_nrep+2);
+   vec maxmutos(max_nrep+2);
+   mutos.zeros();
+   maxmutos.zeros();
+
+
+//end;
+//% ----------------------
+ 
+//% For stopping-rule calculation
+//%so       = 2^(n_data -1-length(belief_in)- length(belief_out)); % do we use this ?
+//% ----------------------
+ 
+ivec all_str = linspace(1,n_data);
+ 
+Array<str_aux> global_best(1); 
+global_best(0) = full;
+
+
+//% MAIN LOOP is here.
+
+str_aux   best;
+for (int n_start = -1; n_start <= max_nrep; n_start++){
+   str_aux  last,best;
+   
+   to = n_start+2;
+   
+   if (n_start == -1){
+      //% start from the full structure
+      last = full;
+   }
+   else {if (n_start == 0)
+      //% start from the empty structure
+      last = empty;     
+       
+      else{
+      //% start from random structure   
+                 
+                ivec last_str = find(to_bvec<int>(::concat<int>(0,floor_i(2*randu(n_data-1)))));// this creates random vector consisting of indexes, and sorted
+                last = sestrremove(full,setdiff(all_str,::concat<int>(::concat<int>(1 ,last_str), empty.strRgr)));
+      
+          }
+   }
+   //% DEBUGging print:
+   //%fprintf('STRUCTURE generated            in loop %2i was %s\n', n_start, strPrintstr(last));
+   
+   //% The loop is repeated until likelihood stops growing (break condition
+   //% used at the end;
+  
+   
+   while (1){
+      //% This structure is going to hold the best elements
+      best = last;
+      //% Nesting by removing elements (enpoorment)
+      ivec  removed_items = setdiff(last.strRgr,belief_in);
+          
+ivec removed_item;
+          str_aux newone;
+          
+          for (int i = 0; i < removed_items.length(); i++){
+        removed_item = vec_1(removed_items(i));
+            newone = sestrremove(last,removed_item);
+                if (nbest>1){
+          add_new(global_best,newone,nbest);
+                }
+                if (newone.loglik>best.loglik){
+           best = newone;
+                }
+          }
+      //% Nesting by adding elements (enrichment)
+      ivec added_items = setdiff(last.strMis,belief_out);
+          ivec added_item;
+          
+          for (int j = 0; j < added_items.length(); j++){
+        added_item = vec_1(added_items(j));
+            newone = sestrinsert(last,added_item);
+                if (nbest>1){
+           add_new(global_best,newone,nbest);
+                }
+                if (newone.loglik>best.loglik){
+           best = newone;
+                }
+     }
+
+   
+
+   
+
+      //% Break condition if likelihood does not change.
+      if (best.loglik <= last.loglik)
+          break;
+      else
+          //% Making best structure last structure.
+          last = best;
+     
+    
+   }
+ 
+
+ 
+
+
+   // % DEBUGging print:
+   //%fprintf('STRUCTURE found (local maxima) in loop %2i was %s randun_seed=%11lu randun_counter=%4lu\n', n_start, strPrintstr(best), randn('seed'), RANDUN_COUNTER);
+   
+   //% Collecting of the best structure in case we don't need the second parameter
+   if (nbest<=1){
+           if (best.loglik > global_best(0).loglik){
+         global_best = best;
+   }
+   }
+ 
+   //% uniqueness of the structure found
+   int append = 1;
+   
+
+   for(int j = 0; j  < local_max.rows() ; j++){
+           if (best.bitstr == local_max.get_row(j)){
+         append = 0;
+                 break;
+           }
+   }
+
+   
+   if (append){
+      local_max.append_row(best.bitstr);
+      muto = muto + 1;
+   } 
+   
+   //% stopping rule:
+   double maxmuto = (to-order_k-1)/lambda-to+1;
+   if (to>2){
+           if (maxmuto>=muto){
+          //% fprintf('*');
+          break;
+           }
+   }      
+ 
+   // do statistics if necessary:
+   //if (nargout>=3){
+      mutos(to-1)    = muto;
+      maxmutos(to-1) = maxmuto;
+   //}
+}
+ 
+//% Aftermath: The best structure was in: global_best
+ 
+//% Updating loglikelihoods: we have to add the constant stuff
+
+
+
+for (int f=0 ; f <global_best.length(); f++){
+   global_best(f).loglik = global_best(f).loglik + seloglik2(global_best(f));
+}
+ 
+/*for f=1:length(global_best);
+   global_best(f).loglik = global_best(f).loglik + seloglik2(global_best(f));
+end;*/
+
+
+//% Making first output parameter:
+
+int max_i = 0;
+for (int j = 1; j < global_best.length(); j++)
+  if  (global_best(max_i).loglik < (global_best(j).loglik)) max_i = j;
+
+best = global_best(max_i);
+
+//% Making the second output parameter
+
+vec logliks(global_best.length());
+for (int j = 0; j < logliks.length(); j++)
+ logliks(j) = global_best(j).loglik;
+
+ivec i = sort_index(logliks);
+rgrsout.set_length(global_best.length());
+
+for (int j = global_best.length() - 1; j >= 0; j--)
+rgrsout(j) = global_best(i(j));
+ 
+//if (nargout>=3);
+   
+   
+str_statistics statistics;
+
+   statistics.allstrs = 2^(n_data -1-length(belief_in) - length(belief_out));
+   statistics.nrand   = to-2;
+   statistics.unique  = muto;
+   statistics.to  = to;
+   statistics.cputime_seconds = timer.get_time();
+   statistics.itemspeed       = statistics.to / statistics.cputime_seconds;
+   statistics.muto = muto;
+   statistics.mutos = mutos;
+   statistics.maxmutos = maxmutos;
+//end;
+
+return best.strRgr;
+
+} 
+
+#ifdef LADIM
+//% randun seed stuff:
+//%randn('seed',SEED);
+ 
+//% --------------------- END of MAIN program --------------------
+ 
+% This is needed for bitstr manipulations
+/*function out = str_bitset(in,ns,nbits)
+   out = in;
+   for n = ns;
+      index = 1+floor((n-2)/nbits);
+      bitindex = 1+rem(n-2,nbits);
+      out(index) = bitset(out(index),bitindex);
+   end;  
+function out = str_bitres(in,ns,nbits)
+   out = in;
+   for n = ns;
+      index = 1+floor((n-2)/nbits);
+      bitindex = 1+rem(n-2,nbits);
+      mask = bitset(0,bitindex);
+      out(index) = bitxor(bitor(out(index),mask),mask);
+   end;*/
+ 
+function out = strPrintstr(in)
+   out = '0';
+   nbits = in.nbits;
+   for f = 2:length(in.d0);
+      index = 1+floor((f-2)/nbits);
+      bitindex = 1+rem(f-2,nbits);
+      if bitget(in.bitstr(index),bitindex);
+          out(f) = '1';
+      else;
+          out(f) = '0';
+      end;
+   end;
+ 
+/*function global_best_out = add_new(global_best,new,nbest)
+% Eventually add to global best, but do not go over nbest values
+% Also avoids repeating things, which makes this function awfully slow
+   if length(global_best)>=nbest;
+      logliks = [global_best.loglik];
+      [loglik i] = min(logliks);
+      global_best_out = global_best;
+      if loglik<new.loglik;
+         %         if ~any(logliks == new.loglik);
+         addit=1;
+         for f = [global_best.bitstr];
+            if f == new.bitstr;
+               addit = 0;
+               break;
+            end;
+         end;         
+         if addit;
+            global_best_out(i) = new;
+            % DEBUGging print:
+            % fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik);
+         end;         
+      end;
+   else;
+      global_best_out = [global_best new];
+   end;*/
+ 
+/*function out = sestrremove(in,removed_elements);
+% Removes elements from regressor
+   n_strL = length(in.strL);
+   out = in;
+   for f=removed_elements;
+      posit1 = find(out.strL==1);
+      positf = find(out.strL==f);
+      for g=(positf-1):-1:posit1;
+         % BEGIN: We are swapping g and g+1 NOW!!!!
+         [out.L, out.d]   = seswapudl(out.L, out.d, g);
+         [out.L0, out.d0]   = seswapudl(out.L0, out.d0, g);
+         out.strL([g g+1]) = [out.strL(g+1) out.strL(g)];
+         % END
+      end;
+   end;
+   out.posit1 = find(out.strL==1);
+   out.strRgr = out.strL((out.posit1+1):n_strL);
+   out.strMis = out.strL(1:(out.posit1-1));
+   out.bitstr = str_bitres(out.bitstr,removed_elements,out.nbits);
+   out.loglik = seloglik1(out);*/
+   
+/*function out = sestrinsert(in,inserted_elements);
+% Moves elements into regressor
+   n_strL = length(in.strL);
+   out = in;
+   for f=inserted_elements;
+      posit1 = find(out.strL==1);
+      positf = find(out.strL==f);
+      for g=positf:(posit1-1);
+          % BEGIN: We are swapping g and g+1 NOW!!!!
+          [out.L,  out.d]   = seswapudl(out.L,  out.d,  g);
+          [out.L0, out.d0]  = seswapudl(out.L0, out.d0, g);
+          out.strL([g g+1]) = [out.strL(g+1) out.strL(g)];
+          % END
+      end;
+   end;         
+   out.posit1 = find(out.strL==1);
+   out.strRgr = out.strL((out.posit1+1):n_strL);
+   out.strMis = out.strL(1:(out.posit1-1));
+   out.bitstr = str_bitset(out.bitstr,inserted_elements,out.nbits);
+   out.loglik = seloglik1(out);*/
+ 
+%
+% seloglik_real = seloglik1 + seloglik2
+%
+ 
+/*function l = seloglik1(in)
+% This is the loglikelihood (non-constant part) - this should be used in
+% frequent computation
+   len = length(in.d);
+   p1  = in.posit1;
+      
+   i1 = -0.5*in.nu *log(in.d (p1)) -0.5*sum(log(in.d ((p1+1):len)));
+   i0 = -0.5*in.nu0*log(in.d0(p1)) -0.5*sum(log(in.d0((p1+1):len)));   
+   l  = i1-i0;
+   
+   % DEBUGGing print:
+   % fprintf('SELOGLIK1: str=%s loglik=%g\n', strPrintstr(in), l);*/
+ 
+ 
+function l = seloglik2(in)
+% This is the loglikelihood (constant part) - this should be added to
+% everything at the end. It needs some computation, so it is useless to
+% make it for all the stuff
+   logpi = log(pi);
+ 
+   i1 = gammaln(in.nu /2) - 0.5*in.nu *logpi;
+   i0 = gammaln(in.nu0/2) - 0.5*in.nu0*logpi;
+   l  = i1-i0;
+ 
+ 
+/*function [Lout, dout] = seswapudl(L,d,i);
+%SESWAPUDL swaps information matrix in decomposition V=L^T diag(d) L
+%
+%  [Lout, dout] = seswapudl(L,d,i);
+%
+% L : lower triangular matrix with 1's on diagonal of the decomposistion
+% d : diagonal vector of diagonal matrix of the decomposition
+% i : index of line to be swapped with the next one 
+% Lout : output lower triangular matrix
+% dout : output diagional vector of diagonal matrix D
+%
+% Description:
+%  Lout' * diag(dout) * Lout = P(i,i+1) * L' * diag(d) * L * P(i,i+1);
+% 
+%  Where permutation matrix P(i,j) permutates columns if applied from the
+%  right and line if applied from the left.
+%   
+% Note: naming:
+%       se = structure estimation
+%       lite = light, simple
+%       udl = U*D*L, or more precisely, L'*D*L, also called as ld
+%   
+% Design  : L. Tesar
+% Updated : Feb 2003
+% Project : post-ProDaCTool
+% Reference: sedydr
+   
+j = i+1;
+ 
+pomd = d(i);
+d(i) = d(j);
+d(j) = pomd;
+ 
+pomL   = L(i,:);
+L(i,:) = L(j,:);
+L(j,:) = pomL;
+ 
+pomL   = L(:,i);
+L(:,i) = L(:,j);
+L(:,j) = pomL;
+ 
+% We must be working with LINES of matrix L !
+ 
+r  = L(i,:)'; 
+f  = L(j,:)';
+Dr = d(i);
+Df = d(j);
+ 
+[r, f, Dr, Df] = sedydr(r, f, Dr, Df, j);
+ 
+r0 = r(i);
+Dr = Dr*r0*r0;
+r  = r/r0;
+ 
+L(i,:) = r'; 
+L(j,:) = f';
+d(i)   = Dr;
+d(j)   = Df;
+ 
+L(i,i) = 1;
+L(j,j) = 1;
+ 
+Lout = L;
+dout = d;*/
+ 
+/*function [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh);
+%SEDYDR dyadic reduction, performs transformation of sum of 2 dyads
+%
+% [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh);
+% [rout, fout, Drout, Dfout] = sedydr(r,f,Dr,Df,R);
+%
+% Description: dyadic reduction, performs transformation of sum of 
+%  2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed by R is zeroed
+%
+%     r    : column vector of reduced dyad
+%     f    : column vector of reducing dyad
+%     Dr   : scalar with weight of reduced dyad
+%     Df   : scalar with weight of reducing dyad
+%     R    : scalar number giving 1 based index to the element of r,
+%            which is to be reduced to 
+%            zero; the corresponding element of f is assumed to be 1.
+%     jl   : lower index of the range within which the dyads are 
+%            modified (can be omitted, then everything is updated)
+%     jh   : upper index of the range within which the dyads are 
+%            modified (can be omitted then everything is updated)
+%     rout,fout,Drout,dfout : resulting two dyads
+%     kr   : coefficient used in the transformation of r
+%            rnew = r + kr*f
+%
+% Description: dyadic reduction, performs transformation of sum of 
+%            2 dyads r*Dr*r'+ f*Df*f' so that the element of r indexed by R is zeroed
+% Remark1: Constant mzero means machine zero and should be modified
+%           according to the precision of particular machine
+% Remark2: jl and jh are, in fact, obsolete. It takes longer time to
+%           compute them compared to plain version. The reason is that we
+%           are doing vector operations in m-file. Other reason is that
+%           we need to copy whole vector anyway. It can save half of time for
+%           c-file, if you use it correctly. (please do tests)
+%
+% Note: naming:
+%       se = structure estimation
+%       dydr = dyadic reduction
+%
+% Original Fortran design: V. Peterka 17-7-89
+% Modified for c-language: probably R. Kulhavy
+% Modified for m-language: L. Tesar 2/2003
+% Updated: Feb 2003
+% Project: post-ProDaCTool
+% Reference: none
+   
+if nargin<6;
+   update_whole=1;
+else
+   update_whole=0;
+end;
+ 
+mzero = 1e-32;
+ 
+if Dr<mzero;
+   Dr=0;
+end;
+ 
+r0   = r(R);
+kD   = Df;
+kr   = r0 * Dr;
+Dfout   = kD + r0 * kr;
+ 
+if Dfout > mzero;
+    kD = kD / Dfout;
+    kr = kr / Dfout;
+else;
+    kD = 1;
+    kr = 0;
+end;
+ 
+Drout = Dr * kD;
+ 
+% Try to uncomment marked stuff (*) if in numerical problems, but I don't
+% think it can make any difference for normal healthy floating-point unit
+if update_whole;
+   rout = r - r0*f;
+%   rout(R) = 0;   % * could be needed for some nonsense cases(or numeric reasons?), normally not
+   fout = f + kr*rout;
+%   fout(R) = 1;   % * could be needed for some nonsense cases(or numeric reasons?), normally not
+else;  
+   rout = r;
+   fout = f;
+   rout(jl:jh) = r(jl:jh) - r0 * f(jl:jh);
+   rout(R) = 0;   
+   fout(jl:jh) = f(jl:jh) + kr * rout(jl:jh);
+end;*/
+ 
+ 
+ 
+#endif
+ 
+ 
+}
+

library/bdm/estim/arx_straux.h

@@ -16 +16 @@
 namespace bdm {
 
+        
+struct str_aux {
+        vec d0;
+        double nu0;
+        mat L0;
+        mat L;
+        vec d;
+        double nu;
+        ivec strL;                 // Current structure of L and d
+        ivec strRgr;               // Structure elements currently inside regressor (after regressand)
+        ivec strMis;               // structure elements, that are currently outside regressor (before regressand)
+        int posit1;                // regressand position
+        int nbits;                                 // number of bits available in double
+        bvec bitstr; 
+        double loglik;          // loglikelihood
+  };
+  
+
+
+struct str_statistics {
+ long long int allstrs;
+ int nrand;
+ int unique;
+ int to;
+ double cputime_seconds;
+ double itemspeed;
+ int muto;
+ ivec mutos;
+ vec maxmutos;
+};
+
+
         //! Rplication of Ludvik Tesar original straux1 from mixtools straux1
-ivec straux1(ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, ivec &rgrsout);
+ivec straux1(ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, Array<str_aux> &rgrsout);
 
 }