function [ figNo ] = calculateEmissionCorr( info, emissionsTrue2D, ...
                            emissionsRecons2D, periods, figNo, ...
                            visualize, saveStuff, oDir, fs)
%calculateEmissionCorr Calculates the correlation between true and
%reconstructed emissions, on a 8-day and 4-periods aggregation basis.
%
%  [ figNo ] = calculateEmissionCorr( info, emissionsTrue2D, ...
%                            emissionsRecons2D, figNo, visualize, saveStuff)
%  calculates the correlation between true and reconstructed emissions.
%  This is done with weekly (strictly, 8-day periods) aggregated emissions
%  and monthly (strictly, 4 periods or 32 days) emissions.
%
%  INPUTS :
%   info            : struct with some info in it
%   emissionsTrue2D : 3D matrix with emissions in 2D, 3rd dim as time axis
%   emissionsRecons2D: ditto, but reconstructed emissions
%   periods         ; start and end of the periods of interest.
%   figNo           : starting figure number
%   visualize       : plot stuff?
%   saveStuff       : save to file?
%   oDir            : output directory
%   fs              : font size
%
% OUTPUTS :
%   figNo    : figNo others should use
%
%  Jaideep Ray, 6/27/2012
% -------------------------------------------------------------------------

periodList = periods(1) : periods(2) ;

%  ====  Loop over each 8-day week and dump out correlation
fprintf('\n\nCorrelation between true and reconstructed emissions per 8-day week\n');
n = length( (info.buffer + 1) : (info.nPeriods - info.buffer) ) ;
weeklyEmissionCorr = zeros(n, 1) ;
count = 1;
for i = (info.buffer + 1) : (info.nPeriods - info.buffer)
    vtrue   = emissionsTrue2D(:, :, i) ;
    vrecons = emissionsRecons2D(:, :, i) ;
    weeklyEmissionCorr(count) = corr(vtrue(:), vrecons(:)) ;
    
    fprintf('Correlation for period %d is %f\n', periodList(count), ...
            weeklyEmissionCorr(count)) ;
    count = count + 1;    
end

% ---- Save correlations to disk
if (saveStuff == true)
    str = sprintf('WeeklyEmissionCorr_P%02d-P%02d_data.dat', periods(1), periods(2));
    save([oDir '/' str], 'weeklyEmissionCorr', '-ascii');
end

% ---- Plot the correlation as a function of weeks. Do so even if visualize
% is false
figNo = figNo + 1 ; figure(figNo) ;
periodList = periods(1) : periods(2) ;
plot((periodList+0.5), weeklyEmissionCorr, '-k') ;
hold on ;

% ====== See correlations over 1 month periods

% ---- Do we have enough data to upscale up to a month?
ndata = info.nPeriods - 2*info.buffer ;
if (ndata >= 4) 
    fprintf('\n\nCorrelation between true and reconstructed emissions per 4-week period\n');
    n = length( (info.buffer + 1) : 4 : (info.nPeriods - info.buffer) ) ;
    monthlyEmissionCorr = zeros(n, 1) ;
    count  = -3;
    count2 = 0;
    for i = (info.buffer + 1) : 4 : (info.nPeriods - info.buffer)
        
        % make sure we're inside good periods
        iend = min( [ (i+3), (info.nPeriods - info.buffer) ] ) ;
        
        % If my start is at the end of the set of good periods, I have run out
        % of multi-week periods to average over. Exit the loop
        if ( i == iend )
            break ;
        end
        
        % OK, I have a multi-week periods to average over before judging
        % correlations
        offset = iend - i;
        vtrue   = mean(emissionsTrue2D  (:, :, i : iend), 3) ;
        vrecons = mean(emissionsRecons2D(:, :, i : iend), 3) ;
        count2 = count2 + 1;
        monthlyEmissionCorr(count2) = corr(vtrue(:), vrecons(:)) ;
        
        count = count + 4 ;
        fprintf('Correlation for period [%d - %d] is %f\n', periodList(count), ...
            periodList(count+offset), monthlyEmissionCorr(count2)) ;
    end

    % ---- Save correlations to disk
    if (saveStuff == true)
        str = sprintf('MonthlyEmissionCorr_P%02d-P%02d_data.dat', periods(1), periods(2));
        save([oDir '/' str], 'monthlyEmissionCorr', '-ascii');
    end
    
    % ---- Plot the monthly correlation as a function of months. This can be
    % done only if we have 2 months or more of data
    if ( length(periodList) > 4 )
        midmonths = 0.5 * (periodList(1 : (end-4)) + periodList( (1+4):end)) ;
    else
        midmonths = 0.5 * (periods(1) + periods(2)) ;
    end
    plot(midmonths(1:4:4*count2), monthlyEmissionCorr(1:count2), 'ob');
end

% ---- Decorate plot
title('Correlation between true and estimated emissions', 'FontSize', fs) ; 
xlabel('Period #', 'FontSize', fs) ; 
ylabel('Correlation', 'FontSize', fs) ;
if (ndata >=4) 
legend('8-day (weekly) emissions', '32-day (monthly) emissions', ...
       'Location', 'Best') ; 
else
    legend('8-day (weekly) emissions', 'Location', 'Best') ; 
end
legend boxoff;
set(gca, 'FontSize', fs);   % Increase fonts
hold off;

% ---- Save correlation figure to disk
if (saveStuff == true)
    str = sprintf('EmissionCorrP%02d-P%02d.jpg', periods(1), periods(2));
    print('-djpeg', [oDir '/' str]);
end

end

