Changeset 648

Timestamp:

10/02/09 12:05:56 (15 years ago)

Author:

mido

Message:

arx_strax astylled

Files:

• library/bdm/estim/arx_straux.cpp (modified) (23 diffs)

library/bdm/estim/arx_straux.cpp

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