function [alphaj,nj,betaj] = patjitDPint(aj,bj,rk,alphak,nk,betak)
% function [alphaj,nj,betaj] = patjitDPint(aj,bj,rk,alphak,nk,betak)
% function [alphaj,nj,betaj] = patjitDPint(aj,bj,struct with r,alpha,n,beta)
%
% used by PatternJitter.m
%
% a single step in the (reverse) dynamic programming algorithm for
% computing the sampling structure

% Copyright (C) 2008 Matthew T Harrison
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    See <http://www.gnu.org/licenses/> for a copy of the license.

%%%%%%%%%%%%%%%%%%%
% INPUT CHECKING
%%%%%%%%%%%%%%%%%%%

% check for struct input
if isstruct(rk)
    % copy
    alphak = rk.alpha;
    nk = rk.n;
    betak = rk.beta;
    rk = rk.r;
end

% change bj to the largest possible value
bj = min(bj,alphak+nk-rk-2);

% check to see if impossible
if isnan(alphak) || bj < aj
    alphaj = NaN; nj = NaN; betaj = NaN;
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPECIAL CASES
%%%%%%%%%%%%%%%%%%%%%%%%%%

% check for fixed output spike
if aj == bj
    alphaj = aj;
    nj = 1;
    betaj = 1;
    return
end

% now the interval is determined
alphaj = aj;
nj = bj-aj+1;

% check to see if no interaction (case E and fixed input spike)
if bj+rk < alphak
    betaj = (nj:-1:1).'/nj;
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAIN ALGORITHM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute betaj

% initialize
betaj = zeros(nj,1);

betasum = 0;
j2k = aj+rk-alphak+1;

% just in case for some reason betak is not normalized
betak1 = betak(1);  % this should be 1

% loop in reverse
for j = nj:-1:1
    % get the index into betak
    k = j + j2k; % should always be <= nk
    if k <= 1
        betasum = betasum + betak1;
    else
        betasum = betasum + betak(k);
    end
    betaj(j) = betasum;
end

% normalize
betaj = betaj ./ betasum;
