function [alphaj,cj,xj,dj,mj,nj,betaj] = patjitDP(aj,bj,rk,alphak,ck,xk,dk,mk,nk)
% function [alphaj,cj,xj,dj,mj,nj,betaj] = patjitDP(aj,bj,rk,alphak,ck,xk,dk,mk,nk)
% function [alphaj,cj,xj,dj,mj,nj,betaj] = patjitDP(aj,bj,struct with r,alpha,c,x,d,m,n)
%
% 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.

% check for struct input
if isstruct(rk)
    % copy
    alphak = rk.alpha;
    ck = rk.c;
    xk = rk.x;
    dk = rk.d;
    mk = rk.m;
    nk = rk.n;
    rk = rk.r;
end

% check to see if impossible
if isnan(alphak(1)) || aj+rk >= alphak(end) || bj < aj
    alphaj = NaN; cj = NaN; xj = NaN; dj = NaN; mj = NaN; nj = NaN; betaj = NaN;
    return
end

% check for fixed output spike
if aj == bj
    alphaj = aj;
    cj = [];
    xj = [];
    dj = [];
    mj = [];
    nj = 0;
    betaj = [];
    return
end
    
% check to see if no interaction (case E)
if bj+rk <= alphak(1)
    alphaj = [aj;bj];
    cj = [1; 0];
    xj = aj;
    dj = bj-aj;
    mj = 2;
    nj = 1;
    betaj = 1;
    return
end

% check for fixed input spike 
if nk == 0
    bj = min(bj,alphak(1)-rk);
    alphaj = [aj;bj];
    cj = [1; 0];
    xj = aj;
    dj = bj-aj;
    mj = 2;
    nj = 1;
    betaj = 1;
    return
end

% transform to r = 0
aj = aj+rk;
bj = min(bj+rk,alphak(nk+1));  % deal with case A
% aj < bj still, otherwise caught by impossible case above

% skip any pieces after bj
for kmax = nk:-1:1
    if alphak(kmax) < bj
        break
        % must exit loop here
    end
    % exit here is caught by case E above 
end
% kmax is largest k such that alphak(k) < bj 

% skip any pieces before aj
for kmin = 1:nk+1
    if alphak(kmin) > aj
        break
        % must exit loop here
    end
    % exit here is caught by impossible above
end
% kmin is smallest k such that alphak(k) > aj 

% allocate
nj = kmax-kmin+2;
alphaj = [aj; alphak(kmin:kmax); bj];
mmax = max(mk(max(1,kmin-1):kmax));
cj = zeros(mmax+1,nj);
xj = zeros(nj,1);
dj = zeros(nj,1);
mj = zeros(nj,1);
betaj = zeros(nj,1);

% iterate
j2k = kmin - 2;
betasum = 0;

for j = nj:-1:1

    k = j + j2k;
    
    % case D
    if k == 0  % j == 1 && kmin == 1
        c = alphaj(2)-aj;
        cj(1,1) = c;
        cj(2,1) = betasum;
        xj(1) = aj;
        dj(1) = c;
        mj(1) = 2;
        betasum = betasum + c;
        betaj(1) = betasum;
        break
    end
    
    % cases B and C
    % grab the relevant constants
    m = mk(k);
    x = xk(k);
    d = dk(k);
    b = 1-(alphaj(j+1)-x)/d;
    a = 1-(alphaj(j)-x)/d;
    % integrate the coefficients
    c = ck(1:m,k)*d./((m:-1:1).');
    % remember the coefficients
    cj(1:m,j) = c;
    % get the two parts of the integral constant term
    bpart = mypolyval(c,b,m)*b;
    apart = mypolyval(c,a,m)*a;
    % get the constant term
    cj(m+1,j) = betasum - bpart;
    % update and record the running integral
    betasum = betasum + (apart - bpart);
    betaj(j) = betasum;
    % record the remaining quantities
    xj(j) = x;
    dj(j) = d;
    mj(j) = m + 1;
end

% normalize
betaj = betaj / betasum;
cj = cj / betasum;

% transform back from r = 0
alphaj = alphaj - rk;
xj = xj - rk;
