Changeset 684

Show
Ignore:
Timestamp:
10/27/09 21:15:03 (14 years ago)
Author:
sarka
Message:

function straux1 converted from straux1.m

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • library/bdm/estim/arx_straux.cpp

    r661 r684  
    1111namespace bdm { 
    1212 
    13  
    14  
    15  
    16 /* 
     13   
     14 
     15         
     16/*       
    1717struct str_aux { 
    1818        vec d0; 
     
    2727        int posit1;                // regressand position 
    2828        int nbits;                                 // number of bits available in double 
    29         bvec bitstr; 
     29        bvec bitstr;  
    3030        double loglik;          // loglikelihood 
    3131  }; 
    32  
     32   
    3333 */ 
    3434 
    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 ) { 
     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 
     64double seloglik1(str_aux in){ 
    6565// This is the loglikelihood (non-constant part) - this should be used in 
    6666// 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  
     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 
     79void 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; 
     126else 
     127   update_whole=0; 
     128end;*/ 
     129  
     130double mzero = 1e-32; 
     131  
     132if (Dr<mzero){ 
     133   Dr=0; 
     134} 
     135  
     136double r0   = r(R,0); 
     137double kD   = Df; 
     138double kr   = r0 * Dr; 
     139 
     140 
     141Df   = kD + r0 * kr; 
     142  
     143if (Df > mzero){ 
     144    kD = kD / Df; 
     145    kr = kr / Df; 
     146}else{ 
     147    kD = 1; 
     148    kr = 0; 
     149} 
     150  
     151Dr = Dr * kD; 
     152  
    153153// Try to uncomment marked stuff (*) if in numerical problems, but I don't 
    154154// think it can make any difference for normal healthy floating-point unit 
    155155//if update_whole; 
    156         r = r - r0 * f; 
     156   r = r - r0*f; 
    157157//   rout(R) = 0;   // * could be needed for some nonsense cases(or numeric reasons?), normally not 
    158         f = f + kr * r; 
     158   f = f + kr*r; 
    159159//   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  
     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); 
     166end;*/ 
     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  
     198int j = i+1; 
     199  
     200double pomd = d(i); 
     201d(i) = d(j); 
     202d(j) = pomd; 
     203  
     204/*vec pomL   = L.get_row(i); 
     205L.set_row(i, L.get_row(j));    
     206L.set_row(j,pomL);*/ 
     207  
     208L.swap_rows(i,j); 
     209L.swap_cols(i,j); 
     210 
     211/*pomL   = L.get_col(i); 
     212L.set_col(i, L.get_col(j)); 
     213L.set_col(j,pomL);*/ 
     214  
    215215//% 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(); 
     216  
     217 
     218 
     219mat r = L.get_row(i); 
     220r = r.transpose(); 
     221r = r.transpose(); 
    222222//???????????????? 
    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 ) { 
     223mat f = L.get_row(j); 
     224f = f.transpose(); 
     225f = f.transpose(); 
     226 
     227 
     228 
     229 
     230double Dr = d(i); 
     231double Df = d(j); 
     232  
     233sedydr(r, f, Dr, Df, j); 
     234  
     235 
     236 
     237double r0 = r(i,0); 
     238Dr = Dr*r0*r0; 
     239r  = r/r0; 
     240 
     241 
     242 
     243mat pom_mat = r.transpose(); 
     244L.set_row(i, pom_mat.get_row(0));  
     245pom_mat = f.transpose(); 
     246L.set_row(j, pom_mat.get_row(0));  
     247 
     248d(i)   = Dr; 
     249d(j)   = Df; 
     250  
     251L(i,i) = 1; 
     252L(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   
     268str_aux sestrremove(str_aux in,ivec removed_elements){ 
    269269//% 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 
     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); 
    288307                } 
    289308        } 
    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  
     309  return a; 
     310  } 
     311 
     312 
     313   
    314314/* 
    315  
     315   
    316316  Array<str_aux> add_new(Array<str_aux> global_best,str_aux newone,int nbest){ 
    317317// Eventually add to global best, but do not go over nbest values 
    318318// Also avoids repeating things, which makes this function awfully slow 
    319  
     319           
    320320          Array<str_aux> global_best_out; 
    321321          if (global_best.length() >= nbest){ 
    322322      //logliks = [global_best.loglik]; 
    323  
    324           vec logliks(1); 
     323       
     324          vec logliks(1);  
    325325          logliks(0) =  global_best(0).loglik; 
    326326          for (int j = 1; j < global_best.length(); j++) 
    327327          logliks = concat(logliks, global_best(j).loglik); 
    328  
     328           
    329329          int i, addit; 
    330330          double loglik = min(logliks, i); 
     
    333333         //         if ~any(logliks == new.loglik); 
    334334          addit=1; 
    335  
    336  
    337  
     335                  
     336                  
     337          
    338338                  if (newone.bitstr.length() == 1) { 
    339339          for (int j = 0; j < global_best.length(); j++){ 
     
    344344               break; 
    345345                          } 
    346                           } 
    347           } 
     346                          }    
     347          }    
    348348                  } 
    349349                 if (addit){ 
     
    351351            // DEBUGging print: 
    352352            // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik); 
    353              } 
    354          } 
     353             }         
     354         }   
    355355} 
    356356   else 
     
    358358 
    359359  return global_best_out; 
    360  
     360   
    361361  } 
    362  
     362   
    363363*/ 
    364  
    365 void add_new ( Array<str_aux> &global_best, str_aux newone, int nbest ) { 
     364    
     365void add_new(Array<str_aux> &global_best,str_aux newone,int nbest){ 
    366366// Eventually add to global best, but do not go over nbest values 
    367367// 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         } else 
    408                 global_best = concat ( global_best, newone ); 
    409  
    410 } 
    411  
    412  
    413  
    414  
    415  
    416  
    417 str_aux sestrinsert ( str_aux in, ivec inserted_elements ) { 
     368           
     369          int  addit, i = 0; 
     370if (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){ 
    418419// Moves elements into regressor 
    419         int n_strL = in.strL.length(); 
    420         str_aux out = in; 
    421         for ( int j = 0; j < inserted_elements.length(); j++ ) { 
    422                 int f = inserted_elements ( j ); 
    423                 int posit1 = ( find ( out.strL == 1 ) ) ( 0 ); 
    424                 int positf = ( find ( out.strL == f ) ) ( 0 ); 
    425                 for ( int g = positf; g <= posit1 - 1; g++ ) { 
    426  
    427                         // BEGIN: We are swapping g and g+1 NOW!!!! 
    428                         seswapudl ( out.L,  out.d,  g ); 
    429                         seswapudl ( out.L0, out.d0, g ); 
    430  
    431  
    432                         int pom_strL = out.strL ( g ); 
    433                         out.strL ( g ) = out.strL ( g + 1 ); 
    434                         out.strL ( g + 1 ) = pom_strL; 
    435  
    436                         // END 
    437                 } 
    438         } 
    439  
    440         out.posit1 = ( find ( out.strL == 1 ) ) ( 0 ) + 1; 
    441         out.strRgr = out.strL.right ( n_strL - out.posit1 ); 
    442         out.strMis = out.strL.left ( out.posit1 - 1 ); 
    443         str_bitset ( out.bitstr, inserted_elements, out.nbits ); 
    444  
    445         out.loglik = seloglik1 ( out ); 
    446  
    447  
    448  
    449         return out; 
    450  
    451 } 
    452  
    453 double seloglik2 ( str_aux in ) { 
     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){ 
    454455// This is the loglikelihood (constant part) - this should be added to 
    455456// everything at the end. It needs some computation, so it is useless to 
    456457// make it for all the stuff 
    457         double logpi = log ( pi ); 
    458  
    459         double i1 = lgamma ( in.nu / 2 ) - 0.5 * in.nu * logpi; 
    460         double i0 = lgamma ( in.nu0 / 2 ) - 0.5 * in.nu0 * logpi; 
    461         return i1 - i0; 
    462 } 
    463  
    464  
    465  
    466  
    467 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*/ ) { 
     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*/){ 
    468469// see utia_legacy/ticket_12/ implementation and str_test.m 
    469470 
    470         const mat &L = Ld._L(); 
    471         const vec &d = Ld._D(); 
    472  
    473         const mat &L0 = Ld0._L(); 
    474         const vec &d0 = Ld0._D(); 
    475  
    476         int n_data = d.length(); 
    477  
    478         ivec belief_out = find ( belief == 4 ) + 2; // we are avoiding to put this into regressor 
    479         ivec belief_in  = find ( belief == 1 ) + 2; // we are instantly keeping this in regressor 
    480  
    481  
    482         str_aux full; 
    483  
    484         full.d0  = d0; 
    485         full.nu0 = nu0; 
    486         full.L0  = L0; 
    487         full.L   = L; 
    488         full.d   = d; 
    489  
    490  
    491  
    492         /*full.d0  = "0.012360650875200       0.975779169502626        1.209840558439000"; 
    493  
    494         full.L0  = "1.0000         0         0;" 
    495                    "0.999427690134298    1.0000         0;" 
    496                    "0.546994424043659   0.534335486953833    1.0000"; 
    497  
    498  
    499         full.L   = "1.0000         0         0;" 
    500                    "0.364376353850780    1.0000         0;" 
    501                    "1.222141096674815   1.286534510946323    1.0000"; 
    502  
    503         full.d   =  "0.001610356837691          3.497566869589465       3.236640487818002"; 
    504         */ 
    505  
    506         fstream F; 
    507  
    508  
    509  
    510  
    511         F.open ( "soubor3.txt", ios::in ); 
    512         F << setiosflags ( ios::scientific ); 
    513         F << setprecision ( 16 ); 
    514 //      F.flags ( 0x1 ); 
    515  
    516  
    517         for ( int i = 0; i < n_data ; i++ ) { 
    518                 F >> full.d0 ( i ); 
    519         } 
    520  
    521  
    522  
    523         for ( int i = 0; i < n_data; i++ ) { 
    524                 for ( int j = 0 ; j < n_data  ; j++ ) { 
    525                         F >> full.L0 ( j, i ); 
    526                 } 
    527         } 
    528  
    529         for ( int i = 0; i < n_data ; i++ ) { 
    530                 F >> full.d ( i ); 
    531         } 
    532  
    533  
    534         for ( int i = 0; i < n_data; i++ ) { 
    535                 for ( int j = 0 ; j < n_data  ; j++ ) { 
    536                         F >> full.L ( j, i ); 
    537                 } 
    538         } 
    539  
    540  
    541         full.nu0 = nu0; 
    542         full.nu  = nu; 
    543         full.strL = linspace ( 1, n_data ); 
    544         full.strRgr = linspace ( 2, n_data ); 
    545         full.strMis = ivec ( 0 ); 
    546         full.posit1 = 1; 
    547         full.bitstr.set_size ( n_data - 1 ); 
    548         full.bitstr.clear(); 
    549         str_bitset ( full.bitstr, full.strRgr, full.nbits ); 
     471const mat &L = Ld._L(); 
     472const vec &d = Ld._D(); 
     473 
     474const mat &L0 = Ld0._L(); 
     475const vec &d0 = Ld0._D(); 
     476 
     477int n_data = d.length(); 
     478 
     479ivec belief_out = find(belief==4)+2; // we are avoiding to put this into regressor 
     480ivec belief_in  = find(belief==1)+2; // we are instantly keeping this in regressor 
     481 
     482 
     483str_aux full; 
     484 
     485full.d0  = d0; 
     486full.nu0 = nu0; 
     487full.L0  = L0; 
     488full.L   = L; 
     489full.d   = d; 
     490 
     491 
     492 
     493/*full.d0  = "0.012360650875200       0.975779169502626        1.209840558439000"; 
     494          
     495full.L0  = "1.0000         0         0;" 
     496           "0.999427690134298    1.0000         0;" 
     497           "0.546994424043659   0.534335486953833    1.0000"; 
     498 
     499 
     500full.L   = "1.0000         0         0;" 
     501           "0.364376353850780    1.0000         0;" 
     502           "1.222141096674815   1.286534510946323    1.0000"; 
     503 
     504full.d   =  "0.001610356837691          3.497566869589465       3.236640487818002"; 
     505*/ 
     506 
     507fstream F; 
     508 
     509 
     510 
     511 
     512F.open("soubor3.txt", ios::in); 
     513F << setiosflags(ios::scientific); 
     514F << setprecision(16); 
     515F.flags(0x1); 
     516 
     517 
     518for (int i = 0; i < n_data ; i++){ 
     519F >> full.d0(i); 
     520} 
     521 
     522 
     523 
     524for (int i =0; i < n_data; i++){ 
     525        for (int j = 0 ; j < n_data  ; j++){ 
     526F >> full.L0(j,i); 
     527} 
     528} 
     529 
     530for (int i = 0; i < n_data ; i++){ 
     531F >> full.d(i); 
     532} 
     533 
     534 
     535for (int i =0; i < n_data; i++){ 
     536        for (int j = 0 ; j < n_data  ; j++){ 
     537F >> full.L(j,i); 
     538} 
     539} 
     540 
     541 
     542full.nu0 = nu0; 
     543full.nu  = nu; 
     544full.strL = linspace(1,n_data); 
     545full.strRgr = linspace(2,n_data); 
     546full.strMis = ivec(0);                      
     547full.posit1 = 1; 
     548full.bitstr.set_size(n_data-1); 
     549full.bitstr.clear(); 
     550str_bitset(full.bitstr,full.strRgr,full.nbits); 
    550551//full.nbits  = std::numeric_lim its<double>::digits-1;  // number of bits available in double 
    551         /*bvec in(n_data-1); 
    552         in.clear(); 
    553         full.bitstr = str_bitset(in,full.strRgr,full.nbits);*/ 
    554  
    555  
    556  
    557  
    558         full.loglik = seloglik1 ( full );    // % loglikelihood 
    559  
     552/*bvec in(n_data-1); 
     553in.clear(); 
     554full.bitstr = str_bitset(in,full.strRgr,full.nbits);*/ 
     555 
     556 
     557 
     558 
     559full.loglik = seloglik1(full);       // % loglikelihood 
     560    
    560561 
    561562 
     
    563564 
    564565//% construct full and empty structure 
    565         full = sestrremove ( full, belief_out ); 
    566         str_aux empty = sestrremove ( full, setdiff ( full.strRgr, belief_in ) ); 
    567  
     566full = sestrremove(full,belief_out); 
     567str_aux empty = sestrremove(full,setdiff(full.strRgr,belief_in)); 
     568  
    568569//% stopping rule calculation: 
    569570 
     
    571572 
    572573 
    573         bmat local_max ( 0, 0 ); 
    574         int to, muto = 0; 
    575  
     574bmat local_max(0,0); 
     575int to, muto = 0; 
     576  
    576577//% statistics: 
    577 //double cputime0 = cputime; 
     578//double cputime0 = cputime; ?????????????????????????????????????????????????????????????? 
    578579//if nargout>=3; 
    579  
    580         CPU_Timer timer; 
    581         timer.start(); 
    582  
    583         ivec mutos ( max_nrep + 2 ); 
    584         vec maxmutos ( max_nrep + 2 ); 
    585         mutos.zeros(); 
    586         maxmutos.zeros(); 
     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(); 
    587588 
    588589 
    589590//end; 
    590591//% ---------------------- 
    591  
     592  
    592593//% For stopping-rule calculation 
    593594//%so       = 2^(n_data -1-length(belief_in)- length(belief_out)); % do we use this ? 
    594595//% ---------------------- 
    595  
    596         ivec all_str = linspace ( 1, n_data ); 
    597  
    598         Array<str_aux> global_best ( 1 ); 
    599         global_best ( 0 ) = full; 
     596  
     597ivec all_str = linspace(1,n_data); 
     598  
     599Array<str_aux> global_best(1); 
     600global_best(0) = full; 
    600601 
    601602 
    602603//% MAIN LOOP is here. 
    603604 
    604         str_aux   best; 
    605         for ( int n_start = -1; n_start <= max_nrep; n_start++ ) { 
    606                 str_aux  last, best; 
    607  
    608                 to = n_start + 2; 
    609  
    610                 if ( n_start == -1 ) { 
    611                         //% start from the full structure 
    612                         last = full; 
    613                 } else { 
    614                         if ( n_start == 0 ) 
    615                                 //% start from the empty structure 
    616                                 last = empty; 
    617  
    618                         else { 
    619                                 //% start from random structure 
    620  
    621                                 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 
    622                                 last = sestrremove ( full, setdiff ( all_str, ::concat<int> ( ::concat<int> ( 1 , last_str ), empty.strRgr ) ) ); 
    623  
    624                         } 
     605str_aux   best; 
     606for (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           
     640ivec 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); 
    625648                } 
    626                 //% DEBUGging print: 
    627                 //%fprintf('STRUCTURE generated            in loop %2i was %s\n', n_start, strPrintstr(last)); 
    628  
    629                 //% The loop is repeated until likelihood stops growing (break condition 
    630                 //% used at the end; 
    631  
    632  
    633                 while ( 1 ) { 
    634                         //% This structure is going to hold the best elements 
    635                         best = last; 
    636                         //% Nesting by removing elements (enpoorment) 
    637                         ivec  removed_items = setdiff ( last.strRgr, belief_in ); 
    638  
    639                         ivec removed_item; 
    640                         str_aux newone; 
    641  
    642                         for ( int i = 0; i < removed_items.length(); i++ ) { 
    643                                 removed_item = vec_1 ( removed_items ( i ) ); 
    644                                 newone = sestrremove ( last, removed_item ); 
    645                                 if ( nbest > 1 ) { 
    646                                         add_new ( global_best, newone, nbest ); 
    647                                 } 
    648                                 if ( newone.loglik > best.loglik ) { 
    649                                         best = newone; 
    650                                 } 
    651                         } 
    652                         //% Nesting by adding elements (enrichment) 
    653                         ivec added_items = setdiff ( last.strMis, belief_out ); 
    654                         ivec added_item; 
    655  
    656                         for ( int j = 0; j < added_items.length(); j++ ) { 
    657                                 added_item = vec_1 ( added_items ( j ) ); 
    658                                 newone = sestrinsert ( last, added_item ); 
    659                                 if ( nbest > 1 ) { 
    660                                         add_new ( global_best, newone, nbest ); 
    661                                 } 
    662                                 if ( newone.loglik > best.loglik ) { 
    663                                         best = newone; 
    664                                 } 
    665                         } 
    666  
    667  
    668  
    669  
    670  
    671                         //% Break condition if likelihood does not change. 
    672                         if ( best.loglik <= last.loglik ) 
    673                                 break; 
    674                         else 
    675                                 //% Making best structure last structure. 
    676                                 last = best; 
    677  
    678  
     649                if (newone.loglik>best.loglik){ 
     650           best = newone; 
    679651                } 
    680  
    681  
    682  
    683  
    684  
    685                 // % DEBUGging print: 
    686                 //%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); 
    687  
    688                 //% Collecting of the best structure in case we don't need the second parameter 
    689                 if ( nbest <= 1 ) { 
    690                         if ( best.loglik > global_best ( 0 ).loglik ) { 
    691                                 global_best = best; 
    692                         } 
     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); 
    693662                } 
    694  
    695                 //% uniqueness of the structure found 
    696                 int append = 1; 
    697  
    698  
    699                 for ( int j = 0; j  < local_max.rows() ; j++ ) { 
    700                         if ( best.bitstr == local_max.get_row ( j ) ) { 
    701                                 append = 0; 
    702                                 break; 
    703                         } 
     663                if (newone.loglik>best.loglik){ 
     664           best = newone; 
    704665                } 
    705  
    706  
    707                 if ( append ) { 
    708                         local_max.append_row ( best.bitstr ); 
    709                         muto = muto + 1; 
    710                 } 
    711  
    712                 //% stopping rule: 
    713                 double maxmuto = ( to - order_k - 1 ) / lambda - to + 1; 
    714                 if ( to > 2 ) { 
    715                         if ( maxmuto >= muto ) { 
    716                                 //% fprintf('*'); 
    717                                 break; 
    718                         } 
    719                 } 
    720  
    721                 // do statistics if necessary: 
    722                 //if (nargout>=3){ 
    723                 mutos ( to - 1 )    = muto; 
    724                 maxmutos ( to - 1 ) = maxmuto; 
    725                 //} 
    726         } 
    727  
     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  
    728729//% Aftermath: The best structure was in: global_best 
    729  
     730  
    730731//% Updating loglikelihoods: we have to add the constant stuff 
    731732 
    732733 
    733734 
    734         for ( int f = 0 ; f < global_best.length(); f++ ) { 
    735                 global_best ( f ).loglik = global_best ( f ).loglik + seloglik2 ( global_best ( f ) ); 
    736         } 
    737  
    738         /*for f=1:length(global_best); 
    739            global_best(f).loglik = global_best(f).loglik + seloglik2(global_best(f)); 
    740         end;*/ 
     735for (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)); 
     741end;*/ 
    741742 
    742743 
    743744//% Making first output parameter: 
    744745 
    745         int max_i = 0; 
    746         for ( int j = 1; j < global_best.length(); j++ ) 
    747                 if ( global_best ( max_i ).loglik < ( global_best ( j ).loglik ) ) max_i = j; 
    748  
    749         best = global_best ( max_i ); 
    750  
     746int max_i = 0; 
     747for (int j = 1; j < global_best.length(); j++) 
     748  if  (global_best(max_i).loglik < (global_best(j).loglik)) max_i = j; 
     749 
     750best = global_best(max_i); 
     751cout << endl << endl << endl << endl << best.strRgr; 
     752cout << endl << endl << endl << endl << best.strRgr;  
    751753//% Making the second output parameter 
    752754 
    753         vec logliks ( global_best.length() ); 
    754         for ( int j = 0; j < logliks.length(); j++ ) 
    755                 logliks ( j ) = global_best ( j ).loglik; 
    756  
    757         ivec i = sort_index ( logliks ); 
    758         rgrsout.set_length ( global_best.length() ); 
    759  
    760         for ( int j = global_best.length() - 1; j >= 0; j-- ) 
    761                 rgrsout ( j ) = global_best ( i ( j ) ); 
    762  
     755vec logliks(global_best.length()); 
     756for (int j = 0; j < logliks.length(); j++) 
     757 logliks(j) = global_best(j).loglik; 
     758 
     759ivec i = sort_index(logliks); 
     760rgrsout.set_length(global_best.length()); 
     761 
     762for (int j = global_best.length() - 1; j >= 0; j--) 
     763rgrsout(j) = global_best(i(j)); 
     764  
    763765//if (nargout>=3); 
    764  
    765  
    766         str_statistics statistics; 
    767  
    768         statistics.allstrs = 2 ^ ( n_data - 1 - length ( belief_in ) - length ( belief_out ) ); 
    769         statistics.nrand   = to - 2; 
    770         statistics.unique  = muto; 
    771         statistics.to  = to; 
    772         statistics.cputime_seconds = timer.get_time(); 
    773         statistics.itemspeed       = statistics.to / statistics.cputime_seconds; 
    774         statistics.muto = muto; 
    775         statistics.mutos = mutos; 
    776         statistics.maxmutos = maxmutos; 
     766    
     767    
     768str_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; 
    777779//end; 
    778780 
    779         return best.strRgr; 
    780  
    781 } 
     781return best.strRgr; 
     782 
     783}  
    782784 
    783785#ifdef LADIM 
    784786//% randun seed stuff: 
    785787//%randn('seed',SEED); 
    786  
     788  
    787789//% --------------------- END of MAIN program -------------------- 
    788  
     790  
    789791% This is needed for bitstr manipulations 
    790792/*function out = str_bitset(in,ns,nbits) 
     
    794796      bitindex = 1+rem(n-2,nbits); 
    795797      out(index) = bitset(out(index),bitindex); 
    796    end; 
     798   end;   
    797799function out = str_bitres(in,ns,nbits) 
    798800   out = in; 
     
    803805      out(index) = bitxor(bitor(out(index),mask),mask); 
    804806   end;*/ 
    805  
    806 function out = strPrintstr ( in ) 
    807                        out = '0'; 
    808 nbits = in.nbits; 
    809 for f = 2: 
    810         length ( in.d0 ); 
    811 index = 1 + floor ( ( f - 2 ) / nbits ); 
    812 bitindex = 1 + rem ( f - 2, nbits ); 
    813 if bitget ( in.bitstr ( index ), bitindex ); 
    814 out ( f ) = '1'; 
    815 else; 
    816 out ( f ) = '0'; 
    817 end; 
    818 end; 
    819  
     807  
     808function 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  
    820821/*function global_best_out = add_new(global_best,new,nbest) 
    821822% Eventually add to global best, but do not go over nbest values 
     
    833834               break; 
    834835            end; 
    835          end; 
     836         end;          
    836837         if addit; 
    837838            global_best_out(i) = new; 
    838839            % DEBUGging print: 
    839840            % fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik); 
    840          end; 
     841         end;          
    841842      end; 
    842843   else; 
    843844      global_best_out = [global_best new]; 
    844845   end;*/ 
    845  
     846  
    846847/*function out = sestrremove(in,removed_elements); 
    847848% Removes elements from regressor 
     
    864865   out.bitstr = str_bitres(out.bitstr,removed_elements,out.nbits); 
    865866   out.loglik = seloglik1(out);*/ 
    866  
     867    
    867868/*function out = sestrinsert(in,inserted_elements); 
    868869% Moves elements into regressor 
     
    879880          % END 
    880881      end; 
    881    end; 
     882   end;          
    882883   out.posit1 = find(out.strL==1); 
    883884   out.strRgr = out.strL((out.posit1+1):n_strL); 
     
    885886   out.bitstr = str_bitset(out.bitstr,inserted_elements,out.nbits); 
    886887   out.loglik = seloglik1(out);*/ 
    887  
     888  
    888889% 
    889890% seloglik_real = seloglik1 + seloglik2 
    890                   % 
    891  
    892                   /*function l = seloglik1(in) 
    893                   % This is the loglikelihood (non-constant part) - this should be used in 
    894                   % frequent computation 
    895                      len = length(in.d); 
    896                      p1  = in.posit1; 
    897  
    898                      i1 = -0.5*in.nu *log(in.d (p1)) -0.5*sum(log(in.d ((p1+1):len))); 
    899                      i0 = -0.5*in.nu0*log(in.d0(p1)) -0.5*sum(log(in.d0((p1+1):len))); 
    900                      l  = i1-i0; 
    901  
    902                      % DEBUGGing print: 
    903                      % fprintf('SELOGLIK1: str=%s loglik=%g\n', strPrintstr(in), l);*/ 
    904  
    905  
    906                   function l = seloglik2 ( in ) 
    907                                % This is the loglikelihood ( constant part ) - this should be added to 
    908                                % everything at the end. It needs some computation, so it is useless to 
    909                                % make it for all the stuff 
    910                                logpi = log ( pi ); 
    911  
    912 i1 = gammaln ( in.nu / 2 ) - 0.5*in.nu *logpi; 
    913 i0 = gammaln ( in.nu0 / 2 ) - 0.5*in.nu0*logpi; 
    914 l  = i1 - i0; 
    915  
    916  
     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  
     907function 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  
    917918/*function [Lout, dout] = seswapudl(L,d,i); 
    918919%SESWAPUDL swaps information matrix in decomposition V=L^T diag(d) L 
     
    922923% L : lower triangular matrix with 1's on diagonal of the decomposistion 
    923924% d : diagonal vector of diagonal matrix of the decomposition 
    924 % i : index of line to be swapped with the next one 
     925% i : index of line to be swapped with the next one  
    925926% Lout : output lower triangular matrix 
    926927% dout : output diagional vector of diagonal matrix D 
     
    928929% Description: 
    929930%  Lout' * diag(dout) * Lout = P(i,i+1) * L' * diag(d) * L * P(i,i+1); 
    930 % 
     931%  
    931932%  Where permutation matrix P(i,j) permutates columns if applied from the 
    932933%  right and line if applied from the left. 
    933 % 
     934%    
    934935% Note: naming: 
    935936%       se = structure estimation 
    936937%       lite = light, simple 
    937938%       udl = U*D*L, or more precisely, L'*D*L, also called as ld 
    938 % 
     939%    
    939940% Design  : L. Tesar 
    940941% Updated : Feb 2003 
    941942% Project : post-ProDaCTool 
    942943% Reference: sedydr 
    943  
     944    
    944945j = i+1; 
    945  
     946  
    946947pomd = d(i); 
    947948d(i) = d(j); 
    948949d(j) = pomd; 
    949  
     950  
    950951pomL   = L(i,:); 
    951952L(i,:) = L(j,:); 
    952953L(j,:) = pomL; 
    953  
     954  
    954955pomL   = L(:,i); 
    955956L(:,i) = L(:,j); 
    956957L(:,j) = pomL; 
    957  
     958  
    958959% We must be working with LINES of matrix L ! 
    959  
    960 r  = L(i,:)'; 
     960  
     961r  = L(i,:)';  
    961962f  = L(j,:)'; 
    962963Dr = d(i); 
    963964Df = d(j); 
    964  
     965  
    965966[r, f, Dr, Df] = sedydr(r, f, Dr, Df, j); 
    966  
     967  
    967968r0 = r(i); 
    968969Dr = Dr*r0*r0; 
    969970r  = r/r0; 
    970  
    971 L(i,:) = r'; 
     971  
     972L(i,:) = r';  
    972973L(j,:) = f'; 
    973974d(i)   = Dr; 
    974975d(j)   = Df; 
    975  
     976  
    976977L(i,i) = 1; 
    977978L(j,j) = 1; 
    978  
     979  
    979980Lout = L; 
    980981dout = d;*/ 
    981  
     982  
    982983/*function [rout, fout, Drout, Dfout, kr] = sedydr(r,f,Dr,Df,R,jl,jh); 
    983984%SEDYDR dyadic reduction, performs transformation of sum of 2 dyads 
     
    986987% [rout, fout, Drout, Dfout] = sedydr(r,f,Dr,Df,R); 
    987988% 
    988 % Description: dyadic reduction, performs transformation of sum of 
     989% Description: dyadic reduction, performs transformation of sum of  
    989990%  2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed by R is zeroed 
    990991% 
     
    994995%     Df   : scalar with weight of reducing dyad 
    995996%     R    : scalar number giving 1 based index to the element of r, 
    996 %            which is to be reduced to 
     997%            which is to be reduced to  
    997998%            zero; the corresponding element of f is assumed to be 1. 
    998 %     jl   : lower index of the range within which the dyads are 
     999%     jl   : lower index of the range within which the dyads are  
    9991000%            modified (can be omitted, then everything is updated) 
    1000 %     jh   : upper index of the range within which the dyads are 
     1001%     jh   : upper index of the range within which the dyads are  
    10011002%            modified (can be omitted then everything is updated) 
    10021003%     rout,fout,Drout,dfout : resulting two dyads 
     
    10041005%            rnew = r + kr*f 
    10051006% 
    1006 % Description: dyadic reduction, performs transformation of sum of 
     1007% Description: dyadic reduction, performs transformation of sum of  
    10071008%            2 dyads r*Dr*r'+ f*Df*f' so that the element of r indexed by R is zeroed 
    10081009% Remark1: Constant mzero means machine zero and should be modified 
     
    10241025% Project: post-ProDaCTool 
    10251026% Reference: none 
    1026  
     1027    
    10271028if nargin<6; 
    10281029   update_whole=1; 
     
    10301031   update_whole=0; 
    10311032end; 
    1032  
     1033  
    10331034mzero = 1e-32; 
    1034  
     1035  
    10351036if Dr<mzero; 
    10361037   Dr=0; 
    10371038end; 
    1038  
     1039  
    10391040r0   = r(R); 
    10401041kD   = Df; 
    10411042kr   = r0 * Dr; 
    10421043Dfout   = kD + r0 * kr; 
    1043  
     1044  
    10441045if Dfout > mzero; 
    10451046    kD = kD / Dfout; 
     
    10491050    kr = 0; 
    10501051end; 
    1051  
     1052  
    10521053Drout = Dr * kD; 
    1053  
     1054  
    10541055% Try to uncomment marked stuff (*) if in numerical problems, but I don't 
    10551056% think it can make any difference for normal healthy floating-point unit 
     
    10591060   fout = f + kr*rout; 
    10601061%   fout(R) = 1;   % * could be needed for some nonsense cases(or numeric reasons?), normally not 
    1061 else; 
     1062else;   
    10621063   rout = r; 
    10631064   fout = f; 
    10641065   rout(jl:jh) = r(jl:jh) - r0 * f(jl:jh); 
    1065    rout(R) = 0; 
     1066   rout(R) = 0;    
    10661067   fout(jl:jh) = f(jl:jh) + kr * rout(jl:jh); 
    10671068end;*/ 
    1068  
    1069  
    1070  
     1069  
     1070  
     1071  
    10711072#endif 
    1072  
    1073  
    1074 } 
    1075  
     1073  
     1074  
     1075} 
     1076