function Parameter_fitting

% Function to find optimal fit between UVic and AMOC-emulator AMOC evolution by tuning free parameters and minimize cost function
% Function uses so-called Very Fast Simulated Annealing algorithm (A.M. Lombardi 2015 % doi:10.1038/srep08417)

clear
% Set random parameter
c=clock;
rand('seed',fix(c(6)));

% Create and load constants
create_constantsmat
load('constants.mat');

% Load GCM-based AMOC timeseries
load('Forcing_data/amoc_gcm.mat');

% Load GCM-based global reference temperature
load('Forcing_data/Tglob_ref_gcm.mat');

disp('##########################')
disp(['Running EBM-AMOC fitting procedure for ',modelname,' model'])
disp(['for ',num2str(numraninit),' iterations of ',num2str(maxlength),' experiments long'])
disp(['and with Tstati=',num2str(Tstati),' and uphti=',num2str(uphti)])
disp('##########################')

iter=1;
while (iter <= numraninit || iter==1) % Start iteration loop until maximum number of iterations (numraninit) is reached
    disp('##########################')
    disp(['Itteration: ',num2str(iter)])
    %% Find initial parameter set
    Par = zeros(numpara,maxlength);

    if raninitpar==0
        % Used prescribed initial guess values of parameters
        Par(1,1)=ki;
        Par(2,1)=thaui;
        Par(3,1)=F1i;
        Par(4,1)=F2i;
        Par(5,1)=V(4);
        Par(6,1)=h1i;
        Par(7,1)=h2i;        
    elseif raninitpar==1
        % Randomly perturbe initial guess values of parameters by 'parinci'%
        Par_temp(1)=ki;
        Par_temp(2)=thaui;
        Par_temp(3)=F1i;
        Par_temp(4)=F2i;
        Par_temp(5)=V(4);
        Par_temp(6)=h1i;
        Par_temp(7)=h2i;        
        for parnum=1:numpara
            Par(parnum,1)=Par_temp(parnum)+Par_temp(parnum)*(rand*parinci*2-parinci);
        end
        % Check if new parameter values are within bounds (parameters k, thau and V(4) have to be >0)
        while (Par(1,1)<=0 || Par(2,1)<=0 || Par(5,1)<=0)
            for parnum=1:numpara
                Par(parnum,1)=Par_temp(parnum)*(rand*parinci-parinci/2); % Randomly perturb by +-500%
            end
        end
        disp('Initial parameters found:')
        for parnum=1:numpara
            disp(['Parameter ',num2str(parnum),' value: ',num2str(Par(parnum,1),5)])
            disp(['Parameter ',num2str(parnum),' ratio: ',num2str(Par(parnum,1)./transpose(Par_temp(parnum)))])
        end
    end

    %% Perform 'maxlength' simulations in fitting procedure
    % Initialize arrays
    R2 = zeros(1,maxlength);
    Tstat = zeros(1,maxlength);
    moutputexp = nan(maxlength,(eq_length+trans_length)/ebm_res*amoc_res);
    Toce = nan(maxlength,(eq_length+trans_length)/ebm_res,4);
    Soce = nan(maxlength,(eq_length+trans_length)/ebm_res,4);

    % Run AMOC-emulator for initial parameters (experiment 1)
    exp=1;
    [moutputexp(exp,:),~,~,~,~]=runmodel_coupledsystem(Par(:,exp));
    % Calculate Cost function value (R-squared between UVic and AMOC-emulator for transient part AMOC evolution for initial parameters (in Sv))
    for time = 1:trans_length/ebm_res*amoc_res;
        R2(exp) = R2(exp)+(moutputexp(exp,time+eq_length/ebm_res*amoc_res)/(10^6*yr2sec)-amoc_gcm(time))^2;
    end
    disp(['Experiment: ',num2str(exp), ' with cost function of: ',num2str(R2(exp),5)])
    R2grad=R2gradStopcrit; % Set initial value of R2grad (cost function gradient) equal to stopping criterion
    
    while (exp <maxlength && moutputexp(exp,eq_length/ebm_res*amoc_res) > 10^6*yr2sec && R2(exp)<1e7 && R2grad<=R2gradStopcrit); % Start exp loop, but only if AMOC is not shut down (>1SV) and stopping criterion is not met.
        exp=exp+1;
        disp('---------------------------')
        disp([modelname,' experiment ',num2str(exp),' output:'])            
        % Calculate 'SA-Temperature' value for simulated annealing
        Tstat(exp)=min(Tstati,max(0.01,2.9*10^-4*log10(R2(exp-1))^3.36));
        disp(['Tstat value: ',num2str(Tstat(exp))])            
        % Calculate randomly perturbed parameters
        Par(:,exp)=Par(:,exp-1)+Par(:,exp-1).*(rand(numpara,1).*Tstat(exp).*2-Tstat(exp));
        % Check if parameters are within limits, otherwise calculate new random parameters (parameters k, thau and V(4) have to be >0)
        while (Par(1,exp)<=0 || Par(2,exp)<=0 || Par(5,exp)<=0)
            Par(:,exp)=Par(:,exp-1)+Par(:,exp-1).*(rand(numpara,1).*Tstat(exp).*2-Tstat(exp));
        end
        % Run AMOC-emulator model
        [moutputexp(exp,:),Toce(exp,:,:),Soce(exp,:,:),~,~]=runmodel_coupledsystem(Par(:,exp));
        % Calculate candidate cost function value
        for time = 1:trans_length/ebm_res*amoc_res;
            R2(exp) = R2(exp)+(moutputexp(exp,time+eq_length/ebm_res*amoc_res)/(10^6*yr2sec)-amoc_gcm(time))^2;
        end
        % Calculate acceptance criterion threshold value (uphill climbing criterion)
        upht=Tstat(exp)*uphti/Tstati;
        randval=rand;
        % Determine if candidate cost function becomes new cost function value
        if (R2(exp)>R2(exp-1) && upht<randval)
            disp('New parameter set not accepted, Cost function new > Cost function old')          
            disp(['Not going uphill, threshold of ',num2str(upht), ' exeeded with ',num2str(randval)])
            R2(exp)=R2(exp-1);
            Par(:,exp)=Par(:,exp-1);
        elseif (R2(exp)>R2(exp-1) && upht>randval)
            disp('New parameter set accepted, Cost function new > Cost function old')            
            disp(['Going uphill, threshold of ',num2str(upht), ' not exeeded with ',num2str(randval)])
        else
            disp('New parameter set accepted, Cost function new < Cost function old')
        end
        disp(['R2new: ',num2str(R2(exp),5)])
        % Calcute new Rgrad (cost function gradient) for stopping criterion
        if (exp>=R2gradlength)
            a=polyfit(exp-R2gradlength+1:exp,R2(exp-R2gradlength+1:exp),1);
            R2grad=a(1)*R2gradlength/a(2);
        end
    end
    % Save arrays that are needed to perform future projections when suitable parameter set is found
    if (R2(exp)<R2exept)
        save(['R2_',num2str(iter),'_',modelname,'.mat'],'R2')
        save(['Par_',num2str(iter),'_',modelname,'.mat'],'Par')
        save(['moutputexp_',num2str(iter),'_',modelname,'.mat'],'moutputexp')
        load('tatmhist.mat')
        save(['tatmhist_',num2str(iter),'_',modelname,'.mat'],'tatmhist')
        load('tbndhist.mat')
        save(['tbndhist_',num2str(iter),'_',modelname,'.mat'],'tbndhist')
        load('Tocehist.mat')
        save(['Tocehist_',num2str(iter),'_',modelname,'.mat'],'Tocehist')
        load('Socehist.mat')
        save(['Socehist_',num2str(iter),'_',modelname,'.mat'],'Socehist')        
        iter=iter+1;
        disp('##########################')    
        disp(['Suitable parameter sets found for ',modelname,' model'])        
    else
        disp('AMOC shutdown or R2 too large, iteration repeated')
    end % End of if AMOC not shutdown loop
end % End iter while loop

