%
% Construct multivariate polynomial response surface for exponential model
%

% Generate 300 training points

%copyfile('PolyfitnTools.mltbx')

% matlab.addons.toolbox.installToolbox('PolyfitnTools.mltbx')
% 
%% copyfile('PolyfitnTools.mltbx');
%unzip('PolyfitnTools.mltbx');

%
%  Compute the solution at n_train = 300 training points and add sigma=0.1 random noise.

  m = 2;
  n_train = 300;                        
  q_train = -1+2*rand(n_train,m);
  f_train = exp(0.7*q_train(:,1)+0.3*q_train(:,2))+normrnd(0,0.01,n_train,1);

%
% Compute multiviariate polynomials of orders 1-5 and compute AIC and BIC.
%

  sigma=.1;
  for i = 1:5
    p = polyfitn(q_train,f_train,i);
    k = length(p.Coefficients);
    ssq = sum((polyvaln(p,q_train)-f_train).^2);                      % Sum of squares
    L = (1/((2*pi*sigma^2)^(n_train/2)))*exp(-ssq/(2*sigma^2));       % Gaussian likelihood
    aic(i) = 2*k-2*log(L);
    bic(i) = -2*log(L)+k*log(n_train);
  end
  order = find(aic == min(aic));                                      % Find minimium AIC value

%
% Construct response surface with appropriate polynomial
%

  p = polyfitn(q_train,f_train,order(1));

  qplot1 = linspace(-1,1,100);
  qplot2 = linspace(-1,1,100);
  [z1 z2] = meshgrid(qplot1,qplot2);
  pairs = [z1(:) z2(:)];

  resp = reshape(polyvaln(p,pairs),100,100);

%
% Generate 300 testing points
%

  n_test= 300;
  q_test= -1+2*rand(n_test,m);
  f_test = exp(0.7*q_test(:,1)+0.3*q_test(:,2))+normrnd(0,sigma,n_test,1);

%
% Compute RMSE
%

  resp_test = polyvaln(p,q_test);
  rmse = sqrt((1/n_test)*sum((f_test-resp_test).^2));

%
% Plot response surface with testing data
%
  figure(3)
  plot3(q_test(:,1),q_test(:,2),f_test,'ob','Linewidth',2)
  hold on
  surf(qplot1, qplot2, resp);
  hold off