% PARSE_HUSBANDRY produces the plots necessary to analyze our husbandry setup.
%
% Blanchoud Group, UNIFR
% Simon Blanchoud
% 01/12/2020
function parse_husbandry(recompute)

  % Checkif we need to recompute everything
  if (nargin==0)
    recompute = false;
  else
    recompute = logical(recompute);
  end

  % Location of the reference files
  %db_file = '/home/blanchou/Documents/Projects/Culture/Long-term/export_2020-11-27_14-55-07_UTC.csv';
  db_file = '/home/blanchou/Documents/Manuscripts/Tunicate_husbandry/data/export_strains.csv';
  %water = '/home/blanchou/Documents/Projects/Culture/Long-term/Water_quality.csv';
  water = '/home/blanchou/Documents/Manuscripts/Tunicate_husbandry/data/Water_quality.csv';
  apex_dir = '/home/blanchou/Code/Apex/Apex/';
  tmpfiles = glob('/home/blanchou/Code/TmpData/*.csv');
  strain_id = 131;

  % Either recompute or gather the existing files
  if (recompute || length(tmpfiles) == 0)
    disp('Recomputing Apex data...')
    files = gather_apex_data(apex_dir);
    disp('done!')
  else
    files = tmpfiles;
  end

  % Parse the data from Apexand dispatch it on several panels
  hfig = figure;hold on;
  nsensors = length(files);
  sensor_label = cell(nsensors, 1);
  plot_indx = [1 1 2 3 3 4 4 4 4 4 5];

  % Read each file, extract its name and plot it accordingly
  for i=1:nsensors
    data = csvread(files{i});
    [x, indxs] = sort(data(:,1));
    [dir, sensor_label{i}, ext] = fileparts(files{i});

    subplot(2,3, plot_indx(i)); hold on;
    plot(x, data(indxs,2));
  end

  % Gather the data of our manual measurements of water quality
  fid = fopen(water, 'r');
  if fid > 0
    line = fgetl(fid);
    headers = strsplit(line, ',');
    data = NaN(500, length(headers));

    % Parse the CSV manually
    i = 1;
    line = fgetl(fid);
    while ischar(line)
      columns = strsplit(line, ',', 'COLLAPSEDELIMITERS', false);

      for j=1:length(columns)
        if (~isempty(columns{j}))
          val = columns{j};
          if (j == 1)
            [curr_time, indx] = strptime(val, '%d/%m/%Y');
            if (indx > length(val))
              data(i, j) = mktime(curr_time);
            end
          else
            data(i, j) = str2double(val);
          end
        end
      end

      line = fgetl(fid);
      i = i + 1;
    end

    fclose(fid);
    data = data(1:i-1, :);
  end

  % Plot the various columns on top of the Apex data
  plot_indx2 = [6 6 6 1 3 5 6];
  for i=2:size(data,2)
    subplot(2,3, plot_indx2(i-1)); hold on;
    plot(data(:,1), data(:,i));
  end

  % Set the axis, legends and labels of all panels
  labels = [sensor_label; headers(2:end).'];
  plot_indxs = [plot_indx plot_indx2];
  limx = xlim;
  limy = [0 300; 0 600; 0 50; 0 30; 7 9; 0 6];

  % Set the limits and legends
  for i=1:6
    subplot(2,3,i);
    xlim(limx);
    ylim(limy(i,:));
    legend(labels(plot_indxs==i));
  end

  % Format the axis into readable format
  t = get(gca, 'XTick');
  tdate = cell(size(t));
  for i=1:numel(t)
    tdate{i} = strftime('%m/%y', localtime(t(i)));
  end
  for i=1:6
    hax = subplot(2,3,i);
    set(hax, 'XTickLabel', tdate);
    xlabel('date (mm/yy)');
  end

  % Parse the strains CSV file
  strains = parse_strains_db(db_file);

  % Estimate the time range of our data
  today = mktime(localtime(time()));
  month = 60*60*24*365.25/12;
  [start, nindx] = min(strains(:,5));
  nmonths = ceil((today - start) / month);
  if (nmonths>200)
    disp(strins(nindx,:))
    error(['Some date issue with entry ' nindx ' !']);
  end
  nstrains = max(strains(:,1));

  % Bin the data to count the number of clones per strain over time
  populations = zeros(nmonths, nstrains);
  death = strains(:,6);
  death(strains(:,4) > 0) = today;

  % Count one "bin" for unknown death times
  lifespan = max(ceil((death - strains(:,5))/month), 1);
  ticks = (start-1)+[0:nmonths]*month;

  % Loop over each clone to add its lifespan over the count of its strain
  for i=1:size(strains, 1)
    id = strains(i, 1);
    new = strains(i, 5);

    % Add one count to the proper strain over the lifespan of the clone
    indx = find(new>ticks(1:end-1) & new<=ticks(2:end));
    populations(indx+[0:lifespan(i)-1],id) = populations(indx+[0:lifespan(i)-1],id) + 1; 
  end
  % Don't show ongoing periods
  populations = populations(1:end-1, :);
  ticks = ticks(1:end-2);

  % Display the statistics using stacked bars
  figure;bar(populations, 'stacked');
  hold on;
  plot([0 nmonths], [128 128], 'k');
  xlim([0 nmonths]);

  % Update the axes, labels and strain index
  pdate = cell(size(ticks));
  for i=1:2:numel(ticks)
    pdate{i} = strftime('%m/%y', localtime(ticks(i)));
  end
  set(gca, 'XTick', [1:length(ticks)], 'XTickLabel', pdate);
  xlabel('date (mm/yy)');
  ylabel('clone number');
  title('Husbandry population')
  hbar = colorbar();
  t = get(hbar, 'YTick');
  set(hbar, 'YTickLabel', cellstr(num2str(round(t.'*(nstrains-1))+1)))

  % Parse one single strain (strain_id) to show its lineage
  clones = strains(strains(:,1)==strain_id, :);

  % Infer death dates (today for living clones)
  missings = isnan(clones(:,end));
  clones(missings, end) = clones(missings, end-1)+month;
  clones(missings & clones(:,4)>0, end) = today;

  % Plot the lifespan of each clone, and link it to its parent clone
  figure;hold on;
  for i=1:size(clones, 1)
    plot(clones(i,[end-1 end]), clones(i, [2 2]), 'r', 'LineWidth', 2);
    mom = find(abs(clones(1:i-1, 4))==clones(i,3), 1, 'last');
    if (~isempty(mom))
      plot([clones(i, end-1)-month/12 clones(i, end-1)], [clones(mom, 2) clones(i,2)], 'k');

    % Here we manually fix a few links to unknown slides
    else
      if clones(i, 3) == 31
        plot([clones(i, end-1)-month/12 clones(i, end-1)], [37 clones(i,2)], 'k');
      elseif clones(i, 3) == 96
        plot([clones(i, end-1)-month/12 clones(i, end-1)], [126 clones(i,2)], 'k');
      elseif ~isnan(clones(i, 3))
        disp(clones(i,[1 3]))
      end
    end
  end

  % Update the axes, labels and strain index
  xval = get(gca, 'XTick');
  xdate = cell(size(xval));
  for i=1:numel(xval)
    xdate{i} = strftime('%m/%y', localtime(xval(i)));
  end
  set(gca, 'XTickLabel', xdate);
  xlabel('date (mm/yy)');
  ylabel('clone ID');

  keyboard

  return;
end