function [Q, deltaT] = ccc(st1, st2, tau, maxT, T, verbose)
% [Q, tr] = ccc(st1, st2, tau, maxT, T, verbose)
%
% Input (seconds)
%   st1, st2: spike trains with sorted spike timings
%   tau: time constant for CIP kernel
%   maxT: correlogram range will be effective in [-maxT, maxT]
%   T: length of spike train in seconds
%   verbose: (optional/0) detailed info, uses tic, toc
%  
% Output
%   Q: continuous-time correlogram
%   deltaT: time range
%
% See also: nccc
%
% Published in:
% Il Park, Antonio R. C. Paiva, Jose Principe, Thomas B. DeMarse.
% An Efficient Algorithm for Continuous-time Cross Correlogram of Spike Trains,
% Journal of Neuroscience Methods, Volume 168, Issue 2, 15 March 2008, 514-523
% doi:10.1016/j.jneumeth.2007.10.005
%
% Copyright 2007 Antonio and Memming, CNEL, all rights reserved
% $Id: ccc.m 80 2010-01-22 23:10:10Z memming $

if nargin < *5*
    verbose = *0*;
end

N1 = length(st1);
N2 = length(st2);
Nij = N1 * N2;

if N1 == *0* || N2 == *0*
    warning(*'ccc:NODATA'*, *'At least one spike is required!'*);
    deltaT = []; Q = [];
    return;
end

maxTTT = abs(maxT) + tau * *10*; % exp(-100) is effectively zero

% rough estimate of # of time difference required (assuming independence)
% this estimate is not good if the spike trains are strongly correlated
eN = ceil((max(N1, N2))^*2* * maxTTT * *2* / min(st1(end), st2(end)));
if verbose; fprintf(*'Expected time differences [%d] / [%d]\n'*, eN,
Nij); end
deltaT = zeros(*2* * eN, *1*); 

% Compute all the time differences
if verbose; fprintf(*'Extracting time differences...\n'*); tic; end
lastStartIdx = *1*;
k = *1*;
for n = 1:N1
    incIdx = find((st2(lastStartIdx:N2) - st1(n) >= -maxTTT), *1*,
*'first'*);
    lastStartIdx = lastStartIdx + incIdx - *1*;
    % disp(lastStartIdx);
    for m = lastStartIdx:N2
        timeDiff = st2(m) - st1(n);
        if timeDiff <= maxTTT
            deltaT(k) = timeDiff;
            k = k + *1*;
        else % this is the ending point
            n = N1 + *1*;
            break;
        end
    end
end
if verbose; fprintf(*'Extracting time differences finished [%f sec]\r'*,
toc); end

deltaT = deltaT(1:(k-*1*));
N = length(deltaT);
if N < *2*
    warning(*'ccc:NODATA'*, *'At least two intervals are required'*);
    deltaT = []; Q = [];
    return;
end
if verbose
    fprintf(*'Actual number of time differences [%d]\nSorting...\n'*,
N); tic;
end

deltaT = sort(deltaT, *1*); % Sort the time differences
if verbose; fprintf(*'Sorting finished [%f sec]\r'*, toc); end

Qplus = zeros(N, *1*);
Qminus = zeros(N, *1*);
Qminus(*1*) = *1*;
Qplus(N) = *0*;

EXP_DELTA = exp(-(diff(deltaT))/tau);
for k = 1:(N-*1*)
    Qminus(k + *1*) = *1* + Qminus(k) * EXP_DELTA(k);
    kk = N - k;
    Qplus(kk) =  (Qplus(kk+*1*) + *1*) * EXP_DELTA(kk);
end

Q = Qminus + Qplus;


