function [fmin,fmax] = plot_surface(func,x,x_min,x_max,options,consts,fig_no)
% [f_min,f_max]=plot_surface(func,x,x_min,x_max,options,consts,fig_no)
%
% Draw a surface plot of J(x) vs. x(i),x(j), where all other values in x
% are held constant.   
%
% INPUT
% ======
%  func    :   the name of the matlab function to be optimized in the form
%               [objective, constraints] = func(x,consts)
%  x       :   the vector of initial parameter values ... a column vector
%  x_min   :   minimum permissible values of the parameters, x
%  x_max   :   maximum permissible values of the parameters, x
%  options :   options(1)  = 3  ... draw a surface
%              options(4)  = tol_g  tolerance on constraint functions
%              options(6)  = penalty  on constraint violations
%              options(7)  = exponent on constraint violations
%              options(11) = 1st index for plotting cost function surface
%              options(12) = 2nd index for plotting cost function surface
%              options(13) = number of 1st index values to plot cost funct surf
%              options(14) = number of 2nd index values to plot cost funct surf
%  consts  :   an optional vector of constats to be passed to func(x,consts)
%  fig_no  :   figure number used to draw the figure
%
% OUTPUT
% ======
%  fmin,fmax :  max and min values of the meshed surface data

%  H.P. Gavin, Civil & Environ. Eng'g, Duke Univ.  

  tol_g   = options(4);
  penalty = options(6);
  q       = options(7);
  i       = options(11);
  j       = options(12);
  Ni      = options(13);
  Nj      = options(14);

  x_init = x;

  x_i = linspace( x_min(i), x_max(i), Ni);
  x_j = linspace( x_min(j), x_max(j), Nj);
  f_mesh = NaN(Ni,Nj);
  for ii=1:Ni
      for jj=1:Nj
          x([i,j]) = [ x_i(ii) ; x_j(jj) ];
          [f,g] = feval(func,x,consts);
%         if g < 0
            f_mesh(ii,jj) = f + penalty*sum(g.*(g>tol_g))^q;
%         end
	  if ( penalty <= 0 &&  max(g) > tol_g ), f_mesh(ii,jj) = NaN; end
      end
  end
  frms = sqrt(sum(sum(f_mesh.^2)));
  fmin = 0.9*min(min(f_mesh));
  fmax = min(max(max(f_mesh)),5*frms);
% cmap = prisim(100);
% cmap = rainbow(100);
% color_index = round(100*(f_mesh-fmin)/(fmax-fmin));  

  [ii,jj] = find(min(min(f_mesh)) == f_mesh);

  figure(fig_no)
    clf
    mesh(x_i,x_j,f_mesh', 'LineWidth',1.5)
%   plot3(x_i'*ones(1,Nj),ones(Ni,1)*x_j,f_mesh,'ok')
      xlabel(sprintf('x_%d',i))
      ylabel(sprintf('x_%d',j))
      zlabel(sprintf('objective   f_{aug}(x_%d,x_%d)',i,j))
      axis([ min(x_i), max(x_i), min(x_j), max(x_j), fmin, fmax]);
    hold on
    [f,g] = feval(func,x_init,consts);
    fa = f + penalty*sum(g.*(g>tol_g))^q;
    plot3(x_init(i),x_init(j),fa+(fmax-fmin)/50,'og', 'MarkerSize',8, 'LineWidth',4);
    plot3(x_i(ii),x_j(jj),f_mesh(ii,jj)+(fmax-fmin)/50,'or', 'MarkerSize',8, 'LineWidth',4);
    title(sprintf('grid opt: x_%d = %f , x_%d = %f , f(%d,%d) = %f \n',  i, x_i(ii), j, x_j(jj), x_i(ii),x_j(jj), f_mesh(ii,jj)))  
 end

 

% plot_surface ================================================================
% updated 2015-09-29  2016-03-23  2016-09-12, 2020-01-20
