function [plower, pupper] = findbounds(NR, n, jconsec)

%Matlab code for "A Set Probability Technique for Detecting
%Relative Time Order Across Multiple Neurons", Neural Computation,
%In Press.

%Takes 3 arguments (NR, n, jconsec)
%  NR is the reference sequence length 
%  n is the word length
%   j is the number of consecutive increasing elements

%Returns plower and pupper = upper and lower bounds on 
%Pr(j or more strictly increasing)
%If an exact solution is possible then pupper=plower

%Anne C Smith April, 2005


res = [];
cnt = 0;

%compute pj = Pr(j in a row or more are increasing)

j      = jconsec;

%Eq. 1
pj     = factorial(NR)/factorial(j)/factorial(NR-j)/(NR^j);

%3 cases depending on relative sizes of n and j
    
if(n==j)

        pfinal = pj;      %subsequence length same as sequence
        gfinal = pj;
        
elseif( 2*j>n )           %case of a very short sequence-exact
      
         pjp1   = factorial(NR)/factorial(j+1)/factorial(NR-j-1)/(NR^(j+1));
  
         pfinal = pj + (n-j)*(pj-pjp1);   % Eq. 8
         gfinal = pfinal;                 % upper and lower bound the same
         
elseif( 2*j<=n )           %deals with longer sequences
        
    pjp1   = factorial(NR)/factorial(j+1)/factorial(NR-j-1)/(NR^(j+1));
     
    f = []; g = [];  % f is upper bound, g is lower bound
    f(1) = pj; 
    g(1) = pj;
    
	for r = 2:j
	    f(r) = pj-pjp1;
        g(r) = f(r);
    end
      
	for r = j+1:n-j+1;
	    xx   = 1:r-j-1;
	    s    = sum(f(xx));
            
        f(r) = (pj-pjp1) * (1-s);           %upper bound, Eq. 15
        g(r) = (pj-pjp1) - (r-j)*pj*pj;     %lower bound, Eq. 16
        if(g(r) < 0)
            g(r) = 0;
        end
                 
    end
 
    pfinal = sum(f);
    gfinal = sum(g);
end 

if(pfinal > 1)  % check pfinal not slightly over 1
    pupper = 1;
else
    pupper = pfinal;
end

plower = gfinal;