%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vertex.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [n0,x,y,x1,x2,f1,f2] = vertex(j,n,u,v,v1,x0,f0,ipar,isplit,
% ichild,z,f,l,L)
% computes the base vertex x and the opposite vertex y of the box # j 
% of MCS and the 'neighboring vertices' x1 and x2 and their function
% values f1 and f2 needed for separable quadratic interpolation
% Input:
% j       label of the box whose vertices etc. are to be computed 
%         (in the main program mcs.m)
% n       dimension of the problem
% [u,v]   original box
% v1      opposite vertex of the original box
% x0(1:n,1:max(L)) coordinates used in the initialization list
%         x^0 = (x0(1,l(1)),...,x0(n,l(n)) base vertex of the original 
%	  box
% f0(1:max(L),1:n) function values appertaining to the initialization 
%	  list
%         f0(m,i) is the function value at x = (x0(1,istar(1)),..., 
%         x0(i-1,istar(i-1)),x0(i,m),x0(i+1,l(i+1)),...,x0(n,l(n)))
% f0(1:max(L),m), m > n  function values appertaining to a later box
%         split according to the initialization list
% ipar(m) label of the parent of box m
% isplit(m) = -(splitting index) of box m  if it is split in the
%             initialization procedure
%           = splitting index of box m  otherwise
% ichild(m) = -(number which child box m is)  if box m was generated
%             by splitting according to the initialization list (in
%             the initialization procedure or later)
%           = number which child box m is  otherwise
% z(1:2,:)  z(1,m) = value of the isplit(m)th coordinate of the base
%           vertex of box m
%           z(2,m) = k  if box m is split according to the 
%                       initialization list and f0(:,k) contains the
%                       function values obtained by splitting the box
%           z(2,m) = splitting value of box m  otherwise
% f(1:2,:)  f(1,m) = function value at the base vertex of box m
%           f(2,m) = function value at the splitting point of box m
%                    (if z(2,m) ~= Inf)
% l(1:n)  pointer defining the initial base vertex x^0 (see above)
% L(1:n)  the init. list contains L(i) values for coordinate i
% Output:
% n0(1:n)  n0(i) indicates that the ith coordinate has been split n0(i)
%          times in the history of box j
% x(1:n)   base vertex of box j
% y(1:n)   opposite vertex of box j
% x1(1:n), x2(1:n), f1(1:n), f2(1:n)
%          x1(i) and x2(i) and the corresponding function values f1(i)
%          and f2(i) are used for quadratic interpolation in the ith
%          coordinate 
% Since it is not always possible to find for each i two vertices that 
% differ from x only in coordinate i, we have to generate 'fictitious'
% vertices and function values using the assumption that the function
% is separable
% Uses the following functions/m-files:
% updtf.m
% vert1.m
% vert2.m
% vert3.m          

function [n0,x,y,x1,x2,f1,f2] = vertex(j,n,u,v,v1,x0,f0,ipar,isplit,ichild,z,f,l,L)
% initialization
% The coordinates of x, y, x1 and x2 are initially set to Inf to 
% indicate that these quantities haven't been found yet in the course of
% pursuing the history of box j
x = Inf*ones(n,1);
y = Inf*ones(n,1);
x1 = Inf*ones(n,1);
x2 = Inf*ones(n,1);
f1 = zeros(n,1);
f2 = zeros(n,1);
n0 = zeros(n,1);
fold = f(1,j);
m = j;
while m > 1
  i = abs(isplit(ipar(m)));
  n0(i) = n0(i) + 1;
  if ichild(m) == 1
    if x(i) == Inf | x(i) == z(1,ipar(m))  
      [x(i),x1(i),x2(i),f1(i),f2(i)] = vert1(2,z(:,ipar(m)),f(:,ipar(m)),x1(i),x2(i),f1(i),f2(i));
    else
      [f1,f2,fold] = updtf(n,i,x1,x2,f1,f2,fold,f(1,ipar(m)));
      [x1(i),x2(i),f1(i),f2(i)] = vert2(1,x(i),z(:,ipar(m)),f(:,ipar(m)),x1(i),x2(i),f1(i),f2(i));
    end
  elseif ichild(m) >= 2
    [f1,f2,fold] = updtf(n,i,x1,x2,f1,f2,fold,f(1,ipar(m)));
    if x(i) == Inf | x(i) == z(2,ipar(m))
      [x(i),x1(i),x2(i),f1(i),f2(i)] = vert1(1,z(:,ipar(m)),f(:,ipar(m)),x1(i),x2(i),f1(i),f2(i));
    else
      [x1(i),x2(i),f1(i),f2(i)] = vert2(2,x(i),z(:,ipar(m)),f(:,ipar(m)),x1(i),x2(i),f1(i),f2(i));
    end 
  end
  if 1 <= ichild(m) & ichild(m) <= 2 & y(i) == Inf
    y(i) = split1(z(1,ipar(m)),z(2,ipar(m)),f(1,ipar(m)),f(2,ipar(m)));
  end
  if ichild(m) < 0 
    % box m was generated by splitting according to the init. list
    % x0(i,j1) = ith coordinate of the base vertex of box m
    % the ith coordinate of the opposite vertex is the golden section 
    % split of x0(i,j2) and x0(i,j2+1) for 1 <= j2 <= L(i) - 1, it is 
    % u(i) for j1 = 1 and v(i) for j2 = L(i) 
    % x0(i,j1+j3) = ith coordinate of the 'neighboring vertex'
    if u(i) < x0(i,1)
      j1 = ceil(abs(ichild(m))/2);  
      j2 = floor(abs(ichild(m))/2);
      if (abs(ichild(m))/2 < j1  & j1 > 1) | j1 == L(i) 
        j3 = -1; 
      else
        j3 = 1;
      end
    else
      j1 = floor(abs(ichild(m))/2) + 1;
      j2 = ceil(abs(ichild(m))/2);
      if abs(ichild(m))/2 + 1 > j1 & j1 < L(i)
        j3 = 1;
      else
        j3 = -1;
      end
    end
    if isplit(ipar(m)) < 0 % box m was generated in the init. procedure
      k = i;
    else
      k = z(2,ipar(m)); 
      % box m was generated by a later split according to the init.
      % list; k points to the corresponding function values  
    end
    if j1 ~= l(i) | (x(i) ~= Inf & x(i) ~= x0(i,l(i)))
      [f1,f2,fold] = updtf(n,i,x1,x2,f1,f2,fold,f0(l(i),k));
    end 
    if x(i) == Inf | x(i) == x0(i,j1)
      x(i) = x0(i,j1);
      if x1(i) == Inf
        [x1(i),x2(i),f1(i),f2(i)] = vert3(j1,x0(i,:),f0(:,k),L(i),x1(i),x2(i),f1(i),f2(i));
      elseif x2(i) == Inf & x1(i) ~= x0(i,j1+j3) 
        x2(i) = x0(i,j1+j3);
        f2(i) = f2(i) + f0(j1+j3,k);
      elseif x2(i) == Inf
        if j1 ~= 1 & j1 ~= L(i)
          x2(i) = x0(i,j1-j3);
          f2(i) = f2(i) + f0(j1-j3,k);
        else
          x2(i) = x0(i,j1+2*j3);
          f2(i) = f2(i) + f0(j1+2*j3,k);
        end
      end
    else
      if x1(i) == Inf
        x1(i) = x0(i,j1);
        f1(i) = f1(i) + f0(j1,k);
        if x(i) ~= x0(i,j1+j3)
          x2(i) = x0(i,j1+j3);
          f2(i) = f2(i) + f0(j1+j3,k);
        end
      elseif x2(i) == Inf  
        if x1(i) ~= x0(i,j1)
          x2(i) = x0(i,j1);
          f2(i) = f2(i) + f0(j1,k);
        elseif x(i) ~= x0(i,j1+j3) 
          x2(i) = x0(i,j1+j3);
          f2(i) = f2(i) + f0(j1+j3,k);
        else
          if j1 ~= 1 & j1 ~= L(i)
            x2(i) = x0(i,j1-j3);
            f2(i) = f2(i) + f0(j1-j3,k);
          else
            x2(i) = x0(i,j1+2*j3);
            f2(i) = f2(i) + f0(j1+2*j3,k);
          end
        end     
      end
    end
    if y(i) == Inf
      if j2 == 0
        y(i) = u(i);
      elseif j2 == L(i) 
        y(i) = v(i);
      else
        y(i) = split1(x0(i,j2),x0(i,j2+1),f0(j2,k),f0(j2+1,k));
      end
    end 
  end
  m = ipar(m);
end
for i = 1:n
  if x(i) == Inf
    x(i) = x0(i,l(i));
    [x1(i),x2(i),f1(i),f2(i)] = vert3(l(i),x0(i,:),f0(:,i),L(i),x1(i),x2(i),f1(i),f2(i));
  end
  if y(i) == Inf
     y(i) = v1(i);
  end
end
