function [logLL, V_hat, transformed_error, pV_hat, transformed_perror]=my_spatial_LL_split_Periphery_v2(EE, pEE, IV, T, N, W, RHO, IsSparse, IndexM)
% computes log likelihood and variance-covariance matrix for given spatial lag rho

p=size(EE,1);
pp=size(pEE,1);
pz=size(IV,1);
if ~exist('IsSparse','var')
    IsSparse=(numel(find(W==0))/numel(W)>0.1 );
elseif isempty(IsSparse)
    IsSparse=(numel(find(W==0))/numel(W)>0.1 );    
end;
NoGroups=size(IndexM,2);

V_hat=cell(NoGroups,1);

GroupedRHO=reshape(RHO(:),p+pz+pp,NoGroups);


try
    if IsSparse %check whether the weighting matrix is sparse
            transformed_error=zeros(p+pz,T,N);
            transformed_perror=zeros(pp,T,N); % keep separate to allow for different samples
            original_error=[EE; IV];
            original_perror=pEE;
            LogLL_det=0;
            p_idx=0;
            while p_idx< (p+pz) && ~isnan(LogLL_det)
                p_idx=p_idx+1;
                if abs(eigs(sparse(W')*spdiags(IndexM*GroupedRHO(p_idx,:)',0,N,N),1)) >= 1-1e-4
                    fprintf(2,'Spatial filter not invertible \n');
                    LogLL_det=NaN;
                    V_hat=[];
                    pV_hat=[];
                else
                    transformed_error(p_idx,:,:)=squeeze(original_error(p_idx,:,:))*(speye(N)-sparse(W')*spdiags(IndexM*GroupedRHO(p_idx,:)',0,N,N)) ;
                    LogLL_det = LogLL_det +  2 * T * sum( log(diag(chol(speye(N)-sparse(W')*spdiags(IndexM*GroupedRHO(p_idx,:)',0,N,N), 'lower'))) );
    %                 [~, U, P] = lu( speye(N)-sparse(W')*RHO(p_idx) );
    %                 du = diag(U);
    % %                     c = det(P) * prod(sign(du));
    % %                     v = log(c) + sum(log(abs(du)));
    %                 LogLL = LogLL +  T * ( log(P) + log( prod(sign(du)) ) + sum(log(abs(du))) )/2;                
                end;
            end;
            p_idx=0;
            while p_idx< pp && ~isnan(LogLL_det)
                p_idx=p_idx+1;
                if abs(eigs(sparse(W')*spdiags(IndexM*GroupedRHO(p+pz+p_idx,:)',0,N,N),1)) >= 1-1e-4
                    fprintf(2,'Spatial filter not invertible \n');
                    LogLL_det=NaN;
                    V_hat=[];
                    pV_hat=[];
                else
                    transformed_perror(p_idx,:,:)=squeeze(original_perror(p_idx,:,:))*(speye(N)-sparse(W')*spdiags(IndexM*GroupedRHO(p+pz+p_idx,:)',0,N,N)) ;
                    LogLL_det = LogLL_det +  2 * T * sum( log(diag(chol(speye(N)-sparse(W')*spdiags(IndexM*GroupedRHO(p+pz+p_idx,:)',0,N,N), 'lower'))) );
    %                 [~, U, P] = lu( speye(N)-sparse(W')*RHO(p_idx) );
    %                 du = diag(U);
    % %                     c = det(P) * prod(sign(du));
    % %                     v = log(c) + sum(log(abs(du)));
    %                 LogLL = LogLL +  T * ( log(P) + log( prod(sign(du)) ) + sum(log(abs(du))) )/2;                
                end;
            end;
            
            if ~isnan(LogLL_det)
                logLL=LogLL_det;
                for gg=1:NoGroups
                    V_hat{gg}=( reshape(transformed_error(:,:,1==IndexM(:,gg)), p+pz,T*sum(IndexM(:,gg)))*...
                        reshape(transformed_error(:,:,1==IndexM(:,gg)), p+pz,T*sum(IndexM(:,gg)))' ) / (sum(IndexM(:,gg))*T);
                    pV_hat{gg}=( reshape(transformed_perror(:,:,1==IndexM(:,gg)), pp,T*sum(IndexM(:,gg)))*...
                            reshape(transformed_perror(:,:,1==IndexM(:,gg)), pp,T*sum(IndexM(:,gg)))' ) / (sum(IndexM(:,gg))*T);

                    logLL = logLL - sum( log(diag(chol(V_hat{gg}, 'lower'))) )* sum(IndexM(:,gg))*T ...
                                    - sum( log(diag(chol(pV_hat{gg}, 'lower'))) )* sum(IndexM(:,gg))*T; % other terms are constant at V_hat
                end;
            else
                logLL=NaN;
            end;
    else
            transformed_error=zeros(p+pz,T,N);
            transformed_perror=zeros(pp,T,N); % keep separate to allow for different samples
            original_error=[EE; IV];
            original_perror=pEE;

            LogLL_det=0;
            p_idx=0;
            while p_idx< (p+pz) && ~isnan(LogLL_det)
                p_idx=p_idx+1;
                if abs(eigs(W'*diag(IndexM*GroupedRHO(p_idx,:)'),1)) >= 1-1e-4
                    fprintf(2,'Spatial filter not invertible \n');
                    LogLL_det=NaN;
                    V_hat=[];
                    pV_hat=[];
                else
%                                                                                         TRUE_U*(eye(N)-D'*diag(IndexM*x(:)))
                    transformed_error(p_idx,:,:)=squeeze(original_error(p_idx,:,:))*(eye(N)-(W')*diag(IndexM*GroupedRHO(p_idx,:)')) ;
                    [~, U, P] = lu( eye(N)-(W')*diag(IndexM*GroupedRHO(p_idx,:)') );
                    du = diag(U);
                    c = det(P) * prod(sign(du));
                    v = log(c) + sum(log(abs(du)));
    %                 LogLL = LogLL +  2 * T * v / 2;
                    LogLL_det = LogLL_det +  T * v ;
                end;
            end;
            p_idx=0;
            while p_idx< pp && ~isnan(LogLL_det)
                p_idx=p_idx+1;
                if abs(eigs(W'*diag(IndexM*GroupedRHO(p+pz+p_idx,:)'),1)) >= 1-1e-4
                    fprintf(2,'Spatial filter not invertible \n');
                    LogLL_det=NaN;
                    V_hat=[];
                    pV_hat=[];
                else
                    transformed_perror(p_idx,:,:)=squeeze(original_perror(p_idx,:,:))*(eye(N)-(W')*diag(IndexM*GroupedRHO(p+pz+p_idx,:)')) ;
                    [~, U, P] = lu( eye(N)-(W')*diag(IndexM*GroupedRHO(p+pz+p_idx,:)') );
                    du = diag(U);
                    c = det(P) * prod(sign(du));
                    v = log(c) + sum(log(abs(du)));
    %                 LogLL = LogLL +  2 * T * v / 2;
                    LogLL_det = LogLL_det +  T * v ;
                end;
            end;
            
            if ~isnan(LogLL_det)
                logLL=LogLL_det;
                for gg=1:NoGroups
%                     V_hat{gg}=zeros(p+pz,p+pz);
%                     pV_hat{gg}=zeros(pp,pp);
%                     for t_idx=1:T
%                         V_hat{gg}=V_hat{gg}+( reshape(transformed_error(:,t_idx,1==IndexM(:,gg)), p+pz,sum(IndexM(:,gg)))*...
%                             reshape(transformed_error(:,t_idx,1==IndexM(:,gg)), p+pz,sum(IndexM(:,gg)))' ) / (sum(IndexM(:,gg))*T);
%                         pV_hat{gg}=pV_hat{gg}+( reshape(transformed_perror(:,t_idx,1==IndexM(:,gg)), pp,sum(IndexM(:,gg)))*...
%                             reshape(transformed_perror(:,t_idx,1==IndexM(:,gg)), pp,sum(IndexM(:,gg)))' ) / (sum(IndexM(:,gg))*T);
%                     end;

                    V_hat{gg}=( reshape(transformed_error(:,:,1==IndexM(:,gg)), p+pz,T*sum(IndexM(:,gg)))*...
                        reshape(transformed_error(:,:,1==IndexM(:,gg)), p+pz,T*sum(IndexM(:,gg)))' ) / (sum(IndexM(:,gg))*T);
                    pV_hat{gg}=( reshape(transformed_perror(:,:,1==IndexM(:,gg)), pp,T*sum(IndexM(:,gg)))*...
                            reshape(transformed_perror(:,:,1==IndexM(:,gg)), pp,T*sum(IndexM(:,gg)))' ) / (sum(IndexM(:,gg))*T);
%                     disp(log(diag(chol(V_hat{gg}, 'lower')))* sum(IndexM(:,gg))*T);
                    logLL = logLL  - sum( log(diag(chol(V_hat{gg}, 'lower'))) )* sum(IndexM(:,gg))*T ...
                                   - sum( log(diag(chol(pV_hat{gg}, 'lower'))) )* sum(IndexM(:,gg))*T; % other terms are constant at V_hat
                 end;
            else
                logLL=NaN;
            end;
    end;
catch
    fprintf(2,'Problem with eigenvalue computation \n');
    logLL=NaN;
    V_hat=[];    
end;

 