function [spd1crits,actual_SC,actual_distances]=IQI_estimator(spd1,spd2,IQI_index_row,crit_indices,actual_crits,truncate_SC_at_final_step,plotting_spd1_crits)
%description:  use simulated_psc_in_noise ('SPIN') with interquantal intervals (IQIs) of 0-N samples between events of fixed size (e.g., first/second event 12/6 pA, or 8/4 pA;  a separate SPIN run for each size) 
    %where N is the length of IQI beyond which the second event can be detected by conventional BC detection;  i.e., the upper limit of our IQI exploration
    %place, say, about 10k two-waveform events at each such IQI, using recorded noise from the experiment being analyzed -- these first SPIN runs are called the 'first-order simulation'
    %the results of each of these SPIN runs is stored and passed in as spd1, a cell array containing results of as many SPIN runs as seem like a good idea
    %for each IQI, average the detection criteria, averaging first within and then, after verifying they are similar (spd1crits is returned for this purpose), across SPIN runs
    %this creates N+1 criterion waveforms, corresponding to IQIs 0-N;  these are called the standard criteria (SC)
    
    %then run SPIN again (the 'second-order simulation'), ~10k events for each IQI 0-N samples, with event amplitude, say, randomly varying from, say, (first/second) 0-12/0-6 pA, or 3-12/3-6 pA (i.e., without or with a minimum size)
    %the results are passed in as spd2, which includes the criterion for each event
    %do B-C detection (at single time point) on that criterion using each criterion waveform in the SC as template -- for each event, save the resulting N+1 values
    %for each IQI, find the mean and SD of these N+1 criterion-detect-criterion ('CDC') values
    
    %then take each recorded waveform from a real experiment, fit with the basic template, and pass in the resulting criterion as actual_crits
    %do criterion-detect-criterion as above -- because we cannot say with certainty when the recorded event begins, CDC this time, unlike above, is done in the usual sliding-across template match way, and we take the maximum 
    %this yields, for each criterion waveform in actual_crits, a vector with one value for every waveform in the SC ('distance' of that actual crit from each waveform of the SC)
    %repmat this vector to have as many columns as there are IQIs
    %take the resultant array and subtract CDC_means, take absolute value ('distance' of that actual crit from each IQI's CDC profile)
    %add CDC_SDs to the result, this penalizes for variability in fits of second order simulations to the SC 
    %perform sum;  the column with the smallest sum is that closest to the recorded event, thus, the best estimate of the IQI of the recorded event;  the results of sum on all the other columns give their 'distance' from the recorded event
    %the distance of all the possible IQIs from our recorded event can be converted into a likelihood that each of them was the actual IQI underlying the observed event
%called by: command line
%calls: multipeak_criterion_analyzer, CDC_creator, template_based_mini_detector
%notes: the amplitude ranges used in spd1 and spd2 should presumably include those seen in the recorded events, but the autoscaling of B-C detection may obviate this
    %in tests the penalty of adding CDC_SDs serves more to offset the curve than to change its shape; both consistently good and consistently bad fits will have low CDC_SD
%code confirmed date: dec 28, 2012
    
close all
tic

%find SC by analyzing all instances of spd contained in spd1, take the within-instance mean criterion by IQI, then average across the instances of spd
spd1crits=[];
for i=1:size(spd1,2)
    [meancriteria,~]=multipeak_criterion_analyzer(spd1{i},IQI_index_row,crit_indices,plotting_spd1_crits);  %variability at this stage is ignored
    spd1crits=cat(3,spd1crits,meancriteria);
end
SC=mean(spd1crits,3);

%get CDC traces by comparing each criterion in spd2 with every column in SC  -- this is the rate-limiting step, hence parfor
num_possible_IQIs=size(SC,2)
SC_template_length=size(SC,1);
second_SC=[spd2(IQI_index_row,:);nan(num_possible_IQIs,size(spd2,2))];
second_SC_rowind=2:(num_possible_IQIs+1);
parfor j=1:size(spd2,2)
    CDCs_this_event=CDC_creator(num_possible_IQIs,spd2(crit_indices,j),SC,SC_template_length);
    second_SC(second_SC_rowind,j)=CDCs_this_event;
end

%find mean CDC values for each IQI
CDC_means=nan(num_possible_IQIs+1,num_possible_IQIs);
CDC_SDs=nan(num_possible_IQIs+1,num_possible_IQIs);
IQIs=second_SC(1,:);
for m=1:num_possible_IQIs
    live=find(IQIs==(m-1));
    CDC_means(:,m)=mean(second_SC(:,live)')'; 
    CDC_SDs(:,m)=std(second_SC(:,live)')';      %CDC_SDs is not normalized by CDC_means
end

if truncate_SC_at_final_step
    SC=SC(1:truncate_SC_at_final_step,:);  %truncating SC at this stage (to appx length(template)+max(IQI)) may improve performance slightly
end

actual_SC=nan(num_possible_IQIs,size(actual_crits,2));
actual_distances=nan(num_possible_IQIs,size(actual_crits,2));
for n=1:size(actual_crits,2)
    for o=1:num_possible_IQIs
        detec=template_based_mini_detector(actual_crits(:,n),SC(:,o));
        actual_SC(o,n)=max(detec);
    end
    a=repmat(actual_SC(:,n),1,num_possible_IQIs);
    b=a-CDC_means(2:end,:);
    c=abs(b)+CDC_SDs(2:end,:);  %small distance means a good fit; large CDC_SD is penalized
    d=sum(c);
    actual_distances(:,n)=d';
end

toc    