%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute log-likelihood in time domain assuming a multivariate gaussian
% distribution
% 
% INPUT
% 
% - res, a vector of AOs containing residuals
% - noise, a vector of AOs containing noise
% - cutbefore, followed by the data samples to cut at the starting of the
% data series
% - cutafter, followed by the data samples to cut at the ending of the
% data series
%
% L Ferraioli 10-10-2010
%
% $Id: loglikelihood_td.m,v 1.1 2011/03/15 16:19:20 miquel Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function loglk = loglikelihood_td(res,noise,varargin)

  cutbefore = [];
  cutafter = [];
  if ~isempty(varargin)
    for j=1:length(varargin)
      if strcmp(varargin{j},'cutbefore')
        cutbefore = varargin{j+1};
      end
      if strcmp(varargin{j},'cutafter')
        cutafter = varargin{j+1};
      end
    end
  end
  
  nres = numel(res);
  loglk = 0;
  
  for ii=1:nres
    
    yres = res(ii).y;
    ynoise = noise(ii).y;
    % willing to work with rows
    if size(yres,2)<size(yres,1)
      yres = yres.';
    end
    if size(ynoise,2)<size(ynoise,1)
      ynoise = ynoise.';
    end
    if ~isempty(cutbefore)
      yres(1:cutbefore) = [];
      ynoise(1:cutbefore) = [];
    end
    if ~isempty(cutafter)
      yres(end-cutafter:end) = [];
      ynoise(end-cutafter:end) = [];
    end
    
%     R = utils.math.Rcovmat(yres);
%     
%     xx = R\yres.';
%     Ndim = numel(yres);
%     
%     loglk = loglk + abs(xx'*xx)./Ndim;

    cmat = utils.math.xCovmat(ynoise);
    
    [L,p] = chol(cmat,'lower');
    if p==0
      xx = L\yres.';
      Ndim = numel(yres);
    else
      q = p-1;
      yyres = yres(1:q);

      xx = L\yyres.';
      Ndim = numel(yyres);
    end
    
    loglk = loglk + abs(xx'*xx)./Ndim;
    
    
  end
  
  
  %%% get cros terms
  if nres>1
  
    for ii=1:nres-1
      for jj=ii+1:nres
        
        yres1 = res(ii).y;
        yres2 = res(jj).y;
        ynoise1 = noise(ii).y;
        ynoise2 = noise(jj).y;
        % willing to work with rows
        if size(yres1,2)<size(yres1,1)
          yres1 = yres1.';
        end
        if size(yres2,2)<size(yres2,1)
          yres2 = yres2.';
        end
        if size(ynoise1,2)<size(ynoise1,1)
          ynoise1 = ynoise1.';
        end
        if size(ynoise2,2)<size(ynoise2,1)
          ynoise2 = ynoise2.';
        end
        if ~isempty(cutbefore)
          yres1(1:cutbefore) = [];
          yres2(1:cutbefore) = [];
          ynoise1(1:cutbefore) = [];
          ynoise2(1:cutbefore) = [];
        end
        if ~isempty(cutafter)
          yres1(end-cutafter:end) = [];
          yres2(end-cutafter:end) = [];
          ynoise1(end-cutafter:end) = [];
          ynoise2(end-cutafter:end) = [];
        end
        
%         Rx = utils.math.Rcovmat(yres1);
%         Ry = utils.math.Rcovmat(yres2);
%         
%         xx = Rx\yres1.';
%         yy = Ry\yres2.';
%         
%         Ndim = numel(yres1);
% 
%         loglk = loglk + 2.*abs(xx'*yy)./Ndim;
        
        cmat = utils.math.xCovmat(ynoise1,ynoise2);
    
        [L,p] = chol(cmat,'lower');
        if p==0
          xx = L\yres1.';
          Ndim = numel(yres1);
          yy = L\yres2.';
        else
          q = p-1;
          yyres1 = yres1(1:q);
          yyres2 = yres2(1:q);

          xx = L\yyres1.';
          Ndim = numel(yyres1);
          yy = L\yyres2.';
        end

        loglk = loglk + 2.*abs(xx'*yy)./Ndim;
        
        
      end
    end
    
  end


end