function [pvalue] = pvt_alpha(n,N,S,METH, CLOCK_MAX, MC_ITER)
% pvt_alpha(n,N,s) -- Computes the p-value of the Poisson variability test.
%	Parameters
%	  1) n = number of trials
%	  2) N = total number of spikes (across all trials)
%	  3) S = sum of squares (across all trials)
%	  4) METH = method of computation
%	    Possible values:
%		  0 = AUTOMATIC (Default)
%		  1 = Dynamic Programming (This is an exact p-value)
%		  2 = Monte Carlo (Approximate p-value)
%		  The AUTOMATIC method tries to compute the p-value by dynamic programming
%		  until CLOCK_MAX seconds elapses, at which point the p-value is computed by MC 
% 		  approximation
%	  5) CLOCK_MAX = maximum number of seconds to solve by dynamic programming
%		  Default value = 20 seconds
%	  6) MC_ITER = number of Monte Carlo samples to use if computed by Monte Carlo
%		  Default value = 10,000
% Two computational methods are available. Dynamic Programming produces an exact p-value but
% can be inefficient if N or S are too large. Monte Carlo samples produce an approximate
% p-value efficiently; the accuracy of the approximation can be increased arbitrarily
% by increasing MC_ITER, the number of Monte Carlo samples.
%
% Only the first three parameters are required.
%
% by Asohan Amarasingham; send comments to asohan@dam.brown.edu

if isempty(who('METH')), METH=0; end
if isempty(who('CLOCK_MAX')), CLOCK_MAX=20; end
if isempty(who('MC_ITER')), MC_ITER=10000; end

if METH==2, ['Monte Carlo'], pvalue=pvtalpha_mc(n,N,S,MC_ITER); return, end
	
if n==1      % if n=1 then the sum of squares=N (w.p.1)
    if S>=N
        prob=1;
    else
        prob=0;
    end
    return
end

if n==2     % if n=2, compute from binomial probability
    
    % for n=2, note x_1^2 + x_2^2 <= k is equivalent to   n/2-k' <= x_1 <= n/2+k', for some k'
    % i.e. find smallest integer c s.t. (x+c)^2 + (N- (x+c) )^2 > S
    x=N/2;
    if round(N/2)~=N/2, j=.5; else, j=0, end
    for c=0:ceil(1+N/2)
        if (x+j+c)^2 + (N - (x+j+c))^2 > S
            break
        end
    end
    
    if c==0
        prob=0;
    else
        c=c-1;
        leftprob=binocdf(x+c+j,N,.5);
        if x-c-j-1<0
            rightprob=0;
        else
            rightprob=binocdf(x-c-j-1,N,.5);        
        end        
        prob=leftprob-rightprob;
    end
    return
    
end
        
% initial loop..

alpha=sparse(N+1,S+1);  % index is shifted by 1 to accomodate 0 indices..

tic

for s1=0:N
    if s1^2 > S
        break
    else
        for s2=s1:N
            if s1^2 + (s2-s1)^2 > S  % ie if t2>S
                break
            else
                t2=s1^2 + (s2-s1)^2;
                alpha(s2+1,t2+1) = alpha(s2+1,t2+1) + coeff( s2,s1,1/n,s2 );
            end
        end
    end
	if toc>CLOCK_MAX, ['Monte Carlo'], pvalue=pvtalpha_mc(n,N,S,MC_ITER); return, end
end

% induction loop

for t=2:n-2

     
	beta=sparse(N+1,S+1);  % again, index shifted by 1.. beta matrix

	                           
	% find non-empty values of alpha matrix..
	% notation in mind: (sm,tm) is (s_m,t_m).. (sm1,tm1) is (s_{n+1},t_{n+1})..
	[i,j]=find(alpha);
	for cc=1:length(i)
        sm=i(cc)-1; tm=j(cc)-1;         % remember the index shift..
        for sm1=sm:N
            tm1=tm + (sm1-sm)^2;
            if tm1>S
                break
            else
                beta( sm1+1,tm1+1 ) = beta(sm1+1,tm1+1) + coeff(sm1,sm,1/n,sm1-sm)*alpha( i(cc),j(cc) );

            end
        end
		if toc>CLOCK_MAX, ['Monte Carlo'], pvalue=pvtalpha_mc(n,N,S,MC_ITER); return, end
	end

    alpha=beta;

    
end

% final loop
beta=sparse(N+1,S+1);
[i,j]=find(alpha);
for cc=1:length(i)
    sm=i(cc)-1; tm=j(cc)-1;         % remember the index shift..
    tm1=tm + (N-sm)^2;
    if tm1 <= S
        beta( N+1,tm1+1 )= beta(N+1,tm1+1) + coeff(N,sm,1/n,N-sm)*alpha( i(cc),j(cc) );
     
    end
end


% sum out t_m and you're done
pvalue=sum(sum(beta));
pvalue=full(pvalue);


return

%--------------------------------

function val = coeff(N,k,p,P)
% computes the coefficient nchoosek(N,k)*(p)^P..
% (using log gamma for more efficient numerics..)
% use gammaln to keep numbers of same order, then exponentiate..
val = exp( gammaln(N+1)-gammaln(k+1)-gammaln(N-k+1) + P*log(p) );
return

%---------------------------------
function [pvalue] = pvtalpha_mc(n,N,S,MC_ITER) 
% Returns the p-value of the Poisson variability test
% by Monte Carlo approximation (using MC_ITER),
% where n is the number of trials, N is the total number of spikes (across trials)
% and S is the sum of squares of the spike counts across trials (see text)

MC_ITER=10000;								     % MC_ITER specifies number of Monte Carlo samples to use
h=hist( ceil( unifrnd(0,1,N,MC_ITER)*n ), n );   % generate multinomial samples
i=sum(h.*h,1);  									 % compute sum of squares
pvalue=sum(i<=S)/MC_ITER;							 % output pvalue as proportion of samples with sum of squares <= k

return