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