function vspar=sparband(x,nspar,nbin)
%function vspar=sparband(x,nspar,nbin)
%% Inputs
%%  x design time points
%%  nspar  number of spars, default 11
%%  nbin    number of binned locations
%% Outputs
%%   
%%   vspar=logspace(log10(sparmin),log10(sparmax),nspar)
%%   
%       
% 
% Copyright 26/12/2001, Jin-Ting Zhang, at NUS.

%% Set the defaults
if nargin<2|length(nspar)==0,nspar=41;end % Default for 41 spars


if nargin<3|length(nbin)==0, 
    nbin=min([length(x),301]);
    %nbin=301;
end  

% Bin the data (allows handling large data sets, and uneven designs)
 xmin = min(x) ;
 xmax = max(x) ;
 [bindat,bincent] = gplbinr([x,ones(size(x))],[xmin,xmax,nbin]) ;
 
 flagn0 = bindat(:,1) > 0 ;
         %  1 where bin is nonempty
 t= bincent(flagn0) ; 
 
 w = bindat(flagn0,1) ;
 n = sum(flagn0) ;
         %  number of nonemty bins
W=diag(w);


%% data transformation
tmin=t(1);tmax=t(n);range=tmax-tmin;
t=(t-tmin)/range; % transfered into [0,1]

%% Compute the related matrices R,Q,K      
   h=t(2:n)-t(1:n-1);
   hinv=h.^(-1);

% Calculate Q:nx(n-2)
   Q=diag([hinv(1:n-2);0],-1)-diag([0;hinv(1:n-2)+hinv(2:n-1);0],0);
   Q=Q+diag([0;hinv(2:n-1)],1);
   Q=Q(2:n-1,:)';
   
% Calculate R:(n-2)x(n-2)
   R=diag(h(2:n-2),-1)/6+diag(h(1:n-2)+h(2:n-1))/3+diag(h(2:n-2),1)/6;

% Calculate K: nxn matrix
   Rinv=inv(R);
   K=Q*Rinv*Q';
 
%% Eigen decomposition
temp=diag(w.^(-.5));
[U,D,V]=svd(temp*K*temp); %% 
[d,temp]=sort(abs(diag(D)));%%

   %% Compute the range of spar
   sparmax=(n/2.2-1)/d(4); %% the largest spar so that tr(S)=2.2,d(n)=max(d)
   sparmin=(n/(.8*n)-1)/d(n-1); %% the smallest spar so that tr(S)=.8*n
   
   vspar=logspace(log10(sparmin),log10(sparmax),nspar)/range^4;