function DP=AndrzejakKreuzDpointprocess(spiketimes,Q,q,s,W,dt,dchoice)

% For further information, please refer to AndrzejakKreuzReadme.pdf and to
% Andrzejak RG, Kreuz T: Characterizing unidirectional couplings between point processes and flows. EPL, 96 (2011) 50012.

N=fix((Q-q)/(s)+1);


num_mass=(mod(dchoice,2)>0)+(mod(dchoice,4)>1);

DP=zeros(num_mass,N,N);


for wc1=W+1:N-1
    winspikes1=spiketimes(spiketimes>0 & spiketimes<Q-wc1*s);
    winspikes2=spiketimes(spiketimes>0+wc1*s & spiketimes<Q)-wc1*s;
    duration=Q-wc1*s;
    
    if dchoice>0
        spikemat=zeros(2,max([length(winspikes1) length(winspikes2)]));
        spikemat(1,1:length(winspikes1))=winspikes1;
        spikemat(2,1:length(winspikes2))=winspikes2;
        [dummy,isi,time_profiles]=f_SPIKE_ISI_distance_new(spikemat,0,duration,dt,mod(dchoice,4));
        % depending on your matlab version, you can replace dummy by the
        % tilde symbol.
        
        for wc2=1:N-wc1
            wmin=(wc2-1)*s;
            wmax=wmin+q;
            
            if dchoice==1 || dchoice==2
                DP(1,wc1+wc2,wc2)=f_windowing(isi,time_profiles(1,:),wmin,wmax,dt);
                DP(1,wc2,wc1+wc2) = DP(1,wc1+wc2,wc2);
            elseif dchoice==3
                DP(1,wc1+wc2,wc2)=f_windowing(isi,time_profiles(2,:),wmin,wmax,dt);
                DP(1,wc2,wc1+wc2)  = DP(1,wc1+wc2,wc2);
                DP(2,wc1+wc2,wc2)=f_windowing(isi,time_profiles(1,:),wmin,wmax,dt);
                DP(2,wc2,wc1+wc2) = DP(2,wc1+wc2,wc2);                
            end            
        end
    end
    
end
DP = squeeze(DP);


function ave=f_windowing(fx,fy,wmin,wmax,dt)

ctfx=cumsum([0 fx]);
pfx=unique([ctfx(1:end) ctfx(2:end-1)-dt/100]);

pfy=reshape([fy; fy],1,length(fy)*2);
first_winspike_ind=find(pfx>=wmin,1,'first');
last_winspike_ind=find(pfx<=wmax,1,'last');
win_pfx=pfx(first_winspike_ind:last_winspike_ind);
win_pfy=pfy(first_winspike_ind:last_winspike_ind);
if wmin<pfx(first_winspike_ind)     % interval to first spike
    win_pfx=[wmin win_pfx];
    win_pfy=[pfy(first_winspike_ind-1) win_pfy];
end
if wmax>pfx(last_winspike_ind)      % interval after last spike
    win_pfx=[win_pfx wmax];
    win_pfy=[win_pfy pfy(last_winspike_ind+1)];
end
win_pfx_isi=diff(win_pfx);
win_pfy_isi=win_pfy(1:end-1);

ave=sum(win_pfy_isi.*repmat(win_pfx_isi,size(win_pfy_isi,1),1),2)/sum(win_pfx_isi);

return

