Include evaluation numbers in return info

Author: lindonroberts

dfols/controller.py

@@ -240,12 +240,12 @@ def initialise_coordinate_directions(self, number_of_samples, num_directions, pa
                 # Handle exit conditions (f < min obj value or maxfun reached)
                 if exit_info is not None:
                     if num_samples_run > 0:
-                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                               x_in_abs_coords=True)
                     return exit_info  # return & quit
 
                 # Otherwise, add new results (increments model.npt_so_far)
-                self.model.change_point(k+1, x - self.model.xbase, rvec_list[0, :])  # expect step, not absolute x
+                self.model.change_point(k+1, x - self.model.xbase, rvec_list[0, :], self.nx)  # expect step, not absolute x
                 for i in range(1, num_samples_run):
                     self.model.add_new_sample(k+1, rvec_extra=rvec_list[i, :])
 
@@ -299,12 +299,12 @@ def initialise_coordinate_directions(self, number_of_samples, num_directions, pa
             # Handle exit conditions (f < min obj value or maxfun reached)
             if exit_info is not None:
                 if num_samples_run > 0:
-                    self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                    self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                           x_in_abs_coords=True)
                 return exit_info  # return & quit
 
             # Otherwise, add new results (increments model.npt_so_far)
-            self.model.change_point(k, x - self.model.xbase, rvec_list[0, :])  # expect step, not absolute x
+            self.model.change_point(k, x - self.model.xbase, rvec_list[0, :], self.nx)  # expect step, not absolute x
             for i in range(1, num_samples_run):
                 self.model.add_new_sample(k, rvec_extra=rvec_list[i, :])
 
@@ -351,13 +351,13 @@ def initialise_random_directions(self, number_of_samples, num_directions, params
                 # Handle exit conditions (f < min obj value or maxfun reached)
                 if exit_info is not None:
                     if num_samples_run > 0:
-                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                               x_in_abs_coords=True)
                     return exit_info  # return & quit
 
                 # Otherwise, add new results (increments model.npt_so_far)
                 self.model.change_point(1 + ndirns, x - self.model.xbase,
-                                        rvec_list[0, :])  # expect step, not absolute x
+                                        rvec_list[0, :], self.nx)  # expect step, not absolute x
                 for i in range(1, num_samples_run):
                     self.model.add_new_sample(1 + ndirns, rvec_extra=rvec_list[i, :])
         else:
@@ -371,12 +371,12 @@ def initialise_random_directions(self, number_of_samples, num_directions, params
                 # Handle exit conditions (f < min obj value or maxfun reached)
                 if exit_info is not None:
                     if num_samples_run > 0:
-                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                               x_in_abs_coords=True)
                     return exit_info  # return & quit
 
                 # Otherwise, add new results (increments model.npt_so_far)
-                self.model.change_point(1 + ndirns, x - self.model.xbase, rvec_list[0, :])  # expect step, not absolute x
+                self.model.change_point(1 + ndirns, x - self.model.xbase, rvec_list[0, :], self.nx)  # expect step, not absolute x
                 for i in range(1, num_samples_run):
                     self.model.add_new_sample(1 + ndirns, rvec_extra=rvec_list[i, :])
 
@@ -408,7 +408,7 @@ def add_new_direction_while_growing(self, number_of_samples, params, min_num_ste
             # Handle exit conditions (f < min obj value or maxfun reached)
             if exit_info is not None:
                 if num_samples_run > 0:
-                    self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                    self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                           x_in_abs_coords=True)
                 return exit_info  # return & quit
 
@@ -422,7 +422,7 @@ def add_new_direction_while_growing(self, number_of_samples, params, min_num_ste
                     return exit_info  # return & quit
 
             # Otherwise, add new results
-            self.model.change_point(kmin, xnew, rvec_list[0, :])  # expect step, not absolute x
+            self.model.change_point(kmin, xnew, rvec_list[0, :], self.nx)  # expect step, not absolute x
             for i in range(1, num_samples_run):
                 self.model.add_new_sample(kmin, rvec_extra=rvec_list[i, :])
 
@@ -548,12 +548,12 @@ def geometry_step(self, knew, adelt, number_of_samples, params):
         # Handle exit conditions (f < min obj value or maxfun reached)
         if exit_info is not None:
             if num_samples_run > 0:
-                self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                       x_in_abs_coords=True)
             return exit_info  # didn't fix geometry - return & quit
 
         # Otherwise, add new results
-        self.model.change_point(knew, xnew, rvec_list[0, :])  # expect step, not absolute x
+        self.model.change_point(knew, xnew, rvec_list[0, :], self.nx)  # expect step, not absolute x
         for i in range(1, num_samples_run):
             self.model.add_new_sample(knew, rvec_extra=rvec_list[i, :])
 
@@ -768,8 +768,8 @@ def soft_restart(self, number_of_samples, nruns_so_far, params, x_in_abs_coords_
         if x_in_abs_coords_to_save is not None:
             assert rvec_to_save is not None, "Soft restart: specified x_to_save but not rvec_to_save"
             assert nsamples_to_save is not None, "Soft restart: specified x_to_save but not nsamples_to_save"
-            self.model.save_point(x_in_abs_coords_to_save, rvec_to_save, nsamples_to_save, x_in_abs_coords=True)
-        self.model.save_point(self.model.xopt(abs_coordinates=True), self.model.ropt(),
+            self.model.save_point(x_in_abs_coords_to_save, rvec_to_save, nsamples_to_save, self.nx, x_in_abs_coords=True)
+        self.model.save_point(self.model.xopt(abs_coordinates=True), self.model.ropt(), self.nx,
                               self.model.nsamples[self.model.kopt], x_in_abs_coords=True)
 
         if self.do_logging:
@@ -820,12 +820,12 @@ def soft_restart(self, number_of_samples, nruns_so_far, params, x_in_abs_coords_
                 # Handle exit conditions (f < min obj value or maxfun reached)
                 if exit_info is not None:
                     if num_samples_run > 0:
-                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                        self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                               x_in_abs_coords=True)
                     return exit_info  # return & quit
 
                 # Otherwise, add new results
-                self.model.add_new_point(xnew, rvec_list[0, :])  # expect step, not absolute x
+                self.model.add_new_point(xnew, rvec_list[0, :], self.nx)  # expect step, not absolute x
                 for i in range(1, num_samples_run):
                     self.model.add_new_sample(self.model.npt() - 1, rvec_extra=rvec_list[i, :])
 
@@ -900,12 +900,12 @@ def move_furthest_points_momentum(self, d, number_of_samples, num_pts_to_move, p
             # Handle exit conditions (f < min obj value or maxfun reached)
             if exit_info is not None:
                 if num_samples_run > 0:
-                    self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                    self.model.save_point(x, np.mean(rvec_list[:num_samples_run, :], axis=0), num_samples_run, self.nx,
                                           x_in_abs_coords=True)
                 return exit_info  # return & quit
 
             # Otherwise, add new results
-            self.model.change_point(knew, xnew, rvec_list[0, :])  # expect step, not absolute x
+            self.model.change_point(knew, xnew, rvec_list[0, :], self.nx)  # expect step, not absolute x
             for i in range(1, num_samples_run):
                 self.model.add_new_sample(knew, rvec_extra=rvec_list[i, :])
         return None

dfols/model.py

@@ -85,6 +85,8 @@ def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, h=None, argsh=
         self.nsamples = np.zeros((npt,), dtype=int)  # number of samples used to evaluate objective at each point
         self.nsamples[0] = r0_nsamples
         self.objbeg = self.objval[0]  # f(x0), saved to check for sufficient reduction
+        self.eval_num = np.zeros((npt,), dtype=int)  # which evaluation number (1-indexed, nx not nf) is currently stored self.points[k,:]
+        self.eval_num[0] = 1
 
         # Termination criteria
         self.abs_tol = abs_tol
@@ -93,13 +95,16 @@ def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, h=None, argsh=
         # Model information
         self.model_const = np.zeros((m, ))  # constant term for model m(s) = c + J*s
         self.model_jac = np.zeros((m, n))  # Jacobian term for model m(s) = c + J*s
+        self.model_jac_eval_nums = np.zeros((npt,), dtype=int)  # which evaluation numbers (1-indexed, nx not nf) were used to build model_jac
 
         # Saved point (in absolute coordinates) - always check this value before quitting solver
         self.xsave = None
         self.rsave = None
         self.objsave = None
         self.jacsave = None
         self.nsamples_save = None
+        self.eval_num_save = None
+        self.jacsave_eval_nums = None
 
         # Factorisation of interpolation matrix
         self.factorisation_current = False
@@ -174,7 +179,7 @@ def distances_to_xopt(self):
             sq_distances[k] = sumsq(self.points[k, :] - xopt)
         return sq_distances
 
-    def change_point(self, k, x, rvec, allow_kopt_update=True):
+    def change_point(self, k, x, rvec, eval_num, allow_kopt_update=True):
         # Update point k to x (w.r.t. xbase), with residual values fvec
         if k >= self.npt_so_far and self.npt_so_far < self.num_pts:
             assert k == self.npt_so_far, "Growing: updating wrong point"
@@ -188,6 +193,7 @@ def change_point(self, k, x, rvec, allow_kopt_update=True):
         if self.h is not None:
             self.objval[k] += self.h(remove_scaling(self.xbase + x, self.scaling_changes), *self.argsh)
         self.nsamples[k] = 1
+        self.eval_num[k] = eval_num
         self.factorisation_current = False
 
         if allow_kopt_update and self.objval[k] < self.objopt():
@@ -198,6 +204,7 @@ def swap_points(self, k1, k2):
         self.points[[k1, k2], :] = self.points[[k2, k1], :]
         self.fval_v[[k1, k2], :] = self.fval_v[[k2, k1], :]
         self.objval[[k1, k2]] = self.objval[[k2, k1]]
+        self.eval_num[[k1, k2]] = self.eval_num[[k2, k1]]
         if self.kopt == k1:
             self.kopt = k2
         elif self.kopt == k2:
@@ -219,14 +226,15 @@ def add_new_sample(self, k, rvec_extra):
         self.kopt = np.argmin(self.objval[:self.npt()])  # make sure kopt is always the best value we have
         return
 
-    def add_new_point(self, x, rvec):
+    def add_new_point(self, x, rvec, eval_num):
         self.points = np.append(self.points, x.reshape((1, self.n())), axis=0)  # append row to xpt
         self.fval_v = np.append(self.fval_v, rvec.reshape((1, self.m())), axis=0)  # append row to fval_v
         obj = sumsq(rvec)
         if self.h is not None:
             obj += self.h(remove_scaling(self.xbase + x, self.scaling_changes), *self.argsh)
         self.objval = np.append(self.objval, obj)  # append entry to fval
         self.nsamples = np.append(self.nsamples, 1)  # add new sample number
+        self.eval_num = np.append(self.eval_num, eval_num)  # add new evaluation number
         self.num_pts += 1  # make sure npt is updated
         self.npt_so_far += 1
 
@@ -249,7 +257,7 @@ def shift_base(self, xbase_shift):
         self.model_const += np.dot(self.model_jac, xbase_shift)
         return
 
-    def save_point(self, x, rvec, nsamples, x_in_abs_coords=True):
+    def save_point(self, x, rvec, nsamples, eval_num, x_in_abs_coords=True):
         xabs = x.copy() if x_in_abs_coords else self.as_absolute_coordinates(x)
         obj = sumsq(rvec)
         if self.h is not None:
@@ -260,16 +268,18 @@ def save_point(self, x, rvec, nsamples, x_in_abs_coords=True):
             self.objsave = obj
             self.jacsave = self.model_jac.copy()
             self.nsamples_save = nsamples
+            self.eval_num_save = eval_num
+            self.jacsave_eval_nums = self.model_jac_eval_nums.copy()
             return True
         else:
             return False  # this value is worse than what we have already - didn't save
 
     def get_final_results(self):
         # Return x and objval for optimal point (either from xsave+objsave or kopt)
         if self.objsave is None or self.objopt() <= self.objsave:  # optimal has changed since xsave+objsave were last set
-            return self.xopt(abs_coordinates=True).copy(), self.ropt().copy(), self.objopt(), self.model_jac.copy(), self.nsamples[self.kopt]
+            return self.xopt(abs_coordinates=True).copy(), self.ropt().copy(), self.objopt(), self.model_jac.copy(), self.nsamples[self.kopt], self.eval_num[self.kopt], self.model_jac_eval_nums
         else:
-            return self.xsave.copy(), self.rsave.copy(), self.objsave, self.jacsave, self.nsamples_save
+            return self.xsave.copy(), self.rsave.copy(), self.objsave, self.jacsave, self.nsamples_save, self.eval_num_save, self.jacsave_eval_nums
 
     def min_objective_value(self):
         # Get termination criterion for f small: f <= abs_tol or f <= rel_tol * f0
@@ -367,6 +377,7 @@ def interpolate_mini_models_svd(self, verbose=False, make_full_rank=False, min_s
         J_old = self.model_jac.copy()
         self.model_jac = dg[1:,:].T
         self.model_const = dg[0,:] - np.dot(self.model_jac, xopt)  # shift base to xbase
+        self.model_jac_eval_nums = self.eval_num.copy()
         if verbose or get_chg_J:
             norm_J_error = np.linalg.norm(self.model_jac - J_old, ord='fro')**2
             linalg_resid = np.linalg.norm(W.dot(dg) - rhs)**2