function sim_psc_data=simulated_psc_in_noise_vardetec(number_reps,recorded_noise,use_noise,splisiterequirements,find_unique_segments,breakpoints_mustbe_consecutive,max_segments,minimal_noise_length,...
    template_type,template_seeds,taper_template,trunc_templ_detec,scaler,two_waveforms,second_waveform_subset,variable_waveform_size,minimum_absolute_peak_amplitude,...
    amp_window,detec_alg,plot_detection,omniscient_start,precede_for_return)
%description: generates psc or gaussian waveforms, embeds them in recombined noise, uses template-match or threshold-based detection to find them, returns true onset time and size, measured onset time and size, ...
    %event confidence (criterion) or difference (if K-T detection),rise slope, plateau length, inter-event interval, and criterion trace in the area of the assigned start time (B-C) or event info (K-T)
    %used to assess event detection, accuracy of onset time and event size measurments of PSCs or symmetric events in noise, and rise time and plateau length with temporally offset PSCs with (or without) noise,...
    %as well as many questions related to the presence, relative timing, and detectability of two events
%called by: command line
%calls: template_generator,noise_segmenter,noise_recombinase,template_based_mini_detector,find_multiple_peaks,kudoh_taguchi_results_sorter,find_riseplateau
%notes:
    %use_noise indicates whether noise should be added or not -- sometimes with double-waveform simulations it may be useful to measure rise time and plateau with no noise at all
    %splisiterequirements (often [3 3 3 3]): required_consec_decreases_belowzero, required_consec_increases_belowzero, required_consec_increases_abovezero, required_consec_decreases_abovezero
    %find_unique_segments, yes or no (find_unique_segments is more restrictive than the name implies, taking the shortest unique segments)
    %breakpoints_mustbe_consecutive, yes or no  (the interaction of find_unique_segments and breakpoints_mustbe_consecutive is slightly complex)
    %set max_segments to zero to have no maximum; if other value, only that many segments will be chosen from the noise for recombination each iteration...
        %this has the advantages of increasing variability of recombined noise from one trial to the next, and possibly of decreasing run times; minimal_noise_length confirms the resulting segment is long enough
    %minimal_noise_length, in samples, generally around 1000; the minimal noise allowed after cold_start and dead_end (below) are removed; 
        %shorter can increase variability in the sense of reducing the number of segments that need to be used on any trial
    %template_type, 1 for psc-like (standard), 2 for gaussian
    %template_seeds has 3 values: for PSC template they are amplitude_scaler,tau_rise,tau_decay (often [23.984085025960670   3.861189308388354  27.526665160068116]); 
        %for gaussian they are peak_height,center_time,bell_width (often [25 20 5])
    %taper_template has two elements, starttaperlength (skip if 0) and endtaperlength (skip if 0); only called if 1 || 2
        %taper_template(1) is negative is a special situation, a kind of non-taper, where abs(taper_template(1)) zeroes are prepended to the front of the template; if used, there is no additional tapering        
    %trunc_templ_detec is set at zero to use the full template for detection, otherwise it, as a vector, specifies the elements of the template that will be used for (bekkers-clements) detection
    %scaler: user inputs desired max amp (or, technically, min amp) -- to examine noise without an event, set scaler to zero
    %two_waveforms is a three-value array (eg 1 50 1):  is there a second waveform, max # samples it can be delayed, amplitude as fraction of first waveform
    %second_waveform_subset is set to zero to have the second waveform use the full waveform (including tapers or prefixes), or to a vector specifying elements of waveform to include
    %variable_waveform_size is a two-value array -- is the first waveform variable, is the second waveform variable (second value meaningless when two_waveforms(1)==0)
        %when variable_waveform_size is [1 0], the first waveform's size will vary, and the second's will be a fixed proportion (two_waveforms(3)) of it every time
        %when variable_waveform_size is [0 1], the first waveform's size will be fixed, and the second's will vary every time, with size determined by scaler, two_waveforms(3), and rand
    %minimum_absolute_peak_amplitude is a two value array, 1st value is gives minimum acceptable magnitude for 1st waveform in pA or zero for no minimum, second element same for 2d waveform
    %amp_window, usually [-3:10], values around the peak amplitude value to average when measuring to estimate the true value
    %detec_alg, 7-element vector, first element is set at 1 to use bekkers-clements detection, 2 to use kudoh-taguchi detection. B-C requires no further settings whereas K-T requires 5:
        %step (usually around 15), the width in samples of a step used in analysis;  roughly corresponds to rise time, though should be set a bit longer than rise time 
        %binwidth (usually around 10), for micro-averaging; peak_time_range (usually around 40), "the temporal period for the peak search";
        %thresh (usually around -10), the initial detection threshold; and noise (usually around -5) baseline noise to ignore, effectively sets a sensitivity value
        %the seventh element is an arbitrarily chosen value (often -1) assigned as 'start time' to allow all the 'no event detected' trials to be readily identified/processed
        %the number of K-T events is capped at two b/c that is the max number of real events
    %plot_detection, yes or no; generally only makes sense when number_reps is very low
    %omniscient_start has two values, binary logicals for using assigned start time versus detected start time for (first value) amplitude and for (second value) rise time and plateau 
    %precede_for_return tells how many samples before the detected event onset time to use as the first point when returning the portion of current trace detected to have an event
%code confirmed date: april 18, 2013

tic
time_vector=1:100;
template=template_generator(template_seeds,time_vector,template_type,taper_template);
templatelength=length(template);templatelengthvec=1:templatelength;
waveform=template*((-1*scaler)/min(template));
detection_template=[0         0         0         0         0   -5.2772   -9.0166  -11.6184  -13.3798  -14.5219  -15.2092  -15.5638  -15.6763  -15.6140  -15.4269  -15.1519  -14.8162  -14.4402  -14.0385  -13.6222  -13.1991  -12.7750  -12.3541  -11.9392-11.5326  -11.1355  -10.7489  -10.3734  -10.0092   -9.6566   -9.3153   -8.9854   -8.6667   -8.3588   -8.0616   -7.7748   -7.4980   -7.2309   -6.9732   -6.7247   -6.4849   -6.2537   -6.0307   -5.8156   -5.6082   -5.4081   -5.2152   -5.0292   -4.8498   -4.6767   -4.5099   -4.3490   -4.1939   -4.0442   -3.9000   -3.7608   -3.6267   -3.4973   -3.3725   -3.2522   -3.1362   -3.0243   -2.9164   -2.8123   -2.7120   -2.6152   -2.5219   -2.4320   -2.3452   -2.2615   -2.1809   -2.1030   -2.0280   -1.9557   -1.8859   -1.8186   -1.7537   -1.6912   -1.6308   -1.5726   -1.5165   -1.4624   -1.4103   -1.3599   -1.3114   -1.2646   -1.2195   -1.1760   -1.1341   -1.0936   -1.0546   -1.0170   -0.9807   -0.9457   -0.9120   -0.8794   -0.8480   -0.8178   -0.7886   -0.7605   -0.7333   -0.7072   -0.6820   -0.6576   -0.6342]

minwindowadder=abs(taper_template(1))+20;
if two_waveforms(1)
    min_window=-5:(minwindowadder+two_waveforms(2));
else
    min_window=-5:minwindowadder; %range around the detected (or assigned) start time within which to search for the minimum value (aka the psc peak)
end

cold_start=ceil(1.5*templatelength);   %the period at the beginning of the noise within which an event is not allowed to begin
if two_waveforms(1)
    dead_end=two_waveforms(2)+ceil(2.5*templatelength); %the period at the end of the noise within which an event is not allowed to begin
else
    dead_end=ceil(2.5*templatelength); %the period at the end of the noise within which an event is not allowed to begin
end

max_plausible_rise_time=50; %in samples -- prevents protracted measurement of 'plateau' when 'rise time' is spuriously high (due to small event --> false criterion peak, far from true start time)
if two_waveforms(1), max_plausible_rise_time=two_waveforms(2)+max_plausible_rise_time;end

if use_noise    %find the possible segmentation sites within the noise
    [below0segs,above0segs]=noise_segmenter(recorded_noise,splisiterequirements,find_unique_segments,breakpoints_mustbe_consecutive)
else
    below0segs=nan;above0segs=nan;
end

kyumsortedseglengths=cumsum(sort([below0segs(2,:)-below0segs(1,:) above0segs(2,:)-above0segs(1,:)])) %reports the length of the shortest noise trace possible for each number of segments

%breaking up variables for use in parfor loop
two_waveforms_1=two_waveforms(1);two_waveforms_2=two_waveforms(2);two_waveforms_3=two_waveforms(3);
variable_waveform_size_1=variable_waveform_size(1);variable_waveform_size_2=variable_waveform_size(2);
omniscient_start_1=omniscient_start(1);omniscient_start_2=omniscient_start(2);
detec_alg_1=detec_alg(1);detec_alg_2=detec_alg(2);detec_alg_3=detec_alg(3);detec_alg_4=detec_alg(4);detec_alg_5=detec_alg(5);detec_alg_6=detec_alg(6);detec_alg_7=detec_alg(7);
double_waveform_excerpt=[-two_waveforms_2:(two_waveforms_2+templatelength)];
minimum_absolute_peak_amplitude_1=minimum_absolute_peak_amplitude(1);minimum_absolute_peak_amplitude_2=minimum_absolute_peak_amplitude(2);

sim_psc_data=[];
parfor i=1:number_reps
    if variable_waveform_size_1
        if minimum_absolute_peak_amplitude_1
            scaled_waveform=0;
            while abs(min(scaled_waveform))<minimum_absolute_peak_amplitude_1
                scaled_waveform=waveform*rand;
            end
        else
            scaled_waveform=waveform*rand;
        end
    else
        scaled_waveform=waveform;
    end
            
    %generate noise
    if use_noise
        noiz=noise_recombinase(below0segs,above0segs,recorded_noise,max_segments);
    else
        noiz=zeros(1,5000);
    end    
        
    viable_noise_length=length(noiz)-(dead_end+cold_start);
    if viable_noise_length>minimal_noise_length
        bb=round(rand*viable_noise_length);
        assigned_startpoint=bb+cold_start;

        padded_event=zeros(1,length(noiz));
        padded_event(assigned_startpoint:(assigned_startpoint+templatelength-1))=scaled_waveform;
        
        if two_waveforms_1
            interevent_interval=randi(two_waveforms_2+1)-1;
            second_start_point=assigned_startpoint+interevent_interval;
            if variable_waveform_size_2
                if minimum_absolute_peak_amplitude_2
                    second_waveform=0;
                    while abs(min(second_waveform))<minimum_absolute_peak_amplitude_2
                        second_waveform=waveform*two_waveforms_3*rand;
                    end
                else
                    second_waveform=waveform*two_waveforms_3*rand;
                end
            else
                second_waveform=waveform*two_waveforms_3;
            end
            if any(second_waveform_subset)
                second_waveform=second_waveform(second_waveform_subset);
            end
            padded_event(second_start_point:(second_start_point+length(second_waveform)-1))=padded_event(second_start_point:(second_start_point+length(second_waveform)-1))+second_waveform;
        else
            interevent_interval=nan;
        end
        
        prenoise_event_size=min(padded_event);
        event_in_noiz=padded_event+noiz;
        
        %detect the event, record the onset time and the size
        %returned_event is based on measured_startpoint rather than assigned_startpoint, and can be returned as zeroes in various failure conditions
        %crit_around_detected_start will also be returned as zeroes if measured_startpoint is too close to the beginning or end
        if detec_alg_1==1
            criterion_trace=template_based_mini_detector(event_in_noiz,detection_template);   if plot_detection, disp('plottting'), figure, plot(event_in_noiz,'k'), hold on, plot(criterion_trace, 'r'), end
            [maxcritamp,maxcritind]=max(criterion_trace);
            measured_startpoint=maxcritind; event_confidence_or_difference=maxcritamp;
            if two_waveforms_1
                if measured_startpoint>two_waveforms_2 && (length(event_in_noiz)-measured_startpoint)>=(two_waveforms_2+templatelength)
                    returned_event=event_in_noiz(measured_startpoint+double_waveform_excerpt)';
                else
                    returned_event=zeros(length(double_waveform_excerpt),1);
                end
                if measured_startpoint>two_waveforms_2 && (length(event_in_noiz)-measured_startpoint)>=(two_waveforms_2+(2*templatelength)) %criterion_trace is shorter than event_in_noiz by templatelength
                    crit_around_detected_start=criterion_trace(measured_startpoint+double_waveform_excerpt);
                else
                    crit_around_detected_start=zeros(1,length(double_waveform_excerpt));
                end
            else
                if measured_startpoint>precede_for_return&&((length(event_in_noiz)-measured_startpoint)>=(templatelength-precede_for_return))
                    returned_event=event_in_noiz(measured_startpoint-precede_for_return+templatelengthvec)';
                else
                    returned_event=zeros(templatelength,1);
                end
                if measured_startpoint>precede_for_return&&(length(event_in_noiz)-measured_startpoint)>=((2*templatelength)-precede_for_return)
                    crit_around_detected_start=criterion_trace(measured_startpoint-precede_for_return+templatelengthvec);
                else
                    crit_around_detected_start=zeros(1,templatelength);
                end
            end
            returned_detection_info=crit_around_detected_start';
            
        elseif detec_alg_1==2
            kt_events=find_multiple_peaks(event_in_noiz,[detec_alg_2 detec_alg_3 detec_alg_4 detec_alg_5 detec_alg_6])
            if length(kt_events)>5  %five values returned per kt-detected event
                kt_events=kudoh_taguchi_result_sorter(kt_events);
                kt_events=kt_events(1:10)
                measured_startpoint=kt_events(1); event_confidence_or_difference=kt_events(10)-kt_events(5);
                if measured_startpoint>precede_for_return&&((length(event_in_noiz)-measured_startpoint)>=(templatelength-precede_for_return))
                    returned_event=event_in_noiz(measured_startpoint-precede_for_return+templatelengthvec)';
                else
                    returned_event=zeros(templatelength,1);
                end
                if plot_detection,close,plot(event_in_noiz,'k'),hold on,
                    plot(kt_events(1,:),kt_events(2,:),'rs'),plot(kt_events(3,:),kt_events(4,:),'ro'),
                    plot(kt_events(6,:),kt_events(7,:),'gs'),plot(kt_events(8,:),kt_events(9,:),'go'),
                end
            elseif length(kt_events)==5  %five values returned per kt-detected event
                kt_events=[kt_events;zeros(5,1)];
                measured_startpoint=kt_events(1); event_confidence_or_difference=kt_events(5);
                if measured_startpoint>precede_for_return&&((length(event_in_noiz)-measured_startpoint)>=(templatelength-precede_for_return))
                    returned_event=event_in_noiz(measured_startpoint-precede_for_return+templatelengthvec)';
                else
                    returned_event=zeros(templatelength,1);
                end                 
                if plot_detection,close,plot(event_in_noiz,'k'),hold on,plot(kt_events(1,:),kt_events(2,:),'rs'),plot(kt_events(3,:),kt_events(4,:),'ro'),end
            else
                kt_events=zeros(10,1); measured_startpoint=detec_alg_7; event_confidence_or_difference=nan;
                returned_event=zeros(templatelength,1);
            end
            returned_detection_info=kt_events;
        end
        
        if omniscient_start_1
            [~,mineventind]=min(event_in_noiz(assigned_startpoint+min_window));  
            overallminindex=assigned_startpoint+min_window(1)+mineventind-1;
            corrected_minamp=mean(event_in_noiz(overallminindex)+amp_window);
        elseif measured_startpoint > abs(min_window(1)) && measured_startpoint <= length(event_in_noiz)-min_window(end)
            [~,mineventind]=min(event_in_noiz(measured_startpoint+min_window));  
            overallminindex=measured_startpoint+min_window(1)+mineventind-1;
            corrected_minamp=mean(event_in_noiz(overallminindex)+amp_window);
        else
            overallminindex=nan;corrected_minamp=nan;
        end
        
        
        if omniscient_start_2 && ~isnan(overallminindex)
            rise_time=overallminindex-assigned_startpoint;startpoint_for_rise_plateau=assigned_startpoint;
        elseif measured_startpoint > abs(min_window(1)) && measured_startpoint<(length(event_in_noiz)-max_plausible_rise_time) && ~isnan(overallminindex)
            rise_time=overallminindex-measured_startpoint;startpoint_for_rise_plateau=measured_startpoint;
        else 
            rise_time=nan;
        end
        rise_slope=corrected_minamp/rise_time;
        if rise_time < max_plausible_rise_time
            [plateau_onset,longest_plateau]=find_riseplateau(startpoint_for_rise_plateau,overallminindex,event_in_noiz);
        else
            plateau_onset=-1;longest_plateau=-1;
        end
        
        results_this_trial=[assigned_startpoint;prenoise_event_size;measured_startpoint;corrected_minamp;event_confidence_or_difference;rise_slope;plateau_onset;longest_plateau;...
            interevent_interval;returned_event;returned_detection_info];
        sim_psc_data=[sim_psc_data results_this_trial];
        
    else
        disp('noise too short')
        disp(num2str(viable_noise_length))
    end       
end
toc