%%% processLabels.m
%%% Nigel Ward, May-June 2023
%%% Creates the tables and results in the Reduction technical report
%%% To be run in the directory: reduction
%%%   so it can access the txt files in annotations-cao/ and annotations-nw/
%%% These are tab-delimited files, created from Elan
%%%    via File->Export->tab-delimited-> {begin,end, ss.msec}
%%%    which gives a .txt file where each line has format:
%%%    tier string1, tier string2, start, end, label


function processLabels()
  global ttlTimes nReduLabels
  pragLabs = {'FI', 'PC', 'UC', 'RE', 'PW', 'DP', 'TC', 'TG', 'PF', ...
	      'PO', 'NE'};
%    pragLabs = {'FI', 'PC', 'UC', 'RE', 'PW', 'DP', 'TC', 'TG', 'PF-PP', 'PF-TH' ...
%	      'PO', 'NE'};
  for k = 1:length(pragLabs)
    nReduLabels = 0;
    fprintf('processing %s\n', pragLabs{k});
    processOnePragLabel(pragLabs{k});
  end
  fprintf('global percents at each reduction level: %.2f %.2f %.2f %.2f\n', ...
	  mean(ttlTimes,1) / sum(ttlTimes));
  fprintf('total reduction-labeled regions: %d\n', nReduLabels);
  fprintf('\n\n');
end


function processOnePragLabel(pragLab)
  nReduLabels = 4;
  likelyMaxPragTimes = 1;   % if exceeded no problem, since will append
  global ttlTimes pragTimes nPragValid;
  global currentPragLabel currentFile currentSide
  ttlTimes = zeros(1,nReduLabels);  % total times at each level  
  nPragValid = 0;
  %% times for each occurance of this prag func at each reduction label
  pragTimes = zeros(likelyMaxPragTimes, nReduLabels);  

  files = dir('annotations-nw/*txt');   % original 
%  files = dir('annotations-nw2/*txt');   % after splitting PF into PF-PP and PF-TH
  
  for k = 1:length(files)
    pragTable = readtable(fullfile(files(k).folder, files(k).name));
    reduTable = readtable(fullfile('annotations-cao', files(k).name));
    currentFile = files(k).name;
    fprintf('   %s %d %d ...', currentFile, size(pragTable,1), size(reduTable,1));
    currentPragLabel = pragLab;   % set the global
    currentSide = 'left';
    processSide(pragLab, currentSide, pragTable, reduTable);
    currentSide = 'right';
    processSide(pragLab, currentSide, pragTable, reduTable); 
  end
  fprintf('\n');
  computeStats(pragLab, pragTimes, ttlTimes);
end
  


function computeStats(pragLab, pragTimes, ttlTimes)
  pragTimesPct = pragTimes ./ sum(pragTimes, 2);   % convert from ms to percents
  avgPercentagePragTimes = mean(pragTimesPct, 1);
  fprintf('for %s, average percent time at each reduction level is:  %.2f %.2f %.2f %.2f\n', ...
	  pragLab, avgPercentagePragTimes);
  %% vector summing up reduced and highly reduced
  pctPragTimesReduced = pragTimesPct(:,addOne(2)) + pragTimesPct(:,addOne(3));

  globalAvgReduced = globalAvgRedu(ttlTimes);
  %%fprintf('globalAvgReduced is %.2f\n', globalAvgReduced);

  [skip,p] = ttest(pctPragTimesReduced, globalAvgReduced, ...
		   'Tail',  'right');
  %%fprintf('  ... p-value of being non-randomly reduced is %.5f\n', p);

  %%fprintf('global avg reduction level is %.2f\n', avgReductionLevel(ttlTimes));
  
  avs = avgReductionLevels(pragTimes);
  avs = avs(~isnan(avs));
  gav = avgReductionLevel(ttlTimes);
  
  [skip,p] = ttest(avs, gav, 'Tail',  'right');
  fprintf('  ... avg reduction=%.2f, (over %d regions) vs global average (%.2f) has p-value %.4f\n', ...
	  mean(avs), length(avs), gav, p);
end



function av = avgReductionLevel(quadruple)
  av = -1 +  sum((1:4) .* quadruple) / sum(quadruple);
end


%% each row of pragtimes is a quadruple
function avs = avgReductionLevels(pragTimes)
  nPragTimes = size(pragTimes, 1);
  avs = zeros(nPragTimes, 1);
  for i = 1:nPragTimes
    avs(i) = avgReductionLevel(pragTimes(i,:));
  end
end



function processSide(pragLab, side, pTable, rTable)
  %%fprintf('    processing %s\n', side);
  pragRegions = parsePrag(pTable, side, pragLab);
  [reduRegions, reduLabels] = parseRedu(rTable, side);
  incrementCounts(pragRegions, reduRegions, reduLabels);
end


%%% for each prag region, save times at 0 reduction, at 1 .. 3
function incrementCounts(pragRegions, reduRegions, reduLabels)
  global ttlTimes nReduLabels pragTimes nPragValid
  global currentPragLabel currentFile currentSide   % for debug output only 
  nPragRegions = size(pragRegions,1);
  nReduRegions = size(reduRegions,1);
  for k = 1:nReduRegions
    labIx = addOne(reduLabels(k,1));
    duration = (reduRegions(k,2) - reduRegions(k,1));
    ttlTimes(labIx) = ttlTimes(labIx) + duration;
  end
  nReduLabels = nReduLabels + nReduRegions;
  for i = 1:nPragRegions
    overlapFound = false;
    onePragTimes = zeros(1,4);
    for j = 1:nReduRegions
      overlap = computeOverlap(pragRegions(i,:), reduRegions(j,:));
      if overlap > 0
	overlapFound = true;
	labIx = addOne(reduLabels(j,1));
	onePragTimes(labIx) = onePragTimes(labIx) + overlap;
	%%fprintf('  increased count for reduction level %d by %.2f\n', labIx - 1, overlap);
      end
    end
    if overlapFound
      nPragValid = nPragValid + 1;
      pragTimes(nPragValid,:) = onePragTimes;
      fprintf('note: for label %s at %.1f~%.1f in %s of %s,  durations at levels are %.1f %.1f %.1f %.1f\n', ...
	      currentPragLabel, pragRegions(i,:), currentSide, currentFile, onePragTimes);
      else
	fprintf('warning: no overlapping reduction-labels  found for %.1f~%.1f in %s of %s\n', ...
		pragRegions(i,:), currentSide, currentFile);
    end
  end
end


function ttlPctReduced = globalAvgRedu(ttlTimes)
  ttlTimesPct = ttlTimes / sum(ttlTimes, 2);
  ttlPctReduced = ttlTimesPct(addOne(2)) + ttlTimesPct(addOne(3));
end


%%% convert reduction labels 0,1,2,3 to indices 1,2,3,4
function index = addOne(label)
  index = 1 + label;
end



%% returns a two-column matrix, listing starts and ends for every
%% region that matches all the specified tracks (string match anywhere
%% in the tier name), and the specified pragmatic-function label
%% (exact match)
function regions = parsePrag(table, targetTrack, targetLabel)
  isMatch = true(size(table,1),1);
  targetTrackIndicator = targetTrack(2:end); % to match {Rr}ight,
					     % {Ll}eft
  for i = 1:size(table,1)
    label = table(i,5).Var5{1};
    tier = table(i,1).Var1{1};
    if ~strcmp(label, targetLabel)
      isMatch(i) = false;
    elseif ~contains(tier, targetTrackIndicator)
      isMatch(i) = false;
    end
  end
  starts = table.Var3;
  ends = table.Var4;
  regions = [starts(isMatch), ends(isMatch)];
  %%fprintf('in %s found %d regions with %s\n', targetTrack, size(regions,1), targetLabel ) ;
end


%%% returns a two-column matrix, listing starts and ends for
%%% every region that matches all the specified track,
%%% and a parallel one-column vector with the specified label:
%%% which will be 0, 1, 2 or 3
function [regions, labels] = parseRedu(table, targetTrack)
  isMatch = true(size(table,1),1);
  targetTrackIndicator = targetTrack(2:end); % to match {Rr}ight,
					     % {Ll}eft
  allLabels = table.Var5;
  
  for i = 1:size(table,1)
    tier = table(i,2).Var2{1};
    if ~contains(tier, targetTrackIndicator)
      isMatch(i) = false;
    elseif isnan(allLabels(i))    % a region with no label 
      isMatch(i) = false;
    end
  end
  starts = table.Var3;
  ends = table.Var4;
  regions = [starts(isMatch), ends(isMatch)];
  labels = allLabels(isMatch);
  %% diagnostic, in case I want to fix up the annotations
  isnanLocations = find(isnan(labels));
  if length(isnanLocations) > 0
    fprintf('!warning: ignoring NaNs in labels');
    regions(isnanLocations,:)
  end
end


function duration = computeOverlap(regionA, regionB)
  duration = min(regionA(2), regionB(2)) - ...
	     max(regionA(1), regionB(1));
  if duration < 0
    duration = 0;
  end
  if duration > 0
%    fprintf('found overlap between %.1f %.1f and %.1f %.1f\n', ...
%	    regionA(1), regionA(2), regionB(1), regionB(2));
  end
end



