%
%                 sparse_int.m
%

  clear all

%
% Input constants
%

%
%  Compute the function on a uniform grid.
%

  N1 = 41;
  N2 = 41;
  xvals1 = linspace(0,1,N1);
  xvals2 = linspace(0,1,N2);
  unit1 = ones(size(xvals1));
  unit2 = ones(size(xvals2));
  unit = ones(N1,N2);
  [X1,X2] = meshgrid(xvals1,xvals2);
  Z = X2.*exp(X1.*X2); 


%
% Compute the interpolant
%

  ell = 3;

  for j = 1:ell
    if j==1;
      n_ell = 1;
      qvals = 0;
    else
      n_ell = 2^(j-1)+1;
      N = n_ell-1; 
      qvals = clencurt(N);
    end
    q = (1/2)*(1-qvals);
    if j==1
      n1 = n_ell; 
      x1 = q;  y1 = q'; 
      f0 = fun_exp(x1,y1,n1,n1);
      F01 = f0*ones(size(xvals1));
      F02 = f0*ones(size(xvals2));
      F0 = f0*ones(N1,N2);
      Iprev = f0*ones(N1,N2);
    elseif j==2
      x2 = q;  y2 = q';
      n2 = n_ell;
      fvals12 = fun_exp(x1,y2,n1,n2);
      fvals21 = fun_exp(x2,y1,n2,n1);
      den1 = (x2(1)-x2(2))*(x2(1)-x2(3));
      den2 = (x2(2)-x2(1))*(x2(2)-x2(3));
      den3 = (x2(3)-x2(1))*(x2(3)-x2(2));
      I12 = fvals12(1)*(xvals2-x2(2)).*(xvals2-x2(3))/den1 ...
          + fvals12(2)*(xvals2-x2(1)).*(xvals2-x2(3))/den2 ...
          + fvals12(3)*(xvals2-x2(1)).*(xvals2-x2(2))/den3 - F02;
      I21 = fvals21(1)*(xvals1-x2(2)).*(xvals1-x2(3))/den1 ...
          + fvals21(2)*(xvals1-x2(1)).*(xvals1-x2(3))/den2 ...
          + fvals21(3)*(xvals1-x2(1)).*(xvals1-x2(2))/den3 - F01;
      I = Iprev + unit1'*I12 + I21'*unit2;
      Iprev = I;
    else
      x3 = q; y3 = q';
      n3 = n_ell;
      fvals13 = fun_exp(x1,y3,n1,n3);
      fvals31 = fun_exp(x3,y1,n3,n1); 
      fvals22 = fun_exp(x2,y2,n2,n2);
      f12 = fvals22(1,2);
      f32 = fvals22(3,2);
      f21 = fvals22(2,1);
      f23 = fvals22(2,3);
      f22 = fvals22(2,2);
      den1 = (x3(1)-x3(2))*(x3(1)-x3(3))*(x3(1)-x3(4))*(x3(1)-x3(5));
      den2 = (x3(2)-x3(1))*(x3(2)-x3(3))*(x3(2)-x3(4))*(x3(2)-x3(5));
      den3 = (x3(3)-x3(1))*(x3(3)-x3(2))*(x3(3)-x3(4))*(x3(3)-x3(5)); 
      den4 = (x3(4)-x3(1))*(x3(4)-x3(2))*(x3(4)-x3(3))*(x3(4)-x3(5));
      den5 = (x3(5)-x3(1))*(x3(5)-x3(2))*(x3(5)-x3(3))*(x3(5)-x3(4));
      I13 = fvals13(1)*(xvals2-x3(2)).*(xvals2-x3(3)).*(xvals2-x3(4)).*(xvals2-x3(5))/den1 ...
          + fvals13(2)*(xvals2-x3(1)).*(xvals2-x3(3)).*(xvals2-x3(4)).*(xvals2-x3(5))/den2 ...
          + fvals13(3)*(xvals2-x3(1)).*(xvals2-x3(2)).*(xvals2-x3(4)).*(xvals2-x3(5))/den3 ...
          + fvals13(4)*(xvals2-x3(1)).*(xvals2-x3(2)).*(xvals2-x3(3)).*(xvals2-x3(5))/den4 ...
          + fvals13(5)*(xvals2-x3(1)).*(xvals2-x3(2)).*(xvals2-x3(3)).*(xvals2-x3(4))/den5;
      I31 = fvals31(1)*(xvals1-x3(2)).*(xvals1-x3(3)).*(xvals1-x3(4)).*(xvals1-x3(5))/den1 ...
          + fvals31(2)*(xvals1-x3(1)).*(xvals1-x3(3)).*(xvals1-x3(4)).*(xvals1-x3(5))/den2 ...
          + fvals31(3)*(xvals1-x3(1)).*(xvals1-x3(2)).*(xvals1-x3(4)).*(xvals1-x3(5))/den3 ...
          + fvals31(4)*(xvals1-x3(1)).*(xvals1-x3(2)).*(xvals1-x3(3)).*(xvals1-x3(5))/den4 ...
          + fvals31(5)*(xvals1-x3(1)).*(xvals1-x3(2)).*(xvals1-x3(3)).*(xvals1-x3(4))/den5;
          den1 = (x2(1)-x2(2))*(x2(1)-x2(3));
          den2 = (x2(2)-x2(1))*(x2(2)-x2(3));
          den3 = (x2(3)-x2(1))*(x2(3)-x2(2));
      I22 = fvals22(1,1)*((xvals1-x2(2)).*(xvals1-x2(3)))'*((xvals2-x2(2)).*(xvals2-x2(3)))/(den1*den1) ...
          + fvals22(2,1)*((xvals1-x2(1)).*(xvals1-x2(3)))'*((xvals2-x2(2)).*(xvals2-x2(3)))/(den2*den1) ...
          + fvals22(3,1)*((xvals1-x2(1)).*(xvals1-x2(2)))'*((xvals2-x2(2)).*(xvals2-x2(3)))/(den3*den1) ...
          + fvals22(1,2)*((xvals1-x2(2)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(3)))/(den1*den2) ...
          + fvals22(2,2)*((xvals1-x2(1)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(3)))/(den2*den2) ...
          + fvals22(3,2)*((xvals1-x2(1)).*(xvals1-x2(2)))'*((xvals2-x2(1)).*(xvals2-x2(3)))/(den3*den2) ...
          + fvals22(1,3)*((xvals1-x2(2)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(2)))/(den1*den3) ...
          + fvals22(2,3)*((xvals1-x2(1)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(2)))/(den2*den3) ...
          + fvals22(3,3)*((xvals1-x2(1)).*(xvals1-x2(2)))'*((xvals2-x2(1)).*(xvals2-x2(2)))/(den3*den3) ...
          - unit1'*(I12+F02) - (I21+F01)'*unit2; - 2*f22*unit - f12*unit - f32*unit - f21*unit - f23*unit;
      %I = Iprev + I31'*unit2 + unit1'*I13 + I22;
      I = I31'*unit2 + unit1'*I13 + I22;
    end
  end

  %I = Iprev;

  xvals = 1/2*ones(1,5);
  nodes1 = [x3' ;
           xvals];
  nodes2 = nodes1'; 
  nodes3 = [0 0.5 1 0 0.5 1 0 0.5 1;
            0 0 0 0.5 0.5 0.5 1 1 1]; 
  nodes = [0 0.5 1 0  0   0  0.5   0   x3(2) 0.5 x3(4) 1 0.5 0 0.5 1;
           0  0  0 0 0.5 1.0 x3(2) 0.5 0.5 0.5 0.5 0.5   x3(4) 1 1 1];
  xnodes = nodes(1,:);
  ynodes = nodes(2,:);
  znodes = ynodes.*exp(xnodes.*ynodes); 

%
% Plot the function and observations
%


  figure(1)
  contour(X1,X2,Z,20,'linewidth',2);
  axis('equal')
  hold on
  scatter(nodes1(1,:),nodes1(2,:),'b','linewidth',4)
  scatter(nodes1(2,:),nodes1(1,:),'b','linewidth',4)
  scatter(nodes3(1,:),nodes3(2,:),'b','linewidth',4)
  hold off
  set(gca,'Fontsize',[22]);
  xlabel('q_1')
  ylabel('q_2')
  %print -depsc full_contour

  figure(2)
  mesh(X1,X2,Z)
  set(gca,'Fontsize',[22]);
  xlabel('q_1')
  ylabel('q_2')



  figure(3)
  contour(X1',X2',I,20,'linewidth',2);
  axis('equal')
  hold on
  scatter(nodes1(1,:),nodes1(2,:),'b','linewidth',4)
  scatter(nodes1(2,:),nodes1(1,:),'b','linewidth',4)
  scatter(nodes3(1,:),nodes3(2,:),'b','linewidth',4)
  hold off
  set(gca,'Fontsize',[22]);
  xlabel('q_1')
  ylabel('q_2')
  %print -depsc sparse_contour


  figure(4)
  mesh(X1',X2',I)
  axis([0 1 0 1 0 3]);
  hold on
  plot3(xnodes,ynodes,znodes,'ok','linewidth',3)
  hold off
  set(gca,'Fontsize',[22]);
  xlabel('q_1')
  ylabel('q_2')
  zlabel('f(q_1,q_2)') 
  %print -depsc function_sparse_pts

  int_points = [I(1,1) I(1,21) I(1,41) I(21,1), I(21,21), I(21,41), I(41,1), I(41,21), I(41,41)]
  Z_points = [Z(1,1) Z(1,21) Z(1,41) Z(21,1), Z(21,21), Z(21,41), Z(41,1), Z(41,21), Z(41,41)]




