Changeset 684

Timestamp:

10/27/09 21:15:03 (15 years ago)

Author:

sarka

Message:

function straux1 converted from straux1.m

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
+   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);
                 }
         }
-        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;
-}
-
-
-
+  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;
-
-
-                        //????????????????????????????????????????????
-                        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
+   
 
 
@@ -563 +564 @@
 
 //% 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:
 
@@ -571 +572 @@
 
 
-        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 ) ) );
-
-                        }
+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);
                 }
-                //% 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;
-
-
+                if (newone.loglik>best.loglik){
+           best = newone;
                 }
-
-
-
-
-
-                // % 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;
-                        }
+          }
+      //% 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);
                 }
-
-                //% 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 (newone.loglik>best.loglik){
+           best = newone;
                 }
-
-
-                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;
-                //}
-        }
-
+     }
+
+   
+
+   
+
+      //% 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);
+cout << endl << endl << endl << endl << best.strRgr;
+cout << endl << endl << endl << endl << best.strRgr; 
 //% 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)
@@ -794 +796 @@
       bitindex = 1+rem(n-2,nbits);
       out(index) = bitset(out(index),bitindex);
-   end;
+   end;  
 function out = str_bitres(in,ns,nbits)
    out = in;
@@ -803 +805 @@
       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
@@ -833 +834 @@
                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
@@ -864 +865 @@
    out.bitstr = str_bitres(out.bitstr,removed_elements,out.nbits);
    out.loglik = seloglik1(out);*/
-
+   
 /*function out = sestrinsert(in,inserted_elements);
 % Moves elements into regressor
@@ -879 +880 @@
           % END
       end;
-   end;
+   end;         
    out.posit1 = find(out.strL==1);
    out.strRgr = out.strL((out.posit1+1):n_strL);
@@ -885 +886 @@
    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
@@ -922 +923 @@
 % 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
@@ -928 +929 @@
 % 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
@@ -986 +987 @@
 % [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
 %
@@ -994 +995 @@
 %     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
@@ -1004 +1005 @@
 %            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
@@ -1024 +1025 @@
 % Project: post-ProDaCTool
 % Reference: none
-
+   
 if nargin<6;
    update_whole=1;
@@ -1030 +1031 @@
    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;
@@ -1049 +1050 @@
     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
@@ -1059 +1060 @@
    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
-
-
-}
-
+ 
+ 
+}
+