% Run automatic calibration on all fibers.

% Clear workspace.
clear;

% Process all fibers.
% Each fiber must have an associated file in both laptop and mt directories.
fibers = {'fiber1', 'fiber2', 'fiber3', 'fiber4', 'fiber5', 'fiber6', 'fiber7', 'fiber8', 'fiber9'}; % list of all fibers to be processed
mt_directory = 'mt_mat_files_new'; % minitweezers data directory
laptop_directory = 'laptop_mat_files'; % laptop (50 kHz) data directory

for fiber = 1:length(fibers)
  fibername = fibers{fiber};
disp(sprintf('Processing fiber "%s"...', fibername));

  % Form filenames.
  mt_filename = sprintf('%s/%s_mt_new.mat', mt_directory, fibername);
  laptop_filename = sprintf('%s/%s.mat', laptop_directory, fibername);

  % Load fibers.
  disp(sprintf('Loading "%s"...', mt_filename));
  mt_file_contents = load(mt_filename);
  mt_fiber = getfield(mt_file_contents, sprintf('%s_mt_new', fibername));
  disp(sprintf('Loading "%s"...', laptop_filename));
  laptop_file_contents = load(laptop_filename);
  laptop_fiber = getfield(laptop_file_contents, fibername);

  % Determine number of traces.
  ntraces = length(mt_fiber);
  disp(sprintf('There are %d traces.', ntraces));

  % Allocate storage for adc-to-force transformation parameters from single fiber.
  % NOTE that these are no longer used.
  a_n = zeros(ntraces,1);
  b_n = zeros(ntraces,1);

  % Allocate storage for mean force and adc readings from aligned regions.
  avgforce_n = zeros(ntraces,1);
  avgadc_n = zeros(ntraces,1);

  % Align minitweezers and laptop traces, extracting average voltage ADC and force readings over aligned regions.
  for trace_index = 1:ntraces
    show_plots = true; % generate plots
%     try
      % Attempt automatic force calibration.
      [a, b, F, avgforce, avgadc] = force_calibration_automatic(mt_fiber{trace_index}, laptop_fiber{trace_index}, show_plots);

      % Store force conversion parameters.
      a_n(trace_index) = a;
      b_n(trace_index) = b;
      
      % Store average force and voltage ADC readings.
      avgforce_n(trace_index) = avgforce;
      avgadc_n(trace_index) = avgadc;

      % Save plot.
      filename = sprintf('plots/%s-%d.eps', fibername, trace_index);
      disp(sprintf('saving plot to "%s"...', filename));
      print('-depsc', filename);
      disp('converting to PDF...');
      command = sprintf('epstopdf %s', filename);
      system(command);
      disp('Done.');
%     catch exception
%       % Automatic force calibration failed for this trace.
%       disp(exception);
%       fprintf('WARNING: Processing of trace %d unsuccessful.\n', trace_index);
%     end

    fprintf('\n');
  end

  % Compute statistics over all traces from this fiber.
  a = mean(a_n);
  da = std(a_n) / sqrt(ntraces);
  b = mean(b_n);
  db = std(b_n) / sqrt(ntraces);  

  fprintf('STASTICS OVER ALL TRACES FOR FIBER %s\n', fibername);
  fprintf('a = %f +- %f\n', a, da);
  fprintf('b = %f +- %f\n', b, db);

  % Perform linear regression over average forces.
  [coeff, coeffint,r,rint,stats] = regress(avgforce_n, [avgadc_n ones(ntraces,1)]);
  a = coeff(1);
  b = coeff(2);

  disp(sprintf('LINEAR REGRESSION OVER ALL TRACES FOR FIBER %s', fibername));
  disp(sprintf('a = %f : 95%% confidence interval [%f,%f]', a, coeffint(1,1), coeffint(1,2)));
  disp(sprintf('b = %f : 95%% confidence interval [%f,%f]', b, coeffint(2,1), coeffint(2,2)));

  % Make plot of average voltages and forces.
  clf;
  plot(avgadc_n, avgforce_n, 'k.');
  axis square;
  oldaxis = axis;
  hold on;
  plot([oldaxis(1) oldaxis(2)], [oldaxis(1) oldaxis(2)] * a + b, 'r-');
  xlabel('ADC voltage reading');
  ylabel('force / pN');
  title(fibername); 

  % Save plot.
  filename = sprintf('plots/%s-tracefits.eps', fibername);
  print('-depsc', filename);
  command = sprintf('epstopdf %s', filename);
  system(command);
  
  % Calibrate forces.
  fiber = cell(ntraces,1);
  for trace_index = 1:ntraces
    fiber{trace_index} = a * laptop_fiber{trace_index} + b;
  end

  % Write data.
  filename = sprintf('calibrated/%s.mat', fibername);
  save(filename, 'fiber');
end


