% Jon Newman and Ming-fai Fong
% GA Tech
% 20 May 2010


% This script modifies data collected on 20100519 where intracellular burst
% input from previously recorded data was provided to cells as input

%% Load intacellular data and begin processing struct
fid = '20100519_3_3';
directory = 'C:\Users\Jon\Desktop\20100519\';
load(strcat([directory fid '.mat']));

synburst = {};
time = [];
vmon = [];
imon = [];

for i = 1:500;

    iid = strcat(['Trace_3_3_' num2str(i) '_1']);
    vid = strcat(['Trace_3_3_' num2str(i) '_2']);


    % if the data matrix exists, put the intracellular current and voltage
    % traces along with a single time trace
    if exist(iid,'var') && exist(vid,'var')
        imat = eval(iid);
        vmat = eval(vid);
        imon = [imon imat(:,2)];
        vmon = [vmon vmat(:,2)];
        if isempty(time)
            time = vmat(:,1);
        end
    end


end

synburst.intrarec.Fs = 20e3;
synburst.intrarec.time = time;
synburst.intrarec.imon = imon;
synburst.intrarec.vmon = vmon;


%% Get stimulation indicies
load('C:\TDTDATA\20100420\20100420-25\20100420-25_cont_320to780');

synburst.burststim.extracellspks = data.extrarec.spiketime;
synburst.burststim.ASDR = data.extrarec.ASDR;
synburst.burststim.bursts = data.extrarec.bursts;
synburst.burststim.stimind = data.burststim.apporder;
synburst.burststim.stimwave = data.burststim.burstinputsnip';

%% Spike detection
%spk detection parameters

maxv = 60e-3;
maxpk = 50e-3;
minpk = -10e-3;
mpd = 40; % 1 ms @ fs = 20khz

spkloc = cell(size(size(vmon,2)));
spkpeak = cell(size(size(vmon,2)));

for i = 1:size(vmon,2)

    tempv = vmon(:,i);
    if max(tempv) < maxv
        dvdt = diff(tempv);

        % peak finding algorithm. (1) look for rapid changes in voltage
        % using estimated time derivative. (2) Find the voltage max within
        % a mpd window after this rapid increase. (3) remove peaks that are
        % not physiological.
        [dvpk dvloc]= findpeaks(dvdt,'MINPEAKHEIGHT',1e-3,'MINPEAKDISTANCE',mpd);
        vloc = zeros(size(dvpk));
        vpk = zeros(size(dvpk));
        for j = 1:length(dvloc)
            snip = tempv(dvloc(j):dvloc(j)+mpd);
            [vp vind] = max(snip);
            vpk(j) = vp;
            vloc(j) = vind+dvloc(j);
        end

        vloc(vpk>maxpk) = [];
        vloc(vpk<minpk) = [];
        vpk(vpk>maxpk) = [];
        vpk(vpk<minpk) = [];

    else
        vloc = [];
        vpk = [];
    end

    spkloc{i} = vloc;
    spkpeak{i} = vpk;

end

synburst.spks.spkloc = spkloc;
synburst.spks.spkpeak = spkpeak;
synburst.spks.param.maxv = maxv;
synburst.spks.param.maxpk = maxpk;
synburst.spks.param.minpk = minpk;
synburst.spks.param.mpd = mpd;


%% Get average firing rate information
meanfr = zeros(size(synburst.burststim.stimwave,2),1); % average firing rate ellicited by each unique stim type
stdfr = zeros(size(synburst.burststim.stimwave,2),1); % variance in firing rate ellicited by each unique stim type

for i = 1:length(meanfr)
    
    Ind = find(synburst.burststim.stimind == i);
    Ind(Ind>size(synburst.spks.spkloc,2)) = [];
    tot = [];
    
    for j = Ind
        tot = [tot length(synburst.spks.spkloc{j})];
    end
    tot = tot/max(synburst.intrarec.time);
    meanfr(i) = mean(tot);
    stdfr(i) = std(tot);

end

synburst.stimresp.meanfr = meanfr;
synburst.stimresp.stdfr = stdfr;

%% Examine Stimuli

% get maximal of each stim waveform
synburst.burststim.meanamp = mean(synburst.burststim.stimwave);

save(strcat([fid '_synburst']),'synburst');

%% Raster across stimuli sorted by increasing stimulus value
numplots = max(synburst.burststim.stimind);
m = 3;
n = ceil(numplots/m);

trialcount = ones(numplots,1);
[Y Ind] = sort(synburst.burststim.meanamp);

h = figure();
for i = 1:length(synburst.spks.spkloc)
    
    subplot(n,m,find(Ind == synburst.burststim.stimind(i)))
    hold on
    plot((1/synburst.intrarec.Fs)*synburst.spks.spkloc{i},...
        trialcount(synburst.burststim.stimind(i))*ones(size(synburst.spks.spkloc{i})),'k.')
    trialcount(synburst.burststim.stimind(i)) = trialcount(synburst.burststim.stimind(i)) + 1;
    xlim([0 max(synburst.intrarec.time)])
    
    
end

for i = 1:numplots
    
    subplot(n,m,i)
    hold on
    ylim([0 max(trialcount)+1])
    if i ~= numplots
        set(gca,'xticklabel',[],'yticklabel',[])
    else
        ylabel('trial')
        xlabel('time (sec)')
    end
    
end

figdir = strcat(strcat([directory fid 'FIG']));
figfid = strcat([fid '_RasterAcrossStim_IncreasingMeanInput']);
if ~exist(figdir,'dir')
    mkdir(figdir);
end
saveas(h, strcat([figdir '\' figfid '.fig']))
saveas(h, strcat([figdir '\' figfid '.jpg']))

%% Raster across stimuli sorted by increasing ASDR value during the input
numplots = max(synburst.burststim.stimind);
m = 3;
n = ceil(numplots/m);

trialcount = ones(numplots,1);
[Y Ind] = sort(synburst.burststim.bursts.bursts(:,2));

h = figure();
for i = 1:length(synburst.spks.spkloc)
    
    subplot(n,m,find(Ind == synburst.burststim.stimind(i)))
    hold on
    plot((1/synburst.intrarec.Fs)*synburst.spks.spkloc{i},...
        trialcount(synburst.burststim.stimind(i))*ones(size(synburst.spks.spkloc{i})),'k.')
    trialcount(synburst.burststim.stimind(i)) = trialcount(synburst.burststim.stimind(i)) + 1;
    xlim([0 max(synburst.intrarec.time)])
end

for i = 1:numplots
    
    subplot(n,m,i)
    hold on
    ylim([0 max(trialcount)+1])
    if i ~= numplots
        set(gca,'xticklabel',[],'yticklabel',[])
    else
        ylabel('trial')
        xlabel('time (sec)')
    end
    
end

figdir = strcat(strcat([directory fid 'FIG']));
figfid = strcat([fid '_RasterAcrossStim_IncreasingASDR']);
if ~exist(figdir,'dir')
    mkdir(figdir);
end
saveas(h, strcat([figdir '\' figfid '.fig']))
saveas(h, strcat([figdir '\' figfid '.jpg']))

%% Look at firing rate versus mean stimulus current
[y Ind] = sort(synburst.burststim.meanamp);

h = figure();
subplot(121)
h1 = errorbar(synburst.burststim.meanamp(Ind),...
    synburst.stimresp.meanfr(Ind),...
    synburst.stimresp.stdfr(Ind),'k');
xlabel('Mean Burst Input Current (pA)')
ylabel('FR of Post Syn. Neuron (Hz+-SD)')
hChildren = get(h1, 'children');
set(hChildren(1), 'linewidth',2)


[y Ind] = sort(synburst.burststim.bursts.bursts(:,2));

subplot(122)
h2 = errorbar(synburst.burststim.bursts.bursts(Ind,2),...
    synburst.stimresp.meanfr(Ind),...
    synburst.stimresp.stdfr(Ind),'k');
xlabel('ASDR During Burst Input Current (Hz)')
ylabel('FR of Post Syn. Neuron (Hz+-SD)')
hChildren = get(h2, 'children');
set(hChildren(1), 'linewidth',2)

figdir = strcat(strcat([directory fid 'FIG']));
figfid = strcat([fid '_InputVersusFR']);
if ~exist(figdir,'dir')
    mkdir(figdir);
end
saveas(h, strcat([figdir '\' figfid '.fig']))
saveas(h, strcat([figdir '\' figfid '.jpg']))

