% Linear analysis of QAR and channel-floodplain characteristics

clear

%% Load QAR and characteristic data

root_dir = '/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/data/';

% QAR data
load(strcat(root_dir,'dist_qar_prop1000m.mat'));
load(strcat(root_dir,'dist_qar_proportional.mat'));
load(strcat(root_dir,'dist_qar_rate1000m.mat'));
load(strcat(root_dir,'dist_qar_rates.mat'));

% Characteristics data
load(strcat(root_dir,'QARw_char_1000m.mat'));
load(strcat(root_dir,'QARw_char_250m.mat'));

% Load WSE minus FP/Bank Z
load(strcat(root_dir,'WSE_minus_bank.mat'));
load(strcat(root_dir,'WSE_minus_fp.mat'));

WSE_minus_bank = abs(WSE_minus_bank);
WSE_minus_fp = abs(WSE_minus_fp);

%% Plot
f1 = figure(1);
f1.Units = 'inches';
%h.Position = [1 1 11 15];
f1.PaperSize = [8.5 11];
set(f1,'color','white');
r_squared_ordinary = zeros(12,2);
r_squared_adjusted = zeros(12,2);
R_value = zeros(12,2);
for i = 1:14
    subplot(5,3,i)
    ax = gca;
    ax.ColorOrderIndex = 1;
    if i < 13
        plot(QARw_char_1000m(:,i),dist_qar_rate1000m(:,1),'o');
        hold on
        plot(QARw_char_1000m(:,i),dist_qar_rate1000m(:,2),'x');
    elseif i == 13
        plot(WSE_minus_bank(:,1),dist_qar_rate1000m(:,1),'o');
        hold on
        plot(WSE_minus_bank(:,2),dist_qar_rate1000m(:,2),'x');
    elseif i == 14
        plot(WSE_minus_fp(:,1),dist_qar_rate1000m(:,1),'o');
        hold on
        plot(WSE_minus_fp(:,2),dist_qar_rate1000m(:,2),'x');
    end
    lsline
    if i < 13
        mdl = fitlm(QARw_char_1000m(:,i),dist_qar_rate1000m(:,1));
        r_squared_ordinary(i,1) = mdl.Rsquared.Ordinary;
        r_squared_adjusted(i,1) = mdl.Rsquared.Adjusted;
        R = corrcoef(QARw_char_1000m(:,i),dist_qar_rate1000m(:,1));
        R_value(i,1) = R(1,2);
        mdl = fitlm(QARw_char_1000m(:,i),dist_qar_rate1000m(:,2));
        r_squared_ordinary(i,2) = mdl.Rsquared.Ordinary;
        r_squared_adjusted(i,2) = mdl.Rsquared.Adjusted;
        R = corrcoef(QARw_char_1000m(:,i),dist_qar_rate1000m(:,2));
        R_value(i,2) = R(1,2);
    elseif i == 13
        mdl = fitlm(WSE_minus_bank(:,1),dist_qar_rate1000m(:,1));
        r_squared_ordinary(i,1) = mdl.Rsquared.Ordinary;
        r_squared_adjusted(i,1) = mdl.Rsquared.Adjusted;
        mdl = fitlm(WSE_minus_bank(:,2),dist_qar_rate1000m(:,2));
        r_squared_ordinary(i,2) = mdl.Rsquared.Ordinary;
        r_squared_adjusted(i,2) = mdl.Rsquared.Adjusted;
    elseif i == 14
        mdl = fitlm(WSE_minus_fp(:,1),dist_qar_rate1000m(:,1));
        r_squared_ordinary(i,1) = mdl.Rsquared.Ordinary;
        r_squared_adjusted(i,1) = mdl.Rsquared.Adjusted;
        mdl = fitlm(WSE_minus_fp(:,2),dist_qar_rate1000m(:,2));
        r_squared_ordinary(i,2) = mdl.Rsquared.Ordinary;
        r_squared_adjusted(i,2) = mdl.Rsquared.Adjusted;
    end
    ylim([0 0.045]);
    if i == 1
        xlabel('Total Area');
    elseif i == 2
        xlabel('FP Area');
    elseif i == 3
        xlabel('Inset Area');
    elseif i == 4
        xlabel('Island Area');
    elseif i == 5
        xlabel('Channel Area');
    elseif i == 6
        xlabel('FP-CH Z');
    elseif i == 7
        xlabel('BA-CH Z');
    elseif i == 8
        xlabel('FP n');
    elseif i == 9
        xlabel('BA n');
    elseif i == 10
        xlabel('FP width');
    elseif i == 11
        xlabel('CH width');
    elseif i == 12
        xlabel('Slope');
    elseif i == 13
        xlabel('Peak Stage - bank z');
    elseif i == 14
        xlabel('Peak Stage - FP z');
    end
end

%%
% print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
%     'influences_chap2/figures/linear_regress8_5x11'],'-dpdf');

%%
% QARw_char_1000m(:,13) = 1;
% 
% [b1,bint1,r1,rint1,stats1] = ...
%     regress(dist_qar_rate1000m(:,1),QARw_char_1000m);
% [b2,bint2,r2,rint2,stats2] = ...
%     regress(dist_qar_rate1000m(:,2),QARw_char_1000m);
% 
% updated_matrix1 = [QARw_char_1000m(:,[1 6:9 12]) WSE_minus_bank(:,1) ...
%     WSE_minus_fp(:,1)];
% updated_matrix2 = [QARw_char_1000m(:,[1 6:9 12]) WSE_minus_bank(:,2) ...
%     WSE_minus_fp(:,2)];
% 
% mdl1 = fitlm(updated_matrix1,dist_qar_rate1000m(:,1))
% mdl2 = fitlm(updated_matrix2,dist_qar_rate1000m(:,2))
% 
% step_mdl1 = step(mdl1,'NSteps',10)
% step_mdl2 = step(mdl2,'NSteps',10)

% [RANKED1,WEIGHT1] = relieff(updated_matrix1,dist_qar_rate1000m(:,1),10);
% [RANKED2,WEIGHT2] = relieff(updated_matrix2,dist_qar_rate1000m(:,2),10);
    
% Column 1 - Total Area
% Column 2 - FP Area
% Column 3 - Inset Area
% Column 4 - Island Area
% Column 5 - Channel Area
% Column 6 - FP - Channel Elevation
% Column 7 - Bank - Channel Elevation
% Column 8 - FP Roughness
% Column 9 - Bank Roughness
% Column 10 - FP Width
% Column 11 - Channel Width
% Column 12 - Slope



% Adjusted matrix of terms
% Column 1 - Total Area
% Column 2 - FP - Channel Elevation
% Column 3 - Bank - Channel Elevation
% Column 4 - FP Roughness
% Column 5 - Bank Roughness
% Column 6 - Slope
% Column 7 - WSE minus bank
% Column 8 - WSE minus fp
    