Changeset 688 for library

Show
Ignore:
Timestamp:
10/30/09 18:04:45 (16 years ago)
Author:
sarka
Message:

function straux1 converted from straux1.m

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • library/bdm/estim/arx_straux.cpp

    r684 r688  
    11 
    2  
    3 #include <limits> 
    42#include "arx.h" 
    53 
    64 
    7 #include <fstream> 
    8 #include<iostream> 
    9 #include<iomanip> 
    10  
    115namespace bdm { 
    126 
    13    
    14  
    15          
    16 /*       
    17 struct str_aux { 
    18         vec d0; 
    19         double nu0; 
    20         mat L0; 
    21         mat L; 
    22         vec d; 
    23         double nu; 
    24         ivec strL;                 // Current structure of L and d 
    25         ivec strRgr;               // Structure elements currently inside regressor (after regressand) 
    26         ivec strMis;               // structure elements, that are currently outside regressor (before regressand) 
    27         int posit1;                // regressand position 
    28         int nbits;                                 // number of bits available in double 
    29         bvec bitstr;  
    30         double loglik;          // loglikelihood 
    31   }; 
    32    
    33  */ 
    34  
    35    
    36  /* bvec str_bitset( bvec in,ivec ns,int nbits){ 
    37       //int index, bitindex,n; 
    38           bvec out = in; 
    39           int n; 
    40            
    41          for (int i = 0; i < ns.length(); i++){ 
    42           n = ns(i); 
    43           out(n-2) = 1; 
    44           cout << out;  
    45           
    46          } 
    47          
    48         return out;  
    49    
    50   }*/ 
    51   
    52    
    53   void str_bitset( bvec &out,ivec ns,int nbits){ 
    54      for (int i = 0; i < ns.length(); i++){ 
    55           out(ns(i)-2) = 1; 
    56          } 
    57 } 
    58  
    59  
    60  
    61  
    62  
    63  
    64 double seloglik1(str_aux in){ 
     7 
     8 
     9         
     10void str_bitset ( bvec &out, ivec ns, int nbits ) { 
     11         
     12for ( int i = 0; i < ns.length(); i++ ) { 
     13                out ( ns ( i ) - 2 ) = 1; 
     14        } 
     15} 
     16 
     17 
     18double seloglik1 ( str_aux in ) { 
    6519// This is the loglikelihood (non-constant part) - this should be used in 
    6620// frequent computation 
    67    int len = length(in.d); 
    68    int p1  = in.posit1 - 1; 
    69     
    70    double i1 = -0.5*in.nu *log(in.d(p1)) -0.5*sum(log(in.d.right(len - p1 -1))); 
    71    double i0 = -0.5*in.nu0*log(in.d0(p1)) -0.5*sum(log(in.d0.right(len - p1 -1)));    
    72    return i1-i0; 
    73  //DEBUGGing print: 
    74  //fprintf('SELOGLIK1: str=%s loglik=%g\n', strPrintstr(in), l);*/ 
    75   
    76   } 
    77  
    78  
    79 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*/){ 
    80 /*SEDYDR dyadic reduction, performs transformation of sum of 2 dyads 
    81 % 
    82 % [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh); 
    83 % [rout, fout, Drout, Dfout] = sedydr(r,f,Dr,Df,R); 
    84 % 
    85 % Description: dyadic reduction, performs transformation of sum of  
    86 %  2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed by R is zeroed 
    87 % 
    88 %     r    : column vector of reduced dyad 
    89 %     f    : column vector of reducing dyad 
    90 %     Dr   : scalar with weight of reduced dyad 
    91 %     Df   : scalar with weight of reducing dyad 
    92 %     R    : scalar number giving 1 based index to the element of r, 
    93 %            which is to be reduced to  
    94 %            zero; the corresponding element of f is assumed to be 1. 
    95 %     jl   : lower index of the range within which the dyads are  
    96 %            modified (can be omitted, then everything is updated) 
    97 %     jh   : upper index of the range within which the dyads are  
    98 %            modified (can be omitted then everything is updated) 
    99 %     rout,fout,Drout,dfout : resulting two dyads 
    100 %     kr   : coefficient used in the transformation of r 
    101 %            rnew = r + kr*f 
    102 % 
    103 % Description: dyadic reduction, performs transformation of sum of  
    104 %            2 dyads r*Dr*r'+ f*Df*f' so that the element of r indexed by R is zeroed 
    105 % Remark1: Constant mzero means machine zero and should be modified 
    106 %           according to the precision of particular machine 
    107 % Remark2: jl and jh are, in fact, obsolete. It takes longer time to 
    108 %           compute them compared to plain version. The reason is that we 
    109 %           are doing vector operations in m-file. Other reason is that 
    110 %           we need to copy whole vector anyway. It can save half of time for 
    111 %           c-file, if you use it correctly. (please do tests) 
    112 % 
    113 % Note: naming: 
    114 %       se = structure estimation 
    115 %       dydr = dyadic reduction 
    116 % 
    117 % Original Fortran design: V. Peterka 17-7-89 
    118 % Modified for c-language: probably R. Kulhavy 
    119 % Modified for m-language: L. Tesar 2/2003 
    120 % Updated: Feb 2003 
    121 % Project: post-ProDaCTool 
    122 % Reference: none*/ 
    123     
    124 /*if nargin<6; 
    125    update_whole=1; 
    126 else 
    127    update_whole=0; 
    128 end;*/ 
    129   
    130 double mzero = 1e-32; 
    131   
    132 if (Dr<mzero){ 
    133    Dr=0; 
    134 } 
    135   
    136 double r0   = r(R,0); 
    137 double kD   = Df; 
    138 double kr   = r0 * Dr; 
    139  
    140  
    141 Df   = kD + r0 * kr; 
    142   
    143 if (Df > mzero){ 
    144     kD = kD / Df; 
    145     kr = kr / Df; 
    146 }else{ 
    147     kD = 1; 
    148     kr = 0; 
    149 } 
    150   
    151 Dr = Dr * kD; 
    152   
     21        int len = length ( in.d ); 
     22        int p1  = in.posit1 - 1; 
     23 
     24        double i1 = -0.5 * in.nu * log ( in.d ( p1 ) ) - 0.5 * sum ( log ( in.d.right ( len - p1 - 1 ) ) ); 
     25        double i0 = -0.5 * in.nu0 * log ( in.d0 ( p1 ) ) - 0.5 * sum ( log ( in.d0.right ( len - p1 - 1 ) ) ); 
     26        return i1 - i0; 
     27//DEBUGGing print: 
     28//fprintf('SELOGLIK1: str=%s loglik=%g\n', strPrintstr(in), l);*/ 
     29} 
     30 
     31 
     32void sedydr ( mat &r, mat &f, double &Dr, double &Df, int R/*,int jl,int jh ,mat &rout, mat &fout, double &Drout, double &Dfoutint &kr*/ ) { 
     33        /*SEDYDR dyadic reduction, performs transformation of sum of 2 dyads 
     34        % 
     35        % [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh); 
     36        % [rout, fout, Drout, Dfout] = sedydr(r,f,Dr,Df,R); 
     37        % 
     38        % Description: dyadic reduction, performs transformation of sum of 
     39        %  2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed by R is zeroed 
     40        % 
     41        %     r    : column vector of reduced dyad 
     42        %     f    : column vector of reducing dyad 
     43        %     Dr   : scalar with weight of reduced dyad 
     44        %     Df   : scalar with weight of reducing dyad 
     45        %     R    : scalar number giving 1 based index to the element of r, 
     46        %            which is to be reduced to 
     47        %            zero; the corresponding element of f is assumed to be 1. 
     48        %     jl   : lower index of the range within which the dyads are 
     49        %            modified (can be omitted, then everything is updated) 
     50        %     jh   : upper index of the range within which the dyads are 
     51        %            modified (can be omitted then everything is updated) 
     52        %     rout,fout,Drout,dfout : resulting two dyads 
     53        %     kr   : coefficient used in the transformation of r 
     54        %            rnew = r + kr*f 
     55        % 
     56        % Description: dyadic reduction, performs transformation of sum of 
     57        %            2 dyads r*Dr*r'+ f*Df*f' so that the element of r indexed by R is zeroed 
     58        % Remark1: Constant mzero means machine zero and should be modified 
     59        %           according to the precision of particular machine 
     60        % Remark2: jl and jh are, in fact, obsolete. It takes longer time to 
     61        %           compute them compared to plain version. The reason is that we 
     62        %           are doing vector operations in m-file. Other reason is that 
     63        %           we need to copy whole vector anyway. It can save half of time for 
     64        %           c-file, if you use it correctly. (please do tests) 
     65        % 
     66        % Note: naming: 
     67        %       se = structure estimation 
     68        %       dydr = dyadic reduction 
     69        % 
     70        % Original Fortran design: V. Peterka 17-7-89 
     71        % Modified for c-language: probably R. Kulhavy 
     72        % Modified for m-language: L. Tesar 2/2003 
     73        % Updated: Feb 2003 
     74        % Project: post-ProDaCTool 
     75        % Reference: none*/ 
     76 
     77        /*if nargin<6; 
     78           update_whole=1; 
     79        else 
     80           update_whole=0; 
     81        end;*/ 
     82 
     83        double mzero = 1e-32; 
     84 
     85        if ( Dr < mzero ) { 
     86                Dr = 0; 
     87        } 
     88 
     89        double r0   = r ( R, 0 ); 
     90        double kD   = Df; 
     91        double kr   = r0 * Dr; 
     92 
     93 
     94        Df   = kD + r0 * kr; 
     95 
     96        if ( Df > mzero ) { 
     97                kD = kD / Df; 
     98                kr = kr / Df; 
     99        } else { 
     100                kD = 1; 
     101                kr = 0; 
     102        } 
     103 
     104        Dr = Dr * kD; 
     105 
    153106// Try to uncomment marked stuff (*) if in numerical problems, but I don't 
    154107// think it can make any difference for normal healthy floating-point unit 
    155108//if update_whole; 
    156    r = r - r0*f; 
     109        r = r - r0 * f; 
    157110//   rout(R) = 0;   // * could be needed for some nonsense cases(or numeric reasons?), normally not 
    158    f = f + kr*r; 
     111        f = f + kr * r; 
    159112//   fout(R) = 1;   // * could be needed for some nonsense cases(or numeric reasons?), normally not 
    160 /*else;   
    161    rout = r; 
    162    fout = f; 
    163    rout(jl:jh) = r(jl:jh) - r0 * f(jl:jh); 
    164    rout(R) = 0;    
    165    fout(jl:jh) = f(jl:jh) + kr * rout(jl:jh); 
    166 end;*/ 
    167   } 
    168  
    169  
    170    
    171  /*mat*/ void seswapudl(mat &L, vec &d , int i/*, vec &dout*/){ 
    172 /*%SESWAPUDL swaps information matrix in decomposition V=L^T diag(d) L 
    173 % 
    174 %   [Lout, dout] = seswapudl(L,d,i); 
    175 % 
    176 % L : lower triangular matrix with 1's on diagonal of the decomposistion 
    177 % d : diagonal vector of diagonal matrix of the decomposition 
    178 % i : index of line to be swapped with the next one  
    179 % Lout : output lower triangular matrix 
    180 % dout : output diagional vector of diagonal matrix D 
    181 % 
    182 % Description: 
    183 %  Lout' * diag(dout) * Lout = P(i,i+1) * L' * diag(d) * L * P(i,i+1); 
    184 %  
    185 %  Where permutation matrix P(i,j) permutates columns if applied from the 
    186 %  right and line if applied from the left. 
    187 %    
    188 % Note: naming: 
    189 %       se = structure estimation 
    190 %       lite = light, simple 
    191 %       udl = U*D*L, or more precisely, L'*D*L, also called as ld 
    192 %    
    193 % Design  : L. Tesar 
    194 % Updated : Feb 2003 
    195 % Project : post-ProDaCTool 
    196 % Reference: sedydr*/ 
    197   
    198 int j = i+1; 
    199   
    200 double pomd = d(i); 
    201 d(i) = d(j); 
    202 d(j) = pomd; 
    203   
    204 /*vec pomL   = L.get_row(i); 
    205 L.set_row(i, L.get_row(j));    
    206 L.set_row(j,pomL);*/ 
    207   
    208 L.swap_rows(i,j); 
    209 L.swap_cols(i,j); 
    210  
    211 /*pomL   = L.get_col(i); 
    212 L.set_col(i, L.get_col(j)); 
    213 L.set_col(j,pomL);*/ 
    214   
     113        /*else; 
     114           rout = r; 
     115           fout = f; 
     116           rout(jl:jh) = r(jl:jh) - r0 * f(jl:jh); 
     117           rout(R) = 0; 
     118           fout(jl:jh) = f(jl:jh) + kr * rout(jl:jh); 
     119        end;*/ 
     120} 
     121 
     122 
     123 
     124/*mat*/ void seswapudl ( mat &L, vec &d , int i/*, vec &dout*/ ) { 
     125        /*%SESWAPUDL swaps information matrix in decomposition V=L^T diag(d) L 
     126        % 
     127        %   [Lout, dout] = seswapudl(L,d,i); 
     128        % 
     129        % L : lower triangular matrix with 1's on diagonal of the decomposistion 
     130        % d : diagonal vector of diagonal matrix of the decomposition 
     131        % i : index of line to be swapped with the next one 
     132        % Lout : output lower triangular matrix 
     133        % dout : output diagional vector of diagonal matrix D 
     134        % 
     135        % Description: 
     136        %  Lout' * diag(dout) * Lout = P(i,i+1) * L' * diag(d) * L * P(i,i+1); 
     137        % 
     138        %  Where permutation matrix P(i,j) permutates columns if applied from the 
     139        %  right and line if applied from the left. 
     140        % 
     141        % Note: naming: 
     142        %       se = structure estimation 
     143        %       lite = light, simple 
     144        %       udl = U*D*L, or more precisely, L'*D*L, also called as ld 
     145        % 
     146        % Design  : L. Tesar 
     147        % Updated : Feb 2003 
     148        % Project : post-ProDaCTool 
     149        % Reference: sedydr*/ 
     150 
     151        int j = i + 1; 
     152 
     153        double pomd = d ( i ); 
     154        d ( i ) = d ( j ); 
     155        d ( j ) = pomd; 
     156 
     157        L.swap_rows ( i, j ); 
     158        L.swap_cols ( i, j ); 
     159 
     160         
    215161//% We must be working with LINES of matrix L ! 
    216   
    217  
    218  
    219 mat r = L.get_row(i); 
    220 r = r.transpose(); 
    221 r = r.transpose(); 
    222 //???????????????? 
    223 mat f = L.get_row(j); 
    224 f = f.transpose(); 
    225 f = f.transpose(); 
    226  
    227  
    228  
    229  
    230 double Dr = d(i); 
    231 double Df = d(j); 
    232   
    233 sedydr(r, f, Dr, Df, j); 
    234   
    235  
    236  
    237 double r0 = r(i,0); 
    238 Dr = Dr*r0*r0; 
    239 r  = r/r0; 
    240  
    241  
    242  
    243 mat pom_mat = r.transpose(); 
    244 L.set_row(i, pom_mat.get_row(0));  
    245 pom_mat = f.transpose(); 
    246 L.set_row(j, pom_mat.get_row(0));  
    247  
    248 d(i)   = Dr; 
    249 d(j)   = Df; 
    250   
    251 L(i,i) = 1; 
    252 L(j,j) = 1; 
    253   
    254    
    255   } 
    256    
    257  
    258  void str_bitres(bvec &out,ivec ns,int nbits){ 
    259      
    260          
    261         for (int i = 0; i < ns.length(); i++){ 
    262           out(ns(i)-2) = 0; 
    263    } 
    264       
    265    
    266   } 
    267    
    268 str_aux sestrremove(str_aux in,ivec removed_elements){ 
     162 
     163    mat r = L.get_row ( i ); 
     164        r.transpose(); 
     165        mat f = L.get_row ( j ); 
     166        f.transpose(); 
     167 
     168    double Dr = d ( i ); 
     169        double Df = d ( j ); 
     170 
     171        sedydr ( r, f, Dr, Df, j ); 
     172 
     173    double r0 = r ( i, 0 ); 
     174        Dr = Dr * r0 * r0; 
     175        r  = r / r0; 
     176 
     177    mat pom_mat = r.transpose(); 
     178        L.set_row ( i, pom_mat.get_row ( 0 ) ); 
     179        pom_mat = f.transpose(); 
     180        L.set_row ( j, pom_mat.get_row ( 0 ) ); 
     181 
     182        d ( i )   = Dr; 
     183        d ( j )   = Df; 
     184 
     185        L ( i, i ) = 1; 
     186        L ( j, j ) = 1; 
     187} 
     188 
     189 
     190void str_bitres ( bvec &out, ivec ns, int nbits ) { 
     191 
     192  for ( int i = 0; i < ns.length(); i++ ) { 
     193                out ( ns ( i ) - 2 ) = 0; 
     194        } 
     195} 
     196 
     197str_aux sestrremove ( str_aux in, ivec removed_elements ) { 
    269198//% Removes elements from regressor 
    270    int n_strL = length(in.strL); 
    271     str_aux out = in; 
    272    for (int i = 0; i < removed_elements.length();i++){ 
    273             
    274           int f = removed_elements(i); 
    275           int posit1 = (find(out.strL==1))(0);  
    276       int positf = (find(out.strL==f))(0); 
    277           int pom_strL; 
    278           for (int g = positf-1; g >posit1 -1; g--) {    
    279              //% BEGIN: We are swapping g and g+1 NOW!!!! 
    280          seswapudl(out.L, out.d, g); 
    281          seswapudl(out.L0, out.d0, g); 
    282         
    283                  pom_strL = out.strL(g); 
    284                  out.strL(g)= out.strL(g+1); 
    285                  out.strL(g+1) = pom_strL; 
    286                           
    287                  //% END 
    288               } 
    289    } 
    290    out.posit1 = (find(out.strL==1))(0)+1; 
    291    out.strRgr = out.strL.right(n_strL - out.posit1); 
    292    out.strMis = out.strL.left(out.posit1-1); 
    293    str_bitres(out.bitstr,removed_elements,out.nbits); 
    294    out.loglik = seloglik1(out); 
    295    
    296   return out; 
    297   } 
    298    
    299  
    300   ivec setdiff(ivec a, ivec b){ 
    301     ivec pos; 
    302            
    303         for (int i = 0; i < b.length(); i++){ 
    304         pos = find(a==b(i)); 
    305                 for (int j = 0; j < pos.length(); j++){ 
    306         a.del(pos(j)-j); 
    307                 } 
    308         } 
    309   return a; 
    310   } 
    311  
    312  
    313    
    314 /* 
    315    
    316   Array<str_aux> add_new(Array<str_aux> global_best,str_aux newone,int nbest){ 
     199         
     200        int n_strL = length ( in.strL ); 
     201        str_aux out = in; 
     202         
     203    for ( int i = 0; i < removed_elements.length(); i++ ) { 
     204 
     205                int f = removed_elements ( i ); 
     206                int posit1 = ( find ( out.strL == 1 ) ) ( 0 ); 
     207                int positf = ( find ( out.strL == f ) ) ( 0 ); 
     208                for ( int g = positf - 1; g > posit1 - 1; g-- ) { 
     209                        //% BEGIN: We are swapping g and g+1 NOW!!!! 
     210                        seswapudl ( out.L, out.d, g ); 
     211                        seswapudl ( out.L0, out.d0, g ); 
     212 
     213                        int pom_strL = out.strL ( g ); 
     214                        out.strL ( g ) = out.strL ( g + 1 ); 
     215                        out.strL ( g + 1 ) = pom_strL; 
     216                        //% END 
     217                } 
     218         
     219         
     220        } 
     221        out.posit1 = ( find ( out.strL == 1 ) ) ( 0 ) + 1; 
     222        out.strRgr = out.strL.right ( n_strL - out.posit1 ); 
     223        out.strMis = out.strL.left ( out.posit1 - 1 ); 
     224        str_bitres ( out.bitstr, removed_elements, out.nbits ); 
     225        out.loglik = seloglik1 ( out ); 
     226 
     227        return out; 
     228} 
     229 
     230 
     231ivec setdiff ( ivec a, ivec b ) { 
     232        ivec pos; 
     233 
     234        for ( int i = 0; i < b.length(); i++ ) { 
     235                pos = find ( a == b ( i ) ); 
     236                for ( int j = 0; j < pos.length(); j++ ) { 
     237                        a.del ( pos ( j ) - j ); 
     238                } 
     239        } 
     240        return a; 
     241} 
     242 
     243 
     244 
     245 
     246void add_new ( Array<str_aux> &global_best, str_aux newone, int nbest ) { 
    317247// Eventually add to global best, but do not go over nbest values 
    318248// Also avoids repeating things, which makes this function awfully slow 
    319            
    320           Array<str_aux> global_best_out; 
    321           if (global_best.length() >= nbest){ 
    322       //logliks = [global_best.loglik]; 
    323        
    324           vec logliks(1);  
    325           logliks(0) =  global_best(0).loglik; 
    326           for (int j = 1; j < global_best.length(); j++) 
    327           logliks = concat(logliks, global_best(j).loglik); 
    328            
    329           int i, addit; 
    330           double loglik = min(logliks, i); 
    331       global_best_out = global_best; 
    332           if (loglik < newone.loglik){ 
    333          //         if ~any(logliks == new.loglik); 
    334           addit=1; 
    335                   
    336                   
    337           
    338                   if (newone.bitstr.length() == 1) { 
    339           for (int j = 0; j < global_best.length(); j++){ 
    340                   for(int i = 0; i < global_best(j).bitstr.length(); i++){ 
    341  
    342                           if (newone.bitstr(0) == global_best(j).bitstr(i)){ 
    343                addit = 0; 
    344                break; 
    345                           } 
    346                           }    
    347           }    
    348                   } 
    349                  if (addit){ 
    350             global_best_out(i) = newone; 
    351             // DEBUGging print: 
    352             // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik); 
    353              }         
    354          }   
    355 } 
    356    else 
    357        global_best_out = concat(global_best, newone); 
    358  
    359   return global_best_out; 
    360    
    361   } 
    362    
    363 */ 
    364     
    365 void add_new(Array<str_aux> &global_best,str_aux newone,int nbest){ 
    366 // Eventually add to global best, but do not go over nbest values 
    367 // Also avoids repeating things, which makes this function awfully slow 
    368            
    369           int  addit, i = 0; 
    370 if (global_best.length() >= nbest){ 
    371       //logliks = [global_best.loglik]; 
    372        
    373            
    374           for (int j = 1; j < global_best.length(); j++){ 
    375                           if (global_best(j).loglik < global_best(i).loglik) { 
    376                                    i = j; 
    377                           } 
    378                   } 
    379           
    380           if (global_best(i).loglik < newone.loglik){ 
    381          //         if ~any(logliks == new.loglik); 
    382           addit=1; 
    383                   
    384                   
    385          //???????????????????????????????????????????? 
    386         if (newone.bitstr.length() == 1) { 
    387           for (int j = 0; j < global_best.length(); j++){ 
    388                   for(int i = 0; i < global_best(j).bitstr.length(); i++){ 
    389  
    390                           if (newone.bitstr(0) == global_best(j).bitstr(i)){ 
    391                addit = 0; 
    392                break; 
    393                           } 
    394                 }    
    395           }    
    396         } 
    397                   
    398                    
    399                 //?????????????????????????????????????????????????   
    400                    
    401                   if (addit){ 
    402             global_best(i) = newone; 
    403             // DEBUGging print: 
    404             // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik); 
    405                   }         
    406          }   
    407 } 
    408    else 
    409        global_best = concat(global_best, newone); 
    410  
    411 } 
    412     
    413    
    414    
    415    
    416    
    417   
    418   str_aux sestrinsert(str_aux in,ivec inserted_elements){ 
     249 
     250        int  addit, i = 0; 
     251        if ( global_best.length() >= nbest ) { 
     252                //logliks = [global_best.loglik]; 
     253 
     254        for ( int j = 1; j < global_best.length(); j++ ) { 
     255                        if ( global_best ( j ).loglik < global_best ( i ).loglik ) { 
     256                                i = j; 
     257                        } 
     258                } 
     259        if ( global_best ( i ).loglik < newone.loglik ) { 
     260                        //         if ~any(logliks == new.loglik); 
     261                        addit = 1; 
     262 
     263 
     264        for (int j = 0; j < global_best.length(); j++){ 
     265                if (newone.bitstr == global_best(j).bitstr){ 
     266         addit = 0; 
     267         break; 
     268                } 
     269        } 
     270        if ( addit ) { 
     271                global_best ( i ) = newone; 
     272                                // DEBUGging print: 
     273                                // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik); 
     274                } 
     275         } 
     276        } else 
     277                global_best = concat ( global_best, newone ); 
     278 
     279} 
     280 
     281 
     282 
     283 
     284 
     285 
     286str_aux sestrinsert ( str_aux in, ivec inserted_elements ) { 
    419287// Moves elements into regressor 
    420    int n_strL = in.strL.length(); 
    421    str_aux out = in; 
    422    for (int j = 0;j < inserted_elements.length(); j++){ 
    423       int f = inserted_elements(j); 
    424       int posit1 = (find(out.strL==1))(0); 
    425       int positf = (find(out.strL==f))(0); 
    426           for (int g = positf; g <= posit1-1; g++ ){ 
    427                    
    428           // BEGIN: We are swapping g and g+1 NOW!!!! 
    429           seswapudl(out.L,  out.d,  g); 
    430           seswapudl(out.L0, out.d0, g); 
    431            
    432            
    433                   int pom_strL = out.strL(g); 
    434                  out.strL(g)= out.strL(g+1); 
    435                  out.strL(g+1) = pom_strL; 
    436                    
    437                   // END 
    438           } 
    439    }          
    440  
    441     out.posit1 = (find(out.strL==1))(0)+1; 
    442    out.strRgr = out.strL.right(n_strL - out.posit1); 
    443    out.strMis = out.strL.left(out.posit1-1); 
    444    str_bitset(out.bitstr,inserted_elements,out.nbits); 
    445     
    446    out.loglik = seloglik1(out); 
    447     
    448     
    449     
    450    return out; 
    451    
    452   } 
    453    
    454   double seloglik2(str_aux in){ 
     288        int n_strL = in.strL.length(); 
     289        str_aux out = in; 
     290        for ( int j = 0; j < inserted_elements.length(); j++ ) { 
     291                int f = inserted_elements ( j ); 
     292                int posit1 = ( find ( out.strL == 1 ) ) ( 0 ); 
     293                int positf = ( find ( out.strL == f ) ) ( 0 ); 
     294                for ( int g = positf; g <= posit1 - 1; g++ ) { 
     295 
     296                        // BEGIN: We are swapping g and g+1 NOW!!!! 
     297                        seswapudl ( out.L,  out.d,  g ); 
     298                        seswapudl ( out.L0, out.d0, g ); 
     299 
     300                        int pom_strL = out.strL ( g ); 
     301                        out.strL ( g ) = out.strL ( g + 1 ); 
     302                        out.strL ( g + 1 ) = pom_strL; 
     303 
     304                        // END 
     305                } 
     306        } 
     307 
     308        out.posit1 = ( find ( out.strL == 1 ) ) ( 0 ) + 1; 
     309        out.strRgr = out.strL.right ( n_strL - out.posit1 ); 
     310        out.strMis = out.strL.left ( out.posit1 - 1 ); 
     311        str_bitset ( out.bitstr, inserted_elements, out.nbits ); 
     312        out.loglik = seloglik1 ( out ); 
     313 
     314        return out; 
     315} 
     316 
     317double seloglik2 ( str_aux in ) { 
    455318// This is the loglikelihood (constant part) - this should be added to 
    456319// everything at the end. It needs some computation, so it is useless to 
    457320// make it for all the stuff 
    458    double logpi = log(pi); 
    459   
    460    double i1 = lgamma(in.nu /2) - 0.5*in.nu *logpi; 
    461    double i0 = lgamma(in.nu0/2) - 0.5*in.nu0*logpi; 
    462    return i1-i0; 
    463   } 
    464  
    465   
    466  
    467    
    468   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*/){ 
     321        double logpi = log ( pi ); 
     322 
     323        double i1 = lgamma ( in.nu / 2 ) - 0.5 * in.nu * logpi; 
     324        double i0 = lgamma ( in.nu0 / 2 ) - 0.5 * in.nu0 * logpi; 
     325        return i1 - i0; 
     326} 
     327 
     328 
     329 
     330 
     331ivec 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*/ ) { 
    469332// see utia_legacy/ticket_12/ implementation and str_test.m 
    470333 
    471 const mat &L = Ld._L(); 
    472 const vec &d = Ld._D(); 
    473  
    474 const mat &L0 = Ld0._L(); 
    475 const vec &d0 = Ld0._D(); 
    476  
    477 int n_data = d.length(); 
    478  
    479 ivec belief_out = find(belief==4)+2; // we are avoiding to put this into regressor 
    480 ivec belief_in  = find(belief==1)+2; // we are instantly keeping this in regressor 
    481  
    482  
    483 str_aux full; 
    484  
    485 full.d0  = d0; 
    486 full.nu0 = nu0; 
    487 full.L0  = L0; 
    488 full.L   = L; 
    489 full.d   = d; 
    490  
    491  
    492  
    493 /*full.d0  = "0.012360650875200       0.975779169502626        1.209840558439000"; 
    494           
    495 full.L0  = "1.0000         0         0;" 
    496            "0.999427690134298    1.0000         0;" 
    497            "0.546994424043659   0.534335486953833    1.0000"; 
    498  
    499  
    500 full.L   = "1.0000         0         0;" 
    501            "0.364376353850780    1.0000         0;" 
    502            "1.222141096674815   1.286534510946323    1.0000"; 
    503  
    504 full.d   =  "0.001610356837691          3.497566869589465       3.236640487818002"; 
    505 */ 
    506  
    507 fstream F; 
    508  
    509  
    510  
    511  
    512 F.open("soubor3.txt", ios::in); 
    513 F << setiosflags(ios::scientific); 
    514 F << setprecision(16); 
    515 F.flags(0x1); 
    516  
    517  
    518 for (int i = 0; i < n_data ; i++){ 
    519 F >> full.d0(i); 
    520 } 
    521  
    522  
    523  
    524 for (int i =0; i < n_data; i++){ 
    525         for (int j = 0 ; j < n_data  ; j++){ 
    526 F >> full.L0(j,i); 
    527 } 
    528 } 
    529  
    530 for (int i = 0; i < n_data ; i++){ 
    531 F >> full.d(i); 
    532 } 
    533  
    534  
    535 for (int i =0; i < n_data; i++){ 
    536         for (int j = 0 ; j < n_data  ; j++){ 
    537 F >> full.L(j,i); 
    538 } 
    539 } 
    540  
    541  
    542 full.nu0 = nu0; 
    543 full.nu  = nu; 
    544 full.strL = linspace(1,n_data); 
    545 full.strRgr = linspace(2,n_data); 
    546 full.strMis = ivec(0);                      
    547 full.posit1 = 1; 
    548 full.bitstr.set_size(n_data-1); 
    549 full.bitstr.clear(); 
    550 str_bitset(full.bitstr,full.strRgr,full.nbits); 
    551 //full.nbits  = std::numeric_lim its<double>::digits-1;  // number of bits available in double 
    552 /*bvec in(n_data-1); 
    553 in.clear(); 
    554 full.bitstr = str_bitset(in,full.strRgr,full.nbits);*/ 
    555  
    556  
    557  
    558  
    559 full.loglik = seloglik1(full);       // % loglikelihood 
    560     
    561  
    562  
     334        const mat &L = Ld._L(); 
     335        const vec &d = Ld._D(); 
     336 
     337        const mat &L0 = Ld0._L(); 
     338        const vec &d0 = Ld0._D(); 
     339 
     340        int n_data = d.length(); 
     341 
     342        ivec belief_out = find ( belief == 4 ) + 2; // we are avoiding to put this into regressor 
     343        ivec belief_in  = find ( belief == 1 ) + 2; // we are instantly keeping this in regressor 
     344 
     345        str_aux full; 
     346 
     347        full.d0  = d0; 
     348        full.nu0 = nu0; 
     349        full.L0  = L0; 
     350        full.L   = L; 
     351        full.d   = d; 
     352        full.nu0 = nu0; 
     353        full.nu  = nu; 
     354        full.strL = linspace ( 1, n_data ); 
     355        full.strRgr = linspace ( 2, n_data ); 
     356        full.strMis = ivec ( 0 ); 
     357        full.posit1 = 1; 
     358        full.bitstr.set_size ( n_data - 1 ); 
     359        full.bitstr.clear(); 
     360        str_bitset ( full.bitstr, full.strRgr, full.nbits ); 
     361        full.loglik = seloglik1 ( full );    // % loglikelihood 
    563362 
    564363 
    565364//% construct full and empty structure 
    566 full = sestrremove(full,belief_out); 
    567 str_aux empty = sestrremove(full,setdiff(full.strRgr,belief_in)); 
    568   
     365        full = sestrremove ( full, belief_out ); 
     366        str_aux empty = sestrremove ( full, setdiff ( full.strRgr, belief_in ) ); 
     367 
    569368//% stopping rule calculation: 
    570369 
    571  
    572  
    573  
    574 bmat local_max(0,0); 
    575 int to, muto = 0; 
    576   
     370        bmat local_max ( 0, 0 ); 
     371        int to, muto = 0; 
     372 
    577373//% statistics: 
    578 //double cputime0 = cputime; ?????????????????????????????????????????????????????????????? 
     374//double cputime0 = cputime; 
    579375//if nargout>=3; 
    580     
    581  CPU_Timer timer; 
    582  timer.start();  
    583   
    584    ivec mutos(max_nrep+2); 
    585    vec maxmutos(max_nrep+2); 
    586    mutos.zeros(); 
    587    maxmutos.zeros(); 
     376 
     377        CPU_Timer timer; 
     378        timer.start(); 
     379        ivec mutos ( max_nrep + 2 ); 
     380        vec maxmutos ( max_nrep + 2 ); 
     381        mutos.zeros(); 
     382        maxmutos.zeros(); 
    588383 
    589384 
    590385//end; 
    591386//% ---------------------- 
    592   
     387 
    593388//% For stopping-rule calculation 
    594389//%so       = 2^(n_data -1-length(belief_in)- length(belief_out)); % do we use this ? 
    595390//% ---------------------- 
    596   
    597 ivec all_str = linspace(1,n_data); 
    598   
    599 Array<str_aux> global_best(1); 
    600 global_best(0) = full; 
     391 
     392        ivec all_str = linspace ( 1, n_data ); 
     393        Array<str_aux> global_best ( 1 ); 
     394        global_best ( 0 ) = full; 
    601395 
    602396 
    603397//% MAIN LOOP is here. 
    604398 
    605 str_aux   best; 
    606 for (int n_start = -1; n_start <= max_nrep; n_start++){ 
    607    str_aux  last,best; 
    608     
    609    to = n_start+2; 
    610     
    611    if (n_start == -1){ 
    612       //% start from the full structure 
    613       last = full; 
    614    } 
    615    else {if (n_start == 0) 
    616       //% start from the empty structure 
    617       last = empty;      
    618         
    619       else{ 
    620       //% start from random structure    
    621                   
    622                 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 
    623                 last = sestrremove(full,setdiff(all_str,::concat<int>(::concat<int>(1 ,last_str), empty.strRgr))); 
    624        
    625           } 
    626    } 
    627    //% DEBUGging print: 
    628    //%fprintf('STRUCTURE generated            in loop %2i was %s\n', n_start, strPrintstr(last)); 
    629     
    630    //% The loop is repeated until likelihood stops growing (break condition 
    631    //% used at the end; 
    632    
    633     
    634    while (1){ 
    635       //% This structure is going to hold the best elements 
    636       best = last; 
    637       //% Nesting by removing elements (enpoorment) 
    638       ivec  removed_items = setdiff(last.strRgr,belief_in); 
    639            
    640 ivec removed_item; 
    641           str_aux newone; 
    642            
    643           for (int i = 0; i < removed_items.length(); i++){ 
    644         removed_item = vec_1(removed_items(i)); 
    645             newone = sestrremove(last,removed_item);//????????????????????????????? 
    646                 if (nbest>1){ 
    647           add_new(global_best,newone,nbest); 
    648                 } 
    649                 if (newone.loglik>best.loglik){ 
    650            best = newone; 
    651                 } 
    652           } 
    653       //% Nesting by adding elements (enrichment) 
    654       ivec added_items = setdiff(last.strMis,belief_out); 
    655           ivec added_item; 
    656            
    657           for (int j = 0; j < added_items.length(); j++){ 
    658         added_item = vec_1(added_items(j)); 
    659             newone = sestrinsert(last,added_item); 
    660                 if (nbest>1){ 
    661            add_new(global_best,newone,nbest); 
    662                 } 
    663                 if (newone.loglik>best.loglik){ 
    664            best = newone; 
    665                 } 
    666      } 
    667  
    668     
    669  
    670     
    671  
    672       //% Break condition if likelihood does not change. 
    673       if (best.loglik <= last.loglik) 
    674           break; 
    675       else 
    676           //% Making best structure last structure. 
    677           last = best; 
    678       
    679      
    680    } 
    681   
    682  
    683   
    684  
    685  
    686    // % DEBUGging print: 
    687    //%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); 
    688     
    689    //% Collecting of the best structure in case we don't need the second parameter 
    690    if (nbest<=1){ 
    691            if (best.loglik > global_best(0).loglik){ 
    692          global_best = best; 
    693    } 
    694    } 
    695   
    696    //% uniqueness of the structure found 
    697    int append = 1; 
    698     
    699  
    700    for(int j = 0; j  < local_max.rows() ; j++){ 
    701            if (best.bitstr == local_max.get_row(j)){ 
    702          append = 0; 
    703                  break; 
    704            } 
    705    } 
    706  
    707     
    708    if (append){ 
    709       local_max.append_row(best.bitstr); 
    710       muto = muto + 1; 
    711    }  
    712     
    713    //% stopping rule: 
    714    double maxmuto = (to-order_k-1)/lambda-to+1; 
    715    if (to>2){ 
    716            if (maxmuto>=muto){ 
    717           //% fprintf('*'); 
    718           break; 
    719            } 
    720    }       
    721   
    722    // do statistics if necessary: 
    723    //if (nargout>=3){ 
    724       mutos(to-1)    = muto; 
    725       maxmutos(to-1) = maxmuto; 
    726    //} 
    727 } 
    728   
     399        str_aux   best; 
     400        for ( int n_start = -1; n_start <= max_nrep; n_start++ ) { 
     401                str_aux  last, best; 
     402 
     403                to = n_start + 2; 
     404 
     405                if ( n_start == -1 ) { 
     406                        //% start from the full structure 
     407                        last = full; 
     408                } else { 
     409                        if ( n_start == 0 ) 
     410                                //% start from the empty structure 
     411                                last = empty; 
     412 
     413                        else { 
     414                                //% start from random structure 
     415             
     416                 
     417                 
     418                          ivec last_str = find ( to_bvec( concat( 0, floor_i ( 2 * randun ( n_data - 1 )) ) ) ) + 1;// this creates random vector consisting of indexes, and sorted 
     419                          last = sestrremove ( full, setdiff ( all_str, concat( concat( 1 , last_str ), empty.strRgr ) ) ); 
     420 
     421                        } 
     422                } 
     423                //% DEBUGging print: 
     424                //%fprintf('STRUCTURE generated            in loop %2i was %s\n', n_start, strPrintstr(last)); 
     425 
     426                //% The loop is repeated until likelihood stops growing (break condition 
     427                //% used at the end; 
     428 
     429 
     430                while ( 1 ) { 
     431                        //% This structure is going to hold the best elements 
     432                        best = last; 
     433                        //% Nesting by removing elements (enpoorment) 
     434                        ivec  removed_items = setdiff ( last.strRgr, belief_in ); 
     435 
     436                        ivec removed_item; 
     437                        str_aux newone; 
     438 
     439                        for ( int i = 0; i < removed_items.length(); i++ ) { 
     440                                removed_item = vec_1 ( removed_items ( i ) ); 
     441                                newone = sestrremove ( last, removed_item ); 
     442                                if ( nbest > 1 ) { 
     443                                        add_new ( global_best, newone, nbest ); 
     444                                } 
     445                                if ( newone.loglik > best.loglik ) { 
     446                                        best = newone; 
     447                                } 
     448                        } 
     449                        //% Nesting by adding elements (enrichment) 
     450                        ivec added_items = setdiff( last.strMis, belief_out ); 
     451                        ivec added_item; 
     452 
     453                        for ( int j = 0; j < added_items.length(); j++ ) { 
     454                                added_item = vec_1 ( added_items ( j ) ); 
     455                                newone = sestrinsert ( last, added_item ); 
     456                                if ( nbest > 1 ) { 
     457                                        add_new ( global_best, newone, nbest ); 
     458                                } 
     459                                if ( newone.loglik > best.loglik ) { 
     460                                        best = newone; 
     461                                } 
     462                        } 
     463 
     464 
     465                        //% Break condition if likelihood does not change. 
     466                        if ( best.loglik <= last.loglik ) 
     467                                break; 
     468                        else 
     469                                //% Making best structure last structure. 
     470                                last = best; 
     471                        } 
     472 
     473 
     474                // % DEBUGging print: 
     475                //%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); 
     476 
     477                //% Collecting of the best structure in case we don't need the second parameter 
     478                if ( nbest <= 1 ) { 
     479                        if ( best.loglik > global_best ( 0 ).loglik ) { 
     480                                global_best = best; 
     481                        } 
     482                } 
     483 
     484                //% uniqueness of the structure found 
     485                int append = 1; 
     486 
     487                for ( int j = 0; j  < local_max.rows() ; j++ ) { 
     488                        if ( best.bitstr == local_max.get_row ( j ) ) { 
     489                                append = 0; 
     490                                break; 
     491                        } 
     492                } 
     493 
     494 
     495                if ( append ) { 
     496                        local_max.append_row ( best.bitstr ); 
     497                        muto = muto + 1; 
     498                } 
     499 
     500                //% stopping rule: 
     501                double maxmuto = ( to - order_k - 1 ) / lambda - to + 1; 
     502                if ( to > 2 ) { 
     503                        if ( maxmuto >= muto ) { 
     504                                //% fprintf('*'); 
     505                                break; 
     506                        } 
     507                } 
     508 
     509                // do statistics if necessary: 
     510                //if (nargout>=3){ 
     511                mutos ( to - 1 )    = muto; 
     512                maxmutos ( to - 1 ) = maxmuto; 
     513                //} 
     514        } 
     515 
    729516//% Aftermath: The best structure was in: global_best 
    730   
     517 
    731518//% Updating loglikelihoods: we have to add the constant stuff 
    732519 
    733  
    734  
    735 for (int f=0 ; f <global_best.length(); f++){ 
    736    global_best(f).loglik = global_best(f).loglik + seloglik2(global_best(f)); 
    737 } 
    738   
    739 /*for f=1:length(global_best); 
    740    global_best(f).loglik = global_best(f).loglik + seloglik2(global_best(f)); 
    741 end;*/ 
    742  
     520        for ( int f = 0 ; f < global_best.length(); f++ ) { 
     521                global_best ( f ).loglik = global_best ( f ).loglik + seloglik2 ( global_best ( f ) ); 
     522        } 
    743523 
    744524//% Making first output parameter: 
    745525 
    746 int max_i = 0; 
    747 for (int j = 1; j < global_best.length(); j++) 
    748   if  (global_best(max_i).loglik < (global_best(j).loglik)) max_i = j; 
    749  
    750 best = global_best(max_i); 
    751 cout << endl << endl << endl << endl << best.strRgr; 
    752 cout << endl << endl << endl << endl << best.strRgr;  
     526        int max_i = 0; 
     527        for ( int j = 1; j < global_best.length(); j++ ) 
     528                if ( global_best ( max_i ).loglik < ( global_best ( j ).loglik ) ) max_i = j; 
     529 
     530        best = global_best ( max_i ); 
     531 
    753532//% Making the second output parameter 
    754533 
    755 vec logliks(global_best.length()); 
    756 for (int j = 0; j < logliks.length(); j++) 
    757  logliks(j) = global_best(j).loglik; 
    758  
    759 ivec i = sort_index(logliks); 
    760 rgrsout.set_length(global_best.length()); 
    761  
    762 for (int j = global_best.length() - 1; j >= 0; j--) 
    763 rgrsout(j) = global_best(i(j)); 
    764   
     534        vec logliks ( global_best.length() ); 
     535        for ( int j = 0; j < logliks.length(); j++ ) 
     536                logliks ( j ) = global_best ( j ).loglik; 
     537 
     538        ivec i = sort_index ( logliks ); 
     539        rgrsout.set_length ( global_best.length() ); 
     540 
     541        for ( int j = global_best.length() - 1; j >= 0; j-- ) 
     542                rgrsout ( j ) = global_best ( i ( j ) ); 
     543 
    765544//if (nargout>=3); 
    766     
    767     
    768 str_statistics statistics; 
    769  
    770    statistics.allstrs = 2^(n_data -1-length(belief_in) - length(belief_out)); 
    771    statistics.nrand   = to-2; 
    772    statistics.unique  = muto; 
    773    statistics.to  = to; 
    774    statistics.cputime_seconds = timer.get_time(); 
    775    statistics.itemspeed       = statistics.to / statistics.cputime_seconds; 
    776    statistics.muto = muto; 
    777    statistics.mutos = mutos; 
    778    statistics.maxmutos = maxmutos; 
     545 
     546        str_statistics statistics; 
     547 
     548        statistics.allstrs = 2 ^ ( n_data - 1 - length ( belief_in ) - length ( belief_out ) ); 
     549        statistics.nrand   = to - 2; 
     550        statistics.unique  = muto; 
     551        statistics.to  = to; 
     552        statistics.cputime_seconds = timer.get_time(); 
     553        statistics.itemspeed       = statistics.to / statistics.cputime_seconds; 
     554        statistics.muto = muto; 
     555        statistics.mutos = mutos; 
     556        statistics.maxmutos = maxmutos; 
    779557//end; 
    780558 
    781 return best.strRgr; 
    782  
    783 }  
    784  
    785 #ifdef LADIM 
    786 //% randun seed stuff: 
    787 //%randn('seed',SEED); 
    788   
    789 //% --------------------- END of MAIN program -------------------- 
    790   
    791 % This is needed for bitstr manipulations 
    792 /*function out = str_bitset(in,ns,nbits) 
    793    out = in; 
    794    for n = ns; 
    795       index = 1+floor((n-2)/nbits); 
    796       bitindex = 1+rem(n-2,nbits); 
    797       out(index) = bitset(out(index),bitindex); 
    798    end;   
    799 function out = str_bitres(in,ns,nbits) 
    800    out = in; 
    801    for n = ns; 
    802       index = 1+floor((n-2)/nbits); 
    803       bitindex = 1+rem(n-2,nbits); 
    804       mask = bitset(0,bitindex); 
    805       out(index) = bitxor(bitor(out(index),mask),mask); 
    806    end;*/ 
    807   
    808 function out = strPrintstr(in) 
    809    out = '0'; 
    810    nbits = in.nbits; 
    811    for f = 2:length(in.d0); 
    812       index = 1+floor((f-2)/nbits); 
    813       bitindex = 1+rem(f-2,nbits); 
    814       if bitget(in.bitstr(index),bitindex); 
    815           out(f) = '1'; 
    816       else; 
    817           out(f) = '0'; 
    818       end; 
    819    end; 
    820   
    821 /*function global_best_out = add_new(global_best,new,nbest) 
    822 % Eventually add to global best, but do not go over nbest values 
    823 % Also avoids repeating things, which makes this function awfully slow 
    824    if length(global_best)>=nbest; 
    825       logliks = [global_best.loglik]; 
    826       [loglik i] = min(logliks); 
    827       global_best_out = global_best; 
    828       if loglik<new.loglik; 
    829          %         if ~any(logliks == new.loglik); 
    830          addit=1; 
    831          for f = [global_best.bitstr]; 
    832             if f == new.bitstr; 
    833                addit = 0; 
    834                break; 
    835             end; 
    836          end;          
    837          if addit; 
    838             global_best_out(i) = new; 
    839             % DEBUGging print: 
    840             % fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik); 
    841          end;          
    842       end; 
    843    else; 
    844       global_best_out = [global_best new]; 
    845    end;*/ 
    846   
    847 /*function out = sestrremove(in,removed_elements); 
    848 % Removes elements from regressor 
    849    n_strL = length(in.strL); 
    850    out = in; 
    851    for f=removed_elements; 
    852       posit1 = find(out.strL==1); 
    853       positf = find(out.strL==f); 
    854       for g=(positf-1):-1:posit1; 
    855          % BEGIN: We are swapping g and g+1 NOW!!!! 
    856          [out.L, out.d]   = seswapudl(out.L, out.d, g); 
    857          [out.L0, out.d0]   = seswapudl(out.L0, out.d0, g); 
    858          out.strL([g g+1]) = [out.strL(g+1) out.strL(g)]; 
    859          % END 
    860       end; 
    861    end; 
    862    out.posit1 = find(out.strL==1); 
    863    out.strRgr = out.strL((out.posit1+1):n_strL); 
    864    out.strMis = out.strL(1:(out.posit1-1)); 
    865    out.bitstr = str_bitres(out.bitstr,removed_elements,out.nbits); 
    866    out.loglik = seloglik1(out);*/ 
    867     
    868 /*function out = sestrinsert(in,inserted_elements); 
    869 % Moves elements into regressor 
    870    n_strL = length(in.strL); 
    871    out = in; 
    872    for f=inserted_elements; 
    873       posit1 = find(out.strL==1); 
    874       positf = find(out.strL==f); 
    875       for g=positf:(posit1-1); 
    876           % BEGIN: We are swapping g and g+1 NOW!!!! 
    877           [out.L,  out.d]   = seswapudl(out.L,  out.d,  g); 
    878           [out.L0, out.d0]  = seswapudl(out.L0, out.d0, g); 
    879           out.strL([g g+1]) = [out.strL(g+1) out.strL(g)]; 
    880           % END 
    881       end; 
    882    end;          
    883    out.posit1 = find(out.strL==1); 
    884    out.strRgr = out.strL((out.posit1+1):n_strL); 
    885    out.strMis = out.strL(1:(out.posit1-1)); 
    886    out.bitstr = str_bitset(out.bitstr,inserted_elements,out.nbits); 
    887    out.loglik = seloglik1(out);*/ 
    888   
    889 % 
    890 % seloglik_real = seloglik1 + seloglik2 
    891 % 
    892   
    893 /*function l = seloglik1(in) 
    894 % This is the loglikelihood (non-constant part) - this should be used in 
    895 % frequent computation 
    896    len = length(in.d); 
    897    p1  = in.posit1; 
    898        
    899    i1 = -0.5*in.nu *log(in.d (p1)) -0.5*sum(log(in.d ((p1+1):len))); 
    900    i0 = -0.5*in.nu0*log(in.d0(p1)) -0.5*sum(log(in.d0((p1+1):len)));    
    901    l  = i1-i0; 
    902     
    903    % DEBUGGing print: 
    904    % fprintf('SELOGLIK1: str=%s loglik=%g\n', strPrintstr(in), l);*/ 
    905   
    906   
    907 function l = seloglik2(in) 
    908 % This is the loglikelihood (constant part) - this should be added to 
    909 % everything at the end. It needs some computation, so it is useless to 
    910 % make it for all the stuff 
    911    logpi = log(pi); 
    912   
    913    i1 = gammaln(in.nu /2) - 0.5*in.nu *logpi; 
    914    i0 = gammaln(in.nu0/2) - 0.5*in.nu0*logpi; 
    915    l  = i1-i0; 
    916   
    917   
    918 /*function [Lout, dout] = seswapudl(L,d,i); 
    919 %SESWAPUDL swaps information matrix in decomposition V=L^T diag(d) L 
    920 % 
    921 %  [Lout, dout] = seswapudl(L,d,i); 
    922 % 
    923 % L : lower triangular matrix with 1's on diagonal of the decomposistion 
    924 % d : diagonal vector of diagonal matrix of the decomposition 
    925 % i : index of line to be swapped with the next one  
    926 % Lout : output lower triangular matrix 
    927 % dout : output diagional vector of diagonal matrix D 
    928 % 
    929 % Description: 
    930 %  Lout' * diag(dout) * Lout = P(i,i+1) * L' * diag(d) * L * P(i,i+1); 
    931 %  
    932 %  Where permutation matrix P(i,j) permutates columns if applied from the 
    933 %  right and line if applied from the left. 
    934 %    
    935 % Note: naming: 
    936 %       se = structure estimation 
    937 %       lite = light, simple 
    938 %       udl = U*D*L, or more precisely, L'*D*L, also called as ld 
    939 %    
    940 % Design  : L. Tesar 
    941 % Updated : Feb 2003 
    942 % Project : post-ProDaCTool 
    943 % Reference: sedydr 
    944     
    945 j = i+1; 
    946   
    947 pomd = d(i); 
    948 d(i) = d(j); 
    949 d(j) = pomd; 
    950   
    951 pomL   = L(i,:); 
    952 L(i,:) = L(j,:); 
    953 L(j,:) = pomL; 
    954   
    955 pomL   = L(:,i); 
    956 L(:,i) = L(:,j); 
    957 L(:,j) = pomL; 
    958   
    959 % We must be working with LINES of matrix L ! 
    960   
    961 r  = L(i,:)';  
    962 f  = L(j,:)'; 
    963 Dr = d(i); 
    964 Df = d(j); 
    965   
    966 [r, f, Dr, Df] = sedydr(r, f, Dr, Df, j); 
    967   
    968 r0 = r(i); 
    969 Dr = Dr*r0*r0; 
    970 r  = r/r0; 
    971   
    972 L(i,:) = r';  
    973 L(j,:) = f'; 
    974 d(i)   = Dr; 
    975 d(j)   = Df; 
    976   
    977 L(i,i) = 1; 
    978 L(j,j) = 1; 
    979   
    980 Lout = L; 
    981 dout = d;*/ 
    982   
    983 /*function [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh); 
    984 %SEDYDR dyadic reduction, performs transformation of sum of 2 dyads 
    985 % 
    986 % [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh); 
    987 % [rout, fout, Drout, Dfout] = sedydr(r,f,Dr,Df,R); 
    988 % 
    989 % Description: dyadic reduction, performs transformation of sum of  
    990 %  2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed by R is zeroed 
    991 % 
    992 %     r    : column vector of reduced dyad 
    993 %     f    : column vector of reducing dyad 
    994 %     Dr   : scalar with weight of reduced dyad 
    995 %     Df   : scalar with weight of reducing dyad 
    996 %     R    : scalar number giving 1 based index to the element of r, 
    997 %            which is to be reduced to  
    998 %            zero; the corresponding element of f is assumed to be 1. 
    999 %     jl   : lower index of the range within which the dyads are  
    1000 %            modified (can be omitted, then everything is updated) 
    1001 %     jh   : upper index of the range within which the dyads are  
    1002 %            modified (can be omitted then everything is updated) 
    1003 %     rout,fout,Drout,dfout : resulting two dyads 
    1004 %     kr   : coefficient used in the transformation of r 
    1005 %            rnew = r + kr*f 
    1006 % 
    1007 % Description: dyadic reduction, performs transformation of sum of  
    1008 %            2 dyads r*Dr*r'+ f*Df*f' so that the element of r indexed by R is zeroed 
    1009 % Remark1: Constant mzero means machine zero and should be modified 
    1010 %           according to the precision of particular machine 
    1011 % Remark2: jl and jh are, in fact, obsolete. It takes longer time to 
    1012 %           compute them compared to plain version. The reason is that we 
    1013 %           are doing vector operations in m-file. Other reason is that 
    1014 %           we need to copy whole vector anyway. It can save half of time for 
    1015 %           c-file, if you use it correctly. (please do tests) 
    1016 % 
    1017 % Note: naming: 
    1018 %       se = structure estimation 
    1019 %       dydr = dyadic reduction 
    1020 % 
    1021 % Original Fortran design: V. Peterka 17-7-89 
    1022 % Modified for c-language: probably R. Kulhavy 
    1023 % Modified for m-language: L. Tesar 2/2003 
    1024 % Updated: Feb 2003 
    1025 % Project: post-ProDaCTool 
    1026 % Reference: none 
    1027     
    1028 if nargin<6; 
    1029    update_whole=1; 
    1030 else 
    1031    update_whole=0; 
    1032 end; 
    1033   
    1034 mzero = 1e-32; 
    1035   
    1036 if Dr<mzero; 
    1037    Dr=0; 
    1038 end; 
    1039   
    1040 r0   = r(R); 
    1041 kD   = Df; 
    1042 kr   = r0 * Dr; 
    1043 Dfout   = kD + r0 * kr; 
    1044   
    1045 if Dfout > mzero; 
    1046     kD = kD / Dfout; 
    1047     kr = kr / Dfout; 
    1048 else; 
    1049     kD = 1; 
    1050     kr = 0; 
    1051 end; 
    1052   
    1053 Drout = Dr * kD; 
    1054   
    1055 % Try to uncomment marked stuff (*) if in numerical problems, but I don't 
    1056 % think it can make any difference for normal healthy floating-point unit 
    1057 if update_whole; 
    1058    rout = r - r0*f; 
    1059 %   rout(R) = 0;   % * could be needed for some nonsense cases(or numeric reasons?), normally not 
    1060    fout = f + kr*rout; 
    1061 %   fout(R) = 1;   % * could be needed for some nonsense cases(or numeric reasons?), normally not 
    1062 else;   
    1063    rout = r; 
    1064    fout = f; 
    1065    rout(jl:jh) = r(jl:jh) - r0 * f(jl:jh); 
    1066    rout(R) = 0;    
    1067    fout(jl:jh) = f(jl:jh) + kr * rout(jl:jh); 
    1068 end;*/ 
    1069   
    1070   
    1071   
    1072 #endif 
    1073   
    1074   
    1075 } 
    1076  
     559        return best.strRgr; 
     560 
     561} 
     562 
     563 
     564 
     565 
     566 
     567 
     568} 
     569