% Jon Newman
% Georgia Tech
% 22 Apr 2010

% Look at instantaneous variance of intracellular current record versus
% culture wide spike rate (analogous to input rate).

% Specify tank and block and plot parameters
% *************************************************************************
tankID = 'C:\TempData\20100420\20100420-34\';
windsize = 200; % ms
recRange = [-100 50]; %pA
zoomrange = [30 40]; % crop range, sec
asdrtimebin = 0.05; % sec
% *************************************************************************

fid = queryfid(tankID);
load(fid);

%% Timeline plot

% Raster and ASDR for Extracellular record
h1 = figure();
st = data.extrarec.spiketime;
st = st-min(st);
if isempty(zoomrange)
    xlimits = [0 max(st)];
else
    xlimits = zoomrange;
end

subplot(511)
plot(st,data.extrarec.channel+0.33*randn(size(st)),'k.','MarkerSize',1);
xlim(xlimits);
set(gca, 'XTickLabelMode', 'Manual')
set(gca, 'XTick', [])
title(fid,'Interpreter','none');
ylabel('Channel');
ylim([0 64]);
subplot(512)
[ASDR bins] = asdr(st,asdrtimebin);
xlim(xlimits);
set(gca, 'XTick', [])

%Intracellular current trace, raw
intracurr = data.intrarec.imon;
time = (1/data.sampfreq)*(1:length(intracurr));

subplot(513)
plot(time,intracurr,'k');
xlim(xlimits);
ylabel('I_{pip} (pA)');
set(gca, 'XTick', [])

windsamp = ceil(data.sampfreq*(windsize/1000));
meaninput = -smooth(intracurr,windsamp);
corrinput = smooth(intracurr.*intracurr,windsamp);
instvar = sqrt(corrinput-meaninput.^2);

% Instanteous current input rate estimation
subplot(514)
plot(time,meaninput,'k')
xlim(xlimits);
ylabel('<I_{mem}> (pA)');
set(gca, 'XTick', [])

% Instanteous SD estimation
subplot(515)
plot(time,instvar,'k')
xlim(xlimits);
ylabel('SD(I_{mem}) (pA)');
xlabel('time (sec)')

findfigs
beep;
reply = input('Save Figure? [y/n]: ','s');
if isempty(reply)
    reply = 'y';
end

if reply == 'y'
    figdir = strcat([tankID fid(1:end-4) '_FIG']);
    figfid = strcat([fid(1:end-4) '_Imem_MEAdata_MomentRelations_FIG']);
    if ~exist(figdir,'dir')
        mkdir(figdir);
    end
    saveas(h1, strcat([figdir '\' figfid '.fig']))
    saveas(h1, strcat([figdir '\' figfid '.jpg']))
end

close all

%% Parametic Plots
meaninput = downsample(meaninput,20);
instvar = downsample(instvar,20);
h2 = figure();
subplot(2,2,1:2)
rainbowplot(meaninput',instvar')
xlabel('<I_{mem}> (pA)');
ylabel('SD(I_{mem}) (pA)');

N = 100; % number of segments
len1 = floor(length(meaninput)/N);
len2 = floor(length(ASDR)/N);
instmeands = [];
instvards = [];
for i = 1:N
    
    temp1 = meaninput((i*len1)-len1+1:(i*len1));
    temp2 = instvar((i*len1)-len1+1:(i*len1));
    instmeands = [instmeands;resample(temp1,len2,len1)];
    instvards = [instvards;resample(temp2,len2,len1)];
    
end

subplot(223)
rainbowplot(ASDR(1:length(instmeands)),instmeands')
xlabel(strcat(['ASDR (spks/',num2str(asdrtimebin),' sec)']));
ylabel('<I_{mem}> (pA)');
subplot(224)
hold on
rainbowplot(ASDR(1:length(instvards)),instvards')
xlabel(strcat(['ASDR (spks/',num2str(asdrtimebin),' sec)']));
ylabel('SD(I_{mem}) (pA)');

findfigs
beep;
reply = input('Save Figure? [y/n]: ','s');
if isempty(reply)
    reply = 'y';
end

if reply == 'y'
    figdir = strcat([tankID fid(1:end-4) '_FIG']);
    figfid = strcat([fid(1:end-4) '_Imem_MEAdata_ParametricMoments_FIG']);
    if ~exist(figdir,'dir')
        mkdir(figdir);
    end
    saveas(h2, strcat([figdir '\' figfid '.fig']))
    saveas(h2, strcat([figdir '\' figfid '.jpg']))
end

close all
