clear all

% Add MCMC codes to path
current_directory = cd;
addpath(strcat(current_directory,'/mcmcstat'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Define model for parameter estimation 
% SS function and Parameters
ssfun = @coupled_func_Bayesian;

model.ssfun = ssfun;

% global sigR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load DFT_cubic
model.N  = length(Pnum);

scaling = 1e3; %all DFT energy density is in units of MPa using the scaling 1e3

%cubic volume
V0 = 4.1476034E+02; %volume of DFT lattice in Bohr^3
V1 = V0*(5.29177249e-11)^3; %conversion of volume to m^3
deltaE = (E_tot - E_tot(1))/V1*1.60217656535e-19/1e9; %Giga-Joules/m^3 conversion from eV/m^3

%input energy paramters from DFT
data.ydata_E_tot = deltaE*scaling;
data.xdata = Pnum; %polarization computed from DFT using the Berry phase approach
%define initial guess for Landau parameters
A = -0.805*scaling;
B = 3.32*scaling;
tmin = [A
        B];
 
%Landau energy model parameter range and initial guess
params = {
    {'A', tmin(1), -5e6*scaling, 5e6*scaling}
    {'B',  tmin(2), -5e3*scaling,5e8*scaling}
    };

% Example of output from this simulation
% MCMC statistics, nsimu = 2000
% 
%                  mean         std      MC_err         tau      geweke
% ---------------------------------------------------------------------
%          A   -0.80461    0.015049  0.00080518      6.8223      0.9996
%          B     3.3196    0.037967   0.0016675      6.3636      0.9997
% ---------------------------------------------------------------------



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%check model initial output
%----------------------------------------------------------------
% THE NEXT LINES OF CODE CHECK THE PARAMETERS PRIOR TO RUNNING A PARAMETER
% ESTIMATION USING MCMC

[Sigma] = coupled_model_Bayesian(tmin,data.xdata);
ss = coupled_func_Bayesian(tmin,data);
SS = ss/length(data.xdata);  %here is the cost calculation

figure(1)
axes1 = axes('FontSize',20);
plot(Pnum,data.ydata_E_tot,'bo--','Linewidth',3)
hold on
plot(Pnum,Sigma,'r-','Linewidth',3)
hold off
ylabel('Energy (MJ/m^3)','FontSize',20)
xlabel('Polarization (C/m^2)','FontSize',20)
legend('Total energy (DFT)','Energy (continuum)')
legend1 = legend(axes1,'show');
set(legend1,'EdgeColor',[1 1 1],'YColor',[1 1 1],'XColor',[1 1 1]);

variance = sum((data.ydata_E_tot-Sigma').^2)/(length(deltaE-1));

%***************************************************************
% BAYESIAN ANALYSIS
model.S20 = variance;%estimate of variance based on prior simulation
model.N0  = 1;
model.N  = length(Pnum);
options.updatesigma=1;
 
options.nsimu = 1000; %run an initial small set of iterations prior to larger run
[results, chain, s2chain] = mcmcrun(model,data,params,options);

options.nsimu = 20000; %run longer simulations using results from smaller run as 
                      %initial conditions
[results, chain, s2chain] = mcmcrun(model,data,params,options,results);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot results
figure(2)
mcmcplot(chain,[],results,'denspanel',2);
figure(3)
mcmcplot(chain,[],results,'pairs');

figure(4); clf
mcmcplot(chain,[],results.names,'chainpanel')
xlabel('Iterations','Fontsize',14)
ylabel('Parameter value','Fontsize',14)

chainstats(chain,results) %summary of parameter means, standard deviations, etc.

%Define a model function that calls the chains for error propagation
modelfun1 = @(d,th)coupled_model_Bayesian(th,d); % NOTE: the order in which d and th appear is important

%The following lines of code propagate the posterior densities of
%parameters to determine prediction and credible intervals of the Landau
%energy with respect to the DFT energy in light of the Bayesian statistics 
nsample = 500;
out = mcmcpred(results,chain,s2chain,data.xdata,modelfun1,nsample);

figure(5)
modelout = mcmcpredplot(out); %,data
hold on
plot(data.xdata,data.ydata_E_tot,'bo-','linewidth',3)
hold off
xlabel('P_{3} (C/m^2)','Fontsize',14);
ylabel('Energy (MPa)','FontSize',20)
legend('95/% Prediction Interval','95% Confidence Interval','Continuum','DFT','Location','Best')