%% J5 Parameterization code
clearvars; close all; clc

% load master tables
data_Beijing        = load('Beijing_master_table_2018-2019_hourly_UTC+8.txt');
data_BeijingSAprx   = load('Beijing_SAproxy_master_table_2019Mar_Apr_hourly_UTC+8.txt');
data_Budapest       = load('Budapest_master_table_2018Mar_Apr_hourly_UTC+2.txt');
data_Cyprus         = load('Cyprus_master_table_2018Feb_Mar_hourly_UTC+2.txt');
data_Hyytiala       = load('Hyytiälä_master_table_2016-2019_hourly_UTC+2.txt');
data_HyytialaSAprx  = load('Hyytiälä_SAproxy_master_table_2016-2019_hourly_UTC+2.txt');
data_Manacapuru     = load('Manacapuru_master_table_2014Aug_Oct_hourly_UTC-4.txt');
data_Varrio         = load('Värriö_master_table_2019Apr_Aug_hourly_UTC+2.txt');

%% Combining datasets
data_all = [data_Hyytiala; data_Beijing; data_Varrio; ...
    data_Budapest; data_Cyprus; data_Manacapuru];

% columns and variables (with NaNs): 
% 1	Time 				(hourly)
% 2	J5			    	(cm^-3 s^-1)
% 3	H2SO4 monomer		(molec cm^-3)
% 4	Temperature		    (oC)
% 5	RH		    		(%)
% 6	CS			    	(s^-1)
% 7					    (marker for data from different sites)

% Apply data filter
idx = data_all(:,2) > 1e-5 & data_all(:,3) > 5e3 & data_all(:,5)<=100 & data_all(:,6) > 0;
data_all = data_all(idx,:);

%% resample data for training
y = randsample(length(data_all),fix(length(data_all)*0.75));
idx1 = ismember(1:length(data_all),y)';
idx0 = find(idx1==0);
train = data_all(idx1,:);
test = data_all(idx0,:);

data_train = train;
data_test = test;

%% Parameterization
variables = data_train;
J = variables(:,2); % compare with the trained result J5fit on L58
SA = variables(:,3);
RH = variables(:,5);
CS = variables(:,6);

var_input = [SA,RH,CS];
gs_function = GlobalSearch('NumTrialPoints',5000,'NumStageOnePoints',500);

% Example for model 1 training
k = string(1);
Y_fit = @(p,var_input) p(1).*SA;   
SSe = @(p) sum((log10(J) - log10(Y_fit(p,var_input))).^2);
problem = createOptimProblem('fmincon','x0',[1e-7],...
    'objective',SSe,'lb',[1e-20],'ub',[2]);
coefficients = run(gs_function,problem);
J5fit = coefficients(1).*SA; % compare with the measurement result J on L43

%%%% model trainig result is analyzed by comparing J (L43) with J5fit (L58) 

% Example for model 1 testing
data_test = test;
data_test = [data_test; data_HyytialaSAprx; data_BeijingSAprx];

% Apply data filter for testing dataset
idx = data_test(:,2) > 1e-5 & data_test(:,3) > 5e3 & data_test(:,5)<=100 & data_test(:,6) > 0;
data_test = data_test(idx,:);

variables2 = data_test;
J_test = variables2(:,2); % compare with J5fit_test on L156
SA = variables2(:,3);
RH = variables2(:,5);
CS = variables2(:,6);

% Testing model 1 with the trained coefficients
J5fit_test = coefficients(1).*SA; % compare with J_test on L150

%%%% model testing result is analyzed by comparing J_test (L150) with J5fit_test (L156) 
