% Calcutes multiregression using the function regress for all possible
% combinations of each input matrix

clear

dir_root = ['/Users/Colin/Documents/MATLAB/DFlow/'...
    'influences_chap2/transform/'];

% Input independent variable matrices
load(strcat(dir_root,'matrix_lny_x_300.mat'));
load(strcat(dir_root,'matrix_lny_x_180.mat'));

% Input dependent variables
%load(strcat(dir_root,'qar_rate1000s180.mat'));
%load(strcat(dir_root,'qar_rate1000s300.mat'));
load(strcat(dir_root,'ln_y300_matrix.mat'));
load(strcat(dir_root,'ln_y180_matrix.mat'));

% Determine all of the combinations for multiregression
[xs_num,ind_var] = size(matrix_lny_x_300);
a = 0;
item_cell = cell(ind_var-1,1);
for i = 2:ind_var
    a = a + 1;
    item_count(a,1) = nchoosek(ind_var,i);
    item_cell(a,1) = {nchoosek(1:ind_var,i)};
end
total_regress = sum(item_count(:,1));

% Calculate regressions
b300lny_x = zeros(total_regress,ind_var+1);
b180lny_x = zeros(total_regress,ind_var+1);
R2_300lny_x = zeros(total_regress,2);
R2_180lny_x = zeros(total_regress,2);
all_combos = zeros(total_regress,ind_var);
a = 0;
for i = 1:ind_var-1
    i_combos = item_cell{i,1};
    [r,c] = size(i_combos);
    for j = 1:r
        a = a + 1;
        ij_ids = i_combos(j,1:c);
        ij_X300lny_x = matrix_lny_x_300(:,ij_ids);
        ij_X180lny_x = matrix_lny_x_180(:,ij_ids);
        mdl300lny_x = fitlm(ij_X300lny_x,ln_y300_matrix);
        mdl180lny_x = fitlm(ij_X180lny_x,ln_y180_matrix);
        R2_300lny_x(a,2) = mdl300lny_x.Rsquared.Ordinary;
        R2_300lny_x(a,1) = mdl300lny_x.Rsquared.Adjusted;
        R2_180lny_x(a,2) = mdl180lny_x.Rsquared.Ordinary;
        R2_180lny_x(a,1) = mdl180lny_x.Rsquared.Adjusted;
        estim = mdl300lny_x.Coefficients.Estimate;
        b300lny_x(a,1:c+1) = mdl300lny_x.Coefficients.Estimate;
        b180lny_x(a,1:c+1) = mdl180lny_x.Coefficients.Estimate;
        all_combos(a,1:c) = ij_ids;
    end
end

% Order results

[R2_300lny_xsort,I_300lny_x] = sortrows(R2_300lny_x,'descend');
[R2_180lny_xsort,I_180lny_x] = sortrows(R2_180lny_x,'descend');
all_combos_300lny_xsort = all_combos(I_300lny_x,:);
all_combos_180lny_xsort = all_combos(I_180lny_x,:);
b300lny_x_sort = b300lny_x(I_300lny_x,:);
b180lny_x_sort = b180lny_x(I_180lny_x,:);

% 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

best_X300lny_xids = ...
    all_combos_300lny_xsort(1,all_combos_300lny_xsort(1,:)>0);
best_X180lny_xids = ...
    all_combos_180lny_xsort(1,all_combos_180lny_xsort(1,:)>0);
best_X300lny_x = matrix_lny_x_300(:,best_X300lny_xids);
best_X180lny_x = matrix_lny_x_180(:,best_X180lny_xids);
best_mdl300lny_x = fitlm(best_X300lny_x,ln_y300_matrix)
best_mdl180lny_x = fitlm(best_X180lny_x,ln_y180_matrix)

figure(1)
set(gcf,'color','white');
%plot(best_mdl180)
plotResiduals(best_mdl300lny_x,'fitted')
%plotDiagnostics(best_mdl180)

figure(2)
set(gcf,'color','white');
% %plot(best_mdl300)
plotResiduals(best_mdl180lny_x,'fitted')
% %plotDiagnostics(best_mdl300)