function info = rate_normalized_coherence(data,lfplabel,interval,window,showplot)
%********* JUDE MITCHELL (jude@salk.edu):  5/20/2008
%******** function info = rate_normalized_coherence(data,lfplabel,interval,window,showplot)
%**********
%    one example:  info = rate_normalized_coherence('mzir4_u1','LFP2',[],200,1);  % use default interval
%
%                  load 'msiv3_u8';
%                  info = rate_normalized_coherence(data,'LFP1',data.sustained,100,1)
%
%                  load 'mnar4_u2';
%                  info = rate_normalized_coherence(data,'LFP1',data.sustained,100,1);
%
%                  load 'msir27_u3';
%                  info = rate_normalized_coherence(data,'LFP2',data.morepause,100,1);  
%
%***** general: computes the coherence in (window size ms, 100 default) hanning windows 
%*****          hopefully the same as in Womelsdorf and Fries
%*****  3/12/08: it runs coherence analysis off the isolated unit
%*****  in the data input 'Su1' with the labeled LFP channel (lpflabel) 
%****** and plots the computed coherence for attended and unattended trials
%*****  NOTE: data has been sorted of LFP1 is same lfp channel as the
%*****        isolated unit was recorded from, and same for MU1 mulit-unit
%** 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
%******
%******
%********** lfplabel - lfp to do spike-lfp coherence, 1 is same lfp
%********** interval - interval of analysis, 1xN array of timepoints
%********** window - duration of LFP window around each spike for FFT
%********** showplot - to plot out results
%*****
%****** showplot - set this to 1 to see the results plotted out
%*****
%******  output:
%******      info.acoho  - attended coherence [1 x F]
%******      info.ascoho - std from permutes of estimate
%******      info.arcoho  - shuffled trials, attended coherence [1 x F]
%******      info.arscoho - shuffled trials, std from permutes of estimate
%******      info.afreq - frequencies [1 x F]
%******      info.ucoho  - attended coherence [1 x F]
%******      info.uscoho - std from jacknife of estimate
%******      info.urcoho  - shuffled trials, attended coherence [1 x F]
%******      info.usrcoho - shuffled trials, std from jacknife 
%******      info.ufreq - frequencies [1 x F]
%******************
%******** recap:
%******** function info = rate_normalized_coherence(data,lfplabel,interval,window,showplot)
%**********
%    one example:  info = rate_normalized_coherence('mzir4_u1','LFP2',[],200,1);  % use default interval
%
%                  load 'msiv3_u8';
%                  info = rate_normalized_coherence(data,'LFP1',data.sustained,100,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;
it_lfp = -1;
for zz = 1:size(data.label,2)
        if (strcmp(data.label{zz},'SU1'))
            it_su = zz;
        end
        if (strcmp(data.label{zz},lfplabel))
            it_lfp = zz;
        end
end
if ( (it_su == -1) | (it_lfp == -1) )
    disp(sprintf('Unable to identify SU1 or lfp label %s',lfplabel));
    info = [];
    return;
end

%********* get the single unit spikes for attended and unattended trials
CNUM = max(data.attend);   % 2 or 3 conditions, or more?
spikes = cell(1,CNUM);   % single unit attended spikes (in sustatined period
lfp = cell(1,CNUM);
for trial = 1:size(data.trials,2)
   cubo = data.attend(trial);
   spikes{cubo} = [spikes{cubo} ; data.trials{trial}(it_su,:)];
   lfp{cubo} = [lfp{cubo} ; data.trials{trial}(it_lfp,:)];
end
%************ subtract out mean LFP's across trials ********
lfpbef = cell(1,CNUM);
for kk = 1:CNUM
  lfpbef{kk} = lfp{kk};
  ameanlfp = mean(lfp{kk});
  for ii = 1:size(lfp{kk},1)
    lfp{kk}(ii,:) = lfp{kk}(ii,:) - ameanlfp;
    meno = mean(lfp{kk}(ii,:));  % subtract out mean inside trial
    lfp{kk}(ii,:) = lfp{kk}(ii,:) - meno;
  end
end
%*****************************************************************

%************ plot the spike rasters and mean firing rates to sanity check
%************ and also show the local fields and mean local fields 
if (showplot == 1)
   
   spiker = spikes;
   
   %****************
   disp('Plotting rastergrams (slow) ...');
   subplot('position',[0.1 0.4 0.4 0.55]);
   make_nice_spike_raster(spiker);
   V = axis;
   axis([0 size(spikes{1},2) V(3) V(4)]);
   grid on;
   ylabel('Trial Number');
   title(sprintf('Unit %s rasters',data.name));
   %*************
   subplot('position',[0.1 0.1 0.4 0.3]);
   smooth_window = 25;  % give sigma of 12.5ms
   make_nice_mean_raster(spiker,smooth_window);
   V = axis;
   axis([0 size(spikes{1},2) 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;
   ylabel('Firing Rate');
   xlabel('Time (ms)');
   
   spiker = lfpbef;
   %****************
   disp('Plotting lfps per trial (slow) ...');
   subplot('position',[0.57 0.4 0.4 0.55]);
   make_nice_lfp_raster(spiker);
   V = axis;
   axis([0 size(spikes{1},2) V(3) V(4)]);
   grid on;
   ylabel('Trial Number');
   title(sprintf('Unit %s lfp %s',data.name,lfplabel));
   %*************
   subplot('position',[0.57 0.1 0.4 0.3]);
   smooth_window = 2;  % give sigma of 12.5ms
   make_nice_mean_raster(spiker,smooth_window);
   V = axis;
   axis([0 size(spikes{1},2) 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;
   ylabel('LFP');
   xlabel('Time (ms)');
   
end


%**************** now call for the coherence estimates ***********
Win = window;  % use specified hanning tapered windows
MinCount = 200;  % use minimum 200 spikes per permuted estimate
%************************************************
for cubo = 1:CNUM
 disp(sprintf('Computing coherence on group %d trials',cubo));
 [coho{cubo},scoho{cubo},phaso{cubo},rcoho{cubo},srcoho{cubo},freq{cubo}] = ...
    fries_coherence_match(lfp{cubo}(:,interval),spikes{cubo}(:,interval),...
                          Win,MinCount);
end
%*********** store results *********
info.coho = coho;
info.scoho = scoho;
info.phaso = phaso;
info.rcoho = rcoho;
info.srcoho = srcoho;
info.freq = freq;
%***********************************

if (showplot == 1)
   figure;
   colo = 'rbgy';
   
   for ii = 1:CNUM
    H = plot(freq{ii},coho{ii},[colo(ii),'-']); hold on;
    set(H,'Linewidth',2);
    H = plot(freq{ii},coho{ii}+scoho{ii},[colo(ii),'-'],...
             freq{ii},coho{ii}-scoho{ii},[colo(ii),'-']); hold on;
    set(H,'Linewidth',1);
    
    H = plot(freq{ii},rcoho{ii},[colo(ii),'--']); hold on;
    set(H,'Linewidth',1.5);
    H = plot(freq{ii},rcoho{ii}+srcoho{ii},[colo(ii),'--'],...
             freq{ii},rcoho{ii}-srcoho{ii},[colo(ii),'--']); hold on;
    set(H,'Linewidth',1);
   end
   
   ylabel('Coherence');
   xlabel('Frequency (Hz)');
   title(sprintf('%s with %s (dashed shuffled)',data.name,lfplabel));
   
end

return;


function [coho,scoho,phaso,urcoho,srcoho,ifcoher] = fries_coherence_match(speco,spiko,Win,MinCount)
%******* [coho,scoho,urcoho,srcoho,ifcoher] = fries_coherence(speco,spiko,Win)
%****** computes the coherence value using spike-triggered LFP similar to 
%****** that of Fries et al, 2001. in 100ms Hanning windowed segments 
%****** ALSO, to rate normalized, total set of spikes is divided into
%****** independent permutations of size MinCount spikes
%******
%****** inputs:
%******    speco - a MxT array where M is trials and T is time in ms
%******    spiko - a MxT array with same dimensions, 0 or 1 for spikes
%******    Win - a window around which to compute LFP segments 
%******    MinCount - number of spikes to use in Fries comp per permutation
%******             - if set to 0, don't rate normalize, run entire set in
%******               of spikes.
%****** outputs:
%******    coho - 1xF array of coherence per frequencies F
%******    scoho - std of the independent permutes in analysis
%******    urcoho - mean coherence of randomly shuffled trials
%******    srcoho - std of coherence from random shuffles
%******    ifcoher - list of frequency values, 1xF

  %**** compute number of spikes, must be at least 200 to run
  if (MinCount == 0)
      yspikes = find( spiko > 0);
      PN = size(yspikes,1);
      PRR = 1:size(yspikes,1);
      NPERM = 1;
      MinCount = PN;   % include every single spike
  else
      sumo = sum(sum(spiko));
      yspikes = find( spiko > 0 );
      NPERM = floor( sumo / MinCount );   % number of random permutations
      PN = size(yspikes,1);  % returns as a column of numbers
      PRR = randperm(PN);
  end
  
  %************* random permutations for shuffle corrected trials
  RB = 20;
  trialperm = [];
  for i = 1:RB
      trialperm = [trialperm ; randperm(size(speco,1))];  %random permutes of trial order
  end
  
  %*********** prep multi-taper fft ***************
  TT = size(speco,2);
  N = Win;
  if (N>(TT/4))
      N = floor(TT/4);
      disp(sprintf('Window size reduced to 1/4 total interval (%d to %d ms)',Win,N));
  end
  Trials = size(speco,1);
  %*******************
  TTa = 1+floor(Win/2);   %tighten up interval of analysis so Win does not fall outside
  TTb = TT-floor(Win/2);
  taper = hanning(Win,'symmetric')';
  %********** compute range of frequencies for FFT analysis
  pad = 0;
  fpass = [0.003 0.088];   % our LFP has a restricted range due to hardware filter
  Fs = 1;
  nfft=2^(nextpow2(N));
  df=Fs/nfft;
  f=0:df:Fs; % all possible frequencies
  f=f(1:nfft);
  findx=find(f>=fpass(1) & f<=fpass(end));
  f=f(findx);
  %**********************
  
  cohopool = [];     % LFP coherences on spike subsets
  urcohopool = [];   % LFP coherences over random permutes of trials (shuffle corrections)
  phaseopool = [];
  %****************
  for uber = 1:NPERM
      
    %********** do analysis many times on different permutes of same data **    
    spika = zeros(size(spiko));
    aa = 1 + ((uber-1)*MinCount);
    bb = (uber*MinCount);
    spika( yspikes(PRR(aa:bb)) ) = 1;   % select random subset of MinCount spikes
    disp(sprintf('Computing Subset %d of %d (contains %d spikes)',uber,NPERM,MinCount));
    
    %******************
    spikecount = 0;
    lfp_pow = [];
    lfp_sta = [];   % time rep of LFP
    rlfp_sta = cell(1,RB);  % random permutations (shuffle correction)
  
    %********* run through all trials, compute for each spike in the
    %********* subset of spikes the LFP segment, its FFT, power etc...
    for ii = 1:Trials
      for tt = TTa:TTb
         if (spika(ii,tt)>0)
            spikecount = spikecount + 1;
            toa = (tt-floor(Win/2));
            tob = toa + Win - 1;
            lfper = speco(ii,toa:tob) .* taper;
            %*************************
            J1=fft(lfper,nfft);
            J1=J1(1,findx);
            %*************************
            if (spikecount == 1)
                lfp_pow = conj(J1) .* J1;  % power spectrum 
                lfp_sta = J1;  
                %****************
                for rr = 1:RB
                   lfper = speco(trialperm(rr,ii),toa:tob);
                   lfper = lfper .* taper;
                   jj1 = fft(lfper,nfft);
                   jj1 = jj1(1,findx);
                   rlfp_sta{1,rr} = jj1;
                end
            else
                lfp_pow = lfp_pow + conj(J1) .* J1;   %Fries 2001 take average of power spectra
                lfp_sta = lfp_sta + J1;
                %****************
                for rr = 1:RB
                   lfper = speco(trialperm(rr,ii),toa:tob);
                   lfper = lfper .* taper;
                   jj1 = fft(lfper,nfft); 
                   jj1 = jj1(1,findx);
                   rlfp_sta{1,rr} = rlfp_sta{1,rr} + jj1;
                end
            end    
                    
         end
      end
    end  % for all trials    
    %**********************************
  
    %********** compute the coherence value
    lfp_sta = (lfp_sta / spikecount);
    S1 = sqrt(lfp_pow / spikecount);
    S12 = abs(lfp_sta);
    phaseo = lfp_sta ./ abs(lfp_sta);
    coho = (S12 ./ S1);
    
    %********* compute null distribution coherences ****
    %****** as average of several randomly permuted trials
    rcoho = [];
    for rr = 1:RB
      rlfp_sta{rr} = (rlfp_sta{rr}/spikecount);
      jja = abs(rlfp_sta{rr}); 
      cobo = (jja ./ S1);
      rcoho = [rcoho ; cobo];
    end
    urcoho = mean(rcoho);
  
    %********* pool the different permutes for later *********
    cohopool = [cohopool ; coho];
    urcohopool = [urcohopool ; urcoho];
    phaseopool = [phaseopool ; phaseo];
  end
  
  ifcoher = f * 1000;
  if (NPERM < 1)  % require min set of samples to estimate cohere
      disp(sprintf('Ret Nan only %d spikes',spikecount));
      coho = NaN * ones(1,size(f,2));
      scoho = NaN * ones(1,size(f,2));
      urcoho = NaN * ones(1,size(f,2));
      srcoho = NaN * ones(1,size(f,2));
      phaso = NaN * ones(1,size(f,2));
  else
    if (NPERM > 1)
      coho = mean(cohopool);
      scoho = std(cohopool) / sqrt( NPERM );
      urcoho = mean(urcohopool);
      srcoho = std(urcohopool) / sqrt( NPERM );
      phaso = mean(phaseopool);
    else
      coho = cohopool;
      phaso = phaseopool;
      urcoho = urcohopool;
      %******** punt on std if too few spikes - could do a real jacknife
      %******** but leave it for later
      scoho = NaN * ones(1,size(f,2));
      srcoho = NaN * ones(1,size(f,2));
    end
  end
  
return;



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

%****************** function to make a nice raster plot ****************
function make_nice_mean_raster(spmat,smooth_window)
%*********** 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

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

return;


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

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}) ];
end
totspike = totspike + 1;
imagesc(totspike);
V = axis;
axis([V(1) V(2) 0 size(totspike,1)]);

return;

%******************* make a rastergram of the actual spikes on each trial
function make_nice_lfp_raster(spmat)

colo = 'rbgy';

totspike = [];
CNUM = size(spmat,2);
for ii = 1:size(spmat,2)
    totspike = [totspike ; spmat{1,ii}];
end
splito = size(totspike,1)/CNUM;

mino = min(min(totspike)) * 0.5;  %allow some partial overlap
maxo = max(max(totspike)) * 0.5;
for ii = 1:size(totspike,1)
   offset = (size(totspike,1)-ii);
   yy = offset + ((totspike(ii,:)-mino)/(maxo-mino));
   cubo = 1 + floor((ii-1)/splito);
   plot(1:size(totspike,2),yy,[colo(cubo),'-']); hold on;
end
V = axis;
axis([V(1) V(2) 0 size(totspike,1)]);

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;