%
%               Helmholtz.DRAM.m
%
%
  clear all
  close all
%

  load Helmholtz.txt;

  P = Helmholtz(:,1);
  udata = Helmholtz(:,2);
  xdata = P;

  alpha_1 = -382.7;
  alpha_11 = 760.3;
  alpha_111 = 70.1;
  p = 3;
  n = length(P);

  sigma = 3.6;
  var = sigma^2;

%
% Construct the solution and residuals.
%

  psi = alpha_1*P.^2 + alpha_11*P.^4 + alpha_111*P.^6;
  R = psi - udata;

%
% Construct the sensitivity matrix, observation error estimator, and covariance estimator.
%

  psi_alpha_1 = P.^2;
  psi_alpha_11 = P.^4;
  psi_alpha_111 = P.^6;

  sigma2 = (1/(n-p))*R'*R;
  sensmat = [psi_alpha_1 psi_alpha_11 psi_alpha_111];
  V = sigma2*inv(sensmat'*sensmat);

%
% Construct the xdata, udata and covariance matrix V employed in DRAM.
% Set the options employed in DRAM.
%

  clear data model options

  data.xdata = xdata;
  data.ydata = udata;
  tcov = V;
  tmin = [alpha_1; alpha_11; alpha_111];

  params = {
    {'q1',tmin(1), -600, -100}
    {'q2',tmin(2), 000, 1200}
    {'q3',tmin(3), -100, 500}
  };
  model.ssfun = @Helmholtz_ss;
  model.sigma2 = sigma2;
  options.qcov = tcov;
  options.nsimu = 20000;
  options.updatesigma = 1;

%
% Run DRAM to construct the chains (stored in chain) and measurement
% variance, which is stored in s2chain.
%

  [res,chain,s2chain] = mcmcrun(model,data,params,options);

%
%
  figure(1)
  mcmcplot(chain,[],res,'chainpanel');

  figure(2)
  mcmcplot(chain,[],res,'pairs');

  alpha_1_val = chain(:,1);
  alpha_11_val = chain(:,2);
  alpha_111_val = chain(:,3);

%
% Compute and plot the prediction intervals.
%

  nsample = 5000;
  out = mcmcpred(res,chain,s2chain,data.xdata,@energy_eval,nsample)
  mcmcpredplot_custom(out)
  hold on
  plot(P,udata,'xb')
  hold off