function explanations=cluster_asynch_based_on_phasic(logfilename,max_number_clusters,plotting,fast_event_thresh_phasic,fast_event_thresh_asynch,jit_range)
%description: Extracts phasic and asynchronous events, removes any underlying phasic current from asynch events, and sorts asynchronous events, assumed to represent a single
% vesicle, from smallest to largest.  For each num_clusters up to max_number_clusters, all possible clusterings of the asynchronous events are tried by
% taking the mean of each cluster to represent that cluster. The clusterings are evaluated by assuming each phasic event represents some combination of
% clusters (release sites), so each possible combination of clusters is compared to each phasic event, over all possible jitters between clusters.  The best
% clustering is taken to be the one producing the lowest total error when used (in the best way found) to fit all phasic events.  When n+1 clusters is 'not much better' 
% at fitting than n clusters, n is taken to be the true number of clusters.
%called by: [command line]
%calls: iterative_curve_fit_singledecay, convert_IPIs_into_timewindows, mini_waveform_extractor, most_basic_colorcoded_plot,
% list_all_nonredundant_permutations_upto_lengthn, find_all_n_possible_contiguous_clusterings, compare_asynch_waveform_clusters_to_phasic_waveform
%notes: it can be difficult to extract waveforms accurately if IPI is very short; a 50 ms minimum is preferred
% this approach potentially involves a very large number of comparisons, so parallelization is desirable;  parallelize prior to any recursive steps
% to test for addition of new release site following induction, run the clustering without additional instructions on both pre- and postinduction, 
% and see if postinduction requires one more cluster than preinduction to be well fit, and how postinduction clusters compare to preinduction clusters
%code confirm date: [begin line 162, recheck variable assignments as final step in code confirmation]

starttime=clock %for display
clearvars -except logfilename max_number_clusters plotting fast_event_thresh_phasic fast_event_thresh_asynch jit_range, close all

%slow_event_thresh=.55;  %slow event threshold and jitter are not currently used, but the slow criterion trace is part of the criterion_vals file -- overly good fit to the slow template (or fit to slow but no fit to fast) is to be used to reject spurious events
    %slow_event_jitter=10;
phasic_window_ends=9.9;  %in ms, the time from the stimulus artifact onset within which the event must be detected (as having *started*) in order to be classed as phasic
asynch_window_begins=10; %in ms
fast_event_jitter=10;    %in samples;  used by mini_waveform_extractor, NOT used in fitting asynch waveforms to phasic waveform
    event_window=-6:57;      %in samples, the membrane current information that will be returned for each event, with time 0 being the time of the peak of the detection criterion
    baseline_offset_window1=1:5;    %the indices of the samples that will be averaged and subtracted from the event for baseline removal
    baseline_offset_window2=47:57;     
    fit_slope_thresh=.1; %used in determining whether the slope of a linear fit to the initial values suggests that the event is riding on some underlying current waveform
latency_search_window=[1 -1*event_window(1)+4]; %in samples, used by iterative_curve_fit_singledecay, NOT used in fitting asynch waveforms to phasic waveform
    return_pareddown_template=1; %used by iterative_curve_fit_singledecay
asynch_event_rise_time_criterion=[1 9.5];   %use/don't use, amount in samples
asynchwindow_afterlastpulse=225; %in ms
%this bunch of variables could vary by recording
stim_pulse_duration=15;  %in samples
samples_per_ms=10;
sample_containing_first_stim_onset=10000; % the criterion_vals file stores each trial with the first AP beginning at the 10000th sample

%----------------

if plotting
    plot_events=1;plot_fits=1;plot_datafitresiduals=1;plot_sloperemovedextractions=1;
else
    plot_events=0;plot_fits=0;plot_datafitresiduals=0;plot_sloperemovedextractions=0;
end

asyn=load(strcat('/Users/tathagatagarbha/dwn_data_analysis_system/extracted_values_dwn/',logfilename,'__fastandslow_criterion_vals'));
pscs_for_IPIs=load(strcat('/Users/tathagatagarbha/dwn_data_analysis_system/extracted_values_dwn/',logfilename,'__pre-postPSC'));
%noi=load(strcat('/Users/tathagatagarbha/dwn_data_analysis_system/extracted_values_dwn/',logfilename,'__noise_traces'));    %noise trace, for random
% noise to add to simulated waveforms of multiple jittered asynch events added together -- not currently implemented
 
im_windows=convert_IPIs_into_timewindows(pscs_for_IPIs(1:2,:),sample_containing_first_stim_onset,samples_per_ms,asynchwindow_afterlastpulse); 

fast_phasic_events_alltrials=[];fast_asynch_events_alltrials=[];
for trial=1:size(asyn,2)/3  %asyn has 3 columns per trial: raw recording, fast criterion, and slow criterion           
    column_adder=(trial-1)*3;
    
    membrane_current_trace=asyn(:,column_adder+1);
    fast_detection_trace=asyn(:,column_adder+2);
    for AP=1:4  %expects a train with 4 APs
        phasewindow=(im_windows(AP,trial)+stim_pulse_duration):(im_windows(AP,trial)+(phasic_window_ends*samples_per_ms));
        fast_phasic_events=mini_waveform_extractor(membrane_current_trace([(phasewindow(1)+event_window(1)):(phasewindow(end)+event_window(end))]),...
            fast_detection_trace(phasewindow),fast_event_thresh_phasic,fast_event_jitter,event_window,baseline_offset_window1,baseline_offset_window2,...
            fit_slope_thresh,plot_sloperemovedextractions);
        fast_phasic_events_alltrials=[fast_phasic_events_alltrials;fast_phasic_events];
        
        asynwindow=(im_windows(AP,trial)+(asynch_window_begins*samples_per_ms)):(im_windows(AP+1,trial)-event_window(end));
        fast_asynch_events=mini_waveform_extractor(membrane_current_trace([(asynwindow(1)+event_window(1)):(asynwindow(end)+event_window(end))]),...
            fast_detection_trace(asynwindow),fast_event_thresh_asynch,fast_event_jitter,event_window,baseline_offset_window1,baseline_offset_window2,...
            fit_slope_thresh,plot_sloperemovedextractions);
        fast_asynch_events_alltrials=[fast_asynch_events_alltrials;fast_asynch_events];
    end
end

if plot_events,clr=colormap(jet(size(fast_phasic_events_alltrials,1)));most_basic_colorcoded_plot(clr,fast_phasic_events_alltrials,'phasic events'),end
if plot_events,clr=colormap(jet(size(fast_asynch_events_alltrials,1)));most_basic_colorcoded_plot(clr,fast_asynch_events_alltrials,'asynch events'),end

phassize=size(fast_phasic_events_alltrials,1)   %for display
asynsize=size(fast_asynch_events_alltrials,1)

fast_phasic_waveform_fits=[];
for present_phase_template=1:size(fast_phasic_events_alltrials,1)
    phase_template=iterative_curve_fit_singledecay(fast_phasic_events_alltrials(present_phase_template,:),latency_search_window,return_pareddown_template,plot_datafitresiduals);
    fast_phasic_waveform_fits=[fast_phasic_waveform_fits;phase_template];
end
if plot_fits,clr=colormap(jet(size(fast_phasic_waveform_fits,1)));most_basic_colorcoded_plot(clr,fast_phasic_waveform_fits(:,1:end-2),'phasic fits before sort'),end
phasic_peak_pscs=fast_phasic_waveform_fits(:,end);
[y_phasic,i_phasic]=sort(phasic_peak_pscs);
fast_phasic_waveform_fits=fast_phasic_waveform_fits(i_phasic,1:end-2); fast_phasic_events_alltrials=fast_phasic_events_alltrials(i_phasic,:);
if plot_fits,most_basic_colorcoded_plot(clr,fast_phasic_waveform_fits,'phasic fits after sort'),end
    
fast_asynch_waveform_fits=[];
for present_asynch_template=1:size(fast_asynch_events_alltrials,1)
    asynch_template=iterative_curve_fit_singledecay(fast_asynch_events_alltrials(present_asynch_template,:),latency_search_window,return_pareddown_template,plot_datafitresiduals);
    fast_asynch_waveform_fits=[fast_asynch_waveform_fits;asynch_template];
end
if asynch_event_rise_time_criterion(1),asynch_rise_times=fast_asynch_waveform_fits(:,end-1);fast_asynch_waveform_fits(find(asynch_rise_times>asynch_event_rise_time_criterion(2)),:)=[];
    fast_asynch_events_alltrials(find(asynch_rise_times>asynch_event_rise_time_criterion(2)),:)=[];end
if plot_fits,clr=colormap(jet(size(fast_asynch_waveform_fits,1)));most_basic_colorcoded_plot(clr,fast_asynch_waveform_fits(:,1:end-2),'asynch fits before sort'),end
asynch_peak_pscs=fast_asynch_waveform_fits(:,end);
[y_asynch,i_asynch]=sort(asynch_peak_pscs);
fast_asynch_waveform_fits=fast_asynch_waveform_fits(i_asynch,1:end-2); fast_asynch_events_alltrials=fast_asynch_events_alltrials(i_asynch,:);
if plot_fits,most_basic_colorcoded_plot(clr,fast_asynch_waveform_fits,'asynch fits after sort'),end

asynpostcullsize=size(fast_asynch_events_alltrials,1) %for display
if plot_events, most_basic_colorcoded_plot(clr,fast_asynch_events_alltrials,'asynch events after cull'),end

list_of_asynch_events=1:size(fast_asynch_waveform_fits,1)
if length(list_of_asynch_events)<2,error('too few asynchronous events'),end

%setting start points and time ranges for waveform comparison
event_prefix=-event_window(1)+1;   %when raw events were extracted, this many samples were returned prior to/simultaneous with what is taken as the true event onset
fast_asynch_latency_prefix=zeros(size(fast_asynch_waveform_fits,1),jit_range(end-2));   %one zero is already prefixed to the fit waveform by iterative_curve_fit_singledecay;  with the length of the zero prefix being one less than jit_range(end), the fitted waveform can begin on the first step straight no prefix
fast_asynch_waveform_fits=[fast_asynch_latency_prefix fast_asynch_waveform_fits]; %prepending zeros to postpone (jitter) the events
phasic_startpoint_forclustering=floor(.5*event_prefix);  %mostly we are looking for delayed exocytosis, so we reduce the duration of pre-onset current we'll consider in the comparison with asynch waveforms to focus more on the later portions
fast_phasic_events_alltrials=fast_phasic_events_alltrials(:,phasic_startpoint_forclustering:end);
asynchfittingsegment=[0:size(fast_phasic_events_alltrials,2)-1];

for num_clusters=1:max_number_clusters
    %find all the ways of parcelling the asynch event waveforms into num_clusters clusters. clusters contain waveforms contiguous in terms of peak
    %size and each cluster will be represented in fitting by the mean of all the waveforms it contains.
	if num_clusters==1
        complete_list_of_clusterings={{list_of_asynch_events}};   %this has to be a cell array containing a cell array, as that is what is expected below
    else
        clear global; global complete_list_of_clusterings_global,complete_list_of_clusterings_global={};
		mobile_dividers=1:num_clusters-1;
		find_all_n_possible_contiguous_clusterings(list_of_asynch_events,[],mobile_dividers);
        complete_list_of_clusterings=complete_list_of_clusterings_global;clear global
    end
        
    %find all the ways that this many clusters could combine to form phasic events (i.e., if 2 clusters of asynch events [2 release sites] are believed to exist, a given phasic event can be from site 1, site 2, or sites 1 and 2)
    clear global, global all_permutations_of_sites_given_numclusters_sites_global,all_permutations_of_sites_given_numclusters_sites_global={};
    for length_sought=1:num_clusters,list_all_nonredundant_permutations_upto_lengthn([1:num_clusters],[],length_sought),end     %it would certainly be preferable to make list_all_nonredundant_permutations_upto_lengthn and find_all_n_possible_contiguous_clusterings use no global variables
    all_permutations_of_sites_given_numclusters_sites=all_permutations_of_sites_given_numclusters_sites_global; clear global
     
    %then consider, for each clustering, how well do all permutations (with all jitters) of the means of those clusters, summed, fit each phasic event, 
    %storing the permutation that gives the best fit and ultimately storing the clustering whos best-fit permutation for each phasic yielded the lowest total error
    %over all phasics -- this is taken to be the best clustering for num_clusters clusters.
    sumsmallesterrorallphasics_foreachpossibleclustering=zeros(1,size(complete_list_of_clusterings,2));
    parfor cluster_counter=1:length(complete_list_of_clusterings)
       current_clustering=complete_list_of_clusterings{cluster_counter};
       asynch_mean_waveform_each_cluster=[];
       for making_means=1:length(current_clustering)
           current_cluster_members=current_clustering{making_means};
           waveforms_this_cluster=fast_asynch_waveform_fits(current_cluster_members,:);
           asynch_mean_waveform_each_cluster=[asynch_mean_waveform_each_cluster;mean(waveforms_this_cluster,1)];
       end
        
       smallest_min_error_all_phasics_thisclustering=inf;     greatest_min_error_all_phasics_thisclustering=-inf;
       for phasic_waveform_under_consideration=1:size(fast_phasic_events_alltrials,1)
           squared_peaksize_present_phasic=(min(fast_phasic_waveform_fits(phasic_waveform_under_consideration,:)))^2;  %used in scaling the fit error
           smallesterrorforthisphasicwaveform_atthisclustering=inf;bestfitsummedasynch_forthisphasic=nan;
           for active_cluster_permutation_counter=1:length(all_permutations_of_sites_given_numclusters_sites)
               smallesterrthisphasicwvform_wthispermutofclustersactive=inf;
               [smallesterrthisphasicwvform_wthispermutofclustersactive,smallesterrorforthisphasicwaveform_atthisclustering,bestfitsummedasynch_forthisphasic]=...
                   compare_asynch_waveform_clusters_to_phasic_waveform(asynch_mean_waveform_each_cluster(all_permutations_of_sites_given_numclusters_sites{active_cluster_permutation_counter},:),[],...
                    fast_phasic_events_alltrials(phasic_waveform_under_consideration,:),squared_peaksize_present_phasic,asynchfittingsegment,jit_range,...
                    smallesterrthisphasicwvform_wthispermutofclustersactive,smallesterrorforthisphasicwaveform_atthisclustering,bestfitsummedasynch_forthisphasic);
               if smallesterrthisphasicwvform_wthispermutofclustersactive<=smallesterrorforthisphasicwaveform_atthisclustering
                    bestactivepermutationofthisclusteringforthisphasic=all_permutations_of_sites_given_numclusters_sites{active_cluster_permutation_counter};
               end
           end
           if smallesterrorforthisphasicwaveform_atthisclustering>greatest_min_error_all_phasics_thisclustering
               greatest_min_error_all_phasics_thisclustering=smallesterrorforthisphasicwaveform_atthisclustering;
               whichphasichadworstbestfitatthisclustering(cluster_counter)=phasic_waveform_under_consideration;
               activepermutationthisclusteringwithworstbestfit{cluster_counter}=bestactivepermutationofthisclusteringforthisphasic;
               summedasynchwithworstbestfitthisclustering{cluster_counter}=bestfitsummedasynch_forthisphasic;
               worstbestfit_allclusterings(cluster_counter)=smallesterrorforthisphasicwaveform_atthisclustering;
           end
           if smallesterrorforthisphasicwaveform_atthisclustering<smallest_min_error_all_phasics_thisclustering
               smallest_min_error_all_phasics_thisclustering=smallesterrorforthisphasicwaveform_atthisclustering;
               whichphasichadbestbestfitatthisclustering(cluster_counter)=phasic_waveform_under_consideration;
               activepermutationthisclusteringwithbestbestfit{cluster_counter}=bestactivepermutationofthisclusteringforthisphasic;
               summedasynchwithbestbestfitthisclustering{cluster_counter}=bestfitsummedasynch_forthisphasic;
           end
           sumsmallesterrorallphasics_foreachpossibleclustering(cluster_counter)=sumsmallesterrorallphasics_foreachpossibleclustering(cluster_counter)+...
                smallesterrorforthisphasicwaveform_atthisclustering;    
        end        
    end
    [sumacrossphasics_smallesterrors_withthismanyclusters(num_clusters) bestind]=min(sumsmallesterrorallphasics_foreachpossibleclustering) %#ok<*AGROW>
    best_clustering_for_this_many_clusters{num_clusters}=complete_list_of_clusterings{bestind};
    sumsmallesterr_allphasics_eachpossibleclustering_allnumclstrs{num_clusters}=sumsmallesterrorallphasics_foreachpossibleclustering;
    worstbestfit_allclusterings_allnumclstrs{num_clusters}=worstbestfit_allclusterings;
    phasicwithbestbestfitatbestclustering(num_clusters)=whichphasichadbestbestfitatthisclustering(bestind);
    bestbestfit_summedasynch_atbestclustering{num_clusters}=summedasynchwithbestbestfitthisclustering{bestind};
    activepermutationthisclusteringwithbestbestfit_allnumclusters{num_clusters}=activepermutationthisclusteringwithbestbestfit{bestind};
    phasicwithworstbestfitatbestclustering(num_clusters)=whichphasichadworstbestfitatthisclustering(bestind);
    worstbestfit_summedasynch_atbestclustering{num_clusters}=summedasynchwithworstbestfitthisclustering{bestind};
    activepermutationthisclusteringwithworstbestfit_allnumclusters{num_clusters}=activepermutationthisclusteringwithworstbestfit{bestind};    
end
explanations={
    sumacrossphasics_smallesterrors_withthismanyclusters 
    best_clustering_for_this_many_clusters
    sumsmallesterr_allphasics_eachpossibleclustering_allnumclstrs
    worstbestfit_allclusterings_allnumclstrs
    fast_phasic_waveform_fits
    fast_phasic_events_alltrials
    fast_asynch_waveform_fits 
    fast_asynch_events_alltrials
    phasicwithbestbestfitatbestclustering
    bestbestfit_summedasynch_atbestclustering
    activepermutationthisclusteringwithbestbestfit_allnumclusters
    phasicwithworstbestfitatbestclustering    
    worstbestfit_summedasynch_atbestclustering
    activepermutationthisclusteringwithworstbestfit_allnumclusters
    };

endtime=clock