Index: /library/bdm/estim/arx_straux.cpp
===================================================================
--- /library/bdm/estim/arx_straux.cpp (revision 646)
+++ /library/bdm/estim/arx_straux.cpp (revision 648)
@@ -11,8 +11,8 @@
 namespace bdm {
 
-  
-
-	
-/*	
+
+
+
+/*
 struct str_aux {
 	vec d0;
@@ -27,304 +27,304 @@
 	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,7 +333,7 @@
          //         if ~any(logliks == new.loglik);
           addit=1;
-		 
-		 
-	 
+
+
+
 		  if (newone.bitstr.length() == 1) {
 	  for (int j = 0; j < global_best.length(); j++){
@@ -344,6 +344,6 @@
                break;
 			  }
-			  }   
-	  }   
+			  }
+	  }
 		  }
 		 if (addit){
@@ -351,6 +351,6 @@
             // DEBUGging print:
             // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik);
-	     }        
-	 }  
+	     }
+	 }
 }
    else
@@ -358,209 +358,204 @@
 
   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,7 +563,7 @@
 
 //% 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,220 +571,220 @@
 
 
-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,5 +794,5 @@
       bitindex = 1+rem(n-2,nbits);
       out(index) = bitset(out(index),bitindex);
-   end;  
+   end;
 function out = str_bitres(in,ns,nbits)
    out = in;
@@ -808,18 +803,19 @@
       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,15 +833,15 @@
                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,5 +864,5 @@
    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,5 +879,5 @@
           % END
       end;
-   end;         
+   end;
    out.posit1 = find(out.strL==1);
    out.strRgr = out.strL((out.posit1+1):n_strL);
@@ -889,34 +885,34 @@
    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,5 +922,5 @@
 % 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,56 +928,56 @@
 % 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,5 +986,5 @@
 % [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,9 +994,9 @@
 %     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,5 +1004,5 @@
 %            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,5 +1024,5 @@
 % Project: post-ProDaCTool
 % Reference: none
-   
+
 if nargin<6;
    update_whole=1;
@@ -1034,16 +1030,16 @@
    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,7 +1049,7 @@
     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,17 +1059,17 @@
    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
- 
- 
-}
-
+
+
+}
+
