function data = dsCalcMetrics(data, varargin)
options=dsCheckOptions(varargin,{...
  'W',1,[],...
  'calc_power_flag',0,[0,1],...
  'auto_gen_test_data_flag',0,[0,1],...
  },false);

%% auto_gen_test_data_flag argin
if options.auto_gen_test_data_flag
  % specific to this function
  varargs = varargin;
  varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
  argin = [{data}, varargs]; % specific to this function
end

if numel(data)>1
  % use dsAnalyze to recursively call dsCalcMetrics on each data set
  data = dsAnalyze(data,@dsCalcMetrics,varargin{:});
  return;
end

time = data.time;
dt = time(2)-time(1);
T = range(time/1000); %[s]
W = options.W;
NW = T*W;

if NW <1.25
  NW = 1.26;
end

if options.calc_power_flag
  data = dsCalcPower(data, 'smooth_factor',1, 'timeBandwidthProduct',NW, varargin{:});
  % data.VARIABLE_Power_SUA.(Pxx,PeakFreq,PeakArea,frequency)
  % data.VARIABLE_Power_MUA.(Pxx,PeakFreq,PeakArea,frequency)
  %   avgs voltage trace then calc power
end

data = dsCalcFRmulti(data, 'bin_size',1, 'bin_shift',1, varargin{:});
% data.VARIABLE_FR_SUA
% data.VARIABLE_FR_MUA
%   sums all spikes in bin then calcs firing rate from bin time length
data = dsCalcFRmulti(data, 'bin_size',1000, 'bin_shift',1000, 'output_suffix','_1sBins' ,varargin{:});
data = dsCalcFRmulti(data, 'bin_size',100, 'bin_shift',100, 'output_suffix','_100msBins' ,varargin{:});

data = dsCalcISI(data, varargin{:});
% data.VARIABLE_ISI_SUA
% data.VARIABLE_ISI_MUA
%   combines all SUA ISIs

data = dsCalcACF(data, varargin{:});
% data.VARIABLE_ACF_SUA
% data.VARIABLE_ACF_MUA
%   adds all spikes for each time point, does gaussian conv, then takes acf

% add acf names to labels to permit power calc
acfTokens = regexp(fieldnames(data),'(.*ACF_MUA)','Tokens');
acfTokens(cellfun(@isempty, acfTokens)) = [];
acfTokens = cellfun(@(x) x{1}, acfTokens);
data.labels = [data.labels, acfTokens'];

%calc power of spike acf for spike power
if options.calc_power_flag
  data = dsCalcPower(data, 'smooth_factor',1, 'timeBandwidthProduct',NW, 'variable','ACF_MUA', 'time_limits',[time(1) min(time(end)-dt, time(min(length(time),1000)-1))], varargin{:});
end

% default variable to process
% if isempty(options.variable)
%   % use first state variable in model
%     options.variable=data.labels{1};
% end

popNames = {data.model.specification.populations.name};
% isiVars = SelectVariables(fieldnames(data),'*_V_ISI*');
% frVars = SelectVariables(fieldnames(data),'*_V_FR*');
% acfVars = SelectVariables(fieldnames(data),'*_V_ACF*');
% sxxVars = SelectVariables(fieldnames(data),'*_V_Power*');

for iPop = 1:length(popNames)
  thisPop = popNames{iPop};
  
  %% NaNs
  thisMetrics.nanV = any(isnan(data.([thisPop '_V'])(:)));
  
  %% Firing Rates
  thisMetrics.suaFR = data.([thisPop '_V_FR_SUA']);
  thisMetrics.muaFR = data.([thisPop '_V_FR_MUA']);
  
  %% Sliding Binned Firing Rates
  thisMetrics.suaFR1sBins = data.([thisPop '_V_FR_SUA_1sBins']);
  thisMetrics.muaFR1sBins = data.([thisPop '_V_FR_MUA_1sBins']);
  thisMetrics.suaFR100msBins = data.([thisPop '_V_FR_SUA_100msBins']);
  thisMetrics.muaFR100msBins = data.([thisPop '_V_FR_MUA_100msBins']);
  
  %% Synchronization Between SUA and Population
  %compare SUA to MUA for synchronization
  thisMetrics.smUnitFRPropDiff = (abs( (thisMetrics.muaFR - thisMetrics.suaFR) / thisMetrics.muaFR )); %spiking
  if options.calc_power_flag
    thisMetrics.smUnitPeakSxxFreqPropDiff = (abs( (data.([thisPop '_V_Power_MUA']).PeakFreq - data.([thisPop '_V_Power_SUA']).PeakFreq) / data.([thisPop '_V_Power_MUA']).PeakFreq )); %subthreshold
  end
  
  %% Bursting
  %check if bimodal ISI (bursting) and if low ISIs
  muaISI = data.([thisPop '_V_ISI_MUA']);
  if ~isempty(muaISI) && length(muaISI)>1
    try
      umD = gmdistribution.fit(muaISI, 1); % single component, unimodal
      bmD = gmdistribution.fit(muaISI, 2); % two components, bimodal
      burstingISIbool(1) = (umD.AIC > bmD.AIC) && any(bmD.mu<=5); % mu <= 5 ms = 200 Hz for bursting
    catch %likely if not enough ISIs
      burstingISIbool(1) = false;
    end
    
    aicOld = inf;
    aicNew = inf;
    nGaussians = 0;
    while aicNew <= aicOld
      try
        aicOld = aicNew;
        nGaussians = nGaussians + 1;
        thisD = gmdistribution.fit(muaISI, nGaussians);
        aicNew = thisD.AIC;
      catch
        break
      end
    end
    
    try
      finalD = thisD;
      burstingISIbool(2) = any(finalD.mu<=5); % mu <= 5 ms = 200 Hz for bursting

      counts = histc(muaISI, [-inf, 5, inf]);
      counts(end) = []; %remove check for count=inf
      burstProp = min(1, counts(1)/counts(2)); % isi < 5ms = 200 Hz for bursting
      burstingISIbool(3) = (burstProp > .3);
    catch
      burstingISIbool = false;
    end
    
  else
    burstingISIbool = false;
  end
  
  thisMetrics.burstingISIbool = any(burstingISIbool); %if any burst measure is true
  
  if ~exist('umD', 'var')
    umD.mu = 0;
    umD.Sigma = 0;
  end
  
  %% Silence
  timeRange = range(time);
  nCells = size(data.([thisPop '_V']),2);
  muaSilentISIs = [];
  suaPropSilence = nan(1, nCells);
  for iCell = 1:length(data.([thisPop '_V_spike_times']))
    thisSpikeTimes = data.([thisPop '_V_spike_times']){iCell};
    thisSpikeTimes = [1; thisSpikeTimes; time(end)];
    %       [~,thisSpikesInd] = intersect(time, thisSpikeTimes);
    %       thisSpikes = zeros(length(time),1);
    %       thisSpikes(thisSpikesInd) = 1;
    %       thisSpikes([1,end]) = 1;
    %       muaSilentISIs = [muaSilentISIs, diff(find(thisSpikes))]*dt;
    thisISIs = diff(thisSpikeTimes);
    suaPropSilence(iCell) = sum(thisISIs(thisISIs > umD.mu + umD.Sigma))/timeRange;
    muaSilentISIs = [muaSilentISIs; thisISIs];
  end
  %     muaSilentISIs = muaISI(muaISI > umD.mu + umD.Sigma);
  muaSilentISIs = muaSilentISIs(muaSilentISIs > umD.mu + umD.Sigma);
  thisMetrics.muaPropSilence = sum(muaSilentISIs) / (timeRange*nCells);
  thisMetrics.suaPropSilence = suaPropSilence;
  
  %% Median Voltage
  thisMetrics.medianV = median(data.([thisPop '_V']));
  
  %% Subthreshold Oscillations
  %subtract normalized power and normalized power of acf
  
  %% Add to data struct
  data.([thisPop '_metrics']) = thisMetrics;
  data.results{end+1} = [thisPop '_metrics'];
  
end %popnames

%% auto_gen_test_data_flag argout
if options.auto_gen_test_data_flag
  argout = {data}; % specific to this function
  
  dsSaveAutoGenTestData(argin, argout);
end

end %function