#include <limits>
#include "arx.h"


#include <fstream>
#include<iostream>
#include<iomanip>

namespace bdm {




/*
struct str_aux {
	vec d0;
	double nu0;
	mat L0;
	mat L;
	vec d;
	double nu;
	ivec strL;                 // Current structure of L and d
	ivec strRgr;               // Structure elements currently inside regressor (after regressand)
	ivec strMis;               // structure elements, that are currently outside regressor (before regressand)
	int posit1;                // regressand position
	int nbits; 				   // number of bits available in double
	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 ) {
// 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;

// 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;
//   rout(R) = 0;   // * could be needed for some nonsense cases(or numeric reasons?), normally not
	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);*/

//% We must be working with LINES of matrix L !



	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 ) {
//% 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;
}



/*

  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);
	  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);
      global_best_out = global_best;
	  if (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_out(i) = newone;
            // DEBUGging print:
            // fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik);
	     }
	 }
}
   else
       global_best_out = concat(global_best, newone);

  return global_best_out;

  }

*/

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 ) {
// 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 ) {
// 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*/ ) {
// 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 );
//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





//% construct full and empty structure
	full = sestrremove ( full, belief_out );
	str_aux empty = sestrremove ( full, setdiff ( full.strRgr, belief_in ) );

//% stopping rule calculation:




	bmat local_max ( 0, 0 );
	int to, muto = 0;

//% statistics:
//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();


//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;


//% 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;
		//}
	}

//% 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;*/


//% 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 );

//% 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 ) );

//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;
//end;

	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)
   out = in;
   for n = ns;
      index = 1+floor((n-2)/nbits);
      bitindex = 1+rem(n-2,nbits);
      out(index) = bitset(out(index),bitindex);
   end;
function out = str_bitres(in,ns,nbits)
   out = in;
   for n = ns;
      index = 1+floor((n-2)/nbits);
      bitindex = 1+rem(n-2,nbits);
      mask = bitset(0,bitindex);
      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 global_best_out = add_new(global_best,new,nbest)
% Eventually add to global best, but do not go over nbest values
% Also avoids repeating things, which makes this function awfully slow
   if length(global_best)>=nbest;
      logliks = [global_best.loglik];
      [loglik i] = min(logliks);
      global_best_out = global_best;
      if loglik<new.loglik;
         %         if ~any(logliks == new.loglik);
         addit=1;
         for f = [global_best.bitstr];
            if f == new.bitstr;
               addit = 0;
               break;
            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;
   else;
      global_best_out = [global_best new];
   end;*/

/*function out = sestrremove(in,removed_elements);
% Removes elements from regressor
   n_strL = length(in.strL);
   out = in;
   for f=removed_elements;
      posit1 = find(out.strL==1);
      positf = find(out.strL==f);
      for g=(positf-1):-1:posit1;
         % BEGIN: We are swapping g and g+1 NOW!!!!
         [out.L, out.d]   = seswapudl(out.L, out.d, g);
         [out.L0, out.d0]   = seswapudl(out.L0, out.d0, g);
         out.strL([g g+1]) = [out.strL(g+1) out.strL(g)];
         % END
      end;
   end;
   out.posit1 = find(out.strL==1);
   out.strRgr = out.strL((out.posit1+1):n_strL);
   out.strMis = out.strL(1:(out.posit1-1));
   out.bitstr = str_bitres(out.bitstr,removed_elements,out.nbits);
   out.loglik = seloglik1(out);*/

/*function out = sestrinsert(in,inserted_elements);
% Moves elements into regressor
   n_strL = length(in.strL);
   out = in;
   for f=inserted_elements;
      posit1 = find(out.strL==1);
      positf = find(out.strL==f);
      for g=positf:(posit1-1);
          % BEGIN: We are swapping g and g+1 NOW!!!!
          [out.L,  out.d]   = seswapudl(out.L,  out.d,  g);
          [out.L0, out.d0]  = seswapudl(out.L0, out.d0, g);
          out.strL([g g+1]) = [out.strL(g+1) out.strL(g)];
          % END
      end;
   end;
   out.posit1 = find(out.strL==1);
   out.strRgr = out.strL((out.posit1+1):n_strL);
   out.strMis = out.strL(1:(out.posit1-1));
   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 [Lout, dout] = seswapudl(L,d,i);
%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

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,:)';
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(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
%
% [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;

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;
    kr = kr / Dfout;
else;
    kD = 1;
    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
if update_whole;
   rout = r - r0*f;
%   rout(R) = 0;   % * could be needed for some nonsense cases(or numeric reasons?), normally not
   fout = f + kr*rout;
%   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;*/



#endif


}