function [result_efir,numcoef] = lyngby_efir_main(pn, X, R, ...
    maxfilterlength, numreshuffle)

% lyngby_efir_main     - Main exhautive FIR function
% 
%       function [result_efir,numcoef] = lyngby_efir_main(PN, X, R,...
%                       maxfilterlength, numreshuffle)
%
%       Input:  PN                Input (the paradigm)
%               X                Target (the datamatrix)
%               R                 Run structure
%               maxfilterlength   Filter length (lag)
%               numreshuffle      Number of crossvalidation
%                                 permutations
%
%       Output: result_efir       Filter coefficients
%               numcoef           Optimal number of filter
%                                 coefficients
%
%       The FIR filter used in combination the SVD pseudo
%       inversion. The optimal model in each voxel is determined by
%       the Generalization error, based on the exhaustive search
%       through a 2D space. The 2D space is span by 'models' and
%       'splitratios'. The 'splitratio' is simply how the dataset
%       is split into a trainset and a testset. The 'models' is the
%       size of the SVD pseudo inversion length.
%
%       See also: lyngby_fir_main

% cvs : $Id: lyngby_efir_main.m,v 1.13 1999/07/09 13:43:53 fnielsen Exp $
%       $Revision: 1.13 $

    % lyngby_global;
        
    lyngby_log('Initializing');

    [rX, cX] = size(X);
    
    
    % Find number of scans
    [a,T] = lyngby_runinfo(R);
    
    if length(T) > 1
       error(['The Exhaustive FIR filter cannot handle runs with ',...
	   'different (unequal) lengths']);
    end    
    if a == 1
      error(sprintf(['The Exhaustive FIR filter cannot work with ',...
	    'just one run.']));
    end
    
    antalmodeller = maxfilterlength*(maxfilterlength+1)/2+1;
    
    mean_testfejl = zeros((a-1),antalmodeller);
    std_testfejl = zeros((a-1),antalmodeller);  
    result_efir = zeros(maxfilterlength,cX);
    opt_model = zeros(cX,2);


    % Generate the resampling matrix
    resample_index = [];
    for i = 1:numreshuffle
      resample_index = [resample_index,randperm(a)'];
    end


    % Generate Input
    p  = [pn;pn(1:T)];
    Xm = zeros((a+1)*T,maxfilterlength);
    for i = 1:maxfilterlength
      Xm(i:((a+1)*T),i) = p(1:((a+1)*T-i+1));
    end
    Xm = Xm(1+T:((a+1)*T),:);


    % Generate model index
    modelindex=zeros(antalmodeller,3);
    modelindex(antalmodeller,:)=[antalmodeller, 0 ,0];
    for fl=1:maxfilterlength
      for eigen=1:fl
	I=(fl-1)*fl/2+eigen;
	modelindex(I,:)=[I,fl,eigen];
      end
    end

    % Find the optimal model
    for v = 1:cX
      lyngby_log(sprintf('Computing Exhaustive FIR-map %i / %i', v, cX));
      y = X(:,v);
      Yresample = zeros(a*T, numreshuffle);
      
      for i=1:numreshuffle
	for j=1:a
	  Yresample(1+(j-1)*T:j*T,i)=...
	      y((1+(resample_index(j,i)-1)*T):(resample_index(j,i)*T));
	end
      end
      
      for fl=1:maxfilterlength
	for i3=1:(a-1)  % number of splitratios
	  Ytrain=Yresample(1:i3*T,:);
	  Ytest =Yresample(1+i3*T:a*T,:);      
	  XT=Xm(1:i3*T,1:fl);
	  XV=Xm(1:(a-i3)*T,1:fl);
	  XTT=XT'*XT;
	  [U,S,V]=svd(XTT);
	  epsilon=S*ones(fl,1);
	  testfejl0 = sum(Ytest.^2)./(T*(a-i3));
	  mean_testfejl(i3,antalmodeller)=mean(testfejl0);
	  std_testfejl(i3,antalmodeller) =std(testfejl0);  
	  for eigen=1:fl  
	    i2=(fl-1)*fl/2+eigen; 
	    SS=(1./(S+eps)).*(S>(epsilon(eigen)-(1e10*eps)));
	    preX=V*SS*U'*XT';
	    h=preX*Ytrain; 
	    Yfit=XV*h;
	    testfejl = sum((Ytest - Yfit).^2)./(T*(a-i3));
	    mean_testfejl(i3,i2)=mean(testfejl);
	    std_testfejl(i3,i2) =std(testfejl);  
	  end % eigen
	end % antalsplit
      end % maxfilterlength  
      [min_fejl split_index]=min(mean_testfejl);
      [E,model_i]=min(min_fejl);
      model(v,1)=model_i;
    end 


    for i = 1:cX
      opt_model(i,:) = [modelindex(model(i),2) modelindex(model(i),3)];
    end
    numcoef = opt_model(:,1)'; 


    %  Calculating the FIR coefficients for the optimal model
    % Iterate over voxels
    for v = 1:cX
      if opt_model(v,1)~=0
	y = X(:,v);
	XT = Xm(:,1:opt_model(v,1)); 
	XTT = XT'*XT;
	[U,S,V] = svd(XTT);
	epsilon = S*ones(opt_model(v,1),1);
	SS = (1./(S+eps)).*(S>(epsilon(opt_model(v,2))-(1e10*eps)));
	preX = V*SS*U'*XT';
	result_efir(1:opt_model(v,1),v)=preX*y; 
      end
    end 













