%
%                               Helmholtz_Example
%
%
  clear all
%  close all
%

  load Helmholtz.txt;

  P = Helmholtz(:,1);
  obs = Helmholtz(:,2);

  alpha_1 = -382.7;
  alpha_11 = 760.3;
  alpha_111 = 155.5;
  p = 3;
  n = length(P);
  Pend = 0.8;


%
% Construct the solution and residuals.
%

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

%
% 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);
  Mu = [alpha_1 alpha_11 alpha_111];

  
  V2 = V(1:2,1:2);
  Mu2 = [alpha_1 alpha_11];

  N = 1e+5; 
  sample_3_param = mvnrnd(Mu,V,N);
  sample_2_param = mvnrnd(Mu2,V2,N);
  s1 = sample_3_param(:,1);
  s2 = sample_3_param(:,2);
  s3 = sample_3_param(:,3);
  mu1 = alpha_1;
  mu2 = alpha_11;
  mu3 = alpha_111;

%
%  Compute constants used to define integrated response.
%

  c1 = (Pend^3)/3;
  c2 = (Pend^5)/5;
  c3 = (Pend^7)/7; 

  f0 = c1*alpha_1*ones(N,1) + c2*alpha_11*ones(N,1) + c3*alpha_111*ones(N,1);

  response_3_param = c1*sample_3_param(:,1) + c2*sample_3_param(:,2) + c3*sample_3_param(:,3);
  response_2_param = c1*sample_2_param(:,1) + c2*sample_2_param(:,2) + c3*alpha_111;

  energy = -20:.01:0;
  dens3 = ksdensity(response_3_param,energy);
  dens2 = ksdensity(response_2_param,energy);

  energy2 = -30:.01:0;
  %energy2 = -80:0.1:30;
  per = 0.10;
  param1 = (alpha_1 - per*abs(alpha_1)) +  2*per*abs(alpha_1)*rand(N,1);
  param2 = (alpha_11 - per*abs(alpha_11)) +  2*per*abs(alpha_11)*rand(N,1);
  param3 = (alpha_111 - per*abs(alpha_111)) +  2*per*abs(alpha_111)*rand(N,1);

  r3 = c1*param1 + c2*param2 + c3*param3;
  r2 = c1*param1 + c2*param2 + c3*alpha_111;

  d3 = ksdensity(r3,energy2);
  d2 = ksdensity(r2,energy2);

%
% Compute Sobol' indices based on the assumption of independent parameters sampled
% within 20% perturbations of the nominal values
%

  per = 0.1;
  a_1 = alpha_1 - per*abs(alpha_1);
  b_1 = alpha_1 + per*abs(alpha_1);
  a_2 = alpha_11 - per*abs(alpha_11);
  b_2 = alpha_11 + per*abs(alpha_11);
  a_3 = alpha_111 - per*abs(alpha_111);
  b_3 = alpha_111 + per*abs(alpha_111);

  D1 = (c1^2)*((b_1 - a_1)^2)/12;
  D2 = (c2^2)*((b_2 - a_2)^2)/12;
  D3 = (c3^2)*((b_3 - a_3)^2)/12;
  D = D1+D2+D3;

  S1 = D1/D;
  S2 = D2/D;
  S3 = D3/D;
  S = [S1 S2 S3]
  sum(S)

%
% Compute Sobol' indices that incorporate the correlation structure associated with V. 
%
 
  f1 = c1*s1 + c2*(mu2 + (V(1,2)/V(1,1))*(s1-mu1*ones(N,1))) ...
             + c3*(mu3 + (V(1,3)/V(1,1))*(s1-mu1*ones(N,1))) - f0;
  f2 = c2*s2 + c1*(mu1 + (V(2,1)/V(2,2))*(s2-mu2*ones(N,1))) ...
             + c3*(mu3 + (V(2,3)/V(2,2))*(s2-mu2*ones(N,1))) - f0;
  f3 = c3*s3 + c1*(mu1 + (V(3,1)/V(3,3))*(s3-mu3*ones(N,1))) ...
             + c2*(mu2 + (V(3,2)/V(3,3))*(s3-mu3*ones(N,1))) - f0;

  f = c1*s1 + c2*s2 + c3*s3;
  
  V12 = V(1:2,1:2);
  V13 = [V(1,1) V(1,3); V(3,1) V(3,3)];
  V23 = V(2:3,2:3);
  V12inv = inv(V12);
  V13inv = inv(V13);
  V23inv = inv(V23);

  for j=1:N
    f12(j,1) = c1*s1(j) + c2*s2(j) - f1(j) - f2(j) - f0(j) ...
               + c3*(mu3 + [V(3,1) V(3,2)]*V12inv*[s1(j)-mu1; s2(j)-mu2]);
    f13(j,1) = c1*s1(j) + c3*s3(j) - f1(j) - f3(j) - f0(j) ...
               + c2*(mu2 + [V(2,1) V(2,3)]*V13inv*[s1(j)-mu1; s3(j)-mu3]);
    f23(j,1) = c2*s2(j) + c3*s3(j) - f2(j) - f3(j) - f0(j) ...
               + c1*(mu1 + [V(1,2) V(1,3)]*V23inv*[s2(j)-mu2; s3(j)-mu3]);
  end
  f123 = f - f0 - f1 - f2 - f3 - f12 - f13 - f23;

  f_center = f - f0;
  var_f = sum(f_center.^2);

  S1_cov = sum(f1.*f_center)/var_f;
  S2_cov = sum(f2.*f_center)/var_f;
  S3_cov = sum(f3.*f_center)/var_f;
  S12_cov = sum(f12.*f_center)/var_f;
  S13_cov = sum(f13.*f_center)/var_f;
  S23_cov = sum(f23.*f_center)/var_f;
  S123_cov = sum(f123.*f_center)/var_f;
  

  S_cov = [S1_cov S2_cov S3_cov S12_cov S13_cov S23_cov S123_cov]
  sum(S_cov)


%
%  Plot the distributions for alpha_1, alpha_11, alpha_111 random and with alpha_111 fixed.
%

  figure(1)
  plot(energy,dens3,'b',energy,dens2,'b--','linewidth',3)
  axis([-16 -6 0 1.5])
  set(gca,'Fontsize',[22]);
  xlabel('Integrated Helmholtz Energy')
  ylabel('PDF')
  legend('3 Random Parameters','2 Random Parameters','Location','NorthEast') 

  
  figure(2)
  plot(energy2,d3,'b',energy2,d2,'b--','linewidth',3)
%  axis([-20 -8 0 1.3])
  set(gca,'Fontsize',[22]);
  xlabel('Integrated Helmholtz Energy')
  ylabel('PDF')
  legend('3 Random Parameters','2 Random Parameters','Location','NorthEast')