First debugging and test script for existing evaluations

Author: Lindon Roberts

dfols/__init__.py

@@ -45,3 +45,5 @@
 from .solver import *
 __all__ = ['solve', 'OptimResults']
 
+from .evaluations_database import *
+__all__ += ['EvaluationDatabase']

dfols/controller.py

@@ -427,8 +427,8 @@ def initialise_from_database(self, eval_database, number_of_samples, params):
                                                                                             dykstra_tol=params("dykstra.d_tol"))
 
         # Add suitable pre-existing evaluations
-        module_logger.debug("Adding %g pre-existing evaluations (aside from x0) to initial model" % len(perturbation_idx))
         for i, idx in enumerate(perturbation_idx):
+            module_logger.info("Adding pre-existing evaluation %g to initial model" % idx)
             x, rx = eval_database.get_eval(idx)
             self.model.change_point(i + 1, x - self.model.xbase, rx, -idx)  # use eval_num = -idx
 

dfols/evaluations_database.py

@@ -89,7 +89,7 @@ def select_starting_evals(self, delta, xl=None, xu=None, projections=[], tol=1e-
         n = len(xbase)
         module_logger.debug("Selecting starting evaluations from existing database")
         module_logger.debug("Have %g evaluations to choose from" % len(self))
-        module_logger.debug("Using base index %g, x0 =" % base_idx, xbase)
+        module_logger.debug("Using base index %g" % base_idx)
 
         # For linear interpolation, we will use the matrix
         # M = [[1, 0], [0, L]] where L has rows (xi-xbase)/delta

dfols/solver.py

@@ -165,12 +165,15 @@ def solve_main(objfun, x0, argsf, xl, xu, projections, npt, rhobeg, rhoend, maxf
     # Evaluate at x0 (keep nf, nx correct and check for f < 1e-12)
     # The hard bit is determining what m = len(r0) should be, and allocating memory appropriately
     if r0_avg_old is None:
+        exit_info = None
         if x0_is_eval_database:
             # We have already got r(x0), so just extract this information
             nf = nf_so_far
             nx = nx_so_far
+            num_samples_run = 1
             r0_avg = x0.get_rx(x0.get_starting_eval_idx())
             m = len(r0_avg)
+            module_logger.info("Using pre-existing evaluation %g as starting point" % (x0.get_starting_eval_idx()))
         else:
             number_of_samples = max(nsamples(rhobeg, rhobeg, 0, nruns_so_far), 1)
             # Evaluate the first time...
@@ -188,7 +191,6 @@ def solve_main(objfun, x0, argsf, xl, xu, projections, npt, rhobeg, rhoend, maxf
             rvec_list[0, :] = r0
             obj_list[0] = obj0
             num_samples_run = 1
-            exit_info = None
 
             for i in range(1, number_of_samples):  # skip first eval - already did this
                 if nf >= maxfun:
@@ -972,6 +974,7 @@ def solve(objfun, x0, h=None, lh=None, prox_uh=None, argsf=(), argsh=(), argspro
         assert len(x0) > 0, "evaluation database x0 cannot be empty"
         assert 0 <= x0.get_starting_eval_idx() < len(x0), "evaluation database must have valid starting index set"
         x0_is_eval_database = True
+        n = len(x0.get_x(x0.get_starting_eval_idx()))
     else:
         x0 = np.array(x0).astype(float)
         n = len(x0)

examples/existing_evaluations.py

@@ -0,0 +1,60 @@
+"""
+Demonstration of using database of existing evaluations to speed up DFO-LS
+
+Test problem is the 'Watson function': for details, see
+J. J. More, B. S. Garbow, K. E. Hillstrom. Testing Unconstrained Optimization Software.
+ACM Transactions on Mathematical Software, 7:1 (1981), pp. 17-41.
+"""
+import numpy as np
+import dfols
+
+# Define the objective function
+def watson(x):
+    n = len(x)
+    m = 31
+    fvec = np.zeros((m,), dtype=float)
+
+    for i in range(1, 30):  # i=1,...,29
+        div = float(i) / 29.0
+        s1 = 0.0
+        dx = 1.0
+        for j in range(2, n + 1):  # j = 2,...,n
+            s1 = s1 + (j - 1) * dx * x[j - 1]
+            dx = div * dx
+        s2 = 0.0
+        dx = 1.0
+        for j in range(1, n + 1):  # j = 1,...,n
+            s2 = s2 + dx * x[j - 1]
+            dx = div * dx
+        fvec[i - 1] = s1 - s2 ** 2 - 1.0
+
+    fvec[29] = x[0]
+    fvec[30] = x[1] - x[0] ** 2 - 1.0
+
+    return fvec
+
+# Define the starting point
+n = 6
+x0 = 0.5 * np.ones((n,), dtype=float)
+
+# When n=6, we expect f(x0) ~ 16.4308 and f(xmin) ~ 0.00228767 at xmin ~ [ -0.0157, 1.0124, 1.2604, -1.5137, 0.992996]
+
+# For optional extra output details
+import logging
+logging.basicConfig(level=logging.INFO, format='%(message)s')
+
+# Now build a database of evaluations
+eval_db = dfols.EvaluationDatabase()
+eval_db.append(x0, watson(x0), make_starting_eval=True)  # make x0 the starting point
+
+# Note: x0, x1 and x2 are colinear, so at least one of x1 and x2 shouldn't be included in the initial model
+x1 = np.ones((n,), dtype=float)
+x2 = np.zeros((n,), dtype=float)
+x3 = np.arange(n).astype(float)
+eval_db.append(x1, watson(x1))
+eval_db.append(x2, watson(x2))
+eval_db.append(x3, watson(x3))
+
+soln = dfols.solve(watson, eval_db)  # replace x0 with eval_db
+
+print(soln)