function [Ampl, Poisson] = lyngby_pois_toptim(x, t, kernelWidth, ...
    convType, spatialType)

% lyngby_pois_toptim   - Possion kernel optimization
%
%      function [Apml, Poisson] = lyngby_pois_toptim(x, t, ...
%          kernelWidth, convType)
%
%      Input:  x         Input to the linear system
%              t         Target output to the linear system
%              kernelwidth   Width of kernel 
%              convType      Convolution type
%              spatialType   
%
%      Output: Ampl      Amplitude of kernel
%              Poisson   Poisson parameter
%
%      Linear system identification with a Poisson kernel optimized in
%      the time domain.
%
%      See also LYNGBY, LYNGBY_PDF_POISSON, LYNGBY_POIS_MAIN,
%               LYNGBY_POIS_ERROR, LYNGBY_POIS_FORWN.
%
% $Id: lyngby_pois_toptim.m,v 1.4 2003/05/23 09:44:26 fnielsen Exp $


    [rX, cX] = size(x);
    [rT, cT] = size(t);

    if rX ~= rT
      error('x and t should have the same number of rows');
    end
    
    method = 1;
    
    % Default parameter and check
    if nargin < 3
      kernelWidth = 30;
    elseif (kernelWidth <= 0) | (kernelWidth > rX)
      error('kernelWidth should be: 0 < kernelWidth < rows of x');
    end
    if nargin < 4
      convType = 1;
    elseif ~any(convType == [1 2 3 4])
      error('convType should be either 1, 2, 3 or 4');
    end
    if nargin < 5
      spatialType = 1;
    elseif ~any(convType == [1 2])
      error('spatialType should be either 1, 2');
    elseif cT == 1
      spatialType = 1;
    end
    
    if convType == 1
      fitIndex = kernelWidth:rT;
    else
      fitIndex = 1:rT;
    end
    
    r = 0.6180;  % (sqrt(5)-1)/2;
    
    if sum(abs(t)) == 0 
      Ampl = 0;
      Poisson = 0;
      return
    end

    maxlag = max([ rT kernelWidth ]);
    nFFT = 2^nextpow2(maxlag);
    if rX < kernelWidth
      x(rT,:) = 0;
    end
    
    if method == 1 
      % Golden ratio method
      
      % First we have to make sure that the curve between a and b is
      % convex. We start from 0 and examine tree values: 0, 1 and 2.

      Y = lyngby_pois_forwn(x, [0 1 2], kernelWidth, convType);
      if spatialType == 1
	% Individual Poisson and amplitude parameter for each voxel
	
	Ampl = sum(Y.^2) .\ sum(Y .* (t(fitIndex) * ones(1,3)));
	Res = lyngby_pois_error(t, ones(length(fitIndex), 1) * Ampl ...
	    .* Y);
	
      else
	% Individual amplitude common poisson parameter for voxels
	
	indexT3 = [1:cT 1:cT 1:cT];
	indexY3 = [ ones(1,cT) 2*ones(1,cT) 3*ones(1, cT) ];
	sumY = sum(Y.^2);
	Ampl = sumY(:, indexY3) .\ sum( Y(:,indexY3) .* ...
	    t(fitIndex, indexT3));
	Res = sum(reshape(lyngby_pois_error(t(:,indexT3), ...
	    ones(length(fitIndex), 1) * Ampl .* Y(:,indexY3)), cT, 3));
	
      end
	
      if Res(1) < Res(2)
	% What should be done here? It is not certain that the curve
	% is convex
	a = 0;
	b = 1;
      elseif Res(2) < Res(3)
	a = 0;
	b = 2;
      elseif Res(1) == Res(2) & Res(1) == Res(3)
	Poisson = 0;
	Ampl = Ampl(1);
	return
      else
	ResOld = Res(3);
	for n = 3:rX

	  Y = lyngby_pois_forwn(x, n, kernelWidth, convType);
	  if spatialType == 1
	    Ampl = sum(Y.^2) .\ sum(Y .* (t(fitIndex) * ones(1,1)));
	    ResNew = lyngby_pois_error(t, ones(length(fitIndex), 1) * ...
		Ampl .* Y);
	  else
	    indexY = [ ones(1,cT)];
	    sumY = sum(Y.^2);
	    Ampl = sumY(:, indexY) .\ sum( Y(:,indexY) .* ...
		t(fitIndex, :));
	    ResNew = sum(lyngby_pois_error(t, ...
		ones(length(fitIndex), 1) * Ampl .* Y(:,indexY)));
	  end

	  if ResNew > ResOld
	    a = n-2;
	    b = n;
	    break
	  else
	    ResOld = ResNew;
	  end
	end
      end

      l1 = a + (1-r)*(b-a);
      l2 = a + r*(b-a);

      Y = lyngby_pois_forwn(x, [l1 l2], kernelWidth, convType);
      if spatialType == 1
	Ampl = sum(Y.^2) .\ sum(Y .* (t(fitIndex) * ones(1,2)));
	Res = lyngby_pois_error(t, ones(length(fitIndex), 1) * Ampl ...
	    .* Y);
      else
	indexT2 = [1:cT 1:cT ];
	indexY2 = [ ones(1,cT) 2*ones(1,cT) ];
	sumY = sum(Y.^2);
	Ampl = sumY(:, indexY2) .\ sum( Y(:,indexY2) .* ...
	    t(fitIndex, indexT2));
	Res = sum(reshape(lyngby_pois_error(t(:,indexT2), ...
	    ones(length(fitIndex), 1) * Ampl .* Y(:,indexY2)), cT, 2));
      end

      cont = 20;
      while cont
	if Res(2) > Res(1)
	  b = l2;
	  l1 = a + (1-r) * (l2-a);
	  l2 = l1;
	  Res(2) = Res(1);
	  l = l1;
	else
	  a = l1;
	  l1 = l2;
	  l2 = b - (1-r)*(b-l1);
	  Res(1) = Res(2);
	  l = l2;
	end
	
	Y = lyngby_pois_forwn(x, l, kernelWidth, convType);
	if spatialType == 1
	  Ampl = sum(Y.^2) .\ sum(Y .* (t(fitIndex) * ones(1,1)));
	  if Res(2) > Res(1)
	    Res(1) = lyngby_pois_error(t, ones(length(fitIndex), 1) * Ampl .* Y);
	  else
	    Res(2) = lyngby_pois_error(t, ones(length(fitIndex), 1) * Ampl .* Y);
	  end
	else
	  indexY = [ ones(1,cT)];
	  sumY = sum(Y.^2);
	  Ampl = sumY(:, indexY) .\ sum( Y(:,indexY) .* ...
	      t(fitIndex, :));
	  if Res(2) > Res(1)
	    Res(1) = sum(lyngby_pois_error(t, ...
		ones(length(fitIndex), 1) * Ampl .* Y(:,indexY)));
	  else
	    Res(2) = sum(lyngby_pois_error(t, ...
		ones(length(fitIndex), 1) * Ampl .* Y(:,indexY)));
	  end
	end
	  
	
	cont = cont - 1;
      end

      if Res(2) > Res(1)
	Poisson = l1;
      else
	Poisson = l2;
      end
    
      Y = lyngby_pois_forwn(x, [Poisson], kernelWidth, convType);
      if spatialType == 1
	Ampl = sum(Y.^2) .\ sum(Y .* (t(fitIndex) * ones(1,1)));
	Res = lyngby_pois_error(t, ones(length(fitIndex), 1) * Ampl .* Y);
      else
	indexY = [ ones(1,cT)];
	Poisson = Poisson(:,indexY);
	sumY = sum(Y.^2);
	Ampl = sumY(:, indexY) .\ sum( Y(:,indexY) .* ...
	    t(fitIndex, :));
	Res = sum(lyngby_pois_error(t, ...
	    ones(length(fitIndex), 1) * Ampl .* Y(:,indexY)));
      end
    end   
 




















