%% Import data from spreadsheet
% Regression models for USMCA
%
clear, clc
%% Setup the Import Options and import the data
opts = spreadsheetImportOptions("NumVariables", 21);


% Specify sheet and range
opts.Sheet = "Regression";
opts.DataRange = "A2:U95322";

% Specify column names and types
opts.VariableNames = ["Wind", "Solar", "Geothermal", "Baseload", "DSM", "SmartElecTransport", "P2H", "Storage", "FossilIndustry", "SynthGas", "GridStability", "VarCFwindOn", "VarCFwindOff","VarCFsolarPV","VarCFchp","VarCFnuclear","VarCFpp", "VarCFpp2","V2Glosses", "RockbedLosses", "RockbedSteam"];
%opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
opts.VariableTypes = repmat("double",size(opts.VariableNames,2),1);
% Import the data
currentFolder = pwd;
pathfile = strcat(currentFolder,'\USMCA.xlsx');
tbl = readtable(pathfile, opts, "UseExcel", false);

%% Clear temporary variables
clear opts


%% Convert to output type
Wind = tbl.Wind;
Solar = tbl.Solar;
Geothermal = tbl.Geothermal;
Baseload = tbl.Baseload;
DSM = tbl.DSM;
SmartElecTransport = tbl.SmartElecTransport;
P2H = tbl.P2H;
Storage = tbl.Storage;
FossilIndustry = tbl.FossilIndustry;
SynthGas = tbl.SynthGas;
GridStability = tbl.GridStability;

VarCFwindOn = tbl.VarCFwindOn;
VarCFwindOff = tbl.VarCFwindOff;
VarCFsolarPV = tbl.VarCFsolarPV;
VarCFchp = tbl.VarCFchp;
VarCFnuclear = tbl.VarCFnuclear;
VarCFpp = tbl.VarCFpp;
VarCFpp2 = tbl.VarCFpp2;
V2Glosses = tbl.V2Glosses;
RockbedLosses = tbl.RockbedLosses;
RockbedSteam = tbl.RockbedSteam;

name_variables = tbl.Properties.VariableNames;

%% Clear temporary variables
clear opts tbl

%% Automatic Linear regression model --> VarCFwind
X_all = [Wind, Solar, Geothermal, Baseload, DSM, SmartElecTransport, P2H, Storage, FossilIndustry, SynthGas, GridStability];
criteria_value = 0.05; % To control the while loop (criteria: when all the correlation values are below 0.05, brake the loop.
iter_y = 1; % to control the name of the sheet in the Excel file

for output = [VarCFwindOn, VarCFwindOff, VarCFsolarPV, VarCFchp, VarCFnuclear, VarCFpp, VarCFpp2, V2Glosses, RockbedLosses, RockbedSteam]
    size_inputs = size(X_all);
    max_val_R = 1;
    A_corr = [Wind, Solar, Geothermal, Baseload, DSM, SmartElecTransport, P2H, Storage, FossilIndustry, SynthGas, GridStability, output];

    Intercept = 0; % Independent term of the equation in the regression model
    dependent_factors = zeros(1,size_inputs(2)); % To save the dependent terms of the regression model, one by input took into account
    n = 1;
    names_result = {}; % To save the names of the inputs selected in correct order

    while max_val_R > criteria_value
        % Matrix of correlation (R)
        R = corrcoef(A_corr);
        max_val_R = max(abs(R(size_inputs(2)+1,1:size_inputs(2))));
        [row_max,c_max] = find(abs(R) == max_val_R);
        names_result(n,1) = name_variables(row_max(2));

        % Regression model
        mdl = fitlm(X_all(:,row_max(2)), output);
        Residuals = table2array(mdl.Residuals(:,1)); % Raw Residuals are selected
        Intercept = Intercept + mdl.Coefficients.Estimate(1);
        dependent_factors(n) = mdl.Coefficients.Estimate(2);

        % New correlation matrix and update variables
        A_corr = [Wind, Solar, Geothermal, Baseload, DSM, SmartElecTransport, P2H, Storage, FossilIndustry, SynthGas, GridStability, Residuals];
        output = Residuals;
        R = corrcoef(A_corr);
        max_val_R = max(abs(R(size_inputs(2)+1,1:size_inputs(2))));
        n = n+1;
    end

    % Write data in the Excel file
    filename = 'USMCA_MatlabResults_Autom.xlsx';

    Sheet_name = char(name_variables(size_inputs(2)+iter_y));
    xlswrite(filename, {'Independent factor'}, Sheet_name,'A1')
    xlswrite(filename, Intercept, Sheet_name,'B1')

    TableResults = table(names_result, transpose(dependent_factors(1:size(names_result))));
    writetable(TableResults, filename, 'Sheet', Sheet_name, 'Range', 'A3')
    iter_y = iter_y + 1;
    
    % Save figure
    fig = plotResiduals(mdl);
    xlabel(strcat('Residuals of' ,{' '}, Sheet_name));
    ylabel('Probability of occurrence, normalized')  % The area of each bar is the relative number of observations. The sum of the bar areas is equal to 1.
    saveas(fig, Sheet_name, 'png');
    
    TableResults = table(names_result, transpose(dependent_factors(1:size(names_result))));
    SummaryTable = varfun(@sum,TableResults,'InputVariables','Var2', 'GroupingVariables','names_result');
    writetable(TableResults, filename, 'Sheet', Sheet_name, 'Range', 'A3')
    writetable(SummaryTable, filename, 'Sheet', Sheet_name, 'Range', 'G3')
end




% % Seventh iteration --> V2Gstorage
% mdl_seventh = fitlm(V2Gstorage,Residuals_sixth)
% Residuals_seventh = table2array(mdl_seventh.Residuals(:,1));
% plotResiduals(mdl_seventh,'probability')
% plotResiduals(mdl_seventh)
% % Show the new correlation matrix
% A_corr_seventh = [Wind, SolarPV, DamHydro, Baseload, V2Gstorage, P2H, Storage, FlexibleDemand, FossilIndustry, SynthGas, VarCFwind, Residuals_seventh];
% R_seventh = corrcoef(A_corr_seventh)


