function info = rate_fano_psth(data,interval,showplot)
%********* JUDE MITCHELL (jude@salk.edu):  5/18/2008
%******** function info = rate_fano_psth(data,interval,showplot)
%**********
%    one example:  
%                  info = rate_fano_psth('mnar4_u2',[],1);  %uses default
%                                                           % interval
%                
%                  load 'msiv3_u8';
%                  info = rate_fano_psth(data,data.sustained,1)
%
%                  load 'mnar4_u2';
%                  info = rate_fano_psth(data,data.sustained,1);
%
%                  load 'msir27_u3';
%                  info = rate_fano_psth(data,data.morepause,1);  
%
%***** general: computes smoothed firing rate and the Fano factor
%*****          at several different smooth windows and count intervals
%*****          and plots those results  
%** inputs:
%********** data - all relevent data fields for a neuron
%******
%****** data.label(nChannels) = {'SU1' ; 'MU1'; 'LFP1'; 'EYE1'};
%****** data.trials{ntrials} = [[Su1 data 1:T];[Mu1 data from 1:T];...];
%****** data.spontaneous{ntrials} = [Su1 data 1:250ms]; %spike data only,
%******                % data from interval when fixation first occurs
%******                % and before the stimuli first onset
%****** data.time{ntrials} = [1:T]
%****** data.attend(ntrials) = 1 if attended in RF, 2 is 2 of 4 unatt, 3 is
%******                             for the 1 of 4 unattended case
%****** data.exactspikes{ntrials} = exact single unit spike times per
%******                              per trial to 0.025 ms precision
%****** data.fsample = 1000;
%****** data.sustained = XA:XB  % to mark the 800 ms sustained period
%****** data.morepause = -250 to +250 after 1000ms pause, 1500 ms
%****** data.pause = exact period of 1000ms pause
%****** data.surround = -950 to -450 ms before pause, stimuli are outside
%******                   the RF (period of non-visual response?)
%****** data.waveform = [1:32];
%****** data.isolation = 1 - well isolated single unit, 
%******                  2 - clear cluster but some overlap
%******                  3 - multi-unit, large overlap
%******                  4 - well isolated, but lost midway in session
%******
%********** interval - interval of specific analysis, 1xN array of timepoints
%**********            - use data from this interval to evaluate
%**********              significant differences in rate and in Fano
%********** showplot - to plot out results
%
%******* input:
%******
%****** data.label(nChannels) = {'SU1' ; 'MU1'; 'LFP1'; 'EYE1'};
%****** data.trials{ntrials} = [[Su1 data 1:T];[Mu1 data from 1:T];...];
%****** data.time{ntrials} = [1:T]
%****** data.attend(ntrials) = 1 if attended in RF
%****** data.fsample = 1000;
%****** data.sustained = XA:XB  % to mark the 800 ms sustained period
%****** data.waveform = [1:32];
%*****
%****** interval - interval of analysis, 1xN array of timepoints
%****** showplot - set this to 1 to see the results plotted out
%*****
%******  output:
%******      info.cintervals = [1 x N] array of counting intervals
%******      info.cinfo{1xN} = detailed fano and rate info (see below)
%******      info.cfano(C,N) = mean fano as function on interval each
%******                         for each of the C attention conditions
%******      info.sfano(C,N) = same as cfano, but sem of the mean
%******
%******  cinfo:
%*****       C - the number of attention conditions 
%***** cinfo.timebin{1xC} - time midpoints of counting intervals
%***** cinfo.smorate{1xC} - gaussian smooth firing rate
%***** cinfo.ratebin{1xC} - mean spike counts in bins
%***** cinfo.varbin{1xC} -  variance spike counts in bins
%***** cinfo.fanobin{1xC} - fano factor, NaN if no spikes in a bin
%*****
%***** cinfo.crange - count windows falling in interval of analysis
%***** cinfo.meanfano(C) - mean fano over range of analysis
%***** cinfo.semfano(C) - sem of fano over range of analysis
%*****
%***** cinfo.airate - AI of rate mod (A-U)/(A+U)
%***** cinfo.ratepval - significance cond 1 vs 2 in firing rate diff
%***** cinfo.aifano - AI of fano mod
%***** cinfo.fanopval - significance cond 1 vs 2 in fano factor diff
%************
%****** recap:
%********
%******** function info = rate_fano_psth(data,interval,showplot)
%**********
%    one example:  
%                  info = rate_fano_psth('mnar4_u2',[],1);  %uses default
%                                                           % interval
%                
%                  load 'msiv3_u8';
%                  info = rate_fano_psth(data,data.sustained,1)

info = [];

%******* if data is just a name, load that file, else it is structure
if (isfield(data,'label') == 0)
    disp(sprintf('Trying to load data file %s',data));
    load(sprintf('exportdata\\%s',data));
    disp(sprintf('Using default sustained interval for analysis'));
    interval = data.sustained;  %must use default interval then
end

%*********** locate su1
it_su = -1;
for zz = 1:size(data.label,2)
        if (strcmp(data.label{zz},'SU1'))
            it_su = zz;
        end
end
if (it_su == -1)
    disp(sprintf('Unable to identify SU1'));
    info = [];
    return;
end

%********* get the single unit spikes for attended and unattended trials
CNUM = max(data.attend);   % 2 or 3 conditions, or more?
                       % condition 1 is attended in RF during 2 of 4 track
                       % condition 2 is unattended in RF during 2 of 4
                       % condition 3 is unattended in RF during 1 of 4
spikes = cell(1,CNUM);   % single unit attended spikes (in sustatined period
sponts = cell(1,CNUM);   % single unit spontaneous before trial
for trial = 1:size(data.trials,2)
   cubo = data.attend(trial);
   spikes{cubo} = [spikes{cubo} ; data.trials{trial}(it_su,:)];
   sponts{cubo} = [sponts{cubo} ; data.spontaneous{trial}(1,:)];
end
%*****************************************************************
spiker = spikes;
sponter = sponts;

%**************************** cintervals *************************
info.cintervals = [12,17,25,35,50,71,100,141,200,283,400];  
choiceint = 7;  % select one countining to plot, default 100ms
for k = 1:size(info.cintervals,2)
  info.cinfo{k} = compute_rate_fano(spikes,interval,info.cintervals(k));
  for kk = 1:size(info.cinfo{k}.meanfano,2)
      info.cfano(kk,k) = info.cinfo{k}.meanfano(kk);
      info.sfano(kk,k) = info.cinfo{k}.semfano(kk);
  end
end
%******************************************************************

%************ plot the spike rasters and mean firing rates, and fano
if (showplot == 1)
   
   %****************
   if (1) % zoom in period of stim in RF
     %if you had to pick a spontaneous period, data.surround might work
     %********* since stimuli are outside the RF during this interval
     xstart = data.surround(1);
     xend = 3100;
   else
     xstart = 1;   %show total period from fixation to final saccade
     xend = size(spiker{1,1},2);
   end
   %****************
   divo = size(sponter{1,1},2)/(size(sponter{1,1},2)+(xend-xstart));
   twid = 0.67;
   awid = divo * twid;
   bwid = (1-divo) * twid;
   
   %****************
   disp('Plotting rastergrams (slow) ...');
   subplot('position',[0.2 0.55 awid 0.4]);
   make_nice_spike_raster(sponter,1,size(sponter{1},2));
   V = axis;
   axis([0 size(sponter{1},2) V(3) V(4)]);
   grid on;
   ylabel('Trial Number');
   title('Spontaneous');
   
   %**************************
   subplot('position',[(0.3+awid) 0.55 bwid 0.4]);
   make_nice_spike_raster(spiker,xstart,xend);
   V = axis;
   axis([xstart xend V(3) V(4)]);
   grid on;
   ylabel('Trial Number');
   title(sprintf('Unit %s rasters',data.name));
   %*************
   subplot('position',[(0.3+awid) 0.35 bwid 0.15]);
   smooth_window = 25;  % give sigma of 12.5ms
   make_nice_mean_raster(spiker,smooth_window,xstart,xend);
   axis tight;
   V = axis;
   maxo = V(4);
   axis([xstart xend V(3) V(4)]);
   plot([interval(1),interval(1)],[V(3),V(4)],'k-'); hold on;
   plot([interval(end),interval(end)],[V(3),V(4)],'k-'); hold on;
   %**** indicate pause period yellow lines
   H = plot([data.pause(1),data.pause(1)],[V(3),V(4)],'y-'); hold on;
   set(H,'Linewidth',2); set(H,'Color',[0.7,0.5,0]);
   H = plot([data.pause(end),data.pause(end)],[V(3),V(4)],'y-'); hold on;
   set(H,'Linewidth',2); set(H,'Color',[0.7,0.5,0]);
   %*******************************
   ylabel('Firing Rate');
   xlabel('Time (ms)');
   
   %**************
   subplot('position',[(0.3+awid) 0.07 bwid 0.2]);
   col = 'rbg';
   y = find( (info.cinfo{choiceint}.timebin{1} >= xstart) & ...
             (info.cinfo{choiceint}.timebin{1} <= xend ) );
   for k = 1:size(spikes,2)
     plot(info.cinfo{choiceint}.timebin{k}(y),info.cinfo{choiceint}.fanobin{k}(y),...
               [col(k),'o-']); hold on;
   end
   plot(info.cinfo{choiceint}.timebin{1},...
         ones(size(info.cinfo{choiceint}.timebin{1})),'k:');
   axis tight;
   V = axis;
   axis([xstart xend V(3) V(4)]);
   plot([interval(1),interval(1)],[V(3),V(4)],'k-'); hold on;
   plot([interval(end),interval(end)],[V(3),V(4)],'k-'); hold on;
   %**** indicate pause period yellow lines
   H = plot([data.pause(1),data.pause(1)],[V(3),V(4)],'y-'); hold on;
   set(H,'Linewidth',2); set(H,'Color',[0.7,0.5,0]);
   H = plot([data.pause(end),data.pause(end)],[V(3),V(4)],'y-'); hold on;
   set(H,'Linewidth',2); set(H,'Color',[0.7,0.5,0]);
   %*******************************
   xlabel(sprintf('Time (%d ms intervals)',info.cintervals(choiceint)));
   ylabel('Fano Factor');
   
   %*************
   subplot('position',[0.2 0.35 awid 0.15]);
   smooth_window = 25;  % give sigma of 12.5ms
   make_nice_mean_raster(sponter,smooth_window,1,size(sponter{1},2));
   V = axis;
   axis([0 size(sponter{1},2) 0 maxo]);
   ylabel('Spontaneous Rate');
   xlabel('Time (ms)');  
   
   %************** plot fanos by counting intervals **************
   subplot('position',[0.07 0.07 (0.3+awid-0.15) 0.2]);
   col = 'rbg';
   for kk = 1:size(info.cinfo{k}.meanfano,2)
      H = semilogx(info.cintervals,info.cfano(kk,:),[col(kk),'.-']); hold on;
      set(H,'Linewidth',2);
      H = semilogx(info.cintervals,(info.cfano(kk,:) + info.sfano(kk,:)),...
                    [col(kk),':']);
      set(H,'Linewidth',1);
      H = semilogx(info.cintervals,(info.cfano(kk,:) - info.sfano(kk,:)),...
                    [col(kk),':']);
      set(H,'Linewidth',1);
   end
   semilogx(info.cintervals,ones(size(info.cintervals)));
   axis tight;
   V = axis;
   axis([ (0.75*info.cintervals(1)) (1.5*info.cintervals(end)) V(3) V(4)]);
   semilogx([info.cintervals(choiceint),info.cintervals(choiceint)],...
            [V(3),V(4)],'k--');   
   xlabel('Log Count Interval');
   ylabel('Fano Factor');
   
end


return;


function info = compute_rate_fano(spmat,interval,cwindow)
%*************
%****** info = compute_rate_fano(spmat,interval,cwindow)
%*** inputs:
%****** spmat - cell matrix {1xC} with spike matrices for each attention
%******         condition, condition 1 - attended in 2 of 4 tracking
%******                    condition 2 - unattended in 2 of 4
%******                    condition 3 - unattended in 1 of 4
%****** interval - the desired interval of analysis to do significance
%******            testing for the rate and fano difference cond 1 vs 2
%****** cwindow - the size of the counting window
%******
%***** outputs:
%***** info.timebin{1xC} - time midpoints of counting intervals
%***** info.smorate{1xC} - gaussian smooth firing rate
%***** info.ratebin{1xC} - mean spike counts in bins
%***** info.varbin{1xC} -  variance spike counts in bins
%***** info.fanobin{1xC} - fano factor, NaN if no spikes in a bin
%*****
%***** info.crange - count windows falling in interval of analysis
%***** info.meanfano(C) - mean fano over range of analysis
%***** info.semfano(C) - sem of fano over range of analysis
%*****
%***** info.airate - AI of rate mod (A-U)/(A+U)
%***** info.ratepval - significance cond 1 vs 2 in firing rate diff
%***** info.aifano - AI of fano mod
%***** info.fanopval - significance cond 1 vs 2 in fano factor diff

  C = size(spmat,2);  % number of attention conditions
  info = [];
  
  %*************** compute time bins for counting windows **********
  timbin = [];
  %align first counting windows so one begins perfectly on interval of
  %***** the analysis specified in input
  tt = interval(1);
  dt = max(1,floor( tt - (cwindow * floor( tt/cwindow )) ));
  bsize = floor( (size(spmat{1,1},2)-dt) / cwindow);  % number of time bins
  begints = [];
  endints = [];
  midints = [];
  ik = 1;
  while (1)       
        base = (ik-1)*cwindow + dt;
        baso = base + cwindow;
        if (baso > size(spmat{1,1},2) )
            break;
        else
            ik = ik + 1;
        end
        midints = [midints (0.5*(base+baso))];
        begints = [begints base];
        endints = [endints baso];
  end
  %*******************************************************************
  %********* store all spike counts in array for later use in Monte-Carlo
  %*********    shuffling to estimate significance of Fano change *******
  trialpercond = size(spmat{1,1},1);  %should be identical # trials per cond
  trialtotal = C * trialpercond;
  cints = size(begints,2);
  storage = zeros(trialtotal,cints);
    
  %********* for each attention condition, compute rate, var, and fano in
  %*********   the specified sized count windows (cwindow) parameter
  for k = 1:C  % for attention condition
     
     startrial = (k-1)*trialpercond; 
     numtrials(k) = size(spmat{1,k},1); 
     rstatbin = [];  %mean and variance per counting interval
     
     for ik = 1:size(begints,2)   % for each counting interval
        
        %***************** 
        base = begints(ik);
        baso = endints(ik);  
        %**** get spike counts of intervals
        for jk = 1:numtrials(k)  
          tval = sum(sum( spmat{1,k}(jk,base:baso) ));  % raw spike counts
          storage( ((k-1)*trialpercond) + jk , ik ) = tval;
        end
        %******* compute mean and variance of counts
        rstatbin(1,ik) = mean( storage((startrial+1):(startrial+numtrials(k)),ik) );
        rstatbin(2,ik) = var( storage((startrial+1):(startrial+numtrials(k)),ik) );
        if (rstatbin(1,ik) > 0)
            rstatbin(3,ik) = (rstatbin(2,ik) / rstatbin(1,ik));  % Fano Factor
        else
            rstatbin(3,ik) = NaN;  % if no spike counts, Fano not defined
        end
     end
     
     %************* store results per attention condition C********
     info.timebin{k} = midints;  % time centers of count intervals in ms
     info.smorate{k} = gauss_smooth(mean( spmat{1,k}),cwindow)*1000;  %smoothed rate
     info.ratebin{k} = rstatbin(1,:);
     info.varbin{k} = rstatbin(2,:);
     info.fanobin{k} = rstatbin(3,:);
  end

  %******* run Monte-Carlo to randomly permute trials between attention
  %******** conditions to assess the significance of the Fano Factor effects
  crange = find( (midints > interval(1)) & (midints < interval(end)) );
  info.crange = crange;
  %*****************
  atto = nanmean( info.fanobin{1}(crange) );
  utto = nanmean( info.fanobin{2}(crange) );
  aival = (atto-utto)/(atto + utto);  % AI value of Fano change in desired 
                                      % interval of analysis
  %*******************************************
  disp(sprintf('Computing Monte-Carlo Its on Counting Interval %d',cwindow));
  %**********
  XA = 1:numtrials(1);
  XB = (numtrials(1)+1):(numtrials(1)+numtrials(2));
  N = 10000;
  aivals = zeros(1,N);
  for ii = 1:N
      %***************
      ax = randperm( (numtrials(1)+numtrials(2)) );
      XXA = ax(XA);  % randomly assign trials as attended cond 1
      XXB = ax(XB);  % randomly assign trials as unattended cond 2
      %************* compute ai value of mean fano for random assign
      amen = mean( storage(XXA,crange) );
      avar = var( storage(XXA,crange) );
      y = find( amen > 0);
      ffa = avar(y) ./ amen(y);  % all non-NAN (var/mean) ratios
      %*********************
      umen = mean( storage(XXB,crange) );
      uvar = var( storage(XXB,crange) );
      y = find( umen > 0);
      ffu = uvar(y) ./ umen(y);
      %*********************
      aa = mean(ffa);
      uu = mean(ffu);
      aivals(ii) = (aa-uu) / (aa+uu);
  end
  y = find( abs(aivals) >= abs(aival));   % fraction of values > observed
  %**************************
  info.aifano = aival;
  info.pfanoval = (size(y,2)/N);

  %********* significant effect of firing rate in analysis window?
  %***** compare 1st and 2nd conditions, attended and unattended 2 of 4
  attrate = (mean(spmat{1,1}(:,interval)'))';
  uttrate = (mean(spmat{1,2}(:,interval)'))';
  %*****************************************
  aa = mean(attrate);
  uu = mean(uttrate);
  info.airate = (aa-uu) / (aa + uu);
  info.prateval = ranksum(attrate,uttrate);  %simple ranksum over trials (rate significance)
  %*******************************************

  %***************************
  for k = 1:C
      info.meanfano(k) = nanmean( info.fanobin{k}(crange) );
      info.semfano(k) = nanstd( info.fanobin{k}(crange) )/sqrt(size(crange,2));
  end
  
  disp(sprintf('CountWin %d, AI rate %5.3f(p=%6.4f), AI fano %5.3f(p=%6.4f)',...
                 cwindow,info.airate,info.prateval,info.aifano,info.pfanoval));

return;

%************** below here are some boring graphics routines to plot things

%****************** function to make a nice raster plot ****************
function maxo = make_nice_mean_raster(spmat,smooth_window,xstart,xend)
%*********** spmat1 and spmat2 are spike matrices of two conditions you wish to compare
%*********** smooth_window ... gaussian smoothing in millisecs
numconds = size(spmat,2);
if (numconds==2)
    colo = [[1,0,0];[0,0,1]];
else
    colo = [[1,0,0];[0,0,1];[0,1,0];[1,1,0]];
end

maxo = 0;
for k = 1:numconds
  spud = spmat{1,k};  
  numtrials(k) = size(spud,1);
  smorate = gauss_smooth(sum( spud(1:numtrials(k),xstart:xend))/....
                      numtrials(k),smooth_window)*1000;
  H = plot(xstart:xend,smorate,'k-'); hold on;
  set(H,'Color',colo(k,:));
  
  if (max(smorate) > maxo)
      maxo = max(smorate);
  end
end


return;


%******************* make a rastergram of the actual spikes on each trial
function make_nice_spike_raster(spmat,xstart,xend)

col = 'rbg';
colo = [[1,1,1];[0.5,0,0];[0,0,0.5];[0,0.5,0]];
colormap(colo);

totspike = [];
for cubo = 1:size(spmat,2)
    totspike = [totspike; (cubo*spmat{1,cubo}(:,xstart:xend)) ];
    totspike = [totspike; zeros(2,size(totspike,2))];
end
totspike = totspike + 1;
imagesc(xstart:xend,(0:(size(totspike,1))+1),totspike); hold on;
%*************
it = 0.5;
for cubo = 1:size(spmat,2)
    plot([xstart,xend],[it,it],[col(cubo),'-']);
    it = it + 1 + size(spmat{1,cubo},1);
    plot([xstart,xend],[it,it],[col(cubo),'-']);
    it = it + 1;
end

return;

%**************************************************************
function output = gauss_smooth(input, window)
% Smoothing function: 
% output = smooth(input, window)
% "Window" is the total kernel width.
% Input array must be one-dimensional.

input_dims = ndims(input);
input_size = size(input);
if input_dims > 2 | min(input_size) > 1,
   disp('Input array is too large.');
   return
end

if input_size(2) > input_size(1),
   input = input';
   toggle_dims = 1;
else
   toggle_dims = 0;
end

if window/2 ~= round(window/2),
   window = window + 1;
end
halfwin = window/2;

input_length = length(input);
%********* gauss window +/- 1 sigma
x = -halfwin:1:halfwin;
kernel = exp(-x.^2/(window/2)^2);
kernel = kernel/sum(kernel);

padded(halfwin+1:input_length+halfwin) = input;
padded(1:halfwin) = ones(halfwin, 1)*input(1);
padded(length(padded)+1:length(padded)+halfwin) = ones(halfwin, 1)*input(input_length);

output = conv(padded, kernel);
output = output(window:input_length+window-1);

if toggle_dims == 1,
   output = output';
end

return;