function [y res] = calcScoreWLSC(path_pred, partition)
% Calculate score for AVEC2012 word-level sub-challenge
% Uses the timings provided in the turn directory to indicate word
% boundaries. Uses the actual label files, not the audio arff files.
%PRE:  path_pred is path to prediction. partition is either 'test' or
%'devel'.
%POST: y is final score (average xcorr over all sequences and dimensions)
%      res is a struct with detailed results and ground truth

    y = 0;
    N = 32;
    
    tdir = '/Users/mfv/Documents/data/AVEC2012/';
    % -- Participants won't have word timings files, but can calculate it from forced
    % aligned word timings in the turn directory
    fn_t = sprintf('%sbaseline/WLSC/data/MLData_WLSC_wordTimings_%s.mat',tdir, partition);
    load(fn_t); % -- Loads  wordTimings
    
    for i = 1:N
        fprintf(1, 'Processing test session %u/%u\n', i, N);
        res(i).truth.arousal = readTruth(i, 'arousal', partition);
        res(i).truth.expectancy = readTruth(i, 'expectancy', partition);
        res(i).truth.power = readTruth(i, 'power', partition);
        res(i).truth.valence = readTruth(i, 'valence', partition);
        
        res(i).pred.arousal = readPred(path_pred, i, 'arousal', partition);
        res(i).pred.expectancy = readPred(path_pred, i, 'expectancy', partition);
        res(i).pred.power = readPred(path_pred, i, 'power', partition);
        res(i).pred.valence = readPred(path_pred, i, 'valence', partition);
        
        if ~checkPredictionFitsTruth(res(i))
            warning('Prediction and truth have different lengths, coping by modifying prediction\n');
            res(i) = modifyPredictionToFitTruth(res(i));
        end
        
        fprintf(1, '\tCalculating correlations...');
        res(i).xcorr.arousal = calccorr(res(i).truth.arousal, ...
            res(i).pred.arousal);
        res(i).xcorr.expectancy = calccorr(res(i).truth.expectancy, ...
            res(i).pred.expectancy);
        res(i).xcorr.power = calccorr(res(i).truth.power, ...
            res(i).pred.power);
        res(i).xcorr.valence = calccorr(res(i).truth.valence, ...
            res(i).pred.valence);
        fprintf(1, '[DONE]\n\n');
        
        fprintf(1, '\tCalculating RMS...');
        res(i).rms.arousal = calcrms(res(i).truth.arousal, ...
            res(i).pred.arousal);
        res(i).rms.expectancy = calcrms(res(i).truth.expectancy, ...
            res(i).pred.expectancy);
        res(i).rms.power = calcrms(res(i).truth.power, ...
            res(i).pred.power);
        res(i).rms.valence = calcrms(res(i).truth.valence, ...
            res(i).pred.valence);
        fprintf(1, '[DONE]\n\n');
        
        y = y + res(i).xcorr.arousal + res(i).xcorr.expectancy + ...
            res(i).xcorr.power + res(i).xcorr.valence;
    end
    y = abs(y/(4*N));
    
    avgres = calcAverageResultsPerDimension(res);
    saveResults(res, path_pred, partition);
    saveGraphics(res, path_pred, partition);
end


function A = calcAverageResultsPerDimension(R)
% -- Calculate the average xcorr and rms per dimension
    n = length(R);
    
    A.arousal.avgxcorr = 0;
    A.expectancy.avgxcorr = 0;
    A.power.avgxcorr = 0;
    A.valence.avgxcorr = 0;
    A.arousal.avgrms = 0;
    A.expectancy.avgrms = 0;
    A.power.avgrms = 0;
    A.valence.avgrms = 0;
    
    for i=1:n
        A.arousal.avgxcorr = A.arousal.avgxcorr + R(i).xcorr.arousal/n;
        A.expectancy.avgxcorr = A.expectancy.avgxcorr + R(i).xcorr.expectancy/n;
        A.power.avgxcorr = A.power.avgxcorr + R(i).xcorr.power/n;
        A.valence.avgxcorr = A.valence.avgxcorr + R(i).xcorr.valence/n;
        A.arousal.avgrms = A.arousal.avgrms + R(i).rms.arousal/n;
        A.expectancy.avgrms = A.expectancy.avgrms + R(i).rms.expectancy/n;
        A.power.avgrms = A.power.avgrms + R(i).rms.power/n;
        A.valence.avgrms = A.valence.avgrms + R(i).rms.valence/n;
    end
    A.arousal.avgxcorr = abs(A.arousal.avgxcorr);
    A.expectancy.avgxcorr = abs(A.expectancy.avgxcorr);
    A.power.avgxcorr = abs(A.power.avgxcorr);
    A.valence.avgxcorr = abs(A.valence.avgxcorr);
end


function y = readTruth(i, dim, partition)
    fprintf(1, '\tLoading ground truth for %s...', dim);
    ddir = '/Users/mfv/Documents/data/AVEC2012/';   % - Partiticipants will have different base directory
    if strcmp(partition, 'test')
        tdir = [ddir 'private/WLSC/labels/test/']; % - WLSC test labels, participants won't have this
    else
        tdir = [ddir 'public/WLSC/labels/devel/']; % - WLSC test labels
    end
    fn = sprintf('%slabels_wordlevel_%s%.3u_%s.dat', tdir, partition, i, dim);
    y = load(fn, '-ASCII');
    fprintf(1, '[DONE]\n');
end

function y = readPred(p, i, dim, partition)
    fprintf(1, '\tLoading predictions for %s...', dim);
    
    % -- Hardcoded Laurens emulator :/
    if ~isempty(regexp(p, 'delft', 'once')) && strcmp(dim, 'expectancy')
        if isempty(regexp(p, 'Final', 'once'))
            dim = 'expectation';
        end
    end
    
    fn = sprintf('%slabels_wordlevel_%s%.3u_%s.dat', p,  partition, i, dim);
    
    % -- Give AVEC2011 style filenaming a chance...
    if ~exist(fn, 'file')
        fn = sprintf('%slabels_wordlevel_%ss%.3u_%s.dat', p,  partition, i, dim);
    end
    if ~exist(fn, 'file')
        error('Participant using incorrect filenaming: %s\n', fn);
    end
    
    y = load(fn, '-ASCII');
    fprintf(1, '[DONE]\n');
end

function B = checkPredictionFitsTruth(R)
    B1 = numel(R.truth.arousal) == numel(R.pred.arousal);
    B2 = numel(R.truth.expectancy) == numel(R.pred.expectancy);
    B3 = numel(R.truth.power) == numel(R.pred.power);
    B4 = numel(R.truth.valence) == numel(R.pred.valence);
    
    B = B1 && B2 && B3 && B4;
    
end

function c = calccorr(x, y)
    C = corrcoef(x, y);
    c = C(1,2);
end

function R = modifyPredictionToFitTruth(R)
    n_t = numel(R.truth.arousal);
    n_p = numel(R.pred.arousal);
    
    if n_t > n_p
        dif = n_t - n_p;
        fprintf(1, 'Prediction has %u fewer instances than ground truth, appending last prediction to fix this\n', dif);
        R.pred.arousal(end+1:end+dif) = R.pred.arousal(end);
        R.pred.expectancy(end+1:end+dif) = R.pred.expectancy(end);
        R.pred.power(end+1:end+dif) = R.pred.power(end);
        R.pred.valence(end+1:end+dif) = R.pred.valence(end);
    else
        if n_p > n_t
            dif = n_p - n_t;
            fprintf(1, 'Prediction has %u more instances than ground truth, appending last prediction to fix this\n', dif);
            R.pred.arousal = R.pred.arousal(1:n_t);
            R.pred.expectancy = R.pred.expectancy(1:n_t);
            R.pred.power = R.pred.power(1:n_t);
            R.pred.valence = R.pred.valence(1:n_t);
        end
    end
end

function saveResults(R, p, part)
% -- Save the results R in a file on path p to email to participants
    pat = regexp(p, '/');
    name = p(pat(end-1)+1:pat(end)-1);
    fn = [p name '_WLSC_' part '.txt'];
    
    A = calcAverageResultsPerDimension(R);
    fid = fopen(fn, 'w');
    
    fprintf(1, 'WLSC results on %s partition for %s\n\n', part, name);
    fprintf(1, '\t\tCorrelation\tRMS\n');
    fprintf(1, 'Arousal\t\t%1.4f\t\t%1.4f\n', A.arousal.avgxcorr, A.arousal.avgrms);
    fprintf(1, 'Expectancy\t%1.4f\t\t%1.4f\n', A.expectancy.avgxcorr, A.expectancy.avgrms);
    fprintf(1, 'Power\t\t%1.4f\t\t%1.4f\n', A.power.avgxcorr, A.power.avgrms);
    fprintf(1, 'Valence\t\t%1.4f\t\t%1.4f\n', A.valence.avgxcorr, A.valence.avgrms);
    fprintf(1, '-------------------------------------------\n');
    fprintf(1, 'Mean\t\t%1.4f\t\t%1.4f\n', mean([A.arousal.avgxcorr, A.expectancy.avgxcorr ...
        A.power.avgxcorr, A.valence.avgxcorr]), mean([A.arousal.avgrms ...
        A.expectancy.avgrms, A.power.avgrms, A.valence.avgrms]));
    
    fprintf(fid, 'WLSC results on %s partition for %s\n\n', part, name);
    fprintf(fid, '\t\tCorrelation\tRMS\n');
    fprintf(fid, 'Arousal\t\t%1.4f\t\t%1.4f\n', A.arousal.avgxcorr, A.arousal.avgrms);
    fprintf(fid, 'Expectancy\t%1.4f\t\t%1.4f\n', A.expectancy.avgxcorr, A.expectancy.avgrms);
    fprintf(fid, 'Power\t\t%1.4f\t\t%1.4f\n', A.power.avgxcorr, A.power.avgrms);
    fprintf(fid, 'Valence\t\t%1.4f\t\t%1.4f\n', A.valence.avgxcorr, A.valence.avgrms);
    fprintf(fid, '------------------------------------------\n');
    fprintf(fid, 'Mean\t\t%1.4f\t\t%1.4f\n', mean([A.arousal.avgxcorr, A.expectancy.avgxcorr ...
        A.power.avgxcorr, A.valence.avgxcorr]), mean([A.arousal.avgrms ...
        A.expectancy.avgrms, A.power.avgrms, A.valence.avgrms]));
    
    fclose(fid);
end


function r = calcrms(x, y)
    r = sqrt(sum((x-y).^2));
end