Initial commit.

+eval/feval.m

@@ -0,0 +1,49 @@
+% Copyright (c) 2018 Adanay Martín & Oliver Schütze.
+% This file is subject to the terms and conditions defined in
+% the file 'LICENSE.txt', which is part of this source code package.
+
+function [fx, fcount, fundef] = feval(f, x, iswarning, opts)
+% Vectorized objective function evaluation.
+% Each row of x is considered an individual to be evaluated.
+% iswarning defines whether to display a warning or an error if the
+% function value is undefined. The function value will be checked only if 
+% opts.FunValCheck is true.
+% opts are the optimization options.
+%
+% The result is always of size (m x nobj) even if m = 1.
+
+fcount = 0;
+fundef = false;
+
+m = size(x, 1); % no of individuals
+
+if m == 1 || opts.UseVectorized
+  fx = feval(f, x);
+  if m == 1
+    fx = fx(:)';
+  else
+    fx = fx(:, :);
+  end
+  fcount = m; 
+  if opts.FunValCheck 
+    fundef = val.checkval(fx, 'f', iswarning);
+  end
+else
+  for i = 1 : m
+    fxi = feval(f, x(i, :));
+    fxi = fxi(:)';
+
+    if i == 1
+      nobj = length(fxi);
+      fx = zeros(m, nobj);
+    end
+
+    fx(i, :) = fxi;
+    fcount = fcount + 1; 
+    if opts.FunValCheck
+      fundef = fundef || val.checkval(fxi, 'f', iswarning);
+    end
+  end
+end
+end
+

+eval/heval.m

@@ -0,0 +1,50 @@
+% Copyright (c) 2018 Adanay Martín & Oliver Schütze.
+% This file is subject to the terms and conditions defined in
+% the file 'LICENSE.txt', which is part of this source code package.
+
+function [Hx, Hcount, Hundef] = heval(H, x, iswarning, opts)
+% Vectorized Hessian function evaluation.
+% Each row of x is considered an individual to be evaluated.
+% iswarning defines whether to display a warning or an error if the
+% Hessian value is undefined. The Hessian value will be checked only if 
+% opts.FunValCheck is true.
+% opts are the optimization options.
+%
+% The result is always of size (m x n x n x nobj) even if m = 1.
+% The result of H(x) is assumed to be (n x n x nobj) if m = 1.
+
+Hcount = 0;
+Hundef = false;
+
+[m, n] = size(x); % m is the number of individuals and n is the number of variables.
+
+if m == 1 || opts.UseVectorized
+  Hx = feval(H, x);
+  if m == 1
+    Hx = Hx(:, :, :);
+    Hx = shiftdim(Hx, -1); % (1 x n x n x nobj)
+  else
+    Hx = Hx(:, :, :, :);
+  end
+  Hcount = m;
+  if opts.FunValCheck 
+    Hundef = val.checkval(Hx, 'H', iswarning);
+  end
+else
+  for i = 1 : m
+    Hxi = feval(H, x(i, :));
+    Hxi = Hxi(:, :, :);
+
+    if i == 1
+      nobj = size(Hxi, 3);
+      Hx = zeros(m, n, n, nobj);
+    end
+
+    Hx(i, :, :, :) = Hxi;
+    Hcount = Hcount + 1;
+    if opts.FunValCheck
+      Hundef = Hundef || val.checkval(Hxi, 'H', iswarning);
+    end
+  end
+end
+end

+eval/hevalandmodchol.m

@@ -0,0 +1,28 @@
+% Copyright (c) 2018 Adanay Martín & Oliver Schütze.
+% This file is subject to the terms and conditions defined in
+% the file 'LICENSE.txt', which is part of this source code package.
+
+function [Hx, Hcount, Hundef, Lx, Lmodif] = hevalandmodchol(H, x, iswarning, opts)
+% Vectorized Hessian function evaluation plus a Modified Cholesky 
+% decomposition.
+% If Lmodif = true, then Hx is not positive definite, thus Lx * Lx' needs  
+% to be computed.
+% Each row of x is considered an individual to be evaluated.
+% iswarning defines whether to display a warning or an error if the
+% Hessian value is undefined. The Hessian value will be checked only if 
+% opts.FunValCheck is true.
+% The Cholesky modification will be performed only if the Hessian is
+% defined.
+% opts are the optimization options.
+%
+% The result is always of size (m x n x n x nobj) even if m = 1.
+% The result of H(x) is assumed to be (n x n x nobj) if m = 1.
+
+Lx = [];
+Lmodif = [];
+
+[Hx, Hcount, Hundef] = eval.heval(H, x, iswarning, opts);
+if ~Hundef && nargout > 3
+  [Lx, Lmodif] = utils.modcholn(Hx, 2, 3);
+end
+end

+eval/hevalorfd.m

@@ -0,0 +1,91 @@
+% Copyright (c) 2018 Adanay Martín & Oliver Schütze.
+% This file is subject to the terms and conditions defined in
+% the file 'LICENSE.txt', which is part of this source code package.
+
+function [Hx, Hcount, Hundef, Happrox, Hident,...
+          fx, fcount, Jx, Jcount, Jundef, nobj] = hevalorfd(H, f, J, x, fx, Jx, lb, ub, iswarning, formident, opts)
+% Vectorized Hessian function evaluation.
+% If H is empty, finite differences (FD) are used.
+% Each row of x is considered an individual to be evaluated.
+% iswarning defines whether to display a warning or an error if the
+% Hessian value is undefined. The Hessian value will be checked only if 
+% opts.FunValCheck is true. I.e., if opts.FunValCheck is true and H(x) is 
+% undef, then
+% + iswarning = true implies that a warning will be displayed and the
+% Hessian will be computed by FD if the Jacobian or the objective functions 
+% are available. Otherwise the Hessian will be approximated to the identity 
+% matrix.
+% + iswarning = false implies that an error will be thrown.
+% In case that FD is applied with the Jacobian, the value of the Jacobian 
+% is the one that will be checked if opts.FunValCheck is true. If also
+% iswarning = true, then FD will be re-computed using the objective 
+% function f if available.
+% formident will force the formation of the identity to be returned in Hx 
+% in case there is no other choice than approximating Hx to the identity. 
+% Otherwise Hx will be empty (but known to be the identity if Hident is true).
+% opts are the optimization options.
+%
+% The result is always of size (m x n x n x nobj) even if m = 1.
+% The result of H(x) is assumed to be (n x n x nobj) if m = 1.
+% fx (if entered) is assumed to be of size (m x obj) even if m = 1 
+% (after being calculated using one of the vec eval functions).
+% Jx (if entered) is assumed to be of size (m x obj x n) even if m = 1 
+% (after being calculated using one of the vec eval functions).
+
+Hcount = 0;
+Hundef = false;
+Happrox = false;
+Hident = false; % the Hessian was approx to the identity matrix
+
+fcount = 0;
+
+Jcount = 0;
+Jundef = false;
+
+m = size(x, 1); % m is the number of individuals
+
+if isempty(H) 
+  hess();
+else
+  [Hx, Hcount, Hundef] = eval.heval(H, x, iswarning, opts);
+  nobj = size(Hx, 4);
+  if Hundef
+    if ~isempty(f) || ~isempty(J)
+      hess();
+    else
+      Happrox = true;
+      Hident = true;
+      Hx = eval.hevalorsd([], x, nobj, iswarning, formident, opts);
+    end
+  end
+end
+
+  function [] = hess()
+    Happrox = true;
+
+    if m == 1 && ~isempty(Jx)
+      Jx = shiftdim(Jx, 1);
+    end
+
+    [Hx, fx, fcount, Jx, Jcount, Jundef] = fd.hessian(f, J, x, fx, Jx, lb, ub, [], opts);
+    nobj = size(Hx, 4);
+    if m == 1
+      Hx = shiftdim(Hx, -1); % (1 x n x n x nobj)
+    end
+    if Jundef
+      if ~isempty(f)
+        [Hx, fx, fcount] = fd.hessian(f, [], x, fx, [], lb, ub, [], opts);
+        if m == 1
+          Hx = shiftdim(Hx, -1); % (1 x n x n x nobj)
+        end
+      else
+        Hident = true;
+        Hx = eval.hevalorsd([], x, nobj, iswarning, formident, opts);
+      end
+    end
+
+    if m == 1 && ~isempty(Jx)
+      Jx = shiftdim(Jx, -1); % (1 x nobj x n)
+    end    
+  end
+end

+eval/hevalorfdandmodchol.m

@@ -0,0 +1,56 @@
+% Copyright (c) 2018 Adanay Martín & Oliver Schütze.
+% This file is subject to the terms and conditions defined in
+% the file 'LICENSE.txt', which is part of this source code package.
+
+function [Hx, Hcount, Hundef, Happrox, Hident, Lx, Lmodif,...
+          fx, fcount, Jx, Jcount, Jundef, nobj] = hevalorfdandmodchol(H, f, J, x, fx, Jx, lb, ub, iswarning, formident, opts)
+% Vectorized Hessian function evaluation plus a modified Cholesky 
+% decomposition.
+% If H is empty, finite differences (FD) are used.
+% A Modified Cholesky decomposition is performed afterwards.
+% If Lmodif = true, then Hx is not positive definite, thus Lx * Lx' needs  
+% to be computed.
+% Each row of x is considered an individual to be evaluated.
+% iswarning defines whether to display a warning or an error if the
+% Hessian value is undefined. The Hessian value will be checked only if 
+% opts.FunValCheck is true. I.e., if opts.FunValCheck is true and H(x) is 
+% undef, then
+% + iswarning = true implies that a warning will be displayed and the
+% Hessian will be computed by FD if the Jacobian or the objective functions 
+% are available. Otherwise the Hessian will be approximated to the identity 
+% matrix.
+% + iswarning = false implies that an error will be thrown.
+% In case that FD is applied with the Jacobian, the value of the Jacobian 
+% is the one that will be checked if opts.FunValCheck is true. If also
+% iswarning = true, then FD will be re-computed using the objective 
+% function f if available.
+% formident will force the formation of the identity to be returned in Hx 
+% (and Lx) in case there is no other choice than approximating Hx to the 
+% identity. Otherwise Hx (and Lx) will be empty (but known to be the 
+% identity if Hident is true).
+% opts are the optimization options.
+%
+% The result is always of size (m x n x n x nobj) even if m = 1.
+% The result of H(x) is assumed to be (n x n x nobj) if m = 1.
+% fx (if entered) is assumed to be of size (m x obj) even if m = 1 
+% (after being calculated using one of the vec eval functions).
+% Jx (if entered) is assumed to be of size (m x obj x n) even if m = 1 
+% (after being calculated using one of the vec eval functions).
+
+Lx = [];
+Lmodif = [];
+
+[Hx, Hcount, Hundef, Happrox, Hident,...
+ fx, fcount, Jx, Jcount, Jundef, nobj] = eval.hevalorfd(H, f, J, x, fx, Jx, lb, ub, iswarning, formident, opts);
+
+if nargout > 5
+  if Hident
+    Lmodif = false;
+    if formident
+      Lx = Hx;
+    end
+  else
+    [Lx, Lmodif] = utils.modcholn(Hx, 2, 3);
+  end
+end
+end

+eval/hevalorqn.m

@@ -0,0 +1,72 @@
+% Copyright (c) 2018 Adanay Martín & Oliver Schütze.
+% This file is subject to the terms and conditions defined in
+% the file 'LICENSE.txt', which is part of this source code package.
+
+function [Hx1, Hcount, Hundef, Happrox, Hident] = hevalorqn(H, x1, Jx1, x0, Jx0, Hx0, iswarning, formident, opts)
+% Vectorized Hessian function evaluation.
+% If H is empty, a QN update will be performed.
+% Each row of x1, (and x0) is considered an individual to be evaluated.
+% iswarning defines whether to display a warning or an error if the
+% Hessian value is undefined. The Hessian value will be checked only if 
+% opts.FunValCheck is true. I.e., if opts.FunValCheck is true and H(x1) is 
+% undef, then
+% + iswarning = true implies that a warning will be displayed and the
+% Hessian will be updated if the previous iteration values are provided.
+% Otherwise it will be approximated to the identity.
+% + iswarning = false implies that an error will be thrown.
+% formident will force the formation of the identity to be returned in Hx1 
+% in case there is no other choice than approximating Hx1 to the identity. 
+% This may happen e.g. if the previous iteration values are not provided. 
+% Otherwise Hx1 will be empty (but known to be the identity if Hident is true).
+% opts are the optimization options.
+%
+% The result is always of size (m x n x n x nobj) even if m = 1.
+% The result of H(x) is assumed to be (n x n x nobj) if m = 1.
+% Jx1 (and Jx0) (if entered) is assumed to be of size (m x obj x n) even  
+% if m = 1 (after being calculated using one of the vec eval functions).
+% Hx0 (if entered) is assumed to be of size (m x n x n x nobj) even if m = 1 
+% (after being calculated using one of the vec eval functions).
+
+Hcount = 0;
+Hundef = false;
+Happrox = false;
+Hident = false; % the Hessian was approx to the identity matrix
+
+m0 = size(x0, 1); 
+m1 = size(x1, 1); 
+m = max(m0, m1); % m is the number of individuals
+
+if isempty(H)                
+  hess();
+else
+  [Hx1, Hcount, Hundef] = eval.heval(H, x1, iswarning, opts);
+  if Hundef
+    if ~isempty(Jx1) && ~isempty(x0) && ~isempty(Jx0)
+      hess();
+    else
+      nobj = size(Hx1, 4);
+      Happrox = true;
+      Hident = true;
+      Hx1 = eval.hevalorsd([], x1, nobj, iswarning, formident, opts);
+    end
+  end
+end
+
+  function [] = hess()
+    Happrox = true;
+
+    if m0 == 1
+      Jx0 = shiftdim(Jx0, 1);
+      Hx0 = shiftdim(Hx0, 1);
+    end
+    if m1 == 1
+      Jx1 = shiftdim(Jx1, 1);
+    end
+
+    Hx1 = qn.bfgs(Hx0, x0, Jx0, x1, Jx1, opts);
+
+    if m == 1
+      Hx1 = shiftdim(Hx1, -1); % (1 x n x n x nobj)
+    end
+  end
+end