function TineAnalysis(varargin)

% By default make the GUI
%from the command line with no parameters
if (nargin == 0)
    fcn = 'makeGUI';
else
    fcn=varargin{1};
    state = get(findobj('tag','AnalysisFigure'),'userdata');
end

%This is the main decision point of the function.
%The switch statement is executed once-per-fuction call
switch fcn

    %This code is executed ONCE when the function enters with no arguments
    case 'makeGUI'
        state.myname = mfilename;
        fig = figure('Position',[100 100 1000 600],...
            'NumberTitle','off',...
            'Name','Tine Analysis 20feb2007',...
            'Interruptible','off',...
            'doublebuffer','on',...
            'HandleVisibility','on',...
            'tag', 'AnalysisFigure', ...
            'toolbar','figure',...
            'Color',[.7 .7 .7]);

        %=== file and channels ===============================
        %         x=0; y=0; w=60; h=20; sp=10;
        %         uicontrol(...
        %             'style', 'pushbutton', ...
        %             'position', [x y w h], ...
        %             'enable','on',...
        %             'tag','textfile',...
        %             'callback', [state.myname,' getfile'] ,...
        %             'TooltipString','Opens a text file containing a single column of voltages with extension *.txt',...
        %             'string','Text File');

        x=0; y=0; w=70; h=17; sp=5;
        uicontrol(...
            'style', 'pushbutton', ...
            'position', [x y w h], ...
            'enable','on',...
            'tag','logfile',...
            'callback', [state.myname,' getfile'] ,...
            'TooltipString','Opens a DAQ log file with *.daq extension',...
            'string','Log File');

        state.ChInvert = uicontrol(...
            'style', 'checkbox', ...
            'position', [x y+h w h], ...
            'enable','on',...
            'TooltipString','Check the box if Channel 1 contains stimulus',...
            'string','Stim in 1');

        %=== sample rate ======================================
        x=x+w+sp;
        w = 55;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','sample rate');
        x=x+w;
        state.FsSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Sample rate in samples/second',...
            'string','20000');

        %         %=== Artifact supression ===============================
        x=x+w+sp;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','Artifact');
        x=x+w;
        state.artifactSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Number of samples to null after a stimulus',...
            'string','20');

        %=== median filtering ===============================
        x=x+w+sp;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','Sample N');
        x=x+w;
        state.medianSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Median filter length in samples-- should be somewhat longer than longest spike',...
            'string','40');

        state.medianUse = uicontrol(...
            'style', 'checkbox', ...
            'position', [x-w h 100 h], ...
            'enable','on',...
            'TooltipString','Use Median filter?',...
            'string','Median filter');
        %=== Threshold =====================================
        x=x+w+sp;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','Threshold');
        x=x+w;
        state.noiseSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Minimum size a spike must be to be counted, can be pos or neg',...
            'string','0.1');

        %=== Kill extra peaks =====================================
        x=x+w+sp;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','Sample N');
        x=x+w;
        state.peakkillSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Number of samples within which two peaks will be counted only as one',...
            'string','25');

        state.peakKillUse = uicontrol(...
            'style', 'checkbox', ...
            'position', [x-w h 100 h], ...
            'enable','on',...
            'TooltipString','Use peak kill?',...
            'string','Peak Kill');

        %=== Histogram control =====================================
        x=x+w+sp;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','NHistBins');
        x=x+w;
        state.histbinsSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Number of bins in the final histogram',...
            'string','50');
        x=x+w+sp;
        uicontrol(...
            'style', 'text', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'string','MaxBin');
        x=x+w;
        state.histLimitSet = uicontrol(...
            'style', 'edit', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'TooltipString','Maximum bin edge in the final histogram',...
            'string','0.020');

        x=x+w+sp;
        uicontrol(...
            'style', 'pushbutton', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'callback', [state.myname,' doanalysis'] ,...
            'string','Analyze');

        x=x+w+sp;
        uicontrol(...
            'style', 'pushbutton', ...
            'position', [x 0 w h], ...
            'enable','on',...
            'callback', [state.myname,' saveISI'] ,...
            'string','Save ISI');

        x=650;w=200;
        state.histNote = uicontrol(...
            'style', 'edit', ...
            'position', [x 20 w h], ...
            'enable','on',...
            'TooltipString','Additional annotation for histograms',...
            'string','Modify this Histogram subtitle');

    case 'saveISI'
        [fname,fpath] = uiputfile('*.txt');
        try
            fid = fopen([fpath,fname],'w');
        catch
            errordlg('You must specify a filename')
            return %bail out of callback
        end
        fprintf(fid,'%s\n', [state.fpath,state.fname]);
        fprintf(fid,'%s\n', get(state.histNote,'string'));
        fprintf(fid,'Analysis region=%f to %f sec\n\n',state.rect(1),state.rect(2));
        %
        fprintf(fid,'burst#, latency, duration, spike#\n')
        for i=1:length(state.duration)
            fprintf(fid,'%4.0f, %8.5f, %8.5f, %4.0f\n',...
                i,state.latency(i),state.duration(i),state.spikeCount(i))
        end
        %
        fprintf(fid,'\n#bursts, latency,  STD,    duration,  STD,    spike#,  STD\n')
        fprintf(fid,'%4.0f, %8.5f, %8.5f, %8.5f, %8.5f, %8.5f, %8.5f\n',...
            length(state.duration), ...
            mean(state.latency), std(state.latency),...
            mean(state.duration), std(state.duration),...
            mean(state.spikeCount), std(state.spikeCount))
        %
        %compute isi sequence number with a burst
        %if ISI>0.1 second set it to sequence #0
        isiSeq = -ones(size(state.isi));
        isiSeq(1) = 1;
        isiIndex = find(state.isi>0.1) ;
        isiSeq(isiIndex) = 0;
        for i = 2:length(state.isi); %for each stimulus burst
            if isiSeq(i) == -1
                isiSeq(i) = isiSeq(i-1) + 1;
            end
        end
        fprintf(fid,'\nList of all ISIs in the region\n')
        for i = 1:length(state.isi)
            fprintf(fid,'%f, %d\n', state.isi(i), isiSeq(i));
        end
        fclose(fid);

    case 'getfile'

        [fname,fpath] = uigetfile('*.daq');
        try
            v = daqread([fpath,fname]);
        catch
            errordlg('You must specify a filename')
            return %bail out of callback
        end
        daqvalues = daqread([fpath,fname],'info');
        nChan = length(daqvalues.ObjInfo.Channel) ;
        str = daqvalues.ObjInfo.Channel(1).ChannelName;
        sFs = (daqvalues.ObjInfo.SampleRate);
        set(state.FsSet,'string',sFs,'enable','off')
        if get(state.ChInvert,'value')==0 %data in channel 1
            y = v(:,daqvalues.ObjInfo.Channel(1).Index);
            stim = v(:,daqvalues.ObjInfo.Channel(2).Index)>0.5 ;
        else %stim in channel 1
            y = v(:,daqvalues.ObjInfo.Channel(2).Index);
            stim = v(:,daqvalues.ObjInfo.Channel(1).Index)>0.5 ;
        end

        %prevent overflow
        amp = max(abs(y)) ;
        state.y = y/amp ;
        state.stim = stim ;
        state.fpath = fpath;
        state.fname = fname;
        Fs = str2num(get(state.FsSet,'string'));
        state.rect(1) = .01;
        state.rect(2) = length(y)/Fs -.01;
        state.newdata = 1;
        state=doAnalysis(state);
        fig = findobj('tag','AnalysisFigure');
        set(fig,'WindowButtonDownFcn',...
            [ 'fig = findobj(''tag'',''AnalysisFigure'');'...
            'state=get(fig,''UserData'');' ...
            'if (gca == state.ax1);' ...
            'rbbox;'...
            'initialpt=get(gca,''currentpoint'');'...
            'state.rect(1)=initialpt(1,1);'...
            'end;'...
            'set(fig,''UserData'',state);' ...
            ]);
        % 'Fs = str2num(get(state.FsSet,''string''));',...
        %                 'if(state.rect(1)>state.rect(2));'...
        %                     'temp=state.rect(1);state.rect(1)=state.rect(2);state.rect(2)=temp;end;'...
        %                  'set(state.rectH,''position'',[state.rect(1), -1, (state.rect(2)-state.rect(1)), 2]);'....
        set(fig,'WindowButtonUpFcn',...
            [ 'fig = findobj(''tag'',''AnalysisFigure'');'...
            'state=get(fig,''UserData'');' ...
            'if (gca == state.ax1);' ...
            'initialpt=get(gca,''currentpoint'');'...
            'state.rect(2)=initialpt(1,1);'...
            'Fs = str2num(get(state.FsSet,''string''));',...
            'if(state.rect(1)>state.rect(2));'...
            'temp=state.rect(1);state.rect(1)=state.rect(2);state.rect(2)=temp;end;'...
            'set(state.rectH,''position'',[state.rect(1), -1, (state.rect(2)-state.rect(1)), 2]);'....
            'end;'...
            'set(fig,''UserData'',state);' ...
            ]);

    case 'doanalysis'
        state=doAnalysis(state);

end %of main case statement
set(findobj('tag','AnalysisFigure'), 'userdata', state);

%==========================================================================
function state=doAnalysis(state)
%get a noise estimate (biased by the spikes, however)
y = state.y;
my = mean(y);
%first estimate of noise level
normY = sum(abs(state.y-mean(state.y)))/length(state.y);
%now get a tighter bound on the noise by zapping spikes
yy=y;
yy(find(y>(2*normY+my))) = my ; %remove spikes
yy(find(y<(-2*normY+my))) = my ;
normY = sum(abs(yy-mean(yy)))/length(yy);
stim = state.stim;

%settable parameters
Fs = str2num(get(state.FsSet,'string'));
% fraction of negative spike to consider artifact (relative to biggest)
StimSupress = str2num(get(state.artifactSet,'string'));
%median filter length in samples-- should be somewhat longer than longest spike
filtLength = str2num(get(state.medianSet,'string'));
medianOn = get(state.medianUse,'value');
% size a spike must be to be counted, in units of the background noise
noiseMultiplier = str2num(get(state.noiseSet,'string'));
%number of samples over which two peaks should be called just one
peakKillRange = str2num(get(state.peakkillSet,'string'));
killOn = get(state.peakKillUse,'value');
% number of bins in the final histogram
numberBins = str2num(get(state.histbinsSet,'string'));
limitBins = str2num(get(state.histLimitSet,'string'));

time = linspace(0, length(y)/Fs, length(y));

figure(findobj('tag','AnalysisFigure'));

np=3;
%==================================
subplot(np ,1,1)
if state.newdata
    cla
    state.newdata=0;
    state.rectH = rectangle...
        ('position',[state.rect(1), -1, (state.rect(2)-state.rect(1)), 2],...
        'facecolor',[1 1 .7]) ;
    hold on
    plot(time,y)
    plot(time,stim,'r')
    set(gca,'ylim',[-1.1 1.1])
    title([state.fpath,state.fname],'interpreter','none');
end

state.ax1 = gca ;
%SET the stimulus artifact threshold
%find very negative spikes
% minAmpStim = min(y);
% set(state.stimLineH,'xdata',[time(1) time(end)],...
%     'ydata',[StimSupressRatio*minAmpStim StimSupressRatio*minAmpStim])
stimLoc1 = find(stim==1);
stimLoc = [stimLoc1; stimLoc1-1; stimLoc1-2; stimLoc1-3];
%bracket the locations found
for i=1:StimSupress
    stimLoc = [stimLoc ; stimLoc1+i];
end
%make sure that all locations actually exist in the vector
stimLoc(find(stimLoc<1)) = 1;
stimLoc(find(stimLoc>length(y))) = length(y);
%replace each location with the mean value
data = y;
data(stimLoc) = mean(y);
%find stim interval > 0.1 second and break into groups of spikes
stimLoc = find(stim==1);
%eliminate stim outside analysis area
stimLoc(find(stimLoc<state.rect(1)*Fs)) = [];
stimLoc(find(stimLoc>state.rect(2)*Fs)) = [];
%get big intervals and the last one
stimEnd = [find(diff(stimLoc)>(0.1*Fs));length(stimLoc)] ;
% and indicate with an x
plot(time(stimLoc(stimEnd)), 1.0, 'bx')

%==================================
subplot(np,1,2)
cla
%prevent underflow
amp = max(abs(data)) ;
data = data/amp ;
%subtract median filtered data
if medianOn
    data = data - medfilt1(data,filtLength);
end
data = data-mean(data);
%peak detection and amplitude threshold
%biggestPeak = max(data);
%threshold = thresholdFraction*(biggestPeak-baseline)+ baseline ;
threshold = noiseMultiplier  ;
plot(time,data, 'b')
hold on
%plot(time(stimLoc(stimEnd)), threshold, 'rx')
line([0,length(data)/Fs],[threshold,threshold],'color','r')
title('Filtered and threshold')

%====================================
subplot(np,1,3)
cla
if threshold>0
    %isolate positive peaks ABOVE threshold
    data = (data>threshold).*data;
    %find peaks where first derivitive changes sign
    dataDiff = diff(data);
    peakPos = find((dataDiff(1:end-1)>0) & (dataDiff(2:end)<=0))+1;
else
    %isolate negative peaks BELOW threshold
    data = (data<threshold).*data;
    %find peaks where first derivitive changes sign
    dataDiff = diff(data);
    peakPos = find((dataDiff(1:end-1)<0) & (dataDiff(2:end)>=0))+1;
end
plot(time,data, 'b')
hold on

%eliminate peaks outside analysis range
peakPos(find(peakPos<state.rect(1)*Fs)) = [];
peakPos(find(peakPos>state.rect(2)*Fs)) = [];

%eliminate all but the biggest peak separated by <peakKillRange samples
if killOn
    peakAmps = data(peakPos);
    peakDiff = diff(peakPos);
    peakKill = 0*length(peakPos);
    for i=1:length(peakDiff)
        if (peakDiff(i)<peakKillRange) & (peakAmps(i)>peakAmps(i+1))
            peakKill(i+1) = i+1 ;
        end
        if (peakDiff(i)<peakKillRange) & (peakAmps(i+1)>peakAmps(i))
            peakKill(i) = i ;
        end
    end
    peakPos(find(peakKill>0)) = [];
end
%get the actual amps for plotting
peakAmps = data(peakPos);
plot(peakPos/Fs,peakAmps,'r.')
title('Marked spikes')
plot(time(stimLoc(stimEnd)), min(peakAmps), 'k.')
xlabel('time (seconds)')
%reset count
state.duration=[];
%find and mark latency
%subplot(np,1,3)
hold on
for i=1:length(stimEnd)
    %find the closest peak to the end of the stim
    state.latency(i) = min(abs(time(stimLoc(stimEnd(i)))-peakPos/Fs));
    line([time(stimLoc(stimEnd(i))),...
        time(stimLoc(stimEnd(i)))+state.latency(i)],...
        [1 1]*min(peakAmps),'color','k','linewidth',2)
end

%find and mark duration use 0.1 seconds as break value to indicate end of
%burst
for i=1:length(stimEnd)
    %eliminate peaks before stim start
    possibleIndex = find(peakPos>stimLoc(stimEnd(i)));
    if length(possibleIndex>0)
        %find ISI>0.1 which marks the end of the burst
        ISIdur = find(diff(peakPos(possibleIndex))/Fs > 0.1);
        if i<length(stimEnd)
            %use the first ISIdur greater than 0.1
            state.duration(i) = peakPos(ISIdur(1)+possibleIndex(1)-1)/Fs -...
                (time(stimLoc(stimEnd(i)))+state.latency(i)) ;
            state.spikeCount(i) = (ISIdur(1)+possibleIndex(1)-1) - possibleIndex(1) + 1 ;
        else %the last spike needs special handling
            length(possibleIndex)
            state.duration(i) = peakPos(end)/Fs - (time(stimLoc(stimEnd(i)))+state.latency(i)) ;
            state.spikeCount(i) = possibleIndex(end) - possibleIndex(1) + 1 ;
        end

        line(...
            [time(stimLoc(stimEnd(i)))+state.latency(i),...
            time(stimLoc(stimEnd(i)))+state.latency(i)+state.duration(i)]...
            ,[1 1]*min(peakAmps),'color','g','linewidth',2)
        % mark number of spikes
        text(time(stimLoc(stimEnd(i)))+state.latency(i)+state.duration(i),...
            min(peakAmps), num2str(state.spikeCount(i),' %d'),'HorizontalAlignment','left')
    end
end
title(['#bursts=',num2str(length(state.duration)) ])
%===========================================
%ISI
figure(2);clf
isi = diff(peakPos)/Fs ;
state.isi = isi ;
hist(isi(find(isi<limitBins)),linspace(0,limitBins,numberBins));
title([state.fname],'interpreter','none');
text(.1,.95,[get(state.histNote,'string'),...
    '  Region=',num2str(state.rect(1)),' to ',num2str(state.rect(2)),' sec'],...
    'units','normalized');
xlabel('isi (seconds)')
ylabel('count')
figure(1)
%==== end doAnylsis ================================================