%
%                      DB_dram.m
%             

  clear all

%
%  Declare the global variables and load the data from Morris and Whitman.  
%

  global Re Pr
  load db_data

  Re = data(:,1);
  Pr = data(:,2);
  Nu = data(:,3);
  n = length(Nu);
  p = 3;

%
% Plot Nu as a function of Re to demonstrate the differing nature of the datasets. 
%

  figure(1) 
  plot(Re,Nu,'x')
  set(gca,'Fontsize',[22]);
  xlabel('Reynolds Number Re')
  ylabel('Nusselt Number Nu')


%
% Input optimal values for this dataset.
%

  q1 = 0.0045;
  q2 = 0.98;
  q3 = 0.4;
  theta0 = [q1;q2;q3];

%
% Construct the sensitivity matrix, Fisher information matrix, covariance matrix, and measurement error estimate.
%

  Nu_q1 = (Re.^(q2)).*(Pr.^(q3));
  Nu_q2 = q1*(Re.^(q2)).*(Pr.^(q3)).*log(Re);
  Nu_q3 = q1*(Re.^(q2)).*(Pr.^(q3)).*log(Pr);
  res = q1*(Re.^(q2)).*(Pr.^(q3)) - Nu;
  sens_mat = [Nu_q1'; Nu_q2'; Nu_q3'];
  Fisher = sens_mat*sens_mat'; 
  sigma2 = (1/(n-p))*res'*res;
  tcov = sigma2*inv(Fisher);

%
% Set the DRAM parameters. We are presently using the inv-chisq with S0 = sigma2 and N0 = 1;  
% We are using the covariance estimate tcov to initiate V and initial estimate sigma2 of the
% variance for the observation error.
%
% The parameter options are: 'name', initial_value, min_value, max_value, prior_mu, prior_sigma, and targetflag 
%

  clear data model options

  data.ydata = Nu;

  params = {
    {'q1',theta0(1),0,2*theta0(1)}
    {'q2',theta0(2),0,2*theta0(2)}
    {'q3',theta0(3),0,2*theta0(3)}};
  model.ssfun = @dbss;
  model.sigma2 = sigma2;
  model.S20 = sigma2;
  model.N0 = 1;
  options.updatesigma = 1;
  options.qcov = tcov;
  
%
% Run DRAM with 1e4 initial burn-in chain elements followed by 1e4
% elements, which are used to construct the posterior distributions. 
%

  options.nsimu = 1e4;
  [results1,chain1,s2chain1] = mcmcrun(model,data,params,options);

  N_total = 1e4;
  options.nsimu = N_total;
  [results,chain,s2chain] = mcmcrun(model,data,params,options,results1);

%
% Plot the results. We plot the chains in Figure 2 to demonstrate that they are burned in
% and pairwise and marginal posterior plots in Figures 3 and 4. The chain for the observation
% errors is in Figure 5.
%

  figure(2)
  mcmcplot(chain,[],results,'chainpanel')

  figure(3)
  mcmcplot(chain,[],results,'pairs')

  figure(4)
  mcmcplot(chain,[],results,'denspanel',2);
  
  figure(5)
  plot(s2chain)

%
% Compute chain statistics
%

  chainstats(chain,results)

%
% Construct chain of length 1000, and compute the uncertainty in Nu at Re = 1000, Pr = 300
%

  clear Nu Re Pr

  i = 1;
  ind = 1;
  N_ind = 10;
  Re = 1000;
  Pr = 300;
  for j = 1:N_total
    if i == N_ind
      param_final(ind,:) = chain(j,:);
      q1 = chain(j,1);
      q2 = chain(j,2);
      q3 = chain(j,3);
      Nu(ind) = q1*(Re^(q2))*(Pr^(q3));
      i = 1;
      ind = ind+1;
    else
      i = i+1;
    end
  end

%
% Plot the histogram and kernel density estimate (KDE) of Nu.
%

  [NU_bandwidth,Nu_density,Nu_mesh,Nu_cdf]=kde(Nu);

  figure(6)
  subplot(2,2,1)
  histogram(Nu) 
  subplot(2,2,2)
  plot(Nu_mesh,Nu_density,'linewidth',3)
  xlabel('Nusselt Number')