% This file relies on the data being imported first: do this by running
% import_all.m, or load the saved result of that script:
load Aimported.mat

%% Calculate means and histograms...
year=1952; % first year with complete data.
for DataTitleC =  {'Alice','CV','WO'}
    DataTitle = DataTitleC{1};
    Data = eval(DataTitle);
    L = length(Data.date_year);
    indices = Data.date_year == year;
    
    dT2 = Data.temp_max(1:(L-1)) - Data.temp_min(2:L);
    cloud = (Data.cloud_cover_9am(1:(L-1)) + Data.cloud_cover_3pm(2:L))/2;
    humid = (Data.humidity_9am(1:(L-1)) + Data.humidity_3pm(2:L))/2;
    
    p = sprintf('hist %s %d',DataTitle,year);          
    figure(1);
    plotHist(1, indices, 'temp min', Data.temp_min);
    plotHist(2, indices, 'temp max', Data.temp_max);
    plotHist(3, indices, 'dT1', Data.temp_max - Data.temp_min);
    plotHist(4, indices, 'dT2', dT2); 
    plotHist(5, indices, 'C', cloud); 
    plotHist(6, indices, 'H', humid); 
    set(gcf,'PaperPositionMode', 'manual', ...
        'PaperType', 'A4', ...
        'PaperUnits','normalized', ...
        'Paperposition',[0  0 1 .5]); % take up half page
    % Save 
    saveas(1, strcat('./images-pdf/',p),'pdf');
    saveas(1, strcat('./images-png/',p),'png');
end;

%% Calculate correlations
fprintf('\nCorrelations:\n\n');
for DataTitleC =  {'Alice','CV','WO'}
    DataTitle = DataTitleC{1};
    Data = eval(DataTitle);
    L = length(Data.date_year);

    % filter to selected year
    indices = Data.date_year == year;
    
    p = sprintf('corr %s %d ',DataTitle,year);
    % calculate derived quantities
    dT2 = Data.temp_max(1:(L-1)) - Data.temp_min(2:L);
    cloud = (Data.cloud_cover_9am(1:(L-1)) + Data.cloud_cover_3pm(2:L))/2;
    humid = (Data.humidity_9am(1:(L-1)) + Data.humidity_3pm(2:L))/2;
    figure(1); 
    % calculate correlations and plot subplots..
    plotCorr([4 2 1], indices, p, 'temp max', 'temp min', Data.temp_max, Data.temp_min);
    plotCorr([4 2 2], indices, p, 'prev temp max ', 'temp min', Data.temp_max(1:(L-1)), Data.temp_min(2:L));
    plotCorr([4 2 3], indices, p, 'dT1', 'hod', Data.temp_max-Data.temp_min, Data.hod);
    plotCorr([4 2 4], indices, p, 'dT2', 'hod', dT2, Data.hod);
    plotCorr([4 2 5], indices, p, 'dT2', 'cloud', dT2, cloud);
    plotCorr([4 2 6], indices, p, 'dT2', 'humid', dT2, humid);
    plotCorr([4 1 4], indices, p, 'humid', 'cloud', humid, cloud);
    % format and save graph
    set(gcf,'PaperPositionMode', 'manual', ...
        'PaperType', 'A4', ...
        'PaperUnits','normalized', ...
        'Paperposition',[0  0 1 1]); % take up whole page
    saveas(1, sprintf('./images-pdf/%s.pdf',p),'pdf');
    saveas(1, sprintf('./images-png/%s.png',p),'png');
end;

%% Monte Carlo errors
subplot(1,1,1);
% For 1- and 5-year blocks, for each data set...
for numYearsInBlock = [1 5]
for DataTitleC =  {'Alice' 'CV' 'WO'}
    DataTitle = DataTitleC{1};
    Data = eval(DataTitle);    
    L = length(Data.date_year);
    years = 1952:numYearsInBlock:2005;
    ds = zeros(length(years), 2);
    tmins = ds;
    tmaxs = ds;    
    i = 0;
    % For each block of years... 
    for year = years
        i = i + 1;
        % filter to selected period
        indices = (Data.date_year >= year) & (Data.date_year < year + numYearsInBlock);        
        % calculate derived quantities
        dT2 = Data.temp_max(1:(L-1)) - Data.temp_min(2:L);        
        % do the Monte Carlo thing for the given data and period
        tmins(i,:) = MonteCarlo(Data.temp_min, indices, sprintf('%d %-5s temp_min', year, DataTitle));
        tmaxs(i,:) = MonteCarlo(Data.temp_max, indices, sprintf('%d %-5s temp_max', year, DataTitle));
        ds(i,:) = MonteCarlo(dT2, indices(1:length(dT2)), sprintf('%d %-5s   dT2   ', year, DataTitle));
    end;
    % Plot everything
    set(gcf,'PaperPositionMode', 'manual', 'PaperType', 'A4',  'PaperUnits','normalized', ...
        'Paperposition',[0  0 .5 .4]); % take up nearly half of page
    
    subplot(3,1,1);
    errorbar(years, tmins(:,1), tmins(:,2), '-');
    plotfit(years, tmins, 'Min temp');        
    subplot(3,1,2);
    errorbar(years, tmaxs(:,1), tmaxs(:,2), 'r-.');
    plotfit(years, tmaxs, 'Max temp');
    subplot(3,1,3);
    errorbar(years, ds(:,1), ds(:,2), 'k:');
    plotfit(years, ds, 'dT2');
    xlabel('year');
    ti = sprintf('Temperature Variation, %d-year periods: %s', numYearsInBlock, DataTitle);
    saveas(1, strcat('./images-pdf/',ti),'pdf');
    saveas(1, strcat('./images-png/',ti),'png');
end;
end;

%% Fit to climate model
set(gcf,'PaperPositionMode', 'manual', 'PaperType', 'A4', 'PaperUnits','normalized', ...
        'Paperposition',[0  0 1 1]); % take up 1/6 of page
for DataTitleC =  {'Alice' 'CV' 'WO' 'BA' 'BP'}
n = 0;
for numYearsInBlock = [1 2 3 5 10 20]
    n = n + 1;
    figure(1);
    subplot(3,2,n);
    DataTitle = DataTitleC{1};
    Data = eval(DataTitle);    
    L = length(Data.date_year);
    years = 1950:numYearsInBlock:2005;
    % sigmas: matrix of fitted sigma values and confidence interval
    sigmas = zeros(length(years), 3);
    i = 0;    
    % loop over blocks of years, generating fitted values for S0 (filling sigmas matrix)
    for year = years
        i = i + 1;
        % Calculate derived quantities
        dT2 = Data.temp_max(1:(L-1)) - Data.temp_min(2:L);            
        cloud = (Data.cloud_cover_9am(1:(L-1)) + Data.cloud_cover_3pm(2:L))/2;
        humid = (Data.humidity_9am(1:(L-1)) + Data.humidity_3pm(2:L))/2;
        nextTmin = Data.temp_min(2:L);
        % Filter to region of interest
        indices = (Data.date_year >= year) & (Data.date_year < year + numYearsInBlock);   
        indices = indices(1:(L-1));  % remove last value so that we can shift by one        
        nextTmin = nextTmin(indices);
        dates = Data.date(indices);
        % build predictor matrix
        X = zeros(sum(indices),4);        
        X(:,1) = cloud(indices);        
        X(:,2) = humid(indices);
        X(:,3) = Data.temp_max(indices);
        X(:,4) = Data.hod(indices) * 3600; % seconds of darkness
        % Filter NaNs and make other cuts if desired
        okrows = all(~isnan(X),2);
        okrows = ~isnan(nextTmin) & okrows;       
        % removing high-humidity data makes the fits worse!
        %okrows = okrows & cloud(indices) < 5;
        %okrows = okrows & humid(indices) < 50;
        nextTmin = nextTmin(okrows);
        X = X(okrows,:);             
        dates = dates(okrows);
        % make sure we have enough data...
        if sum(okrows) < 100
            continue;
        end;
        % Perform fit
        options = statset('Display', 'off', 'Robust', 'off', 'MaxIter', 100);
        [coeffs,r,J,COVB,mse] = nlinfit(X,nextTmin,@modelfun,[ 1 1 1 -1 -1], options);
        % print result to screen
        fprintf('%6s %d coefficients: C1: %-10.4E , C2: %-10.4E ,  S0: %-10.4E ,  S1: %-10.4E ,  S2: %-10.4E\n', ... 
          DataTitle, year, coeffs(1), coeffs(2), coeffs(3), coeffs(4), coeffs(5));      
        % Calculate confidence intervals
        ci = nlparci(coeffs,r,'jacobian',J);
        % Save result for later plotting
        sigmas(i,:) = [coeffs(3) ci(3,1) ci(3,2)];
        % Compute predictions using fitted coefficients
        nextTmin_pred = modelfun(coeffs,X);
        % plot comparisons for each group of years?                
        if 0 == 1
            figure(2);
            subplot(2,1,1);
            plot(datenum(dates), nextTmin, 'k.-'); hold on;
            plot(datenum(dates), nextTmin_pred, 'r.:');             
            plot(datenum(dates), X(:,4)/3600); % plot hours of darkness
            hold off;    
            legend('Actual T\_min', 'Predicted', 'Hours of darkness');
            title(sprintf('%s %d - %d min temp Prediction Error', DataTitle, year, year + numYearsInBlock));
            subplot(2,1,2)
            errors = nextTmin - nextTmin_pred;

            hist(errors, 100)
            title('Error in prediction')        
            xlabel('Size of error (deg)')
            ylabel('Counts')
            text(5, 50, sprintf('Std. Deviation %f', std(errors)));  
            pause;
        end;
    end;    
    % convert confidence interval to error offsets
    sigma0 = sigmas(:,1);
    lowerErrors = sigmas(:,1) - sigmas(:,2);
    upperErrors = sigmas(:,3) - sigmas(:,1);
    % reomve values with very large variance (just for the plot)
    okrows = abs(lowerErrors) < 2 & abs(upperErrors) < 2;
    % plot Sigma0 with error bars
    figure(1);
    errorbar(years(okrows), sigma0(okrows), lowerErrors(okrows), upperErrors(okrows), 'ko');
    ylabel('S0 (95% conf)');    
    xlim([1945 2010]);
    ylim([1 2.5]);
    % filter out yearblocks where we didn't have enough data (or fit didn't work for whatever reason)...
    okrows = sigma0 > 0 & sigma0 < 100;
    % perform a linear (york) fit to see if S0 is changing: S0 = m * YEAR + c, S standard deviation
    [c,m,Sc,Sm] = york_fit(years(okrows), sigma0(okrows)',ones(1,sum(okrows))*.0001,lowerErrors(okrows)'/2,0.0001);
    fprintf('Fitted S0 gradient: %f +- %f (One standard deviation)\n', m, Sm);
    % plot linear fit
    S0fit = m * years + c;
    hold on; plot(years, S0fit); hold off;    
    title(sprintf('%d-year blocks (S0 gradient = %f +- %f)', numYearsInBlock, m, Sm * 2));        
end; % end NumYearsInBlock loop -- have now plotted data for all block sizes for this data set
% save plot
p = sprintf('S0 fit %s', DataTitle);
saveas(1, strcat('./images-pdf/',p),'pdf');
saveas(1, strcat('./images-png/',p),'png');
end;

%% Some other junk...
% Rainfall per year
subplot(1,1,1)
[years netRainfall] = group(Alice.precipitation, Alice.date_year, @sum);
bar(years, netRainfall); xlabel('year'); ylabel('Net Rainfall (mm)');
title('Alice Springs rainfall');
