% Calculates SPIKE- and ISI-Distance from two or more spike trains
%
% Information can be found online under
%
% "http://www.fi.isc.cnr.it/users/thomas.kreuz/Source-Code/Spike-Sync.html" and/or in these papers:
%
% Kreuz T, Chicharro D, Greschner M, Andrzejak RG: Time-resolved and time-scale adaptive measures of spike train synchrony. J Neurosci Methods, 195, 92 (2011)
% 
% Kreuz T, Chicharro D, Andrzejak RG, Haas JS, Abarbanel HDI: Measuring multiple spike train synchrony. J Neurosci Methods 183, 287 (2009).
% 
% Kreuz T, Haas JS, Morelli A, Abarbanel HDI, Politi A: Measuring spike train synchrony. J Neurosci Methods 165, 151 (2007)
%
% function f_compute_gauss_smooth - Copyright: Jude Mitchell, http://www.snl.salk.edu/~jude/sfn2008/index.html
%
% For questions and comment_strings please contact me at "thomaskreuz (at) cnr.it".

function [mean_values,isi,time_profiles]=f_SPIKE_ISI_distance_new(spikes,tmin,tmax,dt,selection)

max_mat_size=5000000;

if nargin<4
    % dt=0.001;           % Define your standard value (has to be smaller than the minimum ISI, ISIs should be multiples of dt)
    dt=f_get_dt(spikes);  % or get it automatically from the precision of the spikes
end

% Choose order and size of the subplots; for multiple figures use a matrix ( in which each row corresponds to one figure)
% Use "subplot_posi" to select the order of the subplots. Use 0 if no subplot is needed.
% Make sure that the range from 1 to the desired number of subplots is covered
% Use "subplot_size" to select the size of the selected subplots.
% Its order refers to the subplots selected in subplot_posi, not to the reference-order
% (i.e., its length should be equal to the number of subplot selected in 'subplot_posi')

%
% The subplot "Stimulus" allows you to add any representation of the stimulus that triggered the spike trains (put 0 if not needed)

% Reference-Order (for position variable only):
% Stimulus  Spikes  ISI <ISI>  Rate <Rate>  I^A I^M  SC^A SC^M  S^A S^M PSTH GPSTH

%subplot_posi=[1 2 3 4 5 6 7 8 9 10 11 12 13 14];
%subplot_posi=[0 1 0 2 0 0 3 0 0 0 4 0 0 0];
if selection==1
    subplot_posi=[0 0 0 0 0 0 0 0 1 0 0 0 0 0];
elseif selection==2
    subplot_posi=[0 0 0 0 0 0 1 0 0 0 0 0 0 0];
else
    subplot_posi=[0 0 0 0 0 0 1 0 2 0 0 0 0 0];
end
subplot_size=ones(1,length(subplot_posi(subplot_posi>0)));

window_mode=1;           % 1-all (recording limits),2:outer spikes,3-inner spikes,4-smaller analysis window

plot_mode=0;             % plotting figure? (0-no,1-yes)

num_fig=1;               % number of (first) figure
comment_string='';       % will be part of the title and the filename
colors='kgkgrg';         % Colors for 1:Stim/STs/ISI,2:Stim/STs/ISI-outside,3:Mean values,4:Mean values-outside,5:Moving averages,6:Moving averages-outside
font_size=10;            % Font size
line_width=1;            % Line width
title_string='';         % will be part of the title
timeunit_string='[s]';   % time unit

publication=0;           % 0-no,1-yes
print_mode=1;            % saving to postscript file? (0-no,1-yes)
imagespath='.\';
filename='SPIKE-ISI-Distance';  % filename of this postscript file

log_mode=0;              % 0-no log-transform,1-log-transform (for the multivariate ISI-Diversity)

ma_mode=0;       % moving average mode: 0-no,1-only,2-both

spike_mao=10;    % order of the moving average
time_mao=45;     % order of the moving average for SPIKE-2 distances
num_bins=100;    % number of bins for the PSTH
gs_window=5;    % window for Gaussian smoothing of PSTH (ceil(num_bins/100))

if ma_mode==0 nma=1; elseif ma_mode==1 nma=2; else nma=[1 2]; end %#ok<SEPEX>


num_trains=size(spikes,1);
if num_trains<2
    error('It does not make much sense to do quantify variability for just one spike train!')
end
num_pairs=(num_trains*num_trains-num_trains)/2;
num_pspikes=zeros(1,num_trains);
for trac=1:num_trains
    if any(spikes(trac,:))
        num_pspikes(trac)=find(spikes(trac,:),1,'last');
    end
end
max_num_pspikes=max(num_pspikes);

pspikes=zeros(num_trains,max_num_pspikes);          % original spikes used for plotting
for trac=1:num_trains
    pspikes(trac,1:num_pspikes(trac))=floor(sort(spikes(trac,1:num_pspikes(trac)))/dt)*dt;                      % ********************************
end

if nargin<3 && exist('tmin','var')==0  && exist('tmax','var')==0
    tmin=min(min(pspikes));
    tmax=max(max(pspikes));
    trange=tmax-tmin;
    tmin=floor((tmin-0.02*trange)/dt)*dt;
    tmax=floor((tmax+0.02*trange)/dt)*dt;
else
    tmin=floor(tmin/dt)*dt;
    tmax=floor(tmax/dt)*dt;
end

itmin=tmin;
itmax=tmax;
if window_mode==1                                                               % all (recording limits)
    wmin=itmin; wmax=itmax;
else
    if window_mode==2                                                           % outer spikes (overall)
        wmin=min(min(min(pspikes(pspikes~=0))));
        wmax=max(max(max(pspikes(pspikes~=0))));
    elseif window_mode==3                                                       % inner spikes (overall)
        wmin=max(pspikes(1:num_trains,1));
        wmax=inf;
        for trac=1:num_trains
            if num_pspikes(trac)>0
                if pspikes(trac,num_pspikes(trac))<itmax
                    wmax=pspikes(trac,num_pspikes(trac));
                end
            end
        end
    elseif window_mode==4                                                          % select
        wmin=tmin+0.2*(tmax-tmin);
        wmax=tmax-0.2*(tmax-tmin);

        %wmin=5;
        %wmax=10;
    end

    wmin=floor(wmin/dt)*dt;
    wmax=floor(wmax/dt)*dt;
    if wmin>=wmax
        disp(' '); disp(' '); error ('********** Error: W_max must be larger than W_min !!!!!')
    end
    if wmin<tmin
        disp(' '); disp(' '); error ('********** Error: W_min must not be smaller than T_min !!!!!')
    end
    if wmax>tmax
        disp(' '); disp(' '); error ('********** Error: W_max must not be larger than T_max !!!!!')
    end
end
itmin=floor(itmin/dt)*dt;
itmax=floor(itmax/dt)*dt;
itrange=itmax-itmin;


result_string='';
num_ispikes=zeros(1,num_trains);
ispikes=zeros(num_trains,max_num_pspikes+2);
for trac=1:num_trains
    dspikes=pspikes(trac,1:num_pspikes(trac));
    dspikes=dspikes(dspikes>=itmin & dspikes<=itmax);
    if any(dspikes)
        num_ispikes(trac)=find(dspikes,1,'last')+2;
        if isempty(find(dspikes(1:num_ispikes(trac)-2)==itmin,1))
            ispikes(trac,1:num_ispikes(trac))=[itmin dspikes(1:num_ispikes(trac)-2) itmax];
        else
            dummy=dspikes(1:num_ispikes(trac)-2);
            dummy2=[itmin dummy(dummy~=itmin) itmax];
            num_ispikes(trac)=length(dummy2);
            ispikes(trac,1:num_ispikes(trac))=dummy2;
        end
    else
        num_ispikes(trac)=2;
        ispikes(trac,1:num_ispikes(trac))=[itmin itmax];
    end
end
max_num_ispikes=max(num_ispikes);
ispikes=ispikes(1:num_trains,1:max_num_ispikes);

num_coins=zeros(num_trains,max_num_ispikes);
num_uspikes=zeros(1,num_trains);
uspikes=zeros(num_trains,max_num_ispikes);
for trac=1:num_trains
    uispikes=unique(ispikes(trac,1:num_ispikes(trac)));
    for uic=1:length(uispikes)-1
        num_coins(trac,uic)=sum(ispikes(trac,1:num_ispikes(trac))==uispikes(uic));
    end
    num_uspikes(trac)=length(uispikes);
    uspikes(trac,1:num_uspikes(trac))=uispikes;
end
max_num_uspikes=max(num_uspikes);
uspikes=uspikes(1:num_trains,1:max_num_uspikes);

for trac=1:num_trains
    uspikes(trac,num_uspikes(trac)+1:max_num_uspikes)=987654321.123456789;
end


mean_values=[];
time_profiles=[];
if any(subplot_posi)                                                        % ISI-SPIKE
    [mall_spikes,mall_indy]=sort(reshape(uspikes',1,numel(uspikes)));
    mall_labels=fix((mall_indy-1)/max_num_uspikes)+1;
    mall_labels=mall_labels(mall_spikes~=987654321.123456789);
    mall_spikes=mall_spikes(mall_spikes~=987654321.123456789);
    mall_labels(1:num_trains)=0;
    mall_labels(end-num_trains+1:end)=0;

    all_spikes=mall_spikes(num_trains:end-num_trains+1);
    all_labels=mall_labels(num_trains:end-num_trains+1);
    all_isi=diff(all_spikes);
    num_all_isi=length(all_isi);
    isi=all_isi(all_isi>0);
    num_isi=length(isi);
    isi_pos=cumsum([itmin isi(1:end-1)])+isi/2;
    cum_isi=all_spikes(1)+[0 cumsum(isi)];
    indy=find(cum_isi<=itmin,1,'last'):find(abs(cum_isi-itmax)<0.0000000001,1,'first')-1;

    isis=zeros(num_trains,max_num_uspikes-1);
    for trac=1:num_trains
        isis(trac,1:num_uspikes(trac)-1)=diff(uspikes(trac,1:num_uspikes(trac)));
        isis(trac,isis(trac,:)~=0)=isis(trac,isis(trac,:)~=0)./num_coins(trac,isis(trac,:)~=0);
    end
    all_ints=zeros(num_trains,num_all_isi);
    all_rates=zeros(num_trains,num_all_isi);
    for trac=1:num_trains
        ivs=[1 find(all_labels==trac)];
        ive=[ivs(2:num_uspikes(trac)-1)-1 num_all_isi];
        for ic=1:num_uspikes(trac)-1
            all_ints(trac,ivs(ic):ive(ic))=isis(trac,ic);
        end
        all_rates(trac,all_ints(trac,:)>0)=1./all_ints(trac,all_ints(trac,:)>0);
    end

    ints=all_ints(:,all_isi>0);
    rates=all_rates(:,all_isi>0);
    mean_ints=mean(ints,1);
    mean_rates=mean(rates,1);

    if window_mode==2
        wmin=max([itmin wmin]);
        wmax=min([itmax wmax]);
    end
end


mat_size=num_pairs*num_all_isi;
if mat_size>max_mat_size
    if max_mat_size<num_pairs
        max_mat_size=num_pairs;
    end
    max_num_all_isi=fix(max_mat_size/num_pairs);
    num_runs=ceil(num_isi/max_num_all_isi);
    run_lengths=[max_num_all_isi*ones(1,num_runs-1) num_isi-max_num_all_isi*(num_runs-1)];
    run_ends=cumsum(run_lengths);
    run_starts=[1 run_ends(1:end-1)+1];
else
    num_runs=1;
    run_lengths=num_isi;
    run_starts=1;
    run_ends=num_isi;
end


if subplot_posi(7)                                         % ISI-Bi
    bi_ints=zeros(1,num_isi);
    run_sum_all_isi=zeros(1,num_runs);
    for ruc=1:num_runs
        run_range=run_starts(ruc):run_ends(ruc);
        run_xy_isi_dist=zeros(num_pairs,run_lengths(ruc));
        pac=0;
        for trac1=1:num_trains-1
            for trac2=trac1+1:num_trains
                pac=pac+1;
                dum=find(ints(trac1,run_range)<ints(trac2,run_range));
                run_xy_isi_dist(pac,dum)=ints(trac1,run_range(dum))./ints(trac2,run_range(dum))-1;
                dum2=find(ints(trac1,run_range)>=ints(trac2,run_range) & ints(trac1,run_range)~=0);
                run_xy_isi_dist(pac,dum2)=-(ints(trac2,run_range(dum2))./ints(trac1,run_range(dum2))-1);
            end
        end
        run_bi_ints=mean(abs(run_xy_isi_dist),1);
        run_indy=indy(indy>=run_starts(ruc) & indy<=run_ends(ruc))-run_starts(ruc)+1;
        run_sum_all_isi(ruc)=sum(all_isi(run_starts(ruc)-1+run_indy));
        bi_ints(run_range)=run_bi_ints;
    end
    clear run_xy_isi_dist
end

if subplot_posi(8)                                                         % ISI-Multi
    if log_mode==0
        multi_ints=std(ints,[],1)./mean_ints;
    else
        log_ints=log(round(ints/dt));
        multi_ints=mean(log_ints,1).*(exp(std(log_ints,[],1)./mean(log_ints,1))-1);
    end
end



if any(subplot_posi(9:12))                                                % SPIKE-C
    previs=zeros(num_trains,max_num_uspikes-1);
    follis=zeros(num_trains,max_num_uspikes-1);
    for trac=1:num_trains
        previs(trac,1:num_uspikes(trac)-1)=uspikes(trac,1:num_uspikes(trac)-1);
        follis(trac,1:num_uspikes(trac)-1)=uspikes(trac,2:num_uspikes(trac));
    end
    all_prev_spikes=zeros(num_trains,num_all_isi);
    all_foll_spikes=zeros(num_trains,num_all_isi);
    for trac=1:num_trains
        ivs=[1 find(all_labels==trac)];
        ive=[ivs(2:num_uspikes(trac)-1)-1 num_all_isi];
        for ic=1:num_uspikes(trac)-1
            all_prev_spikes(trac,ivs(ic):ive(ic))=previs(trac,ic);
            all_foll_spikes(trac,ivs(ic):ive(ic))=follis(trac,ic);
        end
    end
    prev_spikes=all_prev_spikes(:,all_isi>0);
    foll_spikes=all_foll_spikes(:,all_isi>0);
    mean_prev_spikes=mean(prev_spikes,1);
    mean_foll_spikes=mean(foll_spikes,1);

    if subplot_posi(9)                                                  % SPIKE-Bi
        spike_diffs_C=zeros(1,num_isi);
        for ruc=1:num_runs
            run_range=run_starts(ruc):run_ends(ruc);

            run_prev_diffs=zeros(num_pairs,run_lengths(ruc));
            run_foll_diffs=zeros(num_pairs,run_lengths(ruc));
            run_bi_sum_ints=zeros(num_pairs,run_lengths(ruc));
            run_bi_mean_prev_spikes=zeros(num_pairs,run_lengths(ruc));
            run_bi_mean_foll_spikes=zeros(num_pairs,run_lengths(ruc));

            pac=0;
            for trac1=1:num_trains-1
                for trac2=trac1+1:num_trains
                    pac=pac+1;
                    run_prev_diffs(pac,1:run_lengths(ruc))=abs(prev_spikes(trac1,run_range)-prev_spikes(trac2,run_range));
                    run_foll_diffs(pac,1:run_lengths(ruc))=abs(foll_spikes(trac1,run_range)-foll_spikes(trac2,run_range));
                    run_bi_mean_prev_spikes(pac,1:run_lengths(ruc))=(prev_spikes(trac1,run_range)+prev_spikes(trac2,run_range))/2;
                    run_bi_mean_foll_spikes(pac,1:run_lengths(ruc))=(foll_spikes(trac1,run_range)+foll_spikes(trac2,run_range))/2;
                    run_bi_sum_ints(pac,1:run_lengths(ruc))=(ints(trac1,run_range)+ints(trac2,run_range));
                end
            end

            run_bi_pos_prev_diff=zeros(num_pairs,run_lengths(ruc));
            run_bi_pos_foll_diff=zeros(num_pairs,run_lengths(ruc));
            pac=0;
            for trac1=1:num_trains-1
                for trac2=trac1+1:num_trains
                    pac=pac+1;
                    run_bi_pos_prev_diff(pac,1:run_lengths(ruc))=isi_pos(run_range)-run_bi_mean_prev_spikes(pac,1:run_lengths(ruc));
                    run_bi_pos_foll_diff(pac,1:run_lengths(ruc))=run_bi_mean_foll_spikes(pac,1:run_lengths(ruc))-isi_pos(run_range);
                end
            end
            normy=(run_bi_sum_ints/2+(run_bi_sum_ints==0)).^2;
            run_prefo_diffs_C=(run_prev_diffs.*run_bi_pos_foll_diff+run_foll_diffs.*run_bi_pos_prev_diff)./normy;
            run_spike_diffs_C=mean(run_prefo_diffs_C,1);
            spike_diffs_C(run_range)=run_spike_diffs_C;
        end
    end

    if subplot_posi(10)                                                 % SPIKE-Multi
        std_prev_spikes=std(prev_spikes,[],1);
        std_foll_spikes=std(foll_spikes,[],1);

        pos_prev_diff=isi_pos-mean_prev_spikes;
        pos_foll_diff=mean_foll_spikes-isi_pos;
        normy=(mean_ints+(mean_ints==0)).^2;
        spike_stds_C=(std_prev_spikes.*pos_foll_diff+std_foll_spikes.*pos_prev_diff)./normy;
    end
end


if any(subplot_posi([11 12]))                                              % SPIKE
    time=round(itmin/dt)*dt:dt:round(itmax/dt)*dt;
    len=length(time);
    start=zeros(num_trains,max(num_uspikes)-1);
    ende=zeros(num_trains,max(num_uspikes)-1);
    for trac=1:num_trains
        for sc=1:num_uspikes(trac)-1
            start(trac,sc)=fix((uspikes(trac,sc)-itmin)/dt+1);
            ende(trac,sc)=fix((uspikes(trac,sc+1)-itmin)/dt);
        end
    end

    endv=round(cumsum(isi)/dt);
    startv=[1 endv(1:num_isi-1)+1];
    run_length=zeros(1,num_runs+1);
    pf_v_length=zeros(1,num_runs);


    if subplot_posi(11) || subplot_posi(12)                                                  % SPIKE-Bi

        previ=zeros(num_trains,len);
        folli=zeros(num_trains,len);
        for trac=1:num_trains
            for sc=1:num_uspikes(trac)-1
                previ(trac,start(trac,sc):ende(trac,sc))=0:dt:(ende(trac,sc)-start(trac,sc))*dt;            % distance to previous spike
                folli(trac,start(trac,sc):ende(trac,sc))=(ende(trac,sc)-start(trac,sc)+1)*dt:-dt:dt;        % distance to following spike
            end
        end

        if subplot_posi(11)
            spike_diffs=zeros(1,len);
            for ruc=1:num_runs
                if num_runs>1
                    disp([ruc num_runs])
                end
                run_range=run_starts(ruc):run_ends(ruc);        % overall counter
                run_length(ruc+1)=length(run_range);
                pf_v=startv(run_range(1)):endv(run_range(end));
                pf_v_length(ruc)=length(pf_v);

                run_pf_mean_previ=zeros(num_pairs,run_lengths(ruc));
                run_pf_mean_folli=zeros(num_pairs,run_lengths(ruc));
                pac=0;
                for trac1=1:num_trains-1
                    for trac2=trac1+1:num_trains
                        pac=pac+1;
                        run_pf_mean_previ(pac,1:length(pf_v))=(previ(trac1,pf_v)+previ(trac2,pf_v))/2;
                        run_pf_mean_folli(pac,1:length(pf_v))=(folli(trac1,pf_v)+folli(trac2,pf_v))/2;
                    end
                end

                run_prev_diffs=zeros(num_pairs,run_lengths(ruc));
                run_foll_diffs=zeros(num_pairs,run_lengths(ruc));
                pac=0;
                for trac1=1:num_trains-1
                    for trac2=trac1+1:num_trains
                        pac=pac+1;
                        run_prev_diffs(pac,1:run_lengths(ruc))=abs(prev_spikes(trac1,run_range)-prev_spikes(trac2,run_range));
                        run_foll_diffs(pac,1:run_lengths(ruc))=abs(foll_spikes(trac1,run_range)-foll_spikes(trac2,run_range));
                    end
                end

                run_ints=ints(:,run_range);
                run_pf_spike_diffs=zeros(num_pairs,pf_v_length(ruc));
                pac=0;
                for trac1=1:num_trains-1
                    for trac2=trac1+1:num_trains
                        pac=pac+1;
                        for vc=1:run_length(ruc+1)
                            v=startv(sum(run_length(1:ruc))+vc):endv(sum(run_length(1:ruc))+vc);
                            bi_sum_ints=sum(run_ints([trac1 trac2],vc));
                            normy=(bi_sum_ints/2+(bi_sum_ints==0)).^2;
                            run_pf_spike_diffs(pac,v-pf_v(1)+1)=(run_prev_diffs(pac,vc)*run_pf_mean_folli(pac,v-pf_v(1)+1)+...
                                run_foll_diffs(pac,vc)*(run_pf_mean_previ(pac,v-pf_v(1)+1) ))./normy; %
                        end
                    end
                end
                run_prefo_diffs=mean(run_pf_spike_diffs,1);
                spike_diffs(pf_v)=run_prefo_diffs;
            end
        end
    end

    if subplot_posi(12)                                                    % SPIKE-Multi
        pf_std_prevs=zeros(1,len);
        pf_std_folls=zeros(1,len);
        pf_mean_ints=zeros(1,len);
        for vc=1:num_isi
            pf_std_prevs(startv(vc):endv(vc))=std(prev_spikes(:,vc)+1,[],1);
            pf_std_folls(startv(vc):endv(vc))=std(foll_spikes(:,vc)+1,[],1);
            pf_mean_ints(startv(vc):endv(vc))=mean_ints(vc);
        end

        mean_previ=mean(previ,1);
        mean_folli=mean(folli,1);
        normy=(pf_mean_ints+(pf_mean_ints==0)).^2;
        spike_stds=(pf_std_prevs.*mean_folli+pf_std_folls.*mean_previ)./normy;
    end
end

if any(subplot_posi([13 14]))                                              % PSTH
    bin_width=(tmax-tmin)/num_bins;
    bins=tmin:bin_width:tmax;
    bin_centers=bins(1:num_bins)+(bins(2)-bins(1))/2;
    dpsth=zeros(1,num_bins+1);
    for stc=1:num_trains
        dpsth=dpsth+histc(pspikes(stc,1:num_pspikes(stc)),bins);
    end
    psth=dpsth(1:num_bins)/bin_width;
    maxpsth=max(psth);
    if subplot_posi(14)                                                   % Gaussian smoothed PSTH
        gs_psth=f_compute_gauss_smooth(psth,gs_window)';
        maxgspsth=max(gs_psth);
    end
end


singles=zeros(1,max(subplot_posi));
for suc=1:max(subplot_posi)
    singles(suc)=find(subplot_posi==suc,1,'first');                               % new subplot
end
num_subplots=length(singles);

if plot_mode>0
    if window_mode==1                                                              % limits for plotting
        pmin=itmin-0.02*itrange;
        pmax=itmax+0.02*itrange;
    else
        pmax2=-inf;
        for trac=1:num_trains
            if num_pspikes(trac)>0
                pmax2=max([pmax2 pspikes(trac,num_pspikes(trac))]);
            end
        end
        pmin2=min(pspikes(1:num_trains,1));
        prange2=pmax2-pmin2;
        pmin=pmin2-0.02*prange2;
        pmax=pmax2+0.02*prange2;
    end

    maxintval=max(max(ints(:,indy)));
    maxrateval=max(max(rates(:,indy)));


    if any(subplot_posi([4 5]))
        maxmeanintval=max(mean_ints);
    end
    if any(subplot_posi([6 7]))
        maxmeanrateval=max(mean_rates);
    end

    maxbival1=0; maxbival2=0; maxbival3=0;
    if any(subplot_posi([7 9 11]))
        if subplot_posi(7)
            maxbival1=max(bi_ints(indy));
        end
        if subplot_posi(9)
            maxbival2=max(spike_diffs_C(indy));
        end
        if subplot_posi(11)
            maxbival3=max(spike_diffs);
        end
        maxbival=max([maxbival1 maxbival2 maxbival3]);
    end
    maxmultival1=0; maxmultival2=0; maxmultival3=0;
    if any(subplot_posi([8 10 12]))
        if subplot_posi(8)
            maxmultival1=max(multi_ints(indy));
        end
        if subplot_posi(10)
            maxmultival2=max(spike_stds_C(indy));
        end
        if subplot_posi(12)
            maxmultival3=max(spike_stds);
        end
        maxmultival=max([maxmultival1 maxmultival2 maxmultival3]);
    end


    relsubplot_size=subplot_size(1:find(subplot_size>0,1,'last'));
    if length(relsubplot_size)~=max(subplot_posi)
        relsubplot_size=ones(1,max(subplot_posi));
    end

    doubles=intersect(setxor(1:length(subplot_posi),singles),find(subplot_posi>0));     % repeated subplot
    doublesref=singles(subplot_posi(doubles));                                                % corresponding new subplot

    regplotsize=num_subplots*1.1;
    normsubplot_size=relsubplot_size/sum(relsubplot_size);
    dumsubplot_size=normsubplot_size*regplotsize;
    subplot_start2=cumsum(dumsubplot_size);
    subplot_size2=diff([0 subplot_start2]);

    subplot_size=zeros(1,length(subplot_posi));               % normalized size of subplots
    subplot_size(singles)=subplot_size2;
    subplot_size(doubles)=subplot_size(doublesref);

    subplot_start=zeros(1,length(subplot_posi));               % normalized position of subplots
    subplot_start(singles)=subplot_start2;
    subplot_start(doubles)=subplot_start(doublesref);


    figure(num_fig); clf; hold on;
    %set(gcf,'Position',[5 39 1912 1072])
    set(gcf,'Name',['Mean values : ',title_string,'  ',comment_string])
    xlim([pmin pmax])
    ylim([0 regplotsize])
    xl=xlim; yl=ylim;
    for spc=1:length(subplot_start)
        if subplot_posi(spc)>0 && subplot_start(spc)~=regplotsize && subplot_start(spc)~=0
            line(xl,(yl(2)-subplot_start(spc))*ones(1,2),'Color','k','LineStyle','-','LineWidth',2)
        end
    end
    yt=[]; ytl=[];

    if subplot_posi(1)>0                                                                          % Stimulus
        % here you can plot the stimulus. The example below (if uncommented) shows a sine wave.
        % stim=0.5+sin((0:1/((itmax-itmin)/dt):1)*2*pi)/2;
        % xvals=itmin:dt:itmax;
        % if window_mode==1
        %     plot(xvals,yl(2)-subplot_start(1)+0.05+stim/1.1*subplot_size(1),['.-',colors(1)])
        % else
        %     lxvals=intersect(find(xvals>=wmin),find(xvals<=wmax));
        %     plot(xvals,yl(2)-subplot_start(1)+0.05+stim/1.1*subplot_size(1),['.-',colors(2)])
        %     plot(xvals(lxvals),yl(2)-subplot_start(1)+0.05+stim(lxvals)/1.1*subplot_size(1),['.-',colors(1)])
        % end
        % max_stim_val=1;
        % stim_lab=[-1 0 1];
        % stim_pos=[0 0.5 1];
        % yt=[yt yl(2)-subplot_start(1)+(0.05+stim_pos/max_stim_val)/1.1*subplot_size(1)];
        % ytl=[ytl stim_lab];
        % text(xl(1)-0.095*(xl(2)-xl(1)),yl(2)-subplot_start(1)+0.75/1.1*subplot_size(1),'Stimulus','Color','k','FontSize',font_size+1,'FontWeight','bold')
        % line(xl,yl(2)-subplot_start(1)+0.05/1.1*subplot_size(1)*ones(1,2),'Color','k','LineStyle',':')
        % line(xl,yl(2)-subplot_start(1)+1.05/1.1*subplot_size(1)*ones(1,2),'Color','k','LineStyle',':')
    end

    if subplot_posi(2)>0                                                                          % Spikes
        for trac=1:num_trains
            for sc=1:num_pspikes(trac)
                if pspikes(trac,sc)>=wmin && pspikes(trac,sc)<=wmax
                    line(pspikes(trac,sc)*ones(1,2),yl(2)-subplot_start(2)+(0.05+(num_trains-1-(trac-1)+[0.05 0.95])/num_trains)/1.1*subplot_size(2),'Color',colors(1))
                else
                    line(pspikes(trac,sc)*ones(1,2),yl(2)-subplot_start(2)+(0.05+(num_trains-1-(trac-1)+[0.05 0.95])/num_trains)/1.1*subplot_size(2),'Color',colors(2))
                end
            end
        end
        spikesvals=[0.5 1];
        spikeslab=spikesvals*num_trains;
        yt=[yt yl(2)-subplot_start(2)+(1.05-[spikeslab-0.5]/num_trains)/1.1*subplot_size(2)];
        ytl=[ytl spikeslab];
        text(xl(1)-0.12*(xl(2)-xl(1)),yl(2)-subplot_start(2)+0.75/1.1*subplot_size(2),'Spike','Color','k','FontSize',font_size+1,'FontWeight','bold')
        text(xl(1)-0.12*(xl(2)-xl(1)),yl(2)-subplot_start(2)+0.35/1.1*subplot_size(2),'trains','Color','k','FontSize',font_size+1,'FontWeight','bold')
        line(xl,yl(2)-subplot_start(2)+0.05/1.1*subplot_size(2)*ones(1,2),'Color','k','LineStyle',':')
        line(xl,yl(2)-subplot_start(2)+1.05/1.1*subplot_size(2)*ones(1,2),'Color','k','LineStyle',':')
    end
else
    subplot_size=zeros(1,length(subplot_posi));
    subplot_start=zeros(1,length(subplot_posi));
    maxintval=0; maxmeanintval=0; maxrateval=0; maxmeanrateval=0; maxbival=0; maxmultival=0;
    yt=[]; ytl=[]; xl=[]; yl=[];
end

paras=[plot_mode,dt,itmin,itmax,window_mode,wmin,wmax,spike_mao,time_mao,num_subplots,font_size,line_width,nma];
subplot_paras=[subplot_posi',subplot_size',subplot_start'];

if subplot_posi(3)>0                                                                          % ISI
    [dummy,dummy_ave,yt,ytl]=f_subplot(isi,ints,'ISI',subplot_paras(3,:),maxintval,yt,ytl,xl,yl,colors,1,paras);
end
if subplot_posi(4)>0                                                                          % ISI-Mean
    [mean_isi,mean_isi_ave,yt,ytl]=f_subplot(isi,mean_ints,'<ISI>',subplot_paras(4,:),maxmeanintval,yt,ytl,xl,yl,colors,1,paras);
    mean_values=[mean_values mean_isi];
    result_string=[result_string,' <ISI> = ',num2str(mean_isi,3)];
end


if subplot_posi(5)>0                                                                          % Rate
    [dummy,dummy_ave,yt,ytl]=f_subplot(isi,rates,'R',subplot_paras(5,:),maxrateval,yt,ytl,xl,yl,colors,1,paras);
end
if subplot_posi(6)>0                                                                          % Rate-Mean
    [mean_rate,mean_rate_ave,yt,ytl]=f_subplot(isi,mean_rates,'<R>',subplot_paras(6,:),maxmeanrateval,yt,ytl,xl,yl,colors,1,paras);
    mean_values=[mean_values mean_rate];
    result_string=[result_string,'  <R_I> = ',num2str(mean_rate,3)];
end


if subplot_posi(7)>0                                                                          % I^a
    [bi_isi_dist,bi_isi_dist_ave,yt,ytl]=f_subplot(isi,bi_ints,'I^a',subplot_paras(7,:),maxbival,yt,ytl,xl,yl,colors,1,paras);
    mean_values=[mean_values bi_isi_dist];
    time_profiles=[time_profiles bi_ints];
    result_string=[result_string,'  D_I^a = ',num2str(bi_isi_dist,3)];
end
if subplot_posi(8)>0                                                                          % I^m
    [multi_isi_dist,multi_isi_dist_ave,yt,ytl]=f_subplot(isi,multi_ints,'I^m',subplot_paras(8,:),maxmultival,yt,ytl,xl,yl,colors,1,paras);
    mean_values=[mean_values multi_isi_dist];
    result_string=[result_string,'  D_I^m = ',num2str(multi_isi_dist,3)];
end


if subplot_posi(9)>0                                                                          % S_C^a
    [bi_spike_dist_C,bi_spike_dist_C_ave,yt,ytl]=f_subplot(isi,spike_diffs_C,'S_C^a',subplot_paras(9,:),maxbival,yt,ytl,xl,yl,colors,1,paras);
    mean_values=[mean_values bi_spike_dist_C];
    time_profiles=[time_profiles; spike_diffs_C];
    result_string=[result_string,'  D_{SC}^a = ',num2str(bi_spike_dist_C,3)];
end
if subplot_posi(10)>0                                                                          % S_C^m
    [multi_spike_dist_C,multi_spike_dist_C_ave,yt,ytl]=f_subplot(isi,spike_stds_C,'S_C^m',subplot_paras(10,:),maxmultival,yt,ytl,xl,yl,colors,1,paras);
    mean_values=[mean_values multi_spike_dist_C];
    result_string=[result_string,'  D_{SC}^m = ',num2str(multi_spike_dist_C,3)];
end


if subplot_posi(11)>0                                                                          % S^a
    [bi_spike_dist,bi_spike_dist_ave,yt,ytl]=f_subplot(time,spike_diffs,'S^a',subplot_paras(11,:),maxbival,yt,ytl,xl,yl,colors,2,paras);
    mean_values=[mean_values bi_spike_dist];
    result_string=[result_string,'  D_{S}^a = ',num2str(bi_spike_dist,3)];
end
if subplot_posi(12)>0                                                                          % S^m
    [multi_spike_dist,multi_spike_dist_ave,yt,ytl]=f_subplot(time,spike_stds,'S^m',subplot_paras(12,:),maxmultival,yt,ytl,xl,yl,colors,2,paras);
    mean_values=[mean_values multi_spike_dist];
    result_string=[result_string,'  D_{S}^m = ',num2str(multi_spike_dist,3)];
end


if subplot_posi(13)>0                                                                          % PSTH
    [mean_psth,mean_psth_ave,yt,ytl]=f_subplot(bin_centers,psth,'PSTH',subplot_paras(13,:),maxpsth,yt,ytl,xl,yl,colors,3,paras);
    mean_values=[mean_values mean_psth];
    result_string=[result_string,'  PSTH = ',num2str(mean_psth,3)];
end
if subplot_posi(14)>0                                                                          % Gaussian Smoothed PSTH
    [mean_gs_psth,mean_gs_psth_ave,yt,ytl]=f_subplot(bin_centers,gs_psth,'GPSTH',subplot_paras(14,:),maxgspsth,yt,ytl,xl,yl,colors,3,paras);
    mean_values=[mean_values mean_gs_psth];
    result_string=[result_string,'  GPSTH = ',num2str(mean_gs_psth,3)];
end




if plot_mode>0

    if any(subplot_size(9:12))                           % mark perfectly synchronous events (SPIKE-Distances only)
        common=pspikes(1,1:num_pspikes(1));
        for trc=2:num_trains
            spik2=pspikes(trc,1:num_pspikes(trc));
            common=intersect(common,spik2);
        end
        num_commons=length(common);
        for spc=9:12
            if subplot_size(spc)>0
                for cc=1:num_commons
                    if common(cc)>=wmin && common(cc)<=wmax
                        line(common(cc)*ones(1,2),yl(2)-subplot_start(spc)+[0.05 1.05]/1.1*subplot_size(spc),'Color',colors(3),'LineStyle',':')
                    else
                        line(common(cc)*ones(1,2),yl(2)-subplot_start(spc)+[0.05 1.05]/1.1*subplot_size(spc),'Color',colors(4),'LineStyle',':')
                    end
                end
            end
        end
    end

    if window_mode>1
        line(wmin*ones(1,2),yl,'Color',colors(1),'LineStyle','-.')
        line(wmax*ones(1,2),yl,'Color',colors(1),'LineStyle','--')
    end
    line(itmin*ones(1,2),yl,'Color','k','LineStyle',':')
    line(itmax*ones(1,2),yl,'Color','k','LineStyle',':')
    line(tmin*ones(1,2),yl,'Color','k','LineStyle','-.')
    line(tmax*ones(1,2),yl,'Color','k','LineStyle','-.')
    xlabel(['Time  ',timeunit_string],'Color','k','FontSize',font_size+2,'FontWeight','bold')
    if publication==0
        title([title_string,' ',comment_string,result_string],'Color','k','FontSize',font_size+3,'FontWeight','bold')
    end
    [syt,syti]=sort(yt);
    sytl=ytl(syti);
    set(gca,'YTick',syt,'YTickLabel',sytl,'FontSize',font_size);
    set(gcf,'Color','w'); set(gca,'Color','w','Box','on');
    if print_mode                                                                    % Create postscript file
        if publication==1
            set(gcf,'PaperOrientation','Portrait'); set(gcf,'PaperType', 'A4');
            set(gcf,'PaperUnits','Normalized','PaperPosition',[0 0.45 1.0 0.5]);
        else
            set(gcf,'PaperOrientation','Landscape'); set(gcf,'PaperType', 'A4');
            set(gcf,'PaperUnits','Normalized','PaperPosition',[0 0 1.0 1.0]);
        end
        psname=[imagespath,filename,comment_string,'.ps'];
        print(gcf,'-dpsc',psname);
    end
end

%
% *************************************************************************
%

function [ave,ave2,yt,ytl]=f_subplot(fx,fy,fheadline,fsp_paras,fmaxval,fyt,fytl,fxl,fyl,fcols,fdata_type,fparas)

plot_mode=fparas(1);
dt=fparas(2);
itmin=fparas(3);
itmax=fparas(4);
window_mode=fparas(5);
wmin=fparas(6);
wmax=fparas(7);
mao=fparas(8);
time_mao=fparas(9);
num_subplots=fparas(10);
font_size=fparas(11);
line_width=fparas(12);
nma=fparas(13:end);

fsp_posi=fsp_paras(1);
fsp_size=fsp_paras(2);
fsp_start=fsp_paras(3);
if window_mode==1
    dfcols=[fcols([3 3])];
    dmfcols=[fcols([1 1])];
else
    dfcols=[fcols([3 4])];
    dmfcols=[fcols([1 2])];
end

for nmac=1:length(nma)

    if nma(nmac)==2 && size(fy,1)==1                                                       % moving average
        if round(sum(fx)/dt)*dt==round((itmax-itmin)/dt)*dt
            for trac=1:size(fy,1)
                fy(trac,:)=f_moving_average_weighted(fy(trac,:),fx,mao);
            end
        else
            for trac=1:size(fy,1)
                fy(trac,:)=f_moving_average(fy(trac,:),time_mao);
            end
        end
        fheadline=[fheadline,'*'];
        line_width=line_width+1;
        if window_mode==1
            dfcols=[fcols([5 5])];
        else
            dfcols=[fcols([5 6])];
        end
    end

    if fdata_type==1                                            % using vectors of length "num_isi"
        ctfx=cumsum([0 fx]);
        pfx=itmin+unique([ctfx(1:end) ctfx(2:end-1)-dt/100]);

        if size(fy,1)==1
            pfy=reshape([fy; fy],1,length(fy)*2);
            if window_mode==1
                ave=sum(fy.*repmat(fx,size(fy,1),1),2)/sum(fx);
                ave2=mean(fy);
            end
        end
    else                                                        % using vectors of length "len"
        pfx=fx;
        pfy=fy;
        if size(fy,1)==1 && window_mode==1
            ave=mean(fy,2);
            ave2=ave;
        end
    end

    if size(fy,1)==1                                                            % one single profile
        if window_mode>1
            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 fdata_type<3
                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
            end

            win_pfx_isi=diff(win_pfx);
            win_pfy_isi=win_pfy(1:end-1);

            if fdata_type==1                                        % using vectors of length "num_isi"
                ave=sum(win_pfy_isi.*repmat(win_pfx_isi,size(win_pfy_isi,1),1),2)/sum(win_pfx_isi);
                ave2=mean(win_pfy_isi);
            else                                                    % using vectors of length "len"
                ave=mean(win_pfy);
                ave2=ave;
            end
        end
    else
        ave=0;
        ave2=ave;
    end


    if plot_mode>0 && fsp_posi>0                                                              % plotting
        num_vals=2;
        if num_subplots<6
            intvals=[0.5 1];
        elseif num_subplots<10
            intvals=[0.4 0.8];
        else
            intvals=[0.3 0.6];
        end
        intlab=unique([0 f_lab(intvals*fmaxval,num_vals,1,1)]);
        yt=[fyt fyl(2)-fsp_start+(0.05+intlab/fmaxval)/1.1*fsp_size];
        ytl=[fytl intlab];
        text(fxl(1)-0.11*(fxl(2)-fxl(1)),fyl(2)-fsp_start+0.65/1.1*fsp_size,fheadline,'Color','k','FontSize',font_size+1,'FontWeight','bold')
        line(fxl,fyl(2)-fsp_start+0.05/1.1*fsp_size*ones(1,2),'Color','k','LineStyle',':')
        line(fxl,fyl(2)-fsp_start+1.05/1.1*fsp_size*ones(1,2),'Color','k','LineStyle',':')

        if size(fy,1)==1
            plot(pfx,fyl(2)-fsp_start+(0.05+pfy/fmaxval)/1.1*fsp_size,['-',dfcols(2)],'LineWidth',line_width)
            if window_mode>1
                plot(win_pfx,fyl(2)-fsp_start+(0.05+win_pfy/fmaxval)/1.1*fsp_size,['-',dfcols(1)],'LineWidth',line_width)
                %line([wmin wmax],fyl(2)-fsp_start+(0.05+ave/fmaxval)/1.1*fsp_size*ones(1,2),'LineStyle','--','Color',dfcols(1),'LineWidth',line_width)
            end
%             if nmac==length(nma)
%                text(fxl(2)-0.1*(fxl(2)-fxl(1)),fyl(2)-fsp_start+0.5/1.1*fsp_size,num2str(ave,3),'Color',dfcols(1),'FontWeight','bold')
%             end

        elseif nmac<2                                                                   % ISIs and rates for single trials

            if window_mode>1
                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);
                if fdata_type<3
                    if wmin<pfx(first_winspike_ind)
                        win_pfx=[wmin win_pfx];
                    end
                    if wmax>pfx(last_winspike_ind)
                        win_pfx=[win_pfx wmax];
                    end
                end
            end
            for trac=1:size(fy,1)
                fyy=fy(trac,:);
                pfyy=reshape([fyy; fyy],1,length(fyy)*2);
                plot(pfx,fyl(2)-fsp_start+(0.05+pfyy/fmaxval)/1.1*fsp_size,['-',dmfcols(2)],'LineWidth',line_width)
                if window_mode>1
                    win_pfyy=pfyy(first_winspike_ind:last_winspike_ind);
                    if wmin<pfx(first_winspike_ind)
                        win_pfyy=[pfyy(first_winspike_ind-1) win_pfyy];
                    end
                    if wmax>pfx(last_winspike_ind)
                        win_pfyy=[win_pfyy pfyy(last_winspike_ind+1)];
                    end
                    plot(win_pfx,fyl(2)-fsp_start+(0.05+win_pfyy/fmaxval)/1.1*fsp_size,['-',dmfcols(1)],'LineWidth',line_width)
                end
            end
        end
    else
        yt=[];
        ytl=[];
    end
end

%
% *************************************************************************
%

function z=f_moving_average(x,o)    % moving average past+future, inwards averaging at the edges (asymmetric)

if size(x,length(size(x)))==1
    x=x';
end
n=length(x);
if n<2*o+1
    o=fix((n-1)/2);
end
if o~=0
    z=zeros(1,n);
    y=zeros(2*o+1,n-2*o);
    for k=1:o
        z(k)=sum(x(k-(k-1):k+o))/((k-1)+1+o);
    end
    for k=1:2*o+1
        y(k,1:n-2*o)=x(k-1+(1:n-2*o));
    end
    z(o+1:n-o)=mean(y);
    for k=n-o+1:n
        z(k)=sum(x(k-o:k+(n-k)))/(o+1+(n-k));
    end
else
    z=x;
end

return

%
% *************************************************************************
%

function z=f_moving_average_weighted(y,x,o)    % past+future, inwards averaging at the edges (asymmetric), time weighted

if size(y,length(size(y)))==1
    y=y';
end
if size(x,length(size(x)))==1
    x=x';
end

n=length(y);
if n<2*o+1
    o=fix((n-1)/2);
end
if o~=0
    z=zeros(1,n);
    x2=zeros(2*o+1,n-2*o);
    y2=zeros(2*o+1,n-2*o);
    for k=1:o
        z(k)=sum(y(k-(k-1):k+o).*x(k-(k-1):k+o))/sum(x(k-(k-1):k+o));
    end
    for k=1:2*o+1
        y2(k,1:n-2*o)=y(k-1+(1:n-2*o));
        x2(k,1:n-2*o)=x(k-1+(1:n-2*o));
    end
    z(o+1:n-o)=sum(y2.*x2)./sum(x2,1);
    for k=n-o+1:n
        z(k)=sum(y(k-o:k+(n-k)).*x(k-o:k+(n-k)))/sum(x(k-o:k+(n-k)));
    end
else
    z=y;
end

return

%
% *************************************************************************
%

function lab=f_lab(vect,num_values,origin,cut)

if origin==1
    vect2=[0 vect];
else
    vect2=vect;
end
range=max(vect2)-min(vect2);
interval=range/num_values;
min_unit=fix(interval*10^(ceil(-log10(interval))))/10^(ceil(-log10(interval)));

xsep=round((range+min_unit)/num_values/min_unit)*min_unit;
lab=round(min(vect2)/min_unit)*min_unit+(xsep:xsep:(range+min_unit));

if cut==1
    lab=lab(lab>=min(vect2) & lab<=max(vect2));
end

return

%
% % *************************************************************************
%

function dt=f_get_dt(spikes)

vspikes=reshape(spikes,1,numel(spikes));

logs=-10:10;
c=0;
for lc=logs
    c=c+1;
    decades(c)=10^-lc;
    b=vspikes*decades(c);
    yes(c)=all(b==round(b));
    if yes(c)==0
        break;
    end
end

if c==1 || all(yes)
    dt=10^-8;
else
    dt=10^logs(c-1);
end

return
%
% % *************************************************************************
%

function output = f_compute_gauss_smooth(input, window)
% Copyright: Jude Mitchell, http://www.snl.salk.edu/~jude/sfn2008/index.html
% Smoothing function:
% output = smooth(input, window)
% "Window" is the total kernel width.
% Input array must be one-dimensional.

input_dims = ndims(input);
input_size = size(input);
if input_dims > 2 || min(input_size) > 1,
    disp('Input array is too large.');
    return
end

if input_size(2) > input_size(1),
    input = input';
    toggle_dims = 1;
else
    toggle_dims = 0;
end

if window/2 ~= round(window/2),
    window = window + 1;
end
halfwin = window/2;

input_length = length(input);
%********* gauss window +/- 1 sigma
x = -halfwin:1:halfwin;
kernel = exp(-x.^2/(window/2)^2);
kernel = kernel/sum(kernel);

padded(halfwin+1:input_length+halfwin) = input;
padded(1:halfwin) = ones(halfwin, 1)*input(1);
padded(length(padded)+1:length(padded)+halfwin) = ones(halfwin, 1)*input(input_length);

output = conv(padded, kernel);
output = output(window:input_length+window-1);

if toggle_dims == 1,
    output = output';
end

return
