From b904bc2994dcc0f67740e5425615b8de95651158 Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Mon, 11 Oct 2021 17:02:03 -0400
Subject: [PATCH 01/11] Handle arbitrary convex constraints in DFO-LS
---
dfols/__init__.py | 3 +-
dfols/controller.py | 133 +++++++++++++--
dfols/model.py | 11 +-
dfols/params.py | 11 ++
dfols/solver.py | 37 ++++-
dfols/tests/test_model.py | 30 ++--
dfols/tests/test_trust_region.py | 274 ++++++++++++++++++++++++++++++-
dfols/trust_region.py | 159 +++++++++++++++++-
dfols/util.py | 53 +++++-
dfols/version.py | 2 +-
10 files changed, 668 insertions(+), 45 deletions(-)
diff --git a/dfols/__init__.py b/dfols/__init__.py
index 74d7e05..5ae8a0e 100644
--- a/dfols/__init__.py
+++ b/dfols/__init__.py
@@ -7,8 +7,7 @@
It solves the nonlinear least-squares problem:
min_{x} f(x) = r1(x)**2 + ... + rm(x)**2,
-subject to the (optional) bounds
- lb <= x <= ub,
+(optionally) subject to finitely many convex constraints,
where each function ri(x) is differentiable, possibly nonconvex.
Since the derivatives of ri(x) are never required or approximated,
the solver works when the evaluation of ri(x) is noisy.
diff --git a/dfols/controller.py b/dfols/controller.py
index 0796858..5ace32c 100644
--- a/dfols/controller.py
+++ b/dfols/controller.py
@@ -92,13 +92,13 @@ def able_to_do_restart(self):
class Controller(object):
- def __init__(self, objfun, args, x0, r0, r0_nsamples, xl, xu, npt, rhobeg, rhoend, nf, nx, maxfun, params,
+ def __init__(self, objfun, args, x0, r0, r0_nsamples, xl, xu, projections, npt, rhobeg, rhoend, nf, nx, maxfun, params,
scaling_changes, do_logging):
self.do_logging = do_logging
self.objfun = objfun
self.args = args
self.maxfun = maxfun
- self.model = Model(npt, x0, r0, xl, xu, r0_nsamples, precondition=params("interpolation.precondition"),
+ self.model = Model(npt, x0, r0, xl, xu, projections, r0_nsamples, precondition=params("interpolation.precondition"),
abs_tol = params("model.abs_tol"), rel_tol = params("model.rel_tol"), do_logging=do_logging)
self.nf = nf
self.nx = nx
@@ -137,6 +137,107 @@ def initialise_coordinate_directions(self, number_of_samples, num_directions, pa
assert self.model.num_pts <= (self.n() + 1) * (self.n() + 2) // 2, "prelim: must have npt <= (n+1)(n+2)/2"
assert 1 <= num_directions < self.model.num_pts, "Initialisation: must have 1 <= ndirs_initial < npt"
+
+ if self.model.projections:
+ D = np.zeros((self.n(),self.n()))
+ k = 0
+ while k < self.n():
+ ek = np.zeros(self.n())
+ ek[k] = 1
+ p = np.dot(ek,self.delta)
+ yk = dykstra(self.model.projections, self.model.xbase + p, max_iter=params("dykstra.max_iters"), tol=params("dykstra.d_tol"))
+ D[k,:] = yk - self.model.xbase
+
+ k += 1 # move on to next point
+
+ # Have at least one L.D. vector, try negative direction on bad one first
+ k = 0
+ mr_tol = params("matrix_rank.r_tol")
+ D_rank, diag = qr_rank(D,tol=mr_tol)
+ while D_rank != num_directions and k < self.n():
+ if diag[k] < mr_tol:
+ ek = np.zeros(self.n())
+ ek[k] = 1
+ p = -np.dot(ek,self.delta)
+ yk = dykstra(self.model.projections, self.model.xbase + p, max_iter=params("dykstra.max_iters"), tol=params("dykstra.d_tol"))
+ dk = D[k,:].copy()
+ D[k,:] = yk - self.model.xbase
+ D_rank2, _diag2 = qr_rank(D,tol=params("matrix_rank.r_tol"))
+ if D_rank2 <= D_rank:
+ # Did not improve rank, revert change
+ D[k,:] = dk
+ # rank was improved, update D_rank for next comparison
+ D_rank = D_rank2
+ k += 1
+
+ # Try random combination of negatives...
+ k = 0
+ slctr = np.random.randint(0, 1+1, self.n()) # generate rand binary "selector" array
+ D_rank, diag = qr_rank(D,tol=params("matrix_rank.r_tol"))
+ while D_rank != num_directions and k < 100*self.n():
+ if slctr[k%self.n()] == 1: # if selector says make -ve, make -ve
+ ek = np.zeros(self.n())
+ ek[k%self.n()] = 1
+ p = -np.dot(ek,self.delta)
+ yk = dykstra(self.model.projections, self.model.xbase + p, max_iter=params("dykstra.max_iters"), tol=params("dykstra.d_tol"))
+ dk = D[k%self.n(),:].copy()
+ D[k%self.n(),:] = yk - self.model.xbase
+ D_rank2, _diag2 = qr_rank(D,tol=params("matrix_rank.r_tol"))
+ if D_rank2 <= D_rank:
+ # Did not improve rank, revert change
+ D[k%self.n(),:] = dk
+ # rank was improved, update D_rank for next comparison
+ D_rank = D_rank2
+
+ # Go again
+ slctr = np.random.randint(0, 1+1, self.n())
+ k += 1
+
+ # Set still not L.I? Try random directions
+ i = 0
+ D_rank, diag = qr_rank(D,tol=params("matrix_rank.r_tol"))
+ while D_rank != num_directions and i <= 100*num_directions:
+ k = 0
+ while k < self.n():
+ if diag[k] < mr_tol:
+ p = np.random.normal(size=self.n())
+ p = p/np.linalg.norm(p)
+ p = np.dot(p,self.delta)
+ yk = dykstra(self.model.projections, self.model.xbase + p, max_iter=params("dykstra.max_iters"), tol=params("dykstra.d_tol"))
+ dk = D[k,:].copy()
+ D[k,:] = yk - self.model.xbase
+ D_rank2, _diag2 = qr_rank(D,tol=params("matrix_rank.r_tol"))
+ if D_rank2 <= D_rank:
+ # Did not improve rank, revert change
+ D[k,:] = dk
+ # rank was improved, update D_rank for next comparison
+ D_rank = D_rank2
+ k += 1
+ i += 1
+
+ if D_rank != num_directions:
+ raise RuntimeError("Unable to generate suitable initial directions")
+
+ # we have a L.I set of interpolation points
+ for k in range(0,self.n()):
+ # Evaluate objective at this new point
+ x = self.model.as_absolute_coordinates(D[k, :])
+ rvec_list, f_list, num_samples_run, exit_info = self.evaluate_objective(x, 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,
+ 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
+ for i in range(1, num_samples_run):
+ self.model.add_new_sample(k+1, rvec_extra=rvec_list[i, :])
+
+ return None # return & continue
+
at_lower_boundary = (self.model.sl > -0.01 * self.delta) # sl = xl - x0, should be -ve, actually < -rhobeg
at_upper_boundary = (self.model.su < 0.01 * self.delta) # su = xu - x0, should be +ve, actually > rhobeg
@@ -147,17 +248,19 @@ def initialise_coordinate_directions(self, number_of_samples, num_directions, pa
# k = 2n+1, ..., (n+1)(n+2)/2 --> off-diagonal directions
if 1 <= k < self.n() + 1: # first step along coord directions
dirn = k - 1 # direction to move in (0,...,n-1)
- stepa = self.delta if not at_upper_boundary[dirn] else -self.delta
+ stepa = self.delta if not at_upper_boundary[dirn] else -self.delta # take a +delta step if at lower, -delta if at upper
stepb = None
- xpts_added[k, dirn] = stepa
+ xpts_added[k, dirn] = stepa # set new (relative) point to the step since we haven't done any moving, so relative point is all zeros.
elif self.n() + 1 <= k < 2 * self.n() + 1: # second step along coord directions
dirn = k - self.n() - 1 # direction to move in (0,...,n-1)
- stepa = xpts_added[k - self.n(), dirn]
- stepb = -self.delta
+ stepa = xpts_added[k - self.n(), dirn] # previous step
+ stepb = -self.delta # new step
if at_lower_boundary[dirn]:
+ # if at lower boundary, set the second step to be +ve
stepb = min(2.0 * self.delta, self.model.su[dirn]) # su = xu - x0, should be +ve
if at_upper_boundary[dirn]:
+ # if at upper boundary, set the second step to be -ve
stepb = max(-2.0 * self.delta, self.model.sl[dirn]) # sl = xl - x0, should be -ve
xpts_added[k, dirn] = stepb
@@ -325,10 +428,13 @@ def get_new_direction_for_growing(self, step_length):
return dirn * (step_length / LA.norm(dirn))
- def trust_region_step(self):
+ def trust_region_step(self, params):
# Build model for full least squares objectives
gopt, H = self.model.build_full_model()
- d, gnew, crvmin = trsbox(self.model.xopt(), gopt, H, self.model.sl, self.model.su, self.delta)
+ if self.model.projections:
+ d, gnew, crvmin = ctrsbox(self.model.xopt(abs_coordinates=True), gopt, H, self.model.projections, self.delta, d_max_iters=params("dykstra.max_iters"), d_tol=params("dykstra.d_tol"))
+ else:
+ d, gnew, crvmin = trsbox(self.model.xopt(), gopt, H, self.model.sl, self.model.su, self.delta)
return d, gopt, H, gnew, crvmin
def geometry_step(self, knew, adelt, number_of_samples, params):
@@ -337,8 +443,13 @@ def geometry_step(self, knew, adelt, number_of_samples, params):
try:
c, g = self.model.lagrange_gradient(knew)
# c = 1.0 if knew == self.model.kopt else 0.0 # based at xopt, just like d
- # Solve problem: bounds are sl <= xnew <= su, and ||xnew-xopt|| <= adelt
- xnew = trsbox_geometry(self.model.xopt(), c, g, np.minimum(self.model.sl, 0.0), np.maximum(self.model.su, 0.0), adelt)
+ if self.model.projections:
+ # Solve problem: use projection onto arbitrary constraints, and ||xnew-xopt|| <= adelt
+ step = ctrsbox_geometry(self.model.xopt(abs_coordinates=True), c, g, self.model.projections, adelt, d_max_iters=params("dykstra.max_iters"), d_tol=params("dykstra.d_tol"))
+ xnew = self.model.xopt() + step
+ else:
+ # Solve problem: bounds are sl <= xnew <= su, and ||xnew-xopt|| <= adelt
+ xnew = trsbox_geometry(self.model.xopt(), c, g, np.minimum(self.model.sl, 0.0), np.maximum(self.model.su, 0.0), adelt)
except LA.LinAlgError:
exit_info = ExitInformation(EXIT_LINALG_ERROR, "Singular matrix encountered in geometry step")
return exit_info # didn't fix geometry - return & quit
@@ -499,7 +610,7 @@ def reduce_rho(self, current_iter, params):
def calculate_ratio(self, current_iter, rvec_list, d, gopt, H):
exit_info = None
f = sumsq(np.mean(rvec_list, axis=0)) # estimate actual objective value
- pred_reduction = - model_value(gopt, H, d)
+ pred_reduction = - model_value(gopt, H, d) # negative of m since m(0) = 0
actual_reduction = self.model.fopt() - f
self.diffs = [abs(actual_reduction - pred_reduction), self.diffs[0], self.diffs[1]]
if min(sqrt(sumsq(d)), self.delta) > self.rho: # if ||d|| >= rho, successful!
diff --git a/dfols/model.py b/dfols/model.py
index 20373a9..745f2e6 100644
--- a/dfols/model.py
+++ b/dfols/model.py
@@ -36,12 +36,12 @@
import scipy.linalg as LA
from .trust_region import trsbox_geometry
-from .util import sumsq
+from .util import sumsq, dykstra
__all__ = ['Model']
class Model(object):
- def __init__(self, npt, x0, r0, xl, xu, r0_nsamples, n=None, m=None, abs_tol=1e-12, rel_tol=1e-20, precondition=True,
+ def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, n=None, m=None, abs_tol=1e-12, rel_tol=1e-20, precondition=True,
do_logging=True):
if n is None:
n = len(x0)
@@ -63,6 +63,7 @@ def __init__(self, npt, x0, r0, xl, xu, r0_nsamples, n=None, m=None, abs_tol=1e-
self.xbase = x0.copy()
self.sl = xl - self.xbase # lower bound w.r.t. xbase (require xpt >= sl)
self.su = xu - self.xbase # upper bound w.r.t. xbase (require xpt <= su)
+ self.projections = projections
self.points = np.zeros((npt, n)) # interpolation points w.r.t. xbase
# Function values
@@ -123,6 +124,8 @@ def xpt(self, k, abs_coordinates=False):
return np.minimum(np.maximum(self.sl, self.points[k, :].copy()), self.su)
else:
# Apply bounds and convert back to absolute coordinates
+ if self.projections:
+ return dykstra(self.projections, self.xbase + self.points[k,:])
return self.xbase + np.minimum(np.maximum(self.sl, self.points[k, :]), self.su)
def rvec(self, k):
@@ -133,8 +136,10 @@ def fval(self, k):
assert 0 <= k < self.npt(), "Invalid index %g" % k
return self.fval[k]
- def as_absolute_coordinates(self, x):
+ def as_absolute_coordinates(self, x, full_dykstra=False):
# If x were an interpolation point, get the absolute coordinates of x
+ if self.projections:
+ return dykstra(self.projections, self.xbase + x)
return self.xbase + np.minimum(np.maximum(self.sl, x), self.su)
def xpt_directions(self, include_kopt=True):
diff --git a/dfols/params.py b/dfols/params.py
index e3e2622..e168c01 100644
--- a/dfols/params.py
+++ b/dfols/params.py
@@ -109,6 +109,11 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
self.params["growing.full_rank.min_sing_val"] = 1e-6 # absolute floor on singular values
self.params["growing.full_rank.svd_max_jac_cond"] = 1e8 # maximum condition number of Jacobian
self.params["growing.perturb_trust_region_step"] = False # add random direction onto TRS solution?
+ # Dykstra's algorithm
+ self.params["dykstra.d_tol"] = 1e-10
+ self.params["dykstra.max_iters"] = 100
+ # Matrix rank algorithm
+ self.params["matrix_rank.r_tol"] = 1e-16
self.params_changed = {}
for p in self.params:
@@ -257,6 +262,12 @@ def param_type(self, key, npt):
type_str, nonetype_ok, lower, upper = 'float', True, 1.0, None
elif key == "growing.perturb_trust_region_step":
type_str, nonetype_ok, lower, upper = 'bool', False, None, None
+ elif key == "dykstra.d_tol":
+ type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
+ elif key == "dykstra.max_iters":
+ type_str, nonetype_ok, lower, upper = 'int', False, 0, None
+ elif key == "matrix_rank.r_tol":
+ type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
else:
assert False, "ParameterList.param_type() has unknown key: %s" % key
return type_str, nonetype_ok, lower, upper
diff --git a/dfols/solver.py b/dfols/solver.py
index f6279ba..cb487ef 100644
--- a/dfols/solver.py
+++ b/dfols/solver.py
@@ -93,7 +93,7 @@ def __str__(self):
return output
-def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_far, nf_so_far, nx_so_far, nsamples, params,
+def solve_main(objfun, x0, args, xl, xu, projections, npt, rhobeg, rhoend, maxfun, nruns_so_far, nf_so_far, nx_so_far, nsamples, params,
diagnostic_info, scaling_changes, r0_avg_old=None, r0_nsamples_old=None, default_growing_method_set_by_user=None,
do_logging=True, print_progress=False):
# Evaluate at x0 (keep nf, nx correct and check for f < 1e-12)
@@ -160,7 +160,7 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
params('growing.delta_scale_new_dirns', new_value=0.1)
# Initialise controller
- control = Controller(objfun, args, x0, r0_avg, num_samples_run, xl, xu, npt, rhobeg, rhoend, nf, nx, maxfun,
+ control = Controller(objfun, args, x0, r0_avg, num_samples_run, xl, xu, projections, npt, rhobeg, rhoend, nf, nx, maxfun,
params, scaling_changes, do_logging)
# Initialise interpolation set
@@ -271,7 +271,7 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
# Trust region step
- d, gopt, H, gnew, crvmin = control.trust_region_step()
+ d, gopt, H, gnew, crvmin = control.trust_region_step(params)
if do_logging:
logging.debug("Trust region step is d = " + str(d))
xnew = control.model.xopt() + d
@@ -851,7 +851,7 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
return x, rvec, f, jacmin, nsamples, control.nf, control.nx, nruns_so_far, exit_info, diagnostic_info
-def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8, maxfun=None, nsamples=None, user_params=None,
+def solve(objfun, x0, args=(), bounds=None, projections=[], npt=None, rhobeg=None, rhoend=1e-8, maxfun=None, nsamples=None, user_params=None,
objfun_has_noise=False, scaling_within_bounds=False, do_logging=True, print_progress=False):
x0 = x0.astype(float)
n = len(x0)
@@ -869,6 +869,10 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
scaling_within_bounds = False
warnings.warn("Ignoring scaling_within_bounds=True for unconstrained problem/1-sided bounds", RuntimeWarning)
+ if projections and scaling_within_bounds:
+ scaling_within_bounds = False
+ warnings.warn("Ignoring scaling_within_bounds=True for black-box constrained problem", RuntimeWarning)
+
if xl is None:
xl = -1e20 * np.ones((n,)) # unconstrained
if xu is None:
@@ -882,6 +886,18 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
if nsamples is None:
nsamples = lambda delta, rho, iter, nruns: 1 # no averaging
+ # If using arbitrary constraints, create projection from bounds
+ if projections:
+ xlb = xl.copy()
+ xub = xu.copy()
+ bproj = lambda w: pbox(w,xlb,xub)
+ projections = projections.copy()
+ projections.append(bproj)
+
+ # since using arbitrary constraints, don't constrain otherwise
+ xl = -1e20 * np.ones((n,))
+ xu = 1e20 * np.ones((n,))
+
# Set parameters
params = ParameterList(int(n), int(npt), int(maxfun), objfun_has_noise=objfun_has_noise) # make sure int, not np.int
if user_params is not None:
@@ -976,6 +992,13 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
results = OptimResults(None, None, None, None, 0, 0, 0, exit_flag, exit_msg)
return results
+ # Enforce arbitrary constraint bounds on x0
+ if projections:
+ xp = dykstra(projections,x0,max_iter=params("dykstra.max_iters"),tol=params("dykstra.d_tol"))
+ if not np.allclose(xp,x0):
+ warnings.warn("x0 not feasible w.r.t given constraints, adjusting", RuntimeWarning)
+ x0 = xp.copy()
+
# Enforce lower & upper bounds on x0
idx = (x0 <= xl)
if np.any(idx):
@@ -993,7 +1016,7 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
nf = 0
nx = 0
xmin, rmin, fmin, jacmin, nsamples_min, nf, nx, nruns, exit_info, diagnostic_info = \
- solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
+ solve_main(objfun, x0, args, xl, xu, projections, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
diagnostic_info, scaling_changes, default_growing_method_set_by_user=default_growing_method_set_by_user,
do_logging=do_logging, print_progress=print_progress)
@@ -1012,12 +1035,12 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
% (fmin, nf, rhobeg, rhoend))
if params("restarts.hard.use_old_rk"):
xmin2, rmin2, fmin2, jacmin2, nsamples2, nf, nx, nruns, exit_info, diagnostic_info = \
- solve_main(objfun, xmin, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
+ solve_main(objfun, xmin, args, xl, xu, projections, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
diagnostic_info, scaling_changes, r0_avg_old=rmin, r0_nsamples_old=nsamples_min,
do_logging=do_logging, print_progress=print_progress)
else:
xmin2, rmin2, fmin2, jacmin2, nsamples2, nf, nx, nruns, exit_info, diagnostic_info = \
- solve_main(objfun, xmin, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
+ solve_main(objfun, xmin, args, xl, xu, projections, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
diagnostic_info, scaling_changes, do_logging=do_logging, print_progress=print_progress)
if fmin2 < fmin or np.isnan(fmin):
diff --git a/dfols/tests/test_model.py b/dfols/tests/test_model.py
index 9b6a384..66ed89f 100644
--- a/dfols/tests/test_model.py
+++ b/dfols/tests/test_model.py
@@ -46,7 +46,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e20 * np.ones((n,))
xu = 1e20 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
self.assertEqual(model.npt(), 1, 'Wrong npt after initialisation')
self.assertTrue(array_compare(model.xopt(abs_coordinates=True), x0), 'Wrong xopt after initialisation')
self.assertTrue(array_compare(model.ropt(), rosenbrock(x0)), 'Wrong ropt after initialisation')
@@ -97,7 +97,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e20 * np.ones((n,))
xu = 1e20 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
# Now add better point
x1 = np.array([1.0, 0.9])
rvec = rosenbrock(x1)
@@ -124,7 +124,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
self.assertTrue(array_compare(model.sl, xl - x0), 'Wrong sl after initialisation')
self.assertTrue(array_compare(model.su, xu - x0), 'Wrong su after initialisation')
x1 = np.array([1.0, 0.9])
@@ -204,7 +204,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
x1 = np.array([1.0, 0.9])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
x2 = np.array([1.0, 1.0])
@@ -224,17 +224,17 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
x1 = np.array([1.0, 0.9])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
x2 = np.array([2.0, 0.9])
model.change_point(2, x2 - model.xbase, rosenbrock(x2))
self.assertAlmostEqual(model.min_objective_value(), 1e-12, 'Wrong min obj value')
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1, rel_tol=1e-2)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1, rel_tol=1e-2)
self.assertAlmostEqual(model.min_objective_value(), 1e-2 * sumsq(rosenbrock(x0)), 'Wrong min obj value 2')
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1, abs_tol=1.0)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1, abs_tol=1.0)
self.assertAlmostEqual(model.min_objective_value(), 1.0, 'Wrong min obj value 3')
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1, abs_tol=1.0, rel_tol=1e-2)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1, abs_tol=1.0, rel_tol=1e-2)
self.assertAlmostEqual(model.min_objective_value(), 1.0, 'Wrong min obj value 4')
@@ -245,7 +245,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
model.add_new_sample(0, rosenbrock(x0))
x1 = np.array([1.0, 0.9])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
@@ -295,7 +295,7 @@ def runTest(self):
# x0 = np.array([1.0, 2.9])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
model.add_new_sample(0, rosenbrock(x0))
x1 = np.array([1.0, 0.9])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
@@ -326,7 +326,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
model.add_new_sample(0, rosenbrock(x0))
x1 = np.array([-1.2, 0.9])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
@@ -354,7 +354,7 @@ def runTest(self):
delta = 0.5
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
model.add_new_sample(0, rosenbrock(x0))
x1 = x0 + delta * np.array([1.0, 0.0])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
@@ -373,7 +373,7 @@ def runTest(self):
x0 = np.array([-1.2, 1.0])
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+ model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
x1 = np.array([1.0, 0.9])
model.change_point(1, x1 - model.xbase, rosenbrock(x1))
x2 = np.array([2.0, 0.9])
@@ -407,7 +407,7 @@ def runTest(self):
# x0 = np.zeros((n,))
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
- model = Model(npt, A[0,:], np.array([b[0]]), xl, xu, 1)
+ model = Model(npt, A[0,:], np.array([b[0]]), xl, xu, [], 1)
for i in range(1,npt):
xi = A[i,:]
@@ -450,7 +450,7 @@ def runTest(self):
xl = -1e2 * np.ones((n,))
xu = 1e2 * np.ones((n,))
# model = Model(n+1, x0, np.array([0.0]), xl, xu, 1)
- model = Model(n+1, A[0, :], np.array([b[0]]), xl, xu, 1)
+ model = Model(n+1, A[0, :], np.array([b[0]]), xl, xu, [], 1)
for i in range(1,npt):
xi = A[i,:]
diff --git a/dfols/tests/test_trust_region.py b/dfols/tests/test_trust_region.py
index 64524e3..840d1f1 100644
--- a/dfols/tests/test_trust_region.py
+++ b/dfols/tests/test_trust_region.py
@@ -27,7 +27,7 @@
import numpy as np
import unittest
-from dfols.trust_region import trsbox, trsbox_geometry
+from dfols.trust_region import ctrsbox, ctrsbox_geometry, trsbox, trsbox_geometry
from dfols.util import model_value
@@ -70,6 +70,171 @@ def cauchy_pt_box(g, H, delta, lower, upper):
crvmin = -1.0
return s, red, crvmin
+def p_box(x,l,u):
+ return np.minimum(np.maximum(x,l), u)
+
+class TestUncInternalCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 2.0
+ xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
+ sl = -1e20 * np.ones((n,))
+ su = 1e20 * np.ones((n,))
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([-1.0, 0.0, -0.5])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer for internal minimum?
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+ print(crvmin)
+ self.assertAlmostEqual(crvmin, 1.2, 'Wrong crvmin')
+
+class TestUncBdryCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 5.0 / 12.0
+ xopt = np.zeros((n,))
+ sl = -1e20 * np.ones((n,))
+ su = 1e20 * np.ones((n,))
+ pbox = lambda x: p_box(x,sl,su)
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+ self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
+
+class TestUncBdry2CDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+ Delta = 5.0 / 12.0
+ xopt = np.zeros((n,))
+ sl = -1e20 * np.ones((n,))
+ su = 1e20 * np.ones((n,))
+ pbox = lambda x: p_box(x,sl,su)
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.allclose(gnew, g + H.dot(d)), 'Wrong gnew')
+ self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
+
+class TestUncBdry3CDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([0.0, 0.0, 1.0])
+ H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+ Delta = 0.5
+ xopt = np.zeros((n,))
+ sl = -1e20 * np.ones((n,))
+ su = 1e20 * np.ones((n,))
+ pbox = lambda x: p_box(x,sl,su)
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([0.0, 0.0, -0.5])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+ self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
+ # self.assertAlmostEqual(crvmin, crvmin_cauchy, 'Wrong crvmin')
+
+class TestUncHardCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([0.0, 0.0, 1.0])
+ H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+ Delta = sqrt(2.0)
+ xopt = np.zeros((n,))
+ sl = -1e20 * np.ones((n,))
+ su = 1e20 * np.ones((n,))
+ pbox = lambda x: p_box(x,sl,su)
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([1.0, 0.0, -1.0]) # non-unique solution
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+ self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
+
+class TestConInternalCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 2.0
+ xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
+ sl = xopt + np.array([-0.5, -10.0, -10.0])
+ su = xopt + np.array([10.0, 10.0, 10.0])
+ pbox = lambda x: p_box(x,sl,su)
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([-1.0, 0.0, -0.5])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer for internal minimum?
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl-xopt, su-xopt)
+ # print(s_cauchy)
+ # print(d)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+ print(crvmin)
+ self.assertAlmostEqual(crvmin, -1.0, 'Wrong crvmin')
+
+
+# Notes: Gets correct min and gnew. Incorrect crvmin (1.0999999999999999).
+class TestConBdryCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 5.0 / 12.0
+ xopt = np.zeros((n,))
+ sl = xopt + np.array([-0.3, -0.01, -0.1])
+ su = xopt + np.array([10.0, 1.0, 10.0])
+ pbox = lambda x: p_box(x,sl,su)
+ d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
+ true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ # Hope to get actual correct answer
+ # self.assertTrue(np.all(d == true_d), 'Wrong answer')
+ # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
+ s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl - xopt, su - xopt)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.max(np.abs(gnew - g - H.dot(d))) < 1e-10, 'Wrong gnew')
+ print(crvmin)
+ self.assertAlmostEqual(crvmin, -1.0, 'Wrong crvmin')
+ # self.assertAlmostEqual(crvmin, crvmin_cauchy, 'Wrong crvmin')
+
class TestUncInternal(unittest.TestCase):
def runTest(self):
@@ -340,3 +505,110 @@ def runTest(self):
# print(x)
# print(xtrue)
self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeomCDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0])
+ g = np.array([1.0, -1.0])
+ a = np.array([-2.0, -2.0])
+ b = np.array([1.0, 2.0])
+ delta = 2.0
+ c = -1.0
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([-sqrt(2.0), sqrt(2.0)])
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeom2CDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0])
+ g = np.array([1.0, -1.0])
+ a = np.array([-2.0, -2.0])
+ b = np.array([1.0, 2.0])
+ delta = 5.0
+ c = -1.0
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([-2.0, 2.0])
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeom3CDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0]) + 1
+ g = np.array([1.0, -1.0])
+ a = np.array([-2.0, -2.0]) + 1
+ b = np.array([1.0, 2.0]) + 1
+ delta = 5.0
+ c = 3.0 # may want to max instead
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([1.0, -2.0]) + 1
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeomOldBugCDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0])
+ g = np.array([-1.0, -1.0])
+ a = np.array([-2.0, -2.0])
+ b = np.array([0.1, 0.9])
+ delta = sqrt(2.0)
+ c = -1.0 # may want to max instead
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = b
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+ # self.assertFalse(True, "bad")
+
+
+class TestGeomOldBug2CDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0, 0.0])
+ g = np.array([-1.0, -1.0, -1.0])
+ a = np.array([-2.0, -2.0, -2.0])
+ b = np.array([0.9, 0.1, 5.0])
+ delta = sqrt(3.0)
+ c = -1.0 # may want to max instead
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([0.9, 0.1, sqrt(3.0 - 0.81 - 0.01)])
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+ # self.assertFalse(True, "bad")
+
+class TestGeom2WithZerosCDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0])
+ g = np.array([0.0, -1.0])
+ a = np.array([-2.0, -2.0])
+ b = np.array([1.0, 2.0])
+ delta = 5.0
+ c = 0.0
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([0.0, 2.0])
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeom2WithAlmostZerosCDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0])
+ g = np.array([1e-15, -1.0])
+ a = np.array([-2.0, -2.0])
+ b = np.array([1.0, 2.0])
+ delta = 5.0
+ c = 0.0
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([0.0, 2.0])
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeom2WithAlmostZeros2CDFO(unittest.TestCase):
+ def runTest(self):
+ xbase = np.array([0.0, 0.0])
+ g = np.array([1e-15, 0.0])
+ a = np.array([-2.0, -2.0])
+ b = np.array([1.0, 2.0])
+ delta = 5.0
+ c = 0.0
+ x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
+ xtrue = np.array([0.0, 0.0])
+ # print(x)
+ # print(xtrue)
+ self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
diff --git a/dfols/trust_region.py b/dfols/trust_region.py
index ec8c975..8b0a21d 100644
--- a/dfols/trust_region.py
+++ b/dfols/trust_region.py
@@ -11,6 +11,15 @@
The other outputs: gnew is the gradient of the model at d, and crvmin has
information about the curvature of the model at the solution.
+For handling arbitrary constraints, the call is
+ d, gnew, crvmin = ctrsbox(xopt, g, H, projections, delta)
+which produces a new vector d approximately solving the constrained trust region subproblem:
+ min_{d} g'*d + 0.5*d'*H*d
+ s.t. ||d|| <= delta
+ xopt + d is feasible w.r.t. the constraint set C
+The other outputs: gnew is the gradient of the model at d, and crvmin has
+information about the curvature of the model at the solution.
+
We also provide a function for maximising the absolute value of a linear function
inside a similar trust region - this is useful for geometry steps.
The call
@@ -23,6 +32,13 @@
min_s abs(c + g' * d)
s.t. lower <= xbase + d <= upper
||d|| <= delta
+Again, we have a version of this for handling arbitrary constraints
+The call
+ x = ctrsbox_geometry(xbase, c, g, projections, Delta)
+Solves
+ min_s abs(c + g' * d)
+ s.t. xbase + d is feasible w.r.t. the constraint set C
+ ||d|| <= delta
Notes
----
@@ -63,13 +79,77 @@
# Fall back to Python implementation
USE_FORTRAN = False
+from .util import dykstra, pball, pbox, sumsq, model_value
-from .util import sumsq
+__all__ = ['ctrsbox', 'ctrsbox_geometry', 'trsbox', 'trsbox_geometry']
+ZERO_THRESH = 1e-14
-__all__ = ['trsbox', 'trsbox_geometry']
+def ctrsbox(xopt, g, H, projections, delta, d_max_iters=100, d_tol=1e-10, use_fortran=USE_FORTRAN):
+ n = xopt.size
+ assert xopt.shape == (n,), "xopt has wrong shape (should be vector)"
+ assert g.shape == (n,), "g and xopt have incompatible sizes"
+ assert len(H.shape) == 2, "H must be a matrix"
+ assert H.shape == (n,n), "H and xopt have incompatible sizes"
+ assert np.allclose(H, H.T), "H must be symmetric"
+ assert delta > 0.0, "delta must be strictly positive"
-ZERO_THRESH = 1e-14
+ d = np.zeros((n,))
+ gnew = g.copy()
+ gy = g.copy()
+ crvmin = -1.0
+ y = d.copy()
+ eta = 1.2 # L backtrack scaling factor
+ t = 1
+
+ # Initial guess of L is norm(Hessian)
+ L = np.linalg.norm(H, 2)
+
+ # trust region is a ball of radius delta around xopt
+ trproj = lambda w: pball(w, xopt, delta)
+
+ # combine trust region constraints with user-entered constraints
+ P = projections.copy()
+ P.append(trproj)
+ def proj(d0):
+ p = dykstra(P, xopt+d0, max_iter=d_max_iters, tol=d_tol)
+ # we want the step only, so we subtract xopt
+ # from the new point: proj(xk+d) - xk
+ return p - xopt
+
+ MAX_LOOP_ITERS = 100 * n ** 2
+
+ # projected GD loop
+ for ii in range(MAX_LOOP_ITERS):
+ w = y - (1/L)*gy
+ prev_d = d.copy()
+ d = proj(w)
+
+ # size of step taken
+ s = d - prev_d
+ stplen = np.linalg.norm(s)
+
+ # update true gradient
+ gnew += H.dot(s)
+
+ # update CRVMIN
+ crv = s.dot(H).dot(s)/sumsq(s)
+ crvmin = min(crvmin, crv) if crvmin != -1.0 else crv
+
+ # exit condition
+ if stplen <= ZERO_THRESH:
+ break
+
+ # momentum update
+ prev_t = t
+ t = (1 + np.sqrt(1 + 4 * t ** 2))/2
+ prev_y = y.copy()
+ y = d + s*(prev_t - 1)/t
+
+ # update gradient w.r.t y
+ gy += H.dot(y - prev_y)
+
+ return d, gnew, crvmin
def trsbox(xopt, g, H, sl, su, delta, use_fortran=USE_FORTRAN):
@@ -405,8 +485,63 @@ def ball_step(x0, g, Delta):
if sqrt(gsqnorm) < ZERO_THRESH: # Error catching: if g=0, make no step
return 0.0
else:
- return (sqrt(gdotx0**2 + gsqnorm*(Delta**2 - x0sqnorm)) - gdotx0) / gsqnorm
+ # Sqrt had negative input on prob 46 in OG DFOLS with noise
+ # print("Inside of the sqrt:", gdotx0**2 + gsqnorm*(Delta**2 - x0sqnorm))
+ # Got Inside of the sqrt: -3.608971127647144e-42
+ # Added max(0,...) here
+ return (sqrt(np.maximum(0,gdotx0**2 + gsqnorm*(Delta**2 - x0sqnorm))) - gdotx0) / gsqnorm
+
+def ctrsbox_linear(xbase, g, projections, Delta, d_max_iters=100, d_tol=1e-10, use_fortran=USE_FORTRAN):
+ # Solve the convex program:
+ # min_d g' * d
+ # s.t. xbase + d is feasible w.r.t. constraint set C
+ # ||d||^2 <= Delta^2
+
+ n = g.size
+ d = np.zeros((n,))
+ y = d.copy()
+ t = 1
+ dirn = -g
+ cons_dirns = []
+
+ # If g[i] = 0, never step along this direction
+ constant_directions = np.where(np.abs(dirn) < ZERO_THRESH)[0]
+ dirn[constant_directions] = 0.0
+
+ # trust region is a ball of radius delta centered around xbase
+ trproj = lambda w: pball(w, xbase, Delta)
+ # combine trust region constraints with user-entered constraints
+ P = projections.copy()
+ P.append(trproj)
+ def proj(d0):
+ p = dykstra(P, xbase + d0, max_iter=d_max_iters, tol=d_tol)
+ # we want the step only, so we subtract
+ # xbase from the new point: proj(xk + d) - xk
+ return p - xbase
+
+ MAX_LOOP_ITERS = 100 * n ** 2
+
+ # projected GD loop
+ for ii in range(MAX_LOOP_ITERS):
+ w = y + dirn
+ prev_d = d.copy()
+ d = proj(w)
+
+ s = d - prev_d
+ stplen = np.linalg.norm(s)
+
+ # exit condition
+ if stplen <= ZERO_THRESH:
+ break
+
+ # 'momentum' update
+ prev_t = t
+ t = (1 + np.sqrt(1 + 4 * t ** 2))/2
+ prev_y = y.copy()
+ y = d + s*(prev_t - 1)/t
+
+ return d
def trsbox_linear(g, a_in, b_in, Delta, use_fortran=USE_FORTRAN):
# Solve the convex program:
@@ -466,6 +601,22 @@ def trsbox_linear(g, a_in, b_in, Delta, use_fortran=USE_FORTRAN):
dirn[idx_hit] = 0.0 # no more searching this direction
return x
+def ctrsbox_geometry(xbase, c, g, projections, Delta, d_max_iters=100, d_tol=1e-10, use_fortran=USE_FORTRAN):
+ # Given a Lagrange polynomial defined by: L(x) = c + g' * (x - xbase)
+ # Maximise |L(x)| in a box + trust region - that is, solve:
+ # max_x abs(c + g' * (x - xbase))
+ # s.t. x is feasible w.r.t constraint set C
+ # ||x-xbase|| <= Delta
+ # Setting s = x-xbase (or x = xbase + s), this is equivalent to:
+ # max_s abs(c + g' * s)
+ # s.t. xbase + s is is feasible w.r.t constraint set C
+ # ||s|| <= Delta
+ smin = ctrsbox_linear(xbase, g, projections, Delta, d_max_iters=100, d_tol=1e-10, use_fortran=use_fortran) # minimise g' * s
+ smax = ctrsbox_linear(xbase, -g, projections, Delta, d_max_iters=100, d_tol=1e-10, use_fortran=use_fortran) # maximise g' * s
+ if abs(c + np.dot(g, smin)) >= abs(c + np.dot(g, smax)): # choose the one with largest absolute value
+ return smin
+ else:
+ return smax
def trsbox_geometry(xbase, c, g, lower, upper, Delta, use_fortran=USE_FORTRAN):
# Given a Lagrange polynomial defined by: L(x) = c + g' * (x - xbase)
diff --git a/dfols/util.py b/dfols/util.py
index 6107e09..154e041 100644
--- a/dfols/util.py
+++ b/dfols/util.py
@@ -27,11 +27,14 @@
import logging
import numpy as np
+import scipy.linalg as LA
import sys
+import pdb
+
__all__ = ['sumsq', 'eval_least_squares_objective', 'model_value', 'random_orthog_directions_within_bounds',
- 'random_directions_within_bounds', 'apply_scaling', 'remove_scaling']
+ 'random_directions_within_bounds', 'apply_scaling', 'remove_scaling', 'pbox', 'pball', 'dykstra', 'qr_rank']
def sumsq(x):
@@ -207,3 +210,51 @@ def remove_scaling(x_scaled, scaling_changes):
shift, scale = scaling_changes
return shift + x_scaled * scale
+
+def dykstra(P,x0,max_iter=100,tol=1.0e-10):
+ x = x0.copy()
+ p = len(P)
+ y = np.zeros((p,x0.shape[0]))
+
+ n = 0
+ cI = float('inf')
+ while n < max_iter and cI >= tol:
+ cI = 0
+ for i in range(0,p):
+ # Update iterate
+ prev_x = x.copy()
+ x = P[i](prev_x - y[i,:])
+
+ # Update increment
+ # pdb.set_trace()
+ prev_y = y[i,:].copy()
+ y[i,:] = x - (prev_x - prev_y)
+
+ # Stop condition
+ cI += np.linalg.norm(prev_y - y[i,:])**2
+
+ n += 1
+
+ return x
+
+
+def pball(x,c,r):
+ return c + (r/np.max([np.linalg.norm(x-c),r]))*(x-c)
+
+
+def pbox(x,l,u):
+ return np.minimum(np.maximum(x,l), u)
+
+'''
+Calculates rank of square matrix with QR.
+We use the fact that the rank of a square matrix A
+can be given by the number of nonzero diagonal elements of
+R in the QR factorization of A.
+'''
+def qr_rank(A,tol=1e-25):
+ m,n = A.shape
+ assert m == n, "Input matrix must be square"
+ Q,R = LA.qr(A)
+ D = np.abs(np.diag(R))
+ rank = np.sum(D > tol)
+ return rank, D
diff --git a/dfols/version.py b/dfols/version.py
index 94b452c..dc26bd5 100644
--- a/dfols/version.py
+++ b/dfols/version.py
@@ -22,4 +22,4 @@
"""
-__version__ = '1.2.3'
+__version__ = '1.3.0'
From d32b9523e6fd40168707993ae317b402121527cb Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Mon, 11 Oct 2021 17:04:58 -0400
Subject: [PATCH 02/11] Change error message to be more suitable
---
dfols/solver.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/dfols/solver.py b/dfols/solver.py
index cb487ef..0c7ee60 100644
--- a/dfols/solver.py
+++ b/dfols/solver.py
@@ -871,7 +871,7 @@ def solve(objfun, x0, args=(), bounds=None, projections=[], npt=None, rhobeg=Non
if projections and scaling_within_bounds:
scaling_within_bounds = False
- warnings.warn("Ignoring scaling_within_bounds=True for black-box constrained problem", RuntimeWarning)
+ warnings.warn("Ignoring scaling_within_bounds=True for problem with arbitrary constraints", RuntimeWarning)
if xl is None:
xl = -1e20 * np.ones((n,)) # unconstrained
From 0990ef7b9c3f039ba896ae04f0c102e107394659 Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Tue, 12 Oct 2021 15:38:28 -0400
Subject: [PATCH 03/11] Add tests for TR subproblem. Make curvature calc more
robust.
---
dfols/tests/test_trust_region.py | 354 ++++++++++---------------------
dfols/trust_region.py | 2 +-
2 files changed, 112 insertions(+), 244 deletions(-)
diff --git a/dfols/tests/test_trust_region.py b/dfols/tests/test_trust_region.py
index 840d1f1..d2dd2fd 100644
--- a/dfols/tests/test_trust_region.py
+++ b/dfols/tests/test_trust_region.py
@@ -30,6 +30,8 @@
from dfols.trust_region import ctrsbox, ctrsbox_geometry, trsbox, trsbox_geometry
from dfols.util import model_value
+import pdb
+
def cauchy_pt(g, H, delta):
# General expression for the Cauchy point
@@ -73,167 +75,8 @@ def cauchy_pt_box(g, H, delta, lower, upper):
def p_box(x,l,u):
return np.minimum(np.maximum(x,l), u)
-class TestUncInternalCDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([1.0, 0.0, 1.0])
- H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
- Delta = 2.0
- xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
- sl = -1e20 * np.ones((n,))
- su = 1e20 * np.ones((n,))
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([-1.0, 0.0, -0.5])
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer for internal minimum?
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
- print(crvmin)
- self.assertAlmostEqual(crvmin, 1.2, 'Wrong crvmin')
-
-class TestUncBdryCDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([1.0, 0.0, 1.0])
- H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
- Delta = 5.0 / 12.0
- xopt = np.zeros((n,))
- sl = -1e20 * np.ones((n,))
- su = 1e20 * np.ones((n,))
- pbox = lambda x: p_box(x,sl,su)
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
- self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
-
-class TestUncBdry2CDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([1.0, 0.0, 1.0])
- H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
- Delta = 5.0 / 12.0
- xopt = np.zeros((n,))
- sl = -1e20 * np.ones((n,))
- su = 1e20 * np.ones((n,))
- pbox = lambda x: p_box(x,sl,su)
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.allclose(gnew, g + H.dot(d)), 'Wrong gnew')
- self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
-
-class TestUncBdry3CDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([0.0, 0.0, 1.0])
- H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
- Delta = 0.5
- xopt = np.zeros((n,))
- sl = -1e20 * np.ones((n,))
- su = 1e20 * np.ones((n,))
- pbox = lambda x: p_box(x,sl,su)
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([0.0, 0.0, -0.5])
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
- self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
- # self.assertAlmostEqual(crvmin, crvmin_cauchy, 'Wrong crvmin')
-
-class TestUncHardCDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([0.0, 0.0, 1.0])
- H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
- Delta = sqrt(2.0)
- xopt = np.zeros((n,))
- sl = -1e20 * np.ones((n,))
- su = 1e20 * np.ones((n,))
- pbox = lambda x: p_box(x,sl,su)
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([1.0, 0.0, -1.0]) # non-unique solution
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt(g, H, Delta)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
- self.assertAlmostEqual(crvmin, 0.0, 'Wrong crvmin')
-
-class TestConInternalCDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([1.0, 0.0, 1.0])
- H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
- Delta = 2.0
- xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
- sl = xopt + np.array([-0.5, -10.0, -10.0])
- su = xopt + np.array([10.0, 10.0, 10.0])
- pbox = lambda x: p_box(x,sl,su)
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([-1.0, 0.0, -0.5])
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer for internal minimum?
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl-xopt, su-xopt)
- # print(s_cauchy)
- # print(d)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
- print(crvmin)
- self.assertAlmostEqual(crvmin, -1.0, 'Wrong crvmin')
-
-
-# Notes: Gets correct min and gnew. Incorrect crvmin (1.0999999999999999).
-class TestConBdryCDFO(unittest.TestCase):
- def runTest(self):
- n = 3
- g = np.array([1.0, 0.0, 1.0])
- H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
- Delta = 5.0 / 12.0
- xopt = np.zeros((n,))
- sl = xopt + np.array([-0.3, -0.01, -0.1])
- su = xopt + np.array([10.0, 1.0, 10.0])
- pbox = lambda x: p_box(x,sl,su)
- d, gnew, crvmin = ctrsbox(xopt, g, H, sl, su, [], Delta)
- true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
- est_min = model_value(g, H, d)
- true_min = model_value(g, H, true_d)
- # Hope to get actual correct answer
- # self.assertTrue(np.all(d == true_d), 'Wrong answer')
- # self.assertAlmostEqual(est_min, true_min, 'Wrong min value')
- s_cauchy, red_cauchy, crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl - xopt, su - xopt)
- self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
- self.assertTrue(np.max(np.abs(gnew - g - H.dot(d))) < 1e-10, 'Wrong gnew')
- print(crvmin)
- self.assertAlmostEqual(crvmin, -1.0, 'Wrong crvmin')
- # self.assertAlmostEqual(crvmin, crvmin_cauchy, 'Wrong crvmin')
+def p_ball(x,c,r):
+ return c + (r/np.max([np.linalg.norm(x-c),r]))*(x-c)
class TestUncInternal(unittest.TestCase):
@@ -507,108 +350,133 @@ def runTest(self):
self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
-class TestGeomCDFO(unittest.TestCase):
- def runTest(self):
- xbase = np.array([0.0, 0.0])
- g = np.array([1.0, -1.0])
- a = np.array([-2.0, -2.0])
- b = np.array([1.0, 2.0])
- delta = 2.0
- c = -1.0
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([-sqrt(2.0), sqrt(2.0)])
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
-
-
-class TestGeom2CDFO(unittest.TestCase):
- def runTest(self):
- xbase = np.array([0.0, 0.0])
- g = np.array([1.0, -1.0])
- a = np.array([-2.0, -2.0])
- b = np.array([1.0, 2.0])
- delta = 5.0
- c = -1.0
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([-2.0, 2.0])
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
-
-
+# DFO-LS with arbitrary constraints
class TestGeom3CDFO(unittest.TestCase):
def runTest(self):
xbase = np.array([0.0, 0.0]) + 1
g = np.array([1.0, -1.0])
a = np.array([-2.0, -2.0]) + 1
b = np.array([1.0, 2.0]) + 1
+ proj = lambda x: p_box(x,a,b)
delta = 5.0
c = 3.0 # may want to max instead
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([1.0, -2.0]) + 1
+ x = ctrsbox_geometry(xbase, c, g, [proj], delta)
+ xtrue = np.array([1.0, -2.0])
self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+class TestUncInternalCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 2.0
+ xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [], Delta)
+ true_d = np.array([-1.0, 0.0, -0.5])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
-class TestGeomOldBugCDFO(unittest.TestCase):
+class TestUncBdryCDFO(unittest.TestCase):
def runTest(self):
- xbase = np.array([0.0, 0.0])
- g = np.array([-1.0, -1.0])
- a = np.array([-2.0, -2.0])
- b = np.array([0.1, 0.9])
- delta = sqrt(2.0)
- c = -1.0 # may want to max instead
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = b
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
- # self.assertFalse(True, "bad")
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 5.0 / 12.0
+ xopt = np.zeros((n,))
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [], Delta)
+ true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.allclose(gnew, g + H.dot(d)), 'Wrong gnew')
-class TestGeomOldBug2CDFO(unittest.TestCase):
+class TestUncHardCDFO(unittest.TestCase):
def runTest(self):
- xbase = np.array([0.0, 0.0, 0.0])
- g = np.array([-1.0, -1.0, -1.0])
- a = np.array([-2.0, -2.0, -2.0])
- b = np.array([0.9, 0.1, 5.0])
- delta = sqrt(3.0)
- c = -1.0 # may want to max instead
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([0.9, 0.1, sqrt(3.0 - 0.81 - 0.01)])
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
- # self.assertFalse(True, "bad")
+ n = 3
+ g = np.array([0.0, 0.0, 1.0])
+ H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+ Delta = sqrt(2.0)
+ xopt = np.zeros((n,))
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [], Delta)
+ true_d = np.array([1.0, 0.0, -1.0]) # non-unique solution
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt(g, H, Delta)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.allclose(gnew, g + H.dot(d)), 'Wrong gnew')
-class TestGeom2WithZerosCDFO(unittest.TestCase):
+class TestConInternalCDFO(unittest.TestCase):
def runTest(self):
- xbase = np.array([0.0, 0.0])
- g = np.array([0.0, -1.0])
- a = np.array([-2.0, -2.0])
- b = np.array([1.0, 2.0])
- delta = 5.0
- c = 0.0
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([0.0, 2.0])
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 2.0
+ xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
+ sl = xopt + np.array([-0.5, -10.0, -10.0])
+ su = xopt + np.array([10.0, 10.0, 10.0])
+ proj = lambda x: p_box(x,sl,su)
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [proj], Delta)
+ true_d = np.array([-1.0, 0.0, -0.5])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl-xopt, su-xopt)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+class TestConBdryCDFO(unittest.TestCase):
+ def runTest(self):
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 5.0 / 12.0
+ xopt = np.zeros((n,))
+ sl = xopt + np.array([-0.3, -0.01, -0.1])
+ su = xopt + np.array([10.0, 1.0, 10.0])
+ proj = lambda x: p_box(x,sl,su)
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [proj], Delta)
+ true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl - xopt, su - xopt)
+ self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+ self.assertTrue(np.max(np.abs(gnew - g - H.dot(d))) < 1e-10, 'Wrong gnew')
-class TestGeom2WithAlmostZerosCDFO(unittest.TestCase):
+class TestBoxBallInternalCDFO(unittest.TestCase):
def runTest(self):
- xbase = np.array([0.0, 0.0])
- g = np.array([1e-15, -1.0])
- a = np.array([-2.0, -2.0])
- b = np.array([1.0, 2.0])
- delta = 5.0
- c = 0.0
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([0.0, 2.0])
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 2.0
+ xopt = np.ones((n,)) # trying nonzero (since bounds inactive)
+ sl = xopt + np.array([-0.5, -10.0, -10.0])
+ su = xopt + np.array([10.0, 10.0, 10.0])
+ boxproj = lambda x: p_box(x,sl,su)
+ ballproj = lambda x: p_ball(x,xopt,5)
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [boxproj,ballproj], Delta)
+ true_d = np.array([-0.5, 0.0, -0.5])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ self.assertTrue(est_min <= true_min + 1e-3, 'Sufficient decrease not achieved')
-class TestGeom2WithAlmostZeros2CDFO(unittest.TestCase):
+class TestBoxBallBdryCDFO(unittest.TestCase):
def runTest(self):
- xbase = np.array([0.0, 0.0])
- g = np.array([1e-15, 0.0])
- a = np.array([-2.0, -2.0])
- b = np.array([1.0, 2.0])
- delta = 5.0
- c = 0.0
- x = ctrsbox_geometry(xbase, c, g, a, b, [], delta)
- xtrue = np.array([0.0, 0.0])
- # print(x)
- # print(xtrue)
- self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+ n = 3
+ g = np.array([1.0, 0.0, 1.0])
+ H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+ Delta = 5.0 / 12.0
+ xopt = np.zeros((n,))
+ sl = xopt + np.array([-0.3, -0.01, -0.1])
+ su = xopt + np.array([10.0, 1.0, 10.0])
+ boxproj = lambda x: p_box(x,sl,su)
+ ballproj = lambda x: p_ball(x,xopt,0.25)
+ d, gnew, _crvmin = ctrsbox(xopt, g, H, [boxproj,ballproj], Delta)
+ true_d = np.array([-0.22913085, 0.0, -0.09999527])
+ est_min = model_value(g, H, d)
+ true_min = model_value(g, H, true_d)
+ self.assertTrue(est_min <= true_min + 1e-3, 'Sufficient decrease not achieved')
diff --git a/dfols/trust_region.py b/dfols/trust_region.py
index 8b0a21d..74ce269 100644
--- a/dfols/trust_region.py
+++ b/dfols/trust_region.py
@@ -133,7 +133,7 @@ def proj(d0):
gnew += H.dot(s)
# update CRVMIN
- crv = s.dot(H).dot(s)/sumsq(s)
+ crv = s.dot(H).dot(s)/sumsq(s) if sumsq(s) >= ZERO_THRESH else crvmin
crvmin = min(crvmin, crv) if crvmin != -1.0 else crv
# exit condition
From 3ba2e3f965f6c5752912abfe962f1ca32e5cb614 Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Tue, 12 Oct 2021 16:03:28 -0400
Subject: [PATCH 04/11] Add tests for solver and remove pdb in TR tests
---
dfols/tests/test_solver.py | 25 +++++++++++++++++++++++++
dfols/tests/test_trust_region.py | 2 --
2 files changed, 25 insertions(+), 2 deletions(-)
diff --git a/dfols/tests/test_solver.py b/dfols/tests/test_solver.py
index 9c43445..b0601ad 100644
--- a/dfols/tests/test_solver.py
+++ b/dfols/tests/test_solver.py
@@ -41,6 +41,13 @@ def rosenbrock_jacobian(x):
return np.array([[-20.0*x[0], 10.0], [-1.0, 0.0]])
+def p_box(x,l,u):
+ return np.minimum(np.maximum(x,l), u)
+
+def p_ball(x,c,r):
+ return c + (r/np.max([np.linalg.norm(x-c),r]))*(x-c)
+
+
class TestNans(unittest.TestCase):
# Generic objective that only returns NaNs (like optclim code)
# Verify get a sensible termination
@@ -185,3 +192,21 @@ def runTest(self):
self.assertTrue(array_compare(soln.jacobian, jac(soln.x), thresh=1e-1), "Wrong Jacobian")
self.assertTrue(abs(soln.f) < 1e-10, "Wrong fmin")
+
+class TestRosenbrockBoxBall(unittest.TestCase):
+ # Minimise the (2d) Rosenbrock function, where x[1] hits the upper bound
+ def runTest(self):
+ # n, m = 2, 2
+ x0 = np.array([-1.2, 0.7]) # standard start point does not satisfy bounds
+ lower = np.array([0.7, -2.0])
+ upper = np.array([1.0, 2])
+ boxproj = lambda x: p_box(x,lower,upper)
+ ballproj = lambda x: p_ball(x,np.array([0.5,1]),0.25)
+ xmin = np.array([0.70424386, 0.85583188]) # approximate
+ fmin = np.dot(rosenbrock(xmin), rosenbrock(xmin))
+ soln = dfols.solve(rosenbrock, x0, projections=[boxproj,ballproj])
+ print(soln.x)
+ self.assertTrue(array_compare(soln.x, xmin, thresh=1e-2), "Wrong xmin")
+ self.assertTrue(array_compare(soln.resid, rosenbrock(soln.x), thresh=1e-10), "Wrong resid")
+ self.assertTrue(array_compare(soln.jacobian, rosenbrock_jacobian(soln.x), thresh=1e-2), "Wrong Jacobian")
+ self.assertTrue(abs(soln.f - fmin) < 1e-4, "Wrong fmin")
diff --git a/dfols/tests/test_trust_region.py b/dfols/tests/test_trust_region.py
index d2dd2fd..b9561c6 100644
--- a/dfols/tests/test_trust_region.py
+++ b/dfols/tests/test_trust_region.py
@@ -30,8 +30,6 @@
from dfols.trust_region import ctrsbox, ctrsbox_geometry, trsbox, trsbox_geometry
from dfols.util import model_value
-import pdb
-
def cauchy_pt(g, H, delta):
# General expression for the Cauchy point
From 58104c7bf9bfeef278051f1cf70dcdc7f83658c8 Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Tue, 12 Oct 2021 17:16:25 -0400
Subject: [PATCH 05/11] Add more util tests and remove more pdb
---
dfols/params.py | 2 +-
dfols/tests/test_util.py | 133 ++++++++++++++++++++++++++++++++++++++-
dfols/util.py | 7 +--
3 files changed, 135 insertions(+), 7 deletions(-)
diff --git a/dfols/params.py b/dfols/params.py
index e168c01..6d36e2e 100644
--- a/dfols/params.py
+++ b/dfols/params.py
@@ -113,7 +113,7 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
self.params["dykstra.d_tol"] = 1e-10
self.params["dykstra.max_iters"] = 100
# Matrix rank algorithm
- self.params["matrix_rank.r_tol"] = 1e-16
+ self.params["matrix_rank.r_tol"] = 1e-15
self.params_changed = {}
for p in self.params:
diff --git a/dfols/tests/test_util.py b/dfols/tests/test_util.py
index f497904..57ecfab 100644
--- a/dfols/tests/test_util.py
+++ b/dfols/tests/test_util.py
@@ -165,4 +165,135 @@ def runTest(self):
self.assertTrue(np.linalg.norm(dirns[i, :]) <= delta + 1e-10, "Unconstrained: dirn %i too long" % i)
self.assertTrue(np.all(dirns[i, :] >= lower), "Direction %i below lower bound" % i)
self.assertTrue(np.all(dirns[i, :] <= upper), "Direction %i above upper bound" % i)
- # self.assertTrue(False, "bad")
\ No newline at end of file
+ # self.assertTrue(False, "bad")
+
+# Trivial case of full rank
+class TestMatrixRankQR1(unittest.TestCase):
+ def runTest(self):
+ mr_tol = 1e-15
+ A = np.array([
+ [1,0,0,0],
+ [0,1,0,0],
+ [0,0,1,0],
+ [0,0,0,1]])
+ rank, D = qr_rank(A,mr_tol)
+ self.assertTrue(np.all(D > mr_tol), "Incorrect diagonal matrix output")
+ self.assertTrue(rank == 4, "Incorrect rank output")
+
+# Full rank but QR has negative entries for diag(R)
+class TestMatrixRankQR2(unittest.TestCase):
+ def runTest(self):
+ mr_tol = 1e-15
+ A = np.array([
+ [1,2,3,4],
+ [0,6,7,8],
+ [-1,-2,-2,-1],
+ [4,2,2,1]])
+ rank, D = qr_rank(A,mr_tol)
+ self.assertTrue(np.all(D > mr_tol), "Incorrect diagonal matrix output")
+ self.assertTrue(rank == 4, "Incorrect rank output")
+
+
+# Full rank but QR has negative entries for diag(R)
+class TestMatrixRankQR3(unittest.TestCase):
+ def runTest(self):
+ mr_tol = 1e-15
+ A = np.array([
+ [1,2,3,4],
+ [2,6,4,8],
+ [-1,-2,-3,-4],
+ [1,3,2,4]])
+ rank, D = qr_rank(A,mr_tol)
+ self.assertTrue(np.all(D[0:2] > mr_tol), "Incorrect diagonal matrix output (rows 1,2)")
+ self.assertTrue(np.all(D[2:4] <= mr_tol), "Incorrect diagonal matrix output (rows 3,4)")
+ self.assertTrue(rank == 2, "Incorrect rank output")
+
+
+class TestDykstraBoxInt(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([0,0])
+ lower = np.array([-0.01, -0.1])
+ upper = np.array([0.01, 0.5])
+ boxproj = lambda x: pbox(x,lower,upper)
+ P = [boxproj]
+ xproj = dykstra(P,x0)
+ self.assertTrue(np.all(xproj == x0), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxExt(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([-2,5])
+ lower = np.array([-1, -1])
+ upper = np.array([0.5, 0.9])
+ boxproj = lambda x: pbox(x,lower,upper)
+ P = [boxproj]
+ xproj = dykstra(P,x0)
+ xtrue = np.array([-1,0.9])
+ self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+class TestDykstraBallInt(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([0,0])
+ ballproj = lambda x: pball(x,x0+1,2)
+ P = [ballproj]
+ xproj = dykstra(P,x0)
+ self.assertTrue(np.all(xproj == x0), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBallExt(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([-3,5])
+ ballproj = lambda x: pball(x,np.array([-0.5,1]),1)
+ P = [ballproj]
+ xproj = dykstra(P,x0)
+ xtrue = np.array([-1.02999894, 1.8479983])
+ self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxBallInt(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([0.72,1.1])
+ lower = np.array([0.7, -2.0])
+ upper = np.array([1.0, 2])
+ boxproj = lambda x: pbox(x,lower,upper)
+ ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+ P = [boxproj,ballproj]
+ xproj = dykstra(P,x0)
+ self.assertTrue(np.all(xproj == x0), "Incorrect point returned by Dykstra")
+
+class TestDykstraBoxBallExt1(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([0,4])
+ lower = np.array([0.7, -2.0])
+ upper = np.array([1.0, 2])
+ boxproj = lambda x: pbox(x,lower,upper)
+ ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+ P = [boxproj,ballproj]
+ xproj = dykstra(P,x0)
+ xtrue = np.array([0.6940582, 1.1576116])
+ self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxBallExt2(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([0.8,-3])
+ lower = np.array([0.7, -2.0])
+ upper = np.array([1.0, 2])
+ boxproj = lambda x: pbox(x,lower,upper)
+ ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+ P = [boxproj,ballproj]
+ xproj = dykstra(P,x0)
+ xtrue = np.array([0.68976232, 0.8372417])
+ self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxBallBdry(unittest.TestCase):
+ def runTest(self):
+ x0 = np.array([0.7,0.85])
+ lower = np.array([0.7, -2.0])
+ upper = np.array([1.0, 2])
+ boxproj = lambda x: pbox(x,lower,upper)
+ ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+ P = [boxproj,ballproj]
+ xproj = dykstra(P,x0)
+ self.assertTrue(np.allclose(xproj, x0), "Incorrect point returned by Dykstra")
diff --git a/dfols/util.py b/dfols/util.py
index 154e041..02c9884 100644
--- a/dfols/util.py
+++ b/dfols/util.py
@@ -30,8 +30,6 @@
import scipy.linalg as LA
import sys
-import pdb
-
__all__ = ['sumsq', 'eval_least_squares_objective', 'model_value', 'random_orthog_directions_within_bounds',
'random_directions_within_bounds', 'apply_scaling', 'remove_scaling', 'pbox', 'pball', 'dykstra', 'qr_rank']
@@ -211,7 +209,7 @@ def remove_scaling(x_scaled, scaling_changes):
return shift + x_scaled * scale
-def dykstra(P,x0,max_iter=100,tol=1.0e-10):
+def dykstra(P,x0,max_iter=100,tol=1e-10):
x = x0.copy()
p = len(P)
y = np.zeros((p,x0.shape[0]))
@@ -226,7 +224,6 @@ def dykstra(P,x0,max_iter=100,tol=1.0e-10):
x = P[i](prev_x - y[i,:])
# Update increment
- # pdb.set_trace()
prev_y = y[i,:].copy()
y[i,:] = x - (prev_x - prev_y)
@@ -251,7 +248,7 @@ def pbox(x,l,u):
can be given by the number of nonzero diagonal elements of
R in the QR factorization of A.
'''
-def qr_rank(A,tol=1e-25):
+def qr_rank(A,tol=1e-15):
m,n = A.shape
assert m == n, "Input matrix must be square"
Q,R = LA.qr(A)
From 4bf0d53c1cbc628d0cab485306ca025c2e7e420c Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Sat, 16 Oct 2021 21:55:14 -0400
Subject: [PATCH 06/11] Add constrained example
---
examples/rosenbrock_constrained.py | 31 ++++++++++++++++++++++++++++++
1 file changed, 31 insertions(+)
create mode 100644 examples/rosenbrock_constrained.py
diff --git a/examples/rosenbrock_constrained.py b/examples/rosenbrock_constrained.py
new file mode 100644
index 0000000..eeac7bf
--- /dev/null
+++ b/examples/rosenbrock_constrained.py
@@ -0,0 +1,31 @@
+# DFO-LS example: minimize the Rosenbrock function with arbitrary convex constraints
+from __future__ import print_function
+import numpy as np
+import dfols
+
+# Define the objective function
+def rosenbrock(x):
+ return np.array([10.0 * (x[1] - x[0] ** 2), 1.0 - x[0]])
+
+# Define the starting point
+x0 = np.array([-1.2, 1])
+
+def pball(x):
+ c = np.array([0.7,1.5]) # ball centre
+ r = 0.4 # ball radius
+ return c + (r/np.max([np.linalg.norm(x-c),r]))*(x-c)
+
+def pbox(x):
+ l = np.array([-2, 1.1]) # lower bound
+ u = np.array([0.9, 3]) # upper bound
+ return np.minimum(np.maximum(x,l), u)
+
+# For optional extra output details
+import logging
+logging.basicConfig(level=logging.DEBUG, format='%(message)s')
+
+# Call DFO-LS
+soln = dfols.solve(rosenbrock, x0, projections=[pball,pbox])
+
+# Display output
+print(soln)
From 30ee963ae6f04f1f23afdcd9336910cbef6a7084 Mon Sep 17 00:00:00 2001
From: Matt Hough
Date: Sat, 16 Oct 2021 23:24:34 -0400
Subject: [PATCH 07/11] Add examples and update docs
---
dfols/params.py | 2 +-
dfols/tests/test_util.py | 12 +-
docs/advanced.rst | 9 +
docs/build/doctrees/advanced.doctree | Bin 84560 -> 87990 bytes
docs/build/doctrees/diagnostic.doctree | Bin 22026 -> 22021 bytes
docs/build/doctrees/environment.pickle | Bin 34653 -> 36240 bytes
docs/build/doctrees/history.doctree | Bin 15063 -> 15890 bytes
docs/build/doctrees/index.doctree | Bin 13924 -> 14340 bytes
docs/build/doctrees/info.doctree | Bin 34483 -> 34775 bytes
docs/build/doctrees/install.doctree | Bin 17426 -> 17421 bytes
docs/build/doctrees/userguide.doctree | Bin 86188 -> 91374 bytes
docs/build/html/.buildinfo | 2 +-
docs/build/html/_sources/advanced.rst.txt | 9 +
docs/build/html/_sources/history.rst.txt | 4 +
docs/build/html/_sources/index.rst.txt | 6 +-
docs/build/html/_sources/info.rst.txt | 5 +-
docs/build/html/_sources/userguide.rst.txt | 45 +-
.../html/_static/documentation_options.js | 2 +-
docs/build/html/_static/pygments.css | 6 +-
docs/build/html/advanced.html | 21 +-
docs/build/html/diagnostic.html | 4 +-
docs/build/html/genindex.html | 4 +-
docs/build/html/history.html | 11 +-
docs/build/html/index.html | 19 +-
docs/build/html/info.html | 13 +-
docs/build/html/install.html | 4 +-
docs/build/html/objects.inv | Bin 401 -> 399 bytes
docs/build/html/search.html | 4 +-
docs/build/html/searchindex.js | 2 +-
docs/build/html/userguide.html | 52 +-
docs/build/latex/DFOLS.aux | 199 +--
docs/build/latex/DFOLS.fdb_latexmk | 365 ++---
docs/build/latex/DFOLS.fls | 888 ++++++++++-
docs/build/latex/DFOLS.log | 1386 ++++++-----------
docs/build/latex/DFOLS.out | 43 +-
docs/build/latex/DFOLS.pdf | Bin 302355 -> 321015 bytes
docs/build/latex/DFOLS.tex | 104 +-
docs/build/latex/DFOLS.toc | 109 +-
docs/build/latex/sphinxhighlight.sty | 134 +-
docs/history.rst | 4 +
docs/index.rst | 6 +-
docs/info.rst | 5 +-
docs/userguide.rst | 45 +-
examples/rosenbrock_constrained.py | 3 +-
manual.pdf | Bin 302355 -> 321015 bytes
45 files changed, 2029 insertions(+), 1498 deletions(-)
mode change 100644 => 100755 manual.pdf
diff --git a/dfols/params.py b/dfols/params.py
index 6d36e2e..af7b93c 100644
--- a/dfols/params.py
+++ b/dfols/params.py
@@ -113,7 +113,7 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
self.params["dykstra.d_tol"] = 1e-10
self.params["dykstra.max_iters"] = 100
# Matrix rank algorithm
- self.params["matrix_rank.r_tol"] = 1e-15
+ self.params["matrix_rank.r_tol"] = 1e-18
self.params_changed = {}
for p in self.params:
diff --git a/dfols/tests/test_util.py b/dfols/tests/test_util.py
index 57ecfab..51c1b30 100644
--- a/dfols/tests/test_util.py
+++ b/dfols/tests/test_util.py
@@ -170,7 +170,7 @@ def runTest(self):
# Trivial case of full rank
class TestMatrixRankQR1(unittest.TestCase):
def runTest(self):
- mr_tol = 1e-15
+ mr_tol = 1e-18
A = np.array([
[1,0,0,0],
[0,1,0,0],
@@ -183,7 +183,7 @@ def runTest(self):
# Full rank but QR has negative entries for diag(R)
class TestMatrixRankQR2(unittest.TestCase):
def runTest(self):
- mr_tol = 1e-15
+ mr_tol = 1e-18
A = np.array([
[1,2,3,4],
[0,6,7,8],
@@ -194,15 +194,15 @@ def runTest(self):
self.assertTrue(rank == 4, "Incorrect rank output")
-# Full rank but QR has negative entries for diag(R)
+# Not full rank
class TestMatrixRankQR3(unittest.TestCase):
def runTest(self):
- mr_tol = 1e-15
+ mr_tol = 1e-18
A = np.array([
[1,2,3,4],
[2,6,4,8],
- [-1,-2,-3,-4],
- [1,3,2,4]])
+ [0,0,0,0],
+ [0,0,0,0]])
rank, D = qr_rank(A,mr_tol)
self.assertTrue(np.all(D[0:2] > mr_tol), "Incorrect diagonal matrix output (rows 1,2)")
self.assertTrue(np.all(D[2:4] <= mr_tol), "Incorrect diagonal matrix output (rows 3,4)")
diff --git a/docs/advanced.rst b/docs/advanced.rst
index 6073888..1303bee 100644
--- a/docs/advanced.rst
+++ b/docs/advanced.rst
@@ -103,6 +103,15 @@ Dynamically Growing Initial Set
* :code:`growing.gamma_dec` - Trust region decrease parameter during the growing phase. Default is :code:`tr_radius.gamma_dec`.
* :code:`growing.num_new_dirns_each_iter` - Number of new search directions to add with each iteration where we do not have a full set of search directions. Default is 0, as this approach is not recommended.
+Dykstra's Algorithm
+-------------------------------
+* :code:`dykstra.d_tol` - Tolerance on the stopping conditions of Dykstra's algorithm. Default is :math:`10^{-10}`.
+* :code:`dykstra.max_iters` - The maximum number of iterations Dykstra's algorithm is allowed to take before stopping. Default is :math:`100`.
+
+Checking Matrix Rank
+-------------------------------
+* :code:`matrix_rank.r_tol` - Tolerance on what is the smallest posisble diagonal entry value in the QR factorization before being considered zero. Default is :math:`10^{-18}`.
+
References
----------
diff --git a/docs/build/doctrees/advanced.doctree b/docs/build/doctrees/advanced.doctree
index 06e76532206209b518671891b475b7936fe12fd7..678b5b6375c382c4c7a6c1688cab38fbf7044451 100755
GIT binary patch
delta 2877
zcmc&$U2GIp6y^-w)~%qB(wXfK?FA*!;aFPdl~dd|$w?$Xu}9?(AA
zxp&XpbH4NQeVuPC--lOyWv6@>;96Nb)Ng8))s)nl5)GqukEUiK>S$|JjY+MsxO%-0
zuG}}@I$rL-FmrT?^Jn&Ku0H>3!v_c9&sMl;?+d}rIUx)?%AT*k}LWwxjI*!x7|B+sl;jOLQU&Tjk98O(Ti^PJPktGcrYBY!#=3$f0tjfgaB=j
zwWib71U`Ma>P6#vxR*Zw$9moJuFW&ZwaRodSAucr4I+
z;L*Kc{R=)?*d8e^W|Ia>=*&Wtv%KB29%?s+L&vBt69ZnRPk6r+la$Nx!pFRQZ8k-A*miCc-8r
zlj?{RV}{BMaY$kjDXwb8UQ;J|2$`J?Re{`p5-EJ9>S)csMA8dYTj|Y*m#$tcUv3d2
zmhGN1CzW@+IZ2q8Z`JHvSfK}p;bAp0SyvqllB>fDU6oG6!yF;877WcD-kcU{*Mvi%
zLgGvDP)yaCl;3&t-0>#vPJUAG5=}}9Q@&l8l)EQ(y-c%nmY4|g(riSMR}&E%`>b45
zzL-HhSQubSl-HnwJX2yW&Bi=E-hcw~_IN|rs@<{_9U@&gC>mN~l=X_~A#b`DHx}ru
ztQQbZt+GxqX9hzn*3t>l$#rtGTz?>g0#$AP{Utw%B8i^+QI3XMv_%DKL>3*cIYO1f
zsg`tuWmG+(N0JhYCB#8BEhd?iHncHzTuf#qmT+aQr>hU>VD_qLo|{3pr$J&+
z3AGD$=GqzA9d2ic+xC$*X!5`6a`d#z(UTim{_VOp4b
zXAArY+^jHn_C?45_rCkm>aB3Wf5m+XZ69rgb-+Hwb5?FrKEr3OQj_9HIsEvDSqA{K
zn-wl@gL?F2H*AA5+--8S*8hX&!~T979P!VRi->Ynus3aoDED(ATE#NuH?%3AkZ)#x
z-43kUqF@CO&c)9K`~G%NH~&R`@=%oT@Q{Z|eY*>qs{;5U
zA(+OeTq|!?8Uwbl7uHk+@du?jVDH!q!`+RPQ#~!?w2acSj$T+xOC2rbNV!RduV|^D
zV;L^loIfq5Bx5Wi>GVOS1OxU@d!eh1e0Dqjx*cELj-M&V$CPt_Dl0eAjb?NsKPg{T
z5cDhiV7T2uLjPb~Bs7pN5*is7$u#X8Bwx^ukSO=Fh!Kf!9d5?r{FF8l)`*7qhh=v&-(&%5{BbKZTs=v+0Ojs-`<
z_SzA^U_j1m6*(g(%Sx#%&*!wPE{F9@HLqxuvOJfoEM39N-$EleWwQ?L8O6kd1^+f0
z0uH>Q%;=}XSo-v$p&{_yyT0eA_?Iu?HF^NMs{fXhK>Ls7qZEP7kqF;OR5DhmJQcO4dZ`y8!)!oAquvB;Oe4e`lEa#~1_L7wz%mh)cJYLAzqO
z?9y0X7wOd)(rkC{V$keEvuv=jODUD~Qkwl)&*ro_({LsekPoxn)%@LxkscB}pBDRNAr=Mxz04zdR_F41BF#R)$LV?2!Tz)9pQb5FkTAadBt
zLe&to>v^SetDuzmc1cj8z9b@E&cA4JYZk9&acWg{fOn0vs=7m6B+S%{B#J|pQk{>s
pl$sr-b)Io0vXdrK@CuJuK4Ws?y%cWRnaUF~9*i=k`WgATsrtEzB_;YUMfn9u`4#%f`6;RTDQWqea~Y-m
OnM+I3Hg^ZqvjPBPxfHwr
delta 63
zcmZo&!`QWkk)?rU>i>-_9*l}s`nh=}`pNpC#i>Qb`Z<|-Df#*?`N^fZsd**E`stY^
Tn==@t{h32U+%~rb)UyHr&z%-d
diff --git a/docs/build/doctrees/environment.pickle b/docs/build/doctrees/environment.pickle
index 0740510ac7fa914ccc24a531c048c302b294b188..ff76c8c3b763b8e1ecaa39f50624cb97a1f957af 100755
GIT binary patch
literal 36240
zcmcg#35+DyS@vFY-?L}Dj@zEBy*_$oXEt_Bz~ITwuI=n)?W|2qY?|usnyISo?&_^$
zXUAb72*ufz3T%>&5(I=p;U*Fi0&yP+fkZ&WsRVH*1R*W~LI_aA@%{h*>gcZS*{;gA
zt@TvbJO1~-|9kI$|NGyo_ujkjL+|^JyZFCg&TSfIg+n
zCxSlMdg$9)H(JMoJzBHbx)6+2?TY8O|GQ?k$3U!_uy|f`g$<5V&c2hF@-Kn;@hIBQ>qD
z>Ei_=m~?I5spuuw-`FsbYrJCDYI?Z-8b*-Djk^pOFqpBnOaMQ6ZWJe{}&Ri^>TW45R*3Pn58{cyjRsN^ta=uDN>4crY0C?5eh1+Ngu{$gyhS*gFcx
zUIg&^^-9gJ>ZK;I)EyuuL`
zf|i?AD;P7pdadLcdL5xN5?XS+Y8a&AO^6c3t*DYcrYzIAwmBp|t{Q=EjU49jgQKLmNAMGX>K4We!a?n0cD^m*?xlB(E
z8Y2}_EZ6MH734o>x}MqCq!>RKx~A29um*)>6j0GU&|a+z?{67zLA^KV$Nk1au$S8d
zy**$&D8a|ok6<6S|5k8UOClaJ9+n@6F7PuLN4Z)}FF`+%aD!oacUzjoAC@|lLWRJ>
zl7@2dJ&tXAHx}O$#@%=iw9XAx(wmKYrKH2gQ6V+06p7wIjZsJRE#KQHK2da`=_JXN
zO43BbN4G}Bg~lVsqmp2vs&8n1%_})J#Hi|USgTe`8-9az5*n3*zGUbw1b0(%%}|26
zUG-5>;ysGE$8*d&khHyq?dc^)-+(S_kg{uP6}@#!73v!HR!KCFs{)M^#&BkeLZkp#(q>Kz;AY`Y*@KPD-oiwq8PhQdWW@$@v
z8W3T~O)%^F&8DNHF-I{VL`@_(h>ElBEmiD#-A1oNHKEAEg2te&OM26+*_)w2IS8F6
z5j4k=LO%qfUTLGIp=>4HsV0P|Ios%^b10BVh>%Fvyg&S}n{y%i;P}$SLgQVQ&d@pjU*|AM#_lxXH0r_%HIXWS?0$2q+k*`d)mp)eHI1&v8>v?N9pHW;Xx
zuAw^^CX(T?4Ph)u_%xHpFxv3U4HFrHVZU*uVQ+;5rF^%-gSPRHl52!>`DV3+5Ei+x
zULixP^-wZs6RGqtNO9-r(QqP3g^Q+5y+(~Q9AR8Go|h6uY@tQcmM9Qucy>iK;XOiB
zIHzD(d#W2oNnMBn&E-3j<)g!w&`D{yRZ49?*R!sv@Y0r9r9N_(py8K#*Sgh_^ss>HyGg&HaQu%@DzSlE&yb22;;`1&LuS
zaDHROIH}Yk57FH2+V$mJA|lcFNJSfOH6D}PGn@$3a=el!ZBfODeI{0&(zK(T7&_2K
zhhtb+<2Q%mSh;sAT3yOH4q_4F-$Ul^Grq5Fd<`BSg|DHXdrCrFkArU6~Mat6F;Z&
zitUqxl(Q*`HEX7)u)dY7B=x$I4BBXioMS&T!A(aZsHiaVphWVFg2JOok2^S}9o#
zVg(JkAvSEwXe7&}UJ~Y2;0UKlqAEh9!VJYqLNOaPh^=stQV{|*HCWP?p@5}su)&)?
z#>~i$*~fN_;$9je##)#VFbQDKk#c9z2z)Xq(U=oXsdg_GRAwbSHrr!!qP-##+0n*f
zhR_5=gjIqHa7aToSeq-Mlbu%;EdaJ+Or)wb1kkt}%30{~cn6b+=6FsxZgGaUmZ;An
zP_+k)-aM_8b{)Li@Vut`_NAptwPCr7m748WVKq4VA|$hLc>4L(r;+vCxuiP1k4pSexw&vFx%KlUHPfCKRVKp=hk4SpYZr1c`I`soK
zooJRJR_A_6(5V}lA4M}U+U;&Y=M!A>l8vKC*TN*k2ue}Ls__m^a!`^alyW$PbM(04
zOpJ`DIZO$lZ4~1OHQ|PtT>+4!*j0LTSC}c4D;;L<6!l|?ATuhohG*koN5_bZQQ+(@
zrj&^B%sS@67-KpkwruALnT%)mb;f9o3WiRoGmbbMm!cwK?t8MH7zcO3AF2;LBp`3dMzbVU7$C
z=0d@?vJCZ`ZY}8H`#@@Pj3|A=82LWj*}@W
z;oK(JwHzm#ij#?j7ZFI5bEOR7!?p!<4X~d1drnKNk&e`C5V|E
zj4P+!Gk`1>*xhkD@?sI-Xoxi7LfedaNr*p6CN6}RV_FAeT3m71N9+&@rk8xzSt^^2
zCByVRvleGb$SR^Md1j+LVMLs(5Eq2OpZstobUAo|VNZKDy!aTF!}c|UR@LqwMH7W{
zhXdIa_jSfoFtadJu<6EWyBMc&Bq`0%@Kbi(^i_-mNgpgG$YOk}2=-aAM@r?GP!Tf!
z4)!%L8Mqdjkjc1C+sOCJ!K7F-aABi!(#R0Zw}p^291VtGFodxb&TcFM#9sR;wP+Bz
z@rR#N4*h6NB3$nOSD_L~sV7zj4A{!HvrT?M_)27U
zjTIPEgR2K#Q7j@1n+yJ2RqWFk<0O5K)!%51dF#*A+crpQl8
z)BDnO^3szO$RKq6a(n2Qp~02R79*^k;2?KgNCkd73?rD}k#eY_p7TksTj(K&x11gF
z#-8AIe&k>>bZ7w^csPtN1rtoN6ng)Bbc|Y!7CHomWV3&I_4D8Gh2e#k5x0K{dgQMU
zCHo2m{>>+j|NHkRIvIui?cvvc{G{q_^8SH0JK=&H2k5`HIbl7y`>wnA?T24q_ys?J
ztDWp~%yP+oG>!!0%BKyJotFBs~G1zL);Qr)|NrB)MjXOaTgmjlA&ZyL2H__B0TD7TguY`
z)i_cyu*eI`3Z|CzolO2R^d`yW;98`y3eQ5r(J~QPCGTr7A?t&bx>urLNY*b})^{P9
zU)F-Lgjq+LZjXt9oK|XhiU3aWd!Oi6{2GBrFn%Gsi>E9vLw>^^6e%z9sjB&9%S0j8
z>-Zy;WZsO1)j|mB^TttsBz8Ix7zuPgN4H{wsUBG4(~wtPv!Av^7spR
z{1rbG!e77(G<8E1biEzgJwp
z>SK6w3+K%S#dYa+UpYCc`M+`n?tDuk9ef*!Ourx-9tzJWbE!?Qq+*5Ny?fgF;#b!4
zt=G8relS+EG3%g#Hddspw?way0+3!z9Zz}*HEN4nAT=Io0VGTweJo}ZD^zL}WX&wv
zb-e6j*NOO3Bc0r#u2>C^lE<|K&t$zvvg5(Dvb$i(C|y4rx@)0)wM|EB8a*-#Goz6K
z(kmHj&R#k`37GZ0RC(5y#deJP(3dnM%ESdW$`F|H9kV6F*_cclS>GEC>2^VaZz=M4
z<{}lPE&>yZ_cKN7C?4I0l;hTYJz}f}dqq6jD}seQead$*410xVAqjXN+cj1yR=B(X
zlG%o)K1ChOEL4P_B00y2AEol1M3GYX+02D25^E8OXL6FjB5~9^gqhP!>5ii^aoO=*
zuW(8HG06rWv_}aSZsz**i^t@GN1POuXD9_mGjqu)FzR_p)R1wYD9oRSDO-aPil*4K
z;fPL|piKF)h{UX3tGrk5$-nODuv6RXD&a5%b(^Y(TR(md9oL@8FnI3u_3b4vSN1b
z$y=j6dwFKU^Wt7cOwZBscwQr>%N#nX53HYwjKm|%Uo`Q^Z6s0!;acDYy*4A!iVc;6
znafJGta%YWUgZO&YKcvFuC9)IZ1Wf$1<0k
z67EV~5^+CS5`6BoQGQVQjM(BR`0S)`1*3ipJx%2j3`{G2homsyTnR>hSwL7Uki5Pg-KP8Xfmd9u0@jLwJ-nOW2T}{2T
zKk$D$#CZE6rAY-n{efR*if^Dl@GTc?-l?wjzbMdO7HO;Z4+@mv*aCKO&}K2sDbn|g
zUeT51j3PGq#LF3~56^-R;ppDX94QiqBN9h*l1RG~k1(a49VYku=ox#f;MuSw^{$J3
zv*~yd&eNG|)@L6I<`(c#Skbd#Kk7w%yI91&p%gZ2nTt3Xvdl=_y2ROnfXcK`0>0X29D!`>QHYlT(C`dyKu2gCmCQ%Lg6%hsoiCLrLP=<5lqenb}Ig`yY8p3>@AQ
z5xnbm?HZ|G4P7>z7p!mNSduy&5OEq7HrX&8i#^%13)J@x4QR#|oQR-yr9wpwojIH?
zxLl#4y-|eASyd?2pO0o{E(Q3Jyd(w=_}HcJ=!Om`lMw>asmzt=w?B@<<*yEq=rR|0
zDs#Ch`9GPLM5p{?UrXOJ{4VQsGPlo8%Zl0gubd7CYm-GGd&=H8zfI>x%{I@?pTN%|
zP3%SoPPg9|*}aL;zS2HDH?L$aQd(DU&&`%dxR-G6{LIAWn;k_4#lFj(YxgpnT!vdopvPdVVV+@sXS)(jJ-*F{L~1
zp_w~pZ{jmx478Znbeqz2OCaL;Wahf{3)K7qzB(qZc;3=JWr^YW^~~j`;Q5ujBz_qg
zMG6DLRFXSaLaHI|;w@enWIE&8&$f?#2RL@j`)A(?^1$xIup5SoXM#>@~d9-H7ud
z_W9uhv|26U4B*4&!ii$s=^o9@Q3}n6@{$-hnnyd)tYM*{`wqaJlL~t4nJd~adiSKz
z1E;B7fiY^&W#%M>+WEXB29Fxfj7%oXh)HGI*mIA!OX
zj%_XR7lA3w>1F05h1%7;BnFNeYz%0iaKCrY5w3SKLF;=mSFm5S$P`I%xZ@c2k<47A
zkor(w5(7u-ZhW^vTcBx(lTQeGnIk)<6uf>db2anQv!A)i>4pic^RGE06rGb1U;
zejzW3fkQTXo{mbQMZZgM0pnQ74g
z@)vnY3>=hu*S4{@+Qh+Ywe7;wj=e<_WHlixtm_@7{cTeVV8b|!Gx_b>oip(Z3)Pt4
z;->J60gI<+Mf$`)1>fGx%%tGkm6yc8;XAOleT7y&zxdwoc4MbhyKO(O(ma?kiGGpj
zwjLADkeZpxP07EWmqe%h<297qYjfeE+jmpqvCmG+irMkZb=`xI?3zpw?ay6>%w#r>Jg<)_5`TX{*u9p;mQ*q#1>YBpN=>b?!7OT3`+
zubHdVFDi2j@h+2&T_rIn|D3t_6e$0Ym&Cw9c_7?O(z%ty-C-|DznKKS{iAorxKHV;
z;K7A(=Sb(i5x48ckr>n2%sizqoytpM;F#XK#%H@Sf;df4?Jnh{|S
zF+~%82{Z7kcAoDv3wSIuGb#91@{$-hd`WXMX#xbv2m$4#%$4XDl%4I&7c-ZelK%^N
zNp#9ThiyoPod{5jlN=;sE2J@+&_rEZHWl|BTGi4|M?g--yAd
zU0;^R|H$M22d7OlilPb817W@-Ah
z>Y&_%*XyLL@N~dFL=YW~$kiUk+BhF=j^q7tK&|>#vh^yC6qVB>=~v<;S$d5RpQX3t
z;gs|~93DXXDz18)P^V3Xhq)iOeiik#ex3e&g8uv_{rN5W^J)CSoq4!RMcvY}8;3J$
z{BaDD+P(@2L5vxQA6X)Cq
zV>o~^t~0-erC6T@3kt7D;JBQQ@7ZVtC;Gd0aEZ^e<=gv$z4Rqo>+HSK^*>JJ%RM_-
z-p{&bRx5NL56)PXtqvVN>vXhm^9AN6J_DnVrqf^V3GkCr3&d|1ipw(*H7Z>59Zp
zMI`<#CyDfKE&(Ri$+s7UY9(E^b4q{CPRshgsA56RidjQr+`#zw!#7U9rw&VOQ`h5A
zUjq42<}Pk&C0cQD)9zE>6nb}AIb|Wu|
zfdh6V?qTR(%a0e{sWDt%)u9S^NGD+ZXyz*R3)bO8?^-S{KVGynhU%-*%J5^EnMgtP
zgLz2|9I7`Zpdw}XSYfs4P}l0jYNvz(*6(DlWWQh?X$PzLSW#;_)Umpt_36x9q|o|Q
zUJ?UGE74(~-IiC>z_oIscgN~XCs_SU<|_7!RkCBP6qi@jptW)$-Ld{DGZQJO{yr~>
zfkTxTidLv&EuZWK){ZF!uKDqHhxi>cy@RBc0l13GCwBmBCNmo;z$WvOcw+%0+aa{z
zSC%_^R~ulRDFv{%Wv*tw{i_`?SQ()mzp|X}VQwshKro8p~r<&;?btta_CZ`vqOB51|D>V+u1Sx-1+P@@3KK4xxvrw
z=9@Flj{CtBkx5)Ac>}j)Q4sFY6{j*?(h`q}G+B@f&&3qc=a8Bwsm7s(85Vg)eoQ4PE1H
zp<{BB8DI7UW`bQ6;wcz+Y#Y~TZr95;E(|U^^knRFDAS#$FX;^hX
z6<--e*VoB0EG++3Lh>#@7{#|D>V`V1#ISlLh{7ZzmGGWA61{{Xp9*A%M^2qNe~H~&
z+){8?)kfXn_p>j*zvbcrK$V&
zSGy#w-$T#^>$CJX7-gwlbS^d_pD7%-flHP>op-hHrEg&fjt;}S_4}Ol5BTpN;_of%
zkMI{uVt>L9tUb_uQ3)}qd%+Sih}MHBjCv1p)3MBWP=)ZI+TcMo!-IMc9^{k7?O1-W
z+aggMw?2mw(NBi)L4hw9+eEV8iw&JdGSMVHLMp{)ZiFS@xge&pq_#&G_?vEX|d1Pv!^&E4+a7ZTfCQ-f->M
z5hOd_kvEt)&p=?`d@t$`B&ZLm>*B3P5!j35_&VRAj^hPz6IVBWiwTNtOm?yy%srhM8nk;!*Pl1JCa?Rdq;P%HVk#l7cH<-g!3HYoY_oZ^#xct#XGYdwkQ_~j=iYd%TGr0J=kcBIe7|$P^PTUUd(MY;fAZ3cJNTbJ>C`mcyf&+v>$+{3Rn>Is
z{`d(`uar;4f?lrQ_k8`i`i$SB)N1t${$SZExwfj7b-V6AJXbZ{xzgN|j%qt|72Pab
zb4ymqBVx|nn(oe(SFMUOhdj<)kx0`cWq0a+uVa^(%9~P9LAH`*uIg)be?(Rwc<((a
zB0`P_taa74^)d>3uI}$~RM+DgsMlQ8bacyf>Q_*)P!QML*{b4dhH|a$?+>JdK-5F?
ziZx{ogjB!3qL^zQULcGS$MWowT5!D8RUNs8N>-(!mRzRo4_&m&s;!on^paZ#SC{Jk
z?pfP$>lYkFZ=2eT%Y`KW%*`<}jifh>_*!HNaZo8*o;yqNb$`rM*7Y^cBSS_}aejA#07q5ql2){?$$%jURdH4E
zraJxvk52Xp6E7_(83Q?C^uPOdbyurU!7u*|JXyUQQoOtVN02xaY-
z8mgro@q4R^=_!@ETUY!6JVW}_+|WKjVGsABx_0_XT|1&I@1#T%pr@Ow7PTPJ--5(k
z91?3sw7dP>MN`e8W2!lHgOaV+NRBycH7A?eEM$6KcS>pn&04iQhdCHkYhJNXa-71d
zUST%;QI>o|kO0GoB^oikD*9$&)w0o(SZ(w&Ng9b}hS1rcWgQMG<#J)wGg&90QQ6SU
zn(9Ds*F~;zN>H`R9x6(_hY@$Wwq6C2mTOwBTCmks=rWU(T}>&e^&7HKN3%8xLIXK6
zP&=m0YdVDPPioLy71CazkRV3Mr>xgRl*Z!L8CK07b=7K(I+a7pNCO6f>{t+9C_=rH
zCKll%S5>uE*idW}A`H0kCmgR18LGI
zT)%LE_!$>ru>tBi#2{iADuFdZ!Lus?EB&VK*7VZV&H4>aG%DW<`kF~P`HINfD>`L>
z*p%#}=$CjZl+_v}yS$SHqP<6ZuaL|bjf4dxz#vmlDwQDj1c$3nLry_24|+?G1HB@s
z{-9fN3K)d+)y*Jih=NLJK{R&Btb-m}bO8zvCP+Zdf!N%qJv6;gn3J)vg2FAj?)N;s
zgwc%A+#8RVqZlg)e<0|F+>dt|3KVER4`|er$O!0xfH6Ensk(L56ZA${3x?WXyIVlP
zp~1*{=Fxs>Pwl2~(m?(orx4?v_Nw+Og-u3bg;+D^IGLWuKkaUjeu&e{ksTV16AFfc
zx}Y(sW?jSxg$)KOx}&K!hKYE1WJ74vB7BU=V;D7EeN{&Wzt1zTn$|`@P?YbKc+l2<
zBy#n0E>ADl5yB!Dte42ps^1q6T0<%|2$I~{YB-#Tq=Ji@MZHFjGaR8^)}9n43b6$i
zNkbxENW-;CvI+MHqJnb@hP5ZVVVKl`D9~I!7B3$izJN|j!>v?md%2!fMTQqP^fL94
z9U_;g4eLF0MjFd>lSbGEDA1w{!m7J2YBJ6=V}mfskWm#euwtP`iXN<~FlK`DY}h&t
zd!(c>Iu63r)1N_N7z>x4TAtIT?$HH9k_&+Wp!CBKJ5aLbYtS;EJ{=
zV}yMstU9S_2RJcwpoI?0s(zS5R$+IcCNfZar}h-bYVXqC-B5}PTddT`o+60Sz{*C|7$>dyAkfl3_H1v|!0qK!hO0FvtCF
zIhk1skInYToM^0wh-_-(FhghpB7{|d3b08-nyk&G(8m5>MSHP
zr(phxsdUkqley*=X5TXVz#Q8_?yTWF0?A{HefwqgUHQfHEBW*1&lMiKeEyte<4MYj
zhi?#yk2Ngj(SQ}HvuNzw^}ZRj_V?XaR|OPd=Cba<{MzIJa@YIbONp)oi9$kxn?MbW
zt*MwbUUSKyB9QCX{XQAQW8-)bPqsAA>8=~k?A)<~#>n9y&Z8mbMi@RBgkw0g^?F$Nz9_et_7_}LsbVH~Zt=pICl}8k>$!7}R`glRUW;WL7IQ(>T%(b}@2?TU
z5Dc4S+DUDRTi>8aCG;qVjEfNXTF|gqM@cL0u!x|ZFP^F}^`UDkUe%QLpPYYE5=;%t
zfmwJ+q<5`aMO{~^AF$~}vkYN%?iC4|btCg5X~ssoT_$us!8IkaaTMuVn1mQXDN0+^
zPH>X_B1uds2SYeV4@u61k#UB@qyQR5F^-TEZkX9607;5nrbl;$nNqsiWcH3yKNb;W
zMg`XJL=vFyVesT}7|J2jaA8{YH*
z*&|0diK%WtN|3+5+cPog61{bb<7g}n$AWMIK;<%NIWK7#BZNJA4T>IO-pYBe#YCufFS&;%!FNp~+rwWWN>0VA~n=JoJ4!2<&C
zQt2ublA{II=RQCQ*xB>snPoXIufxYN5XP`BoeE!}6geT(<^|&?Ql0LNqRZ+!Ch%+Q
zLffO7dANDPWCx>Ul1ebQ@psNg$ws4O!omv)#LBr^gz(`@AQ$GaaCPN{=P^t@9??N+
zp3R1FF-o>tC9gwc{DCO1n6(F;?nX%BP?V&qo6yyEwNQZ;h60Er6BC>RwTMn#cIw*q
zLLxC9h07-%ITn>5%w&H^I`tj{WRbuQN9o9mMS!Cr(u50bGvXyC{xF$vAzX=Q?T=`2
z#bF<@L&P7O^BjAwsGD<|?zws;$`X@RNLS>U2=fFHQLbE65C(tp!xhlw;01;~?OJf<
zV^|K_*Em{Lqkj}l6wVzEWLMlH+M_VDFjTPV#%Ws^r%@y+&A{-JcHP+17zyG&SWJ+G
z@hv0RXT=^VnPXT+i21j_$HZh{Juo5Tam}`o=N0`CVa>pWjm}9UgFn>}Lej7mH16EY
zy|SC!dmmhR5!sGoH71x()pMSmURWb9=JuZiqNOQp0?~w@*TpTB4*$>QK@hl0KU?|
z|2`{pKdjIC@-F|NX4L>HN}#xO5yvO!fTsY@Id*y66MqCN{e8}Eu(a8e)+AQzgOZ$Q)*x(FA$m)i_54OB=wV)zPrGPPv
zyp9-guKU9^I3{Tt!5Yoa1v7w+jg47R8?qng(0S&N(wW0z4&?})uGF-K?bT9@#vK6-
z?w7N=>aH>kvKJ^nrEbVSs+$c{GwnD0QF2<+T)%LQT=padvJM@u*cdvfDex<^`3Td;
z-_N}mQh^f>4FR(~(iCOXb1LqF3*6;!owIjd>+!d9D*GdWcMI6S-JyNYA7+w;z!m8E
z{pbO->cKD=Z`-wF$C<@f^Xub3zyFx=Nci?YKfY)@8ovGNn}22eTJSc}`}6yaUk~4g
z>@ON`58uW<(_<`#Z--|->-li3lQoVxEm@1&A%94EuwkIn!k%y?qGiC8>7h>Za{VLq
zOLZngJ$W2X3k4o7ZalL!*xffBzh5_bnXOzmuA(@@#2?D#k6<>bDzGdn8_K3r|2W83
zVPg>unY4nW{%%<~=jmm*B7rclG#)s!CXGE<0OS<33Mp1|2W(|S`X)5f?=NXsn}sa|
z|B9giB3|(JcdwB*Qi0p>HvA+_P+PLq7$rFR
z?rY=7^9L&y?06cjB1Ok~MyLYnsjTr0>Km->4QeZ_=?%h)wY))+VbyFfM_BWT4<$nT
z4HZ4?;}9OoP-ix27*>-5+mY{nmx#tX~=sf#B92Pl-KvSuqatxaaZ%s@Db
zw@b##u2~nywTST&RhkuAZ%t)T9nr}Bs5M}L5
zU6v&AtB~FQ#*=T&!PB~yQ_OPiBCkO-V$~{vKtFL(?!=VjSkJ)qG8Gm_W)%
z%|!yKwX7t%j?{t6w&%FHOY9}gomWg4`?RyIu^Y8bC2;yg>gsih)4^$%_{$3~W*+|X
zfWHJxA4|+9`T>CeQGU0~)AiVa=Iz
zXS)%ueKLX8-MhOqfAYx`iE+5)829?rTqKY>oRvh^kva?`gcKmd1ovg7Em&Ko6nLFV
zUCmB=*B#R^MkJqk+MHL}60#Gi8A(9)NLCVEhiu|J?MA_}H4I2BQPJMshN*Ikcma=|
zx=NkGGdayWWWW=QM#)n`wpBKBXKZVOCdhI^mRrF(
zMFT|JrWSy`l)Ac|0+yRDNBrj31iu)tczTwnb(I8sA5G0n0=^fslIS{ogO`F&E}A-1
zJVJo-O6p2<3d#_aoc!B%0-~reCA^Y3P;`ZpI8GA};(kg#wYeV!u@a0m`RSSE^G`rlw;n
zivd>&Tt1b$`~)s9XC={$O9lfjvH^w7kby+R#lK{=T%--iyM@}Vb+;WPGuwcPy+!dR
zXS@;`e}9`ApVm!`s<)l-_cy7_N;IdhGl|=1{2dH-aLvyi%gsEPJFOOTIIH8W_cac7
z?URX~{FBsG?@-MCfoZw9H9wy}mVYq6suuIaU~;eRuWx@>UA-wNEtqmTddjCwVJWo
zDxm=D`P7x{6s$vyVC5gkD>a)sRtvN~n3{_OTF+)B(RH+99R}KMenAdg3&+}btmbqA
zt3OIz#ZIw`cdUi{{DK^`7LFx5)*qy1A_3L!XC={fsA5CW0(GqUScLov>7hlUdeUv^F{<7)z`A}Z_XH5-Gm%j*|E{~@Mw6!83
zd2Hu&Ig`4qMEhG}5?eT3{+IaZ!+tb6{8z1I@18wjXHonE=
z52atgceCPKX?lu%7fr9xm(}zZ|CXBG!_TQ-z~QFD@@viD_>vjpY1G$n=#NK#Ht5eL
z{dpFDa7rvr50mGT?ZU^73ZHC1BLeil-4-N*K^D
z4!XdnwadUzeEVdX6`bti58?yPsLuQvlwv#&79>C=JmMr)K2)miKak?2-8uefo6oED
zchf-~_2s*S&inw8FSd0@W?i!zUo6oHboh{1>^yJU^^lAf+`P=(r1?U#!+zpVOydL8
zMvkOZqa6@8QQSfGmwb5Sj<`xsk!;f%-g5Aqt7`TH1Ib<
z4gB@g^=aL*s9xJ?;IE}FE76X=$|SbXz_;3+I{B2QnmO0X;kb#GpMkbYAUf@jQdhS_
z(fU)<8{jJMS~%oJ&Ew$YJPs{MblV@KEM-9Gh}~*jCc+OEbEn0j1C8qt
zE#d_#V>rt+-4>=!QJI_$KTXWbgAkIS45uzW0m?vD5?u%7jVo4Voi<=zzyOZXRIzg_
zn>Yf#>d;4Vr||haZA)(&^E%BW;5(kW@}0tW@3bSh$5#ZB`@AE#rw``Tko!bOJC>TK
z1g80{B)X32T}yoU86$|}hNKRRz8r4scH5>FC_a_Cx}Bnk4Yq<`!VLUM?HGMOoa~a9
zQ!|r*?;Tl5bRE9T(nE_p*{`mVXA@FfuQNagu>V
z#Kl*#4!F|1%;8HI
zV04@aNE7Eqo!S^B`s%k*SFuBJ`Xks3ALB$yCYlU>Gj+)c;eI14iKw6awZP|ATb29e
zQ8eQ7IgihNsp+*wJlZ7^$owpIp!7}d~CR7zVI3gQ*!vK(j9^iLOJ_cNrg{
zp;b1*5rzr~9!br9ryZmpA8?EE6FuT>sS8br{h_QRx-Rz7Cvg`{tkK1XV81z1fKW~_%EmCC;`%!vXY3p+7|@Qx7uW!cnrsxY~den5ifB1
z=hRi|6sO5)+*h`pclaMu7oPy-+gV9;9h6BrzLPG<$t~iH=oWLTRtW?wJ%hK#LQH5c
zc_}aH0+(O@zr|&IF}I@Qx;r&@>Y75ojh13ET3~W-YMFG}Rj~Yt
z>j#OGyec_LAaYOY(i4c>os~q_d&-@QxKsun`xE#et?cm8g6DB}P|NcLw@E9oT25Wn
zPO*CZG>#XgIWI0j1+RHJaIgwq+nzf3SZY=hz@5rUqU*p-%O3PFAJ4QIJ09&6oifHg
zg2vv}Q4Rs3mAcxUf;c@b`x&3rv>7`aO^w)2biqn$eiBIPSxIyq$)h2X(VFHd%wxqt
z-g8Sh6utEtXPS%x*OyZZp;KIshPXzHpfnj5mB`b{?(G*vZ>DA|0q;v$Npv0FVZp=fAKhYkG<|jf^`E5{NT*OA3!x4cOqIlv={75-Urfzh0`OO|lIS|%@s+js6uHHc
z+jmk|qEnEzUUGXib-4-o|9w^x&GOISKZ#s0V)CAWM8w6*S*>N^g7N$@?DwQwgmFcP
z*q%vuDA1K1%OiZ?&2ldv*}go;=d~~I;=|ezNIqoz&_rQ*(DO%$Ozb|e=Ws6t1>qV*
zu|twhg!jkg=2BeT#aC-kAk)R3&~5f>I{r9s?ZhQVxHLsvbx*9)E>B$H!BlSYcF*XI
z_I*Zgv;j1Fqn)798*K><-#mYWHh!Wc-U-V$d`txUJR3p>#oaS}^9Pvmca~708Gneb
z7B4uP)uM&_o|ecH)Vm_yg6ln4&EduvN5_T0=GtNr7xZFFYw&(T#p#`R-^2G&;Q<$#
z;64tV6N8B%ZGDv&5#e!6cP^>%*c<*pA$D7cKT=TO`8XNfA3_%niBnx!{wD~+9iBgc
zb13DVLR5)9`AQH4N%~8{%|ax45sG{=kRk3rdHVb%{&Jby1}-jHt=jy4{-1WfSa)!N
zlT1!`-!x?*2-qpe8x5%y-gG>_2iLWcT#b8%uzU9mLX9JMY80h$AA&9zN9nIWz+$^-
zU#vkuqrqL|Y@qNxJwXWO*LXMH%vtmNcNTwd7<2fGOGM`R2lA2^pu|KWFMo>3A>80h
zVZ=KK&|ymWp)ZK>L;IldLzTl1^&tFEeQ|w?=kGE|8Z*Xm6pFsmha*sM<$D7}yZOGQ
zSfUXea$1Ns_XH(h=Yp?uxi`4j>s;z}F7yVMdEEkpr2AaE{UJLu&RG)A$3$!3=QsVDfVty}|>FXB2L-28x@U?Bu7cPacA>c;aZ
z!1yq;h&SR@GK{yI^!5oE7wf>XYa-$WEDg$L;Nfy@aW%I;Bq8N%?)>r7xaFAcoxu9e
z(@WwWCfxR5*l6!TKflRuFVr0hpop7LP34N%!*6^RCHW&uxbhr#c9n^oBfd)TCvoos
lj_t>lcU(5E?9$Q9`~m9M1*wgph1FHMKa8%VM66ew{Xd_Ei&Fpq
diff --git a/docs/build/doctrees/history.doctree b/docs/build/doctrees/history.doctree
index e0c1d81cfcd442ce22cad1306ae563bcdeda9d5b..c8a9010aa79add5e4ab2cec1f902bb2bf2faf423 100755
GIT binary patch
delta 621
zcmca!I;n=Wfpsdo-9}bRMj0dhjQreG{oKTo5`CAV{DP$X3jO5#lvMqcwEWHSjK_2t
zvnH?B&z#&WCOBDDzhLqT{c^^l%?buWjFXiGWW=PyQj3Z+^Yau8^^Ell6f_LY6#SD*
zCf8|8sv_|fj0}tnHK(*rX`E6!MI(bXYdSLnLxxn=ROZS30tNz;Ks?3_<_zf!tqd7&
zqshIjvh`6OiFqkGsS0VCd6^}tIh6{ziFuU@iA717B}Iuvl?uuEd1a{;AgUNBkeOFf
ztdN3Xu@2x0%Si(NVErTD=n-=
z_#bSeRMsly$pu_lxHX9fgFfo3cywF&TxwIs0^Dg5yX0{7JBeKp;b~6uB?FTV?
zK*R|UaTG|D_HgDU=BAcZ7Niz~q9*IY!E`p+Hi8a)gx`WANmDD-QsjBgzy2
delta 256
zcmbPKbG?+cfo19$tBtIdjEd&^xp^h}$@-zisYS*5IhlDW`T8#T$)&lec_qdA>6s;)
zqZp6rGGwEoI;
iig|K}&0!`v7NGSiEK*t0AcoxJ$+l*Ul9Nx_dH?{%uvSz6
diff --git a/docs/build/doctrees/index.doctree b/docs/build/doctrees/index.doctree
index c4b2eb245a662a7ad6f86ca75772ff06ff04b509..a87d31bf8267aeadba8b94582694ef8ff316add2 100755
GIT binary patch
delta 1699
zcmb_cZD?C%6wZ5lvnJh2leTHIBrdPrntm*6lDaC)AuXwjo3_;XQM6+7aqHW=BsU~C
zr75_s^DE9^%iil?IuZA);AgqYu;5UC38MI?AOjVtDE<)!RS;qGoHwcM+V0=}Jonsl
z-sd^bIrruMntCy7J#BsCpL>?91ClQsVp>`WrDM7t8q8_gxHcV1Xh|iMoY2hYq&3of
zS5Me@$(S>qH#(WMY>o6lf@|_>xFK>k9Zu+4T*=YiaBpPOscj|ZVhR3^yVv~CdQP&A
zT8PEMkJ*d65?W?TndV>EN;@Ub@MCprRtdSjSKy&Majay&K%6)!+-L4%wHTQh`hftzmNNBZG5cpEb*AXH+o6v
zFyssD0Q-ntIZj|t9(>2`NI|->Lh*)G9%(k@hbmfryZRtRB7Lp-$x1h#aE;1il>1j{(kEK9QdlxrM*Tv%VtKF1q_n{~eDLauv77
z=l2!m)vAtIUfBOAlwRJL;5=E8{{O3vaDT^JRi3%h(Lm%=8$=u5)3r!O%ul;)#5uo>
z`CWHb(e42zi7H2NsQ-y{E9gaSei
z;R3>WgwGL{0SaY1-blTeRr2`bR+x53C&W0
zxlu(J?`D_TQpI&2S`x0)7zx)47zx*3VI*Atg%P@@YCmAbb_ye5djTV~t$t$e>wAFw
E2a~H4tpET3
delta 1352
zcmb7EUrbY180TD01H~#7L1;xhNEuk5m)7Wm^@W5mBZei@s>?L4r8hghE4M&<(UK?>
z|5a_4tzUgI>Vqa*qDI|~Vc|iRh4^MB@47v(B|h$9ahhz|&UbDhhZp-XN7P7Kn@Fq4gb`1v1E~>nR5LPZ_47o=iYb4P
zhEQNv=|U-2I#cT6j)Ef?5DC!fvftGgWCu(`W1+rKa81j*dEUrJAVovg14mI&3je@h
zb=DC}nMNFzoU_%u>O_gxtxa1?&PhTzJlt^1!=mhXQX6Nd?MJLsZkN@%^;M{K<+^7N
zbB&o}W7-6Z8gZ6Pof46pmu4x0Wkhu{I+2LRk{TS?H1J`{%%l@>O&097oq?4x|5MBN
zpIWQC<8Q<+Uc>J-A?tGOONCH4+v2ry^@j-ENGd9TRQceq%~uGuE;Vc=YOBbZ+{YjA
zCqslc>J-m0u`8vAt=Tr7unJ3!CA!PjX>kdn9ux;z`FhJydl*zU9ixMGPrFke#)kU^
z2DMEL?lw6^)8nR$`-n)tGZ^+Ze^D0rp?QQR$^wc$4hKE`qVlCD=gz*zV&1LdF8jRw
z3Mq=KUZ&SqiMOF|+knWL+pf6pt|Jzqt~E=G*o&XdLTiMOMO+oAEK8r+8a{Z=J{IZr
zwxavzbpngW9jdXBzwg}yhFdER~-D-_r6aRC9J@7iw6RoQ*;2!q^bj!IuDXC
zw03>vj@as}FER`ZUDwMT%iLK@D7Q&Ixaa?tgsn(-HF0O(V;Mc8E0bRLWI?f&D^V4}bGcV`m<*dA%lb18{Qh!)#i>lz{rDopb&8Q-c`}rJyWV_yumZa;W
o7)jSt7)jUnFp{o+!3bS5-7nC#eFP(EdkQ18t>3hs>=`8g0+CpQ&Hw-a
diff --git a/docs/build/doctrees/info.doctree b/docs/build/doctrees/info.doctree
index 8f40b31685afa109faa54d8fe7e594b92c847d04..1302d15f3a1035ebf2397664b61b643aac070245 100755
GIT binary patch
delta 3503
zcma)8du&tJ8RuO4`V|OCFqqg$>;xyaLxR%~W1tBEoJUE(ETv!yNo*%}a$}O%&Vzu#5A>HoBD!kYMD?e6WdA|)4p@=bpk4X
z^q*sYzwh-s-}%mWuJ{sn^&rPT%O9C~VUk}=a-2>f+7ov61P2D3HGR=oTXe|T9u0+^
zp^hjtkq=1>a=GT*_Gs_U@Q?+6`Um=gUA+VSjwMM!-I5gMCzyr)Q5L`{I*A|VJeCc2
z6xA6+mgcTriznqbtx(Gf@F;;ChpWo&3f&XozOZE_(h^n%JON8rza?ntjrNv?dtw8-
zEWzH8rM)}aAI2YvS$}vyT+2UANCo^{wIoaOPmU!=lfy~7z$X}-Cc}bJXot~abB-sF
zB11eg#)1kBk`PxFw-)`V%u8l`6F=62akx-cENiZ@{ygPR(;UYX7Gt|SfESr|0
zloR|Z(+=-tE|ru<$S`Q&Sg{G7n`e{b`LjNiJ-KxRg6<+EOY_FO9Tfhl#YOatf?wZ9
z^4n-mA}9otSIF@WwngS5;7A?W;^Jb+Utog;3)a9=!y;{pVjDcWV34#SJ4#>?d}Xks
z)NJrDtR9Xpq;kf>P181vB=(o
z`aJB;$tUM&VYrW^aOYum?iGYS$-T+lATYVP5N_q|(|jeB6CKQdt-6Q8px}u77k?AF
zzHkRNcW-+*lx!qe^7B+n51hk0HPf+FqFTmtAzZSY(O5ZlMg?%FqEw%bN;3sm0fk~A
zr_m(7pI@&aWsq4IhrLC1{P#xTMkX6!eo-gvw%K)M=)=;LL*}A#_@Zd7Y*1yZ;`|&v
zv%F~=vho2&C}7Sq4h~r=_#B!Z?%MjfAr5zKhL}|8l59ez9?JAF
z%S5KcGfoBJrljb_EX;Rhx43T}4pYJuwdF^&oPyPNmIBgCEBG=NK4)nU(Pl`BbV*Up
z96qcF5)+iVMzjt|Yopo{b}m+m%$$o${Ff`85EmS;F+!eu4Os;Nw+-cSx07=zSUs0v
zP%Fe6^`Kfl%J~!z)pR*i1~=-(yUQmOqzeNTz#l5jFlUt=rP)>Q%GftmW^hzfY}u+w
z8Pipp$x)h6(}lFyH9lG02(MTBxGpJo^>p}#)#V(SWZ{3g)(q=Cc9ha;7ctBVPERiA
zYbZA8aWl*a?|C-Bb%s^dv@>iCJYPHG=dATXK^>J()w*OKIycmHCb-G)#yUU43*e7+
zW>NwdytcwgnQDS(!+xd0X9EY~O*Q0Qv@4jNZ)e5Nn%s)4
zvoz6z!|`jdZ5}6oX56l9)~No529fXK)6I=@>BFGn()o$74WO(+LVo=Wdq!4vE%!&tEDOh
z(ngsRFu_}`RFT`dE~OZ4bwRJ+3df3+sQA2f>7xi){S9Cb$b{f${!J-wIV?E{O)WO`
zZwjZa-_UM`$3u290f*azk|G~g>?nl%utZ`4LZPR#CVp6z
zE1@UhREI=jLTm`9kyPARM>RatNn`${V?zpD-RXuwfr7_6M^j*>uma!n$yfoDt>RJj
z3!#y!DBq_fzJ*=43FlB@wIy1u-FD)D7j_(!z&bGYFrd3h01U8EUBKQtJIPh#6p;f^~p_E5;g
z9;L0QibE`(t+_)xnhG=|g%vzZ#rKjX)_k3P4((|vlBlB1XQVIO?9SV@|;pHJm>b}T^
z_r*&?22y#Oy8Yn7NH+9!F0mvZc%4x~db0Ly$lkq^Q>yXtEMV<_72etHnW>(W%OFslGaBa|{Ptug~{Y~+K>w97b8yY@>4s`7`DXAeijds2jIb4L(dzTj;
zMktt2hobGh!JcqZXraDG5tNUig%%UCh67}~*gxFH8HX``B0b#K7wz-o*GM$f)!Ug&
z@Ph*#;;oU(imY=qe_K$qs75X5k%)+^56qCS7jO(d98J
zF2Q9mp{4sA*cA)+(@*?}Zv;^#tIKOW|eV$64rc7Iu74Fw;PE
z2SOU7FD#b0{|}E#;!Z!fB<}Qym$-jPH*i9=3|BfTBlIL9w0hDuXa`H%xH!Fyf6{++
g8xcD+jX5;$;9sopIe2bth+KqwV}5b~wm-l1|8XWLKmY&$
delta 3275
zcmZuzdr(y875CibF1tXGH!KS*%X1gRRRn=3Dsce`vS3jGc`OSoTy9{&1q5F~Qzgs<
z$Ck6i)KseV5jCmCxqS?2C&rG`OebTLPLp(O(ubW`jft(M?M$Z8w%_^gu4v3Z-@WJg
z`$Ghea(jnteZ8Kp0h@iGbEMzx86L9r
z^bX^Ccn3WAU7^L;F)-rkO6?r*3=Ml-z07X)Nx`i?DRzq6qaL>-XnH`x-)b5Ssjk$P
zKKIttG4;_A2<6s@uUm0R5JDeb4(e0fVk5v}d{etP%(vgS*SFhek;Eyk5Pedj)QK0&
z`N|^QT7bhisBeXnxVJbp^f)se9_e7@+VO5vJVu1fF?(b}JYEa2;uXDF>?03c5}M-I
zcZcR-NN6UL=0e|RD@zLtl}ll>ApVWHeu9yje2BzLVKK_9;b%dRzQ8RJeS&nK<9)?k
zzTkKV@i*ZI@rB4)@*u$O&_(#GhzvYh97L+wD5q*M;jF=k-l%L6rRYCC?7t$~LgE`S
z&(MRYM04b=v0jL(g#`;FkTjCfAtgB&=so;0ZWtP;4e27|A=lzd@gqVl;H6RnI)Zh$
zB*BP%@#&0#sH0%wiQi91KupL8_z}ye6Rc2+7DIhNNx?+&e+d?-RsL#d0A@dDE`(Ul
zW}mT$_2C}AmGl~%3!tCHmy*rGSrwhKK9fwfi*OcmQr>0je^NddJ_J0rD-)kc{kHCX
zswPopVX0S5H;J8Ol|N(aJGfLU9uR3WC72c3$#YYBTT@^l4#7eI>F6Sny@7(1Z~S~Vlo+(n^py2
zl%x6^2zv-QWHf-X4O??Y@wr9&g;AbWG*cixm3K}5vRcvi1cZVDZF&GL26rsp1Ft-S
zmdVi?Z3ck3N6>yrXu3H+8n!PT!-l2%ggLbWYanMd=H!0_A3c&YmlNq+HlzQCO5`>t
zBHB3ClLh^tzpom;QVmOaz^8>SNW}i36zpB@(#L7|uEc2gu0)bFzbgewbWu1!8opJU
zh%1ZBAQzz`o5`n&Z9&g=G6QyIobPWD=@reBztF&D7GDNV9=jQHR$b(p
zbbMlU3c@N&$YzGhcXT~oSP_SVE0*J#)h&J}>zWkIU2S2~;F>eSW)*Z5Prk9NLV#J8
zJRRR$o5DIy~fAXB#;aS|cu81|T
ztdcnXc-@dNtJboH*V0in41zuA!5m;$8EUw*T`8@O2i9G4W1+a#-|ZMYzJ3O}nZta?
z%en2vnjzheAik+%80jb!t9bru#|984D8Csm)|%iT±`b)_)DQ&)4p@w#w<=EU^3
z?C=`%QODG2MW`>*M4hI*4^1=q1FYUW2^Tqyi@S8%x0!={1EHaMp?_9T$Vi8L3!58o
z*}Ek8M~3_*n;U#os4i^$!{Yh2k-!L@=G1C7a#`?%tAPAy;%(r4P+~fY1Yt%k!_#FF
zW~tgNn8BZRxm3LYH}&P=fi4SEPP&_tW*%#kHSqm5;JFS1#&z4!=~k)DD4y<62qg0L
ztj6+g&Wk0x*HaTPbV)_{yB_+j~p!Wyzu=K;HGg_a@=>UJKc&f_k-rD)}x5
zaaw*I7*oNd
At each step, we compute a trial step \(s_k\) designed to make our approximation \(m_k(s)\) small (this task is called the trust region subproblem). We evaluate the objective at this new point, and if this provided a good decrease in the objective, we take the step (\(x_{k+1}=x_k+s_k\)), otherwise we stay put (\(x_{k+1}=x_k\)). Based on this information, we choose a new value \(\Delta_{k+1}\), and repeat the process.
In DFO-LS, we construct our approximation \(m_k(s)\) by interpolating a linear approximation for each residual \(r_i(x)\) at several points close to \(x_k\). To make sure our interpolated model is accurate, we need to regularly check that the points are well-spaced, and move them if they aren’t (i.e. improve the geometry of our interpolation points).
-A complete description of the DFO-LS algorithm is given in our paper [CFMR2018].
+A complete description of the DFO-LS algorithm is given in our paper [CFMR2018].
References
diff --git a/docs/build/html/install.html b/docs/build/html/install.html
index d74c88f..b9a3e9a 100644
--- a/docs/build/html/install.html
+++ b/docs/build/html/install.html
@@ -8,7 +8,7 @@
-
Installing DFO-LS — DFO-LS v1.2.3 documentation
+
Installing DFO-LS — DFO-LS v1.3.0 documentation
@@ -61,7 +61,7 @@
- 1.2.3
+ 1.3.0
diff --git a/docs/build/html/objects.inv b/docs/build/html/objects.inv
index 0a40ac70446e2f79c64912950d2f8f543bdf4bdf..540a6f6429ef5e575695ecea53587d35cac46bb8 100644
GIT binary patch
delta 292
zcmV+<0o(qO1CIlcJOeW>Fp)nzf6r^fAQZm)S3KCRj&-Ne+aQz`T2ctz-No1BfzeJ)
zw(PH;iDor0*iFA*zW1S!gU~ffa9(busex}mBXz)TAYTtshv>Xx&dof7X)o#*X
D2V**ZULzuJoDvXGAe
delta 294
zcmV+>0oneK1CaxeJOeT=Gm$?$e@jcmFc7}yR|NDLaXrc2gtEA>OBagXV>4|B^4Lw%
zYW?-5O=}Ynyd~c&Gn3(B;JOAuofmsyYRYz$FW5upm837L;?40ENZly=^tgMT4U&I=
zoGiTq+x946_T1f5Mg;O!8|Sf3vo|!R7n|y+8)FssXf$a-qZaU$zp|pye|nsbq(kS8
z9n&`!D5~DUx{yd;BS^N`MV%xThmxJDJuSfE5DFQ99V`GXuO1bCgQcNDdIIh|t8P8F
z;OH6k#Kp2Pw2)^vCPtS7*a6`)NSx~*O9G6G@gy7#2NM
- Search — DFO-LS v1.2.3 documentation
+ Search — DFO-LS v1.3.0 documentation
@@ -60,7 +60,7 @@
- 1.2.3
+ 1.3.0
diff --git a/docs/build/html/searchindex.js b/docs/build/html/searchindex.js
index fd398c6..9520167 100644
--- a/docs/build/html/searchindex.js
+++ b/docs/build/html/searchindex.js
@@ -1 +1 @@
-Search.setIndex({docnames:["advanced","diagnostic","history","index","info","install","userguide"],envversion:{"sphinx.domains.c":2,"sphinx.domains.changeset":1,"sphinx.domains.citation":1,"sphinx.domains.cpp":2,"sphinx.domains.index":1,"sphinx.domains.javascript":2,"sphinx.domains.math":2,"sphinx.domains.python":2,"sphinx.domains.rst":2,"sphinx.domains.std":1,sphinx:56},filenames:["advanced.rst","diagnostic.rst","history.rst","index.rst","info.rst","install.rst","userguide.rst"],objects:{},objnames:{},objtypes:{},terms:{"00000000e":6,"00180000e":6,"00670067e":6,"00e":6,"01256863e":6,"01941978e":6,"04752848e":6,"06227943e":6,"07858081e":6,"09280752e":6,"09843514e":6,"10862447e":6,"10g":6,"12897463e":6,"19971362e":6,"1e20":6,"20e":6,"22992989e":6,"24129962e":6,"24991776e":6,"33017181e":6,"34676365e":6,"35e":6,"41166495e":6,"42808544e":6,"43e":6,"45394186e":6,"47252555e":6,"50e":6,"51525603e":6,"550476524e":6,"59085679e":6,"59634974e":6,"61e":6,"62398179e":6,"63036198e":6,"650274685e":6,"65956876e":6,"70205419e":6,"70749826e":6,"71355033e":6,"75304364e":6,"77e":6,"79999999e":6,"80e":6,"827884295e":6,"83184867e":6,"90335675e":6,"95045269e":6,"95108811e":6,"96161167e":6,"97239623e":6,"98196347e":6,"98830861e":6,"99950530e":6,"99999998e":6,"case":[0,4,6],"const":0,"default":[0,2,5,6],"float":6,"function":[0,4,6],"import":6,"long":2,"new":[0,4],"public":3,"return":[1,2,6],"short":6,"throw":0,"true":[0,1,4,6],"try":[0,3,4],"while":[0,4],Adding:3,And:6,For:[0,3,4,5,6],Such:4,That:[3,4],The:[0,1,3,6],Then:4,There:[0,4],These:6,Use:2,Using:[0,1,3],__future__:6,a_i:6,abl:6,about:[1,6],abov:[4,5,6],abs_tol:0,absolut:0,accept:0,access:6,accur:[4,6],accuraci:6,achiev:6,acm:[0,3,4,6],actual:1,add:[0,6],added:0,adding:0,addit:0,addition:5,additive_noise_level:0,adjust:[2,6],admin:5,advanc:[3,6],after:0,algebra:6,algorithm:[2,3,6],all:[0,1,4,6],allow:[0,6],along:6,alpha1:0,alpha2:0,alpha_1:0,alpha_2:0,also:[0,4,5,6],altern:[0,3,5,6],alwai:[0,4],amount:[0,6],ani:[0,1,3,6],anoth:0,appli:[0,6],applic:6,approach:0,approx:[4,6],approxim:[4,6],apr:3,aren:4,arg:[2,6],argument:[2,3],arrai:6,ask:[0,3,6],attempt:4,author:3,auto:2,auto_detect:0,automat:[0,5,6],avail:[0,6],averag:6,avoid:[0,2],awai:6,axes:6,b_i:6,base:[0,2,3,4,6],bash:5,basicconfig:6,becom:6,been:[1,6],befor:[0,6],behavior:2,below:[0,6],benjamin:[0,4,6],best:[0,1,4,6],between:2,block:5,bobyqa:3,both:6,bound:[0,1,2,3,4],bug:2,bugfix:2,build:[0,6],calcul:[2,4,6],calibr:[4,6],call:[0,1,4,6],can:[0,4,5,6],cannot:[0,4],cap:0,carlo:4,carti:[0,3,4,6],categori:4,caus:[0,6],cdl:6,cdot:[0,6],centr:3,certainli:4,cfmr2018:[0,4,6],chang:[0,1,2],check:4,check_objfun_for_overflow:0,choic:4,choos:4,clone:5,close:[2,4],closest:4,code:[2,5,6],coeffici:0,collabor:3,column:1,com:[5,6],common:4,commonli:6,compar:6,compil:5,complet:4,compon:0,comput:[1,4],computation:0,condit:[0,1,6],consecut:0,consid:4,consist:0,constrain:6,constraint:[3,4,6],construct:[4,6],contact:3,contain:[1,5,6],content:3,coordin:[0,6],coralia:[0,3,4,6],correct:[2,6],correctli:[2,6],correl:0,cost:[0,2,6],could:6,count:3,creat:2,criteria:0,criterion:0,crvmin:2,csv:1,current:[0,3,4,6],customis:2,data:[0,3,4],datafram:[1,6],date:3,deactiv:6,dec:0,decai:6,decreas:[0,4,6],def:6,defin:[1,6],deflat:4,delta:[1,6],delta_:4,delta_k:[0,1,4],delta_scale_new_dirn:0,demonstr:6,depend:[0,1,2,4,6],deprec:2,depth:6,deriv:[0,4,6],describ:[0,4,6],descript:[1,4,6],design:[4,6],detail:[0,3,6],detect:2,determin:[0,4,6],determinist:2,develop:[3,4],dfo:[0,1,2],dfol:[0,1,5,6],diagnost:[0,3,6],diagnostic_info:[1,6],dictionari:[0,6],did:[0,6],differ:[0,2,4,6],differenc:4,difficult:4,dimension:6,direct:[0,6],directli:4,directori:[5,6],disciplin:4,displai:6,distanc:1,divid:2,do_geom_step:0,do_log:6,do_safety_step:0,doctor:3,document:[5,6],doe:[4,6],doi:2,doing:4,don:2,download:5,dure:[0,2],dynam:3,each:[0,1,4,6],easi:[5,6],easy_instal:5,effect:0,either:6,enabl:2,encount:0,end:0,enough:[0,6],ensur:2,entri:6,epsrc:3,equal:1,equat:3,error:[0,1,4,6],estim:[3,4],eta1:0,eta2:0,eta_1:0,eta_2:0,etc:[0,1],eval:6,evalu:[0,1,2,3,4],even:4,everi:[4,6],exactli:[1,6],exampl:[3,4,5],exclud:0,exist:4,exit:[2,6],exit_false_success_warn:6,exit_input_error:6,exit_linalg_error:6,exit_maxfun_warn:6,exit_slow_warn:6,exit_success:6,exit_tr_increase_error:6,exp:6,expect:6,expens:[0,1,4],experi:4,explain:1,exponenti:6,extend:6,extra:[0,2,6],factor:0,fals:[0,6],far:[1,6],fast:5,faster:[2,6],feasibl:2,featur:3,feb:3,februari:[],fewest:4,fiala:[0,3,4,6],figur:6,file:[0,1,5,6],filemod:6,filenam:6,find:[3,4,6],finish:6,finit:4,first:6,fit:[0,3],fix:[0,2,6],flag:[0,6],flexibl:[0,3,4,6],floor:0,focus:3,follow:[5,6],form:4,format:6,formul:4,fortran:[2,5],found:[0,1,6],framework:6,free:[0,4,6],frobeniu:1,from:[1,2,4,5,6],full:[0,2,3,6],full_geom_step:0,full_rank:0,fulli:1,further:0,g_k:1,gamma:0,gamma_:0,gamma_dec:0,gamma_inc:0,gamma_inc_overlin:0,gaussian:6,gca:6,gener:[3,6],geometri:[0,2,4],geq:0,get:[4,6],gfortran:5,git:5,github:[5,6],give:[0,6],given:[0,3,4,6],gnu:3,good:[4,6],gracefulli:0,grad:6,gradient:1,grid:6,group:3,grow:3,guarante:4,had:1,hand:[2,4,5],handl:2,happen:4,hard:0,has:[0,1,2,3,4,6],have:[0,4,5,6],help:6,here:6,hessian:2,high:6,higher:5,histori:[0,3],history_for_slow:0,home:5,how:[0,1,3,5],howev:[4,6],htm:6,html:[5,6],http:[5,6],idea:4,imlug:6,imlug_genstatexpls_sect004:6,impact:[2,6],imposs:4,impract:0,improv:[0,2,3,4,6],inaccur:4,inc:0,includ:[0,1,6],increas:[0,6],increase_ndirs_initial_amt:0,increase_npt:0,increase_npt_amt:0,increase_num_extra_steps_with_restart:0,indic:6,industri:3,infinit:4,info:[1,6],inform:[3,4,6],infti:6,init:0,initi:[2,3,6],initialis:[2,6],input:[0,2,4,6],inspect:0,instal:[2,3],instanc:[0,4],instead:[4,5],integ:6,interest:[3,4],interfac:6,intermedi:0,intern:[0,4,6],interpol:[2,3,4,6],interpolation_change_j_norm:1,interpolation_condition_numb:1,interpolation_error:1,interpolation_total_residu:1,introduc:0,invers:2,involv:[0,4],iter:[0,2,3,6],iter_this_run:1,iter_typ:1,iters_tot:1,its:[0,6],j_k:[0,1],jacobian:[0,1,2,6],jan:[0,3,4,6],jun:3,june:3,just:2,keep:4,known:4,l015803:3,label:6,lambda:[0,1],larg:[0,4],larger:0,largest:0,last:[0,1],lastli:6,latest:5,ldot:4,least:4,least_squar:6,left:4,legend:6,len:[0,6],length:0,leq:[0,3,4,6],less:0,let:6,level:[0,5,6],librari:4,licens:3,like:[4,6],linalg:0,linalgerror:0,lindon:[0,3,4,6],line:[0,6],linear:[0,1,4,6],link:2,linspac:6,list:[2,6],loc:6,local:[3,6],locat:5,log:[1,2,3,6],longer:2,lower:[0,1,3,6],lsqcurvefit:6,m_k:4,magnitud:[0,6],mai:[1,3,4,6],main:[4,6],maintain:4,major:4,make:[0,2,4,5,6],manag:3,mani:[0,3,4,6],manual:3,marteau:[0,3,4,6],math:6,mathbb:[3,4,6],mathemat:[0,3,4,6],mathrm:4,mathwork:6,matplotlib:6,matrix:[1,6],max:6,max_distance_xk:1,max_fake_successful_step:0,max_npt:0,max_slow_it:0,max_unsuccessful_restart:0,maxfun:[0,6],maximum:[0,1,6],mean:[1,4],measur:[0,4],messag:6,method:[0,2,4],min:6,min_:[3,4,6],min_chgj_slop:0,min_correl:0,min_sing_v:0,minim:[2,4],minimum:[0,6],minor:2,model:[2,3,4,6],modifi:6,modul:6,momentum:0,momentum_extra_step:0,mont:4,more:[0,3,4],most:[0,1,4,6],move:[0,4,6],move_xk:0,msg:6,multipl:[2,3,4,6],multiplicative_noise_level:0,must:6,myfil:[1,6],n_to_print_whole_x_vector:0,nag:3,nan:[0,2],navig:5,ndirs_initi:0,necessari:2,need:[4,6],never:[3,6],next:6,nfev:6,nois:[3,4,6],noisi:[0,3,4],non:[3,6],nonconvex:4,none:[0,6],nonlinear:3,nonlinear_system:6,nonzero:0,nor:4,norm:1,norm_gk:1,norm_sk:1,normal:6,note:[5,6],now:[0,6],npt:[0,1,6],nrestart:6,nrun:[1,6],nsampl:[1,6],num_extra_step:0,num_geom_step:0,num_new_dirns_each_it:0,number:[0,1,4,6],numer:3,numericalalgorithmsgroup:5,numpi:[0,2,5,6],obj:6,object:[1,2,3,4],objfun:[0,1,2,6],objfun_has_nois:[0,6],obs:4,observ:[4,6],occur:6,often:4,older:5,one:[0,1,4,5,6],onli:[0,1,4,6],oper:2,oppos:0,opposit:0,opt:6,optim:[0,1,4,6],option:[0,2,3,4,5],order:6,org:5,origin:6,orthogon:0,other:0,otherwis:[0,4],our:[3,4,6],out:6,output:[1,2,3,4],outsid:[3,6],over:[0,6],overdetermin:4,overflow:2,overflowerror:0,overlin:0,overrid:0,overridden:6,overview:[3,6],packag:[2,3,5],page:6,panda:[1,5,6],paper:[0,3,4,6],param1:6,param2:6,paramet:[2,3],part:[2,6],partial:6,particularli:[4,6],pass:[2,6],past:0,per:[0,1,2,6],perform:[0,2,4,6],perturb:0,perturb_trust_region_step:0,phase:[0,6],physic:4,piec:[0,1],pip:3,pleas:3,plot:6,plt:6,point:[1,2,3,4,6],pois:1,poised:[0,1],popul:1,possibl:[4,6],post:0,practic:6,precondit:0,predict:[1,4],prediction_error:6,prefer:4,preprint:[0,3,4,6],present:5,preserv:0,previou:[0,6],previous:1,print:[0,2,6],print_funct:6,print_progress:6,privat:5,privileg:5,probabl:4,problem:[0,1,2,3,4,6],proccess:4,process:[0,4],produc:[4,6],progress:[3,6],prohibit:4,provid:[4,5,6],pull:5,pure:5,purpos:6,put:4,pydata:5,pyplot:6,python:[0,5,6],quad:[3,4,6],quantit:4,quantiti:4,quit:0,quit_on_noise_level:0,r_1:[0,4,6],r_2:4,r_i:[0,4,6],r_m:[0,4,6],radii:0,radiu:[0,1,2,4,6],random:[0,2,6],random_directions_make_orthogon:0,random_initial_direct:0,rang:6,rank:0,rate:0,rather:2,ratio:[0,1],reach:[0,6],recommend:0,recycl:0,reduc:[0,2],reduce_delta:0,reduct:1,refer:3,regim:2,region:[2,3,4,5,6],regress:3,regular:0,regularli:4,rel:0,rel_tol:0,relationship:4,relax:[3,6],releas:[2,3],remov:[2,5],repeat:4,replac:6,reproduc:[2,6],request:6,requir:[0,3,4,6],rerun:5,reset:0,reset_delta:0,reset_rho:0,resid:6,residu:[1,4,6],respect:0,restart:[1,2,3,6],result:[0,1,2,4,6],retriev:2,rho:[1,6],rho_:0,rho_k:[0,1],rhobeg:6,rhoend:6,rhoend_scal:0,right:[2,6],risk:0,robert:[0,3,4,6],robust:[0,3,4,6],root:5,rosenbrock:6,rosenbrock_noisi:6,rounding_error_const:0,roundoff:0,row:1,run:[0,5,6],run_in_parallel:0,runtim:2,runtimewarn:6,s_k:[0,1,4],safeti:[0,1,2],safety_step_thresh:0,sai:4,same:[0,3,4,6],sampl:6,sas:6,satisfi:[4,6],save:[0,1,6],save_diagnostic_info:[0,1],save_poised:[0,1],save_rk:[0,1],save_xk:[0,1],saw:1,scale:[0,2,6],scale_factor:0,scale_factor_for_quit:0,scaling_within_bound:6,scipi:[4,5,6],screen:0,script:6,search:[0,2,4],section:[0,1,2,6],see:[1,5,6],seed:[2,6],sens:4,sensibl:6,set:[1,2,3,4,6],set_xlabel:6,set_ylabel:6,setup:[0,5],sever:[0,4,6],shape:6,shift:[0,6],should:[0,5,6],show:6,side:[2,4],similar:4,similarli:4,simpl:[3,5],simpli:4,simul:4,sin:4,sinc:[0,1],singular:[0,6],site:5,situat:4,size:[4,6],slope:0,slow:[1,3,6],slow_it:1,small:[3,4,6],smaller:0,smallest:[0,1],smooth:0,soft:0,softwar:[0,3,4,5,6],soln:[1,6],solut:[2,3,4,5,6],solv:[0,1,3],solver:[0,2,3,4,6],some:[0,4,5,6],sourc:5,space:[0,2,4],specifi:[0,6],squar:[1,4],stai:4,stand:3,start:[2,4,6],statu:6,step:[0,1,2,4,6],still:0,stochast:[3,6],stop:2,store:[0,2],str:6,string:6,strongli:4,structur:3,subproblem:[2,4,5,6],success:[0,1,6],successfulli:6,sudo:5,suffici:6,suitabl:4,sum:[1,4],sum_:[3,4,6],supervis:3,support:[3,5,6],suppos:[4,6],sure:[4,5],svd_max_jac_cond:0,svd_scale_factor:0,system:[0,1,3,5],t_i:6,tabl:6,take:[4,6],taken:6,task:4,tdata:6,technic:6,techniqu:4,termin:[1,3,6],test:[3,6],text:[0,1,3,4,6],than:[0,2,6],thei:[0,4],them:[2,4],therefor:6,thi:[0,1,2,3,4,5,6],thresh_for_slow:0,threshold:0,through:[0,2],throw_error_on_nan:0,time:[0,1,4,6],to_csv:1,toler:0,top:5,total:[1,6],tr_radiu:0,train:3,transact:[0,3,4,6],trend:0,trial:4,triangular:2,trigger:0,troubl:6,trust:[2,3,4,5,6],trustregion:[2,5],tupl:6,turn:1,two:[0,4],type:[1,4],typic:4,unabl:6,under:3,underdetermin:4,uninstal:3,unknown:4,unless:6,unpack:5,unsuccess:0,updat:2,upgrad:5,upper:[2,3,6],usag:[3,6],use:[0,1,2,3,5],use_full_rank_interp:0,use_old_rk:0,use_restart:0,use_soft_restart:0,used:[0,1,4,6],useful:[0,4,6],user:[0,4,5,6],user_param:[0,1,6],uses:0,using:[0,2,3,4,6],usual:0,val1:6,val2:6,valu:[1,2,3,4,6],vari:0,variabl:[3,6],vast:4,vdot:4,vector:[0,1,4,6],veri:[0,4,6],version:[0,3,5],via:6,viewer:6,visibl:6,wai:[0,4,6],want:[5,6],warn:[2,6],well:4,what:1,when:[0,2,3,6],where:[0,4,6],whether:[0,4,6],which:[0,1,2,3,4,5,6],whole:2,why:6,wish:[3,4,6],within:[0,2,6],without:[0,3],work:5,wors:6,write:6,written:5,www:5,x_0:[0,6],x_1:[4,6],x_2:[4,6],x_b:0,x_j:6,x_k:[0,4],x_n:4,xmin:6,xtol:6,y_1:4,y_i:[4,6],y_k:[0,1],y_m:4,y_t:0,ydata:6,yet:0,you:[3,4,5,6],your:[4,5,6],zenodo:2,zero:2},titles:["Advanced Usage","Diagnostic Information","Version History","DFO-LS: Derivative-Free Optimizer for Least-Squares Minimization","Overview","Installing DFO-LS","Using DFO-LS"],titleterms:{Adding:6,Using:6,acknowledg:3,advanc:0,algorithm:[0,1,4],apr:2,argument:6,bound:6,count:1,current:1,data:6,deriv:3,detail:4,dfo:[3,4,5,6],diagnost:1,dynam:0,equat:[4,6],estim:6,evalu:6,exampl:6,feb:2,fit:[4,6],free:3,gener:0,grow:0,histori:2,how:6,inform:[0,1],initi:0,instal:5,interpol:[0,1],iter:1,jan:2,jun:2,least:[3,6],log:0,manag:0,manual:5,minim:[3,6],model:[0,1],more:6,multipl:0,nois:0,noisi:6,nonlinear:[4,6],object:[0,6],optim:3,option:6,output:[0,6],overview:4,paramet:[0,4,6],pip:5,point:0,progress:[0,1],refer:[0,4,6],region:[0,1],regress:0,requir:5,restart:0,set:0,simpl:6,slow:0,small:0,solv:[4,6],squar:[3,6],stochast:0,system:[4,6],termin:0,test:5,trust:[0,1],uninstal:5,usag:0,use:[4,6],using:5,valu:0,version:2,when:4}})
\ No newline at end of file
+Search.setIndex({docnames:["advanced","diagnostic","history","index","info","install","userguide"],envversion:{"sphinx.domains.c":2,"sphinx.domains.changeset":1,"sphinx.domains.citation":1,"sphinx.domains.cpp":2,"sphinx.domains.index":1,"sphinx.domains.javascript":2,"sphinx.domains.math":2,"sphinx.domains.python":2,"sphinx.domains.rst":2,"sphinx.domains.std":1,sphinx:56},filenames:["advanced.rst","diagnostic.rst","history.rst","index.rst","info.rst","install.rst","userguide.rst"],objects:{},objnames:{},objtypes:{},terms:{"00000000e":6,"00004412e":6,"00180000e":6,"00670067e":6,"00e":6,"01256863e":6,"01941978e":6,"04752848e":6,"06227943e":6,"07858081e":6,"09280752e":6,"09843514e":6,"10862447e":6,"10g":6,"12897463e":6,"19971362e":6,"1e20":6,"20e":6,"22992989e":6,"24129962e":6,"24991776e":6,"33017181e":6,"34676365e":6,"35e":6,"41166495e":6,"42808544e":6,"43e":6,"45394186e":6,"47252555e":6,"50e":6,"51525603e":6,"550476524e":6,"59085679e":6,"59634974e":6,"61e":6,"62398179e":6,"63036198e":6,"650274685e":6,"65956876e":6,"70205419e":6,"70749826e":6,"71355033e":6,"75304364e":6,"77e":6,"79826221e":6,"79999999e":6,"80e":6,"81262102e":6,"827884295e":6,"83184867e":6,"90335675e":6,"95045269e":6,"95108811e":6,"96161167e":6,"97239623e":6,"98196347e":6,"98830861e":6,"99950530e":6,"99999998e":6,"case":[0,4,6],"const":0,"default":[0,2,5,6],"float":6,"function":[0,4,6],"import":6,"long":2,"new":[0,4,6],"public":3,"return":[1,2,6],"short":6,"throw":0,"true":[0,1,4,6],"try":[0,3,4],"while":[0,4],Adding:3,And:6,For:[0,3,4,5,6],Such:4,That:[3,4],The:[0,1,3,6],Then:4,There:[0,4],These:6,Use:2,Using:[0,1,3],__future__:6,a_i:6,abl:6,about:[1,6],abov:[4,5,6],abs_tol:0,absolut:0,accept:0,access:6,accur:[4,6],accuraci:6,achiev:6,acm:[0,3,4,6],actual:1,add:[0,6],added:0,adding:0,addit:[0,2],addition:5,additive_noise_level:0,adjust:[2,6],admin:5,advanc:[3,6],after:0,algebra:6,algorithm:[2,3,6],all:[0,1,4,6],allow:[0,6],alon:6,along:6,alpha1:0,alpha2:0,alpha_1:0,alpha_2:0,also:[0,4,5,6],altern:[0,3,5,6],alwai:[0,4],amount:[0,6],ani:[0,1,3,6],anoth:0,appli:[0,6],applic:6,approach:0,approx:[4,6],approxim:[4,6],apr:3,arbitrari:2,aren:4,arg:[2,6],argument:[2,3],arrai:6,ask:[0,3,6],assum:3,attempt:4,author:3,auto:2,auto_detect:0,automat:[0,5,6],avail:[0,6],averag:6,avoid:[0,2],awai:6,axes:6,b_i:6,ball:6,base:[0,2,3,4,6],bash:5,basicconfig:6,becom:6,been:[1,6],befor:[0,6],behavior:2,being:0,below:[0,6],benjamin:[0,4,6],best:[0,1,4,6],between:2,block:5,bobyqa:3,both:6,bound:[0,1,2,6],box:6,bug:2,bugfix:2,build:[0,6],calcul:[2,4,6],calibr:[4,6],call:[0,1,4,6],can:[0,4,5,6],cannot:[0,4],cap:0,carlo:4,carti:[0,3,4,6],categori:4,caus:[0,6],cdl:6,cdot:[0,6],centr:[3,6],certainli:4,cfmr2018:[0,4,6],chang:[0,1,2],check:[3,4],check_objfun_for_overflow:0,choic:4,choos:4,clone:5,close:[2,3,4,6],closest:4,code:[2,5,6],coeffici:0,collabor:3,column:1,com:[5,6],common:4,commonli:6,compar:6,compil:5,complet:4,complex:6,compon:0,comput:[1,4],computation:0,condit:[0,1,6],consecut:0,consid:[0,4],consist:0,constrain:6,constraint:[2,3,4],construct:[4,6],contact:3,contain:[1,5,6],content:3,convex:[2,3,4,6],coordin:[0,6],coralia:[0,3,4,6],correct:[2,6],correctli:[2,6],correl:0,correspond:6,cost:[0,2,6],could:6,count:3,creat:2,criteria:0,criterion:0,crvmin:2,csv:1,current:[0,3,4,6],customis:2,d_tol:0,data:[0,3,4],datafram:[1,6],date:3,deactiv:6,dec:0,decai:6,decreas:[0,4,6],def:6,defin:[1,6],deflat:4,delta:[1,6],delta_:4,delta_k:[0,1,4],delta_scale_new_dirn:0,demonstr:6,depend:[0,1,2,4,6],deprec:2,depth:6,deriv:[0,4,6],describ:[0,4,6],descript:[1,4,6],design:[4,6],detail:[0,3,6],detect:2,determin:[0,4,6],determinist:2,develop:[3,4],dfo:[0,1,2],dfol:[0,1,5,6],diagnost:[0,3,6],diagnostic_info:[1,6],diagon:0,dictionari:[0,6],did:[0,6],differ:[0,2,4,6],differenc:4,difficult:4,dimension:6,direct:[0,6],directli:4,directori:[5,6],disciplin:4,displai:6,distanc:1,divid:2,do_geom_step:0,do_log:6,do_safety_step:0,doctor:3,document:[5,6],doe:[4,6],doi:2,doing:4,don:2,download:5,dure:[0,2],dykstra:3,dynam:3,each:[0,1,4,6],easi:[5,6],easy_instal:5,effect:0,either:6,empti:[3,4,6],enabl:2,encount:0,end:0,enough:[0,6],ensur:2,entri:[0,6],epsrc:3,equal:1,equat:3,error:[0,1,4,6],estim:[3,4],eta1:0,eta2:0,eta_1:0,eta_2:0,etc:[0,1],eval:6,evalu:[0,1,2,3,4],even:4,everi:[4,6],exactli:[1,6],exampl:[3,4,5],exclud:0,exist:4,exit:[2,6],exit_false_success_warn:6,exit_input_error:6,exit_linalg_error:6,exit_maxfun_warn:6,exit_slow_warn:6,exit_success:6,exit_tr_increase_error:6,exp:6,expect:6,expens:[0,1,4],experi:4,explain:1,exponenti:6,extend:6,extra:[0,2,6],factor:0,fals:[0,6],far:[1,6],fast:5,faster:[2,6],feasibl:[2,3,6],featur:3,feb:3,februari:[],fewest:4,fiala:[0,3,4,6],figur:6,file:[0,1,5,6],filemod:6,filenam:6,find:[3,4,6],finish:6,finit:[2,4],first:6,fit:[0,3],fix:[0,2,6],flag:[0,6],flexibl:[0,3,4,6],floor:0,focus:3,follow:[5,6],form:4,format:6,formul:4,fortran:[2,5],found:[0,1,6],framework:6,free:[0,4,6],frobeniu:1,from:[1,2,4,5,6],full:[0,2,3,6],full_geom_step:0,full_rank:0,fulli:1,further:0,g_k:1,gamma:0,gamma_:0,gamma_dec:0,gamma_inc:0,gamma_inc_overlin:0,gaussian:6,gave:[],gca:6,gener:[3,6],geometri:[0,2,4],geq:0,get:[4,6],gfortran:5,git:5,github:[5,6],give:[0,6],given:[0,3,4,6],gnu:3,good:[4,6],gracefulli:0,grad:6,gradient:1,grid:6,group:3,grow:3,guarante:4,had:1,hand:[2,4,5],handl:2,happen:4,hard:0,has:[0,1,2,3,4,6],have:[0,4,5,6],help:6,here:[4,6],hessian:2,high:6,higher:5,histori:[0,3],history_for_slow:0,home:5,how:[0,1,3,5],howev:[4,6],htm:6,html:[5,6],http:[5,6],idea:4,imlug:6,imlug_genstatexpls_sect004:6,impact:[2,6],imposs:4,impract:0,improv:[0,2,3,4,6],inaccur:4,inc:0,includ:[0,1,6],increas:[0,6],increase_ndirs_initial_amt:0,increase_npt:0,increase_npt_amt:0,increase_num_extra_steps_with_restart:0,indic:6,industri:3,infinit:4,info:[1,6],inform:[3,4,6],infti:6,init:0,initi:[2,3,6],initialis:[2,6],input:[0,2,4,6],inspect:0,instal:[2,3],instanc:[0,4],instead:[4,5,6],integ:6,interest:[3,4],interfac:6,intermedi:0,intern:[0,4,6],interpol:[2,3,4,6],interpolation_change_j_norm:1,interpolation_condition_numb:1,interpolation_error:1,interpolation_total_residu:1,intersect:6,introduc:0,invers:2,involv:[0,4],iter:[0,2,3,6],iter_this_run:1,iter_typ:1,iters_tot:1,its:[0,6],j_k:[0,1],jacobian:[0,1,2,6],jan:[0,3,4,6],jun:3,june:[],just:2,keep:4,known:4,l015803:3,label:6,lambda:[0,1],larg:[0,4],larger:0,largest:0,last:[0,1],lastli:6,latest:5,ldot:4,least:4,least_squar:6,left:4,legend:6,len:[0,6],length:0,leq:[0,6],less:0,let:6,level:[0,5,6],librari:4,licens:3,like:[4,6],linalg:[0,6],linalgerror:0,lindon:[0,3,4,6],line:[0,6],linear:[0,1,4,6],link:2,linspac:6,list:[2,6],loc:6,local:[3,6],locat:5,log:[1,2,3,6],longer:2,lower:[0,1,6],lsqcurvefit:6,m_k:4,magnitud:[0,6],mai:[1,3,4,6],main:[4,6],maintain:4,major:4,make:[0,2,4,5,6],manag:3,mani:[0,2,3,4,6],manual:3,marteau:[0,3,4,6],math:6,mathbb:[3,4,6],mathemat:[0,3,4,6],mathrm:4,mathwork:6,matplotlib:6,matrix:[1,3,6],matrix_rank:0,max:6,max_distance_xk:1,max_fake_successful_step:0,max_it:0,max_npt:0,max_slow_it:0,max_unsuccessful_restart:0,maxfun:[0,6],maximum:[0,1,6],mean:[1,4],measur:[0,4],messag:6,method:[0,2,4],min:6,min_:[3,4,6],min_chgj_slop:0,min_correl:0,min_sing_v:0,minim:[2,4],minimum:[0,6],minor:2,model:[2,3,4,6],modifi:6,modul:6,momentum:0,momentum_extra_step:0,mont:4,more:[0,3,4],moreov:3,most:[0,1,4,6],move:[0,4,6],move_xk:0,msg:6,multipl:[2,3,4,6],multiplicative_noise_level:0,must:6,myfil:[1,6],n_to_print_whole_x_vector:0,nag:3,nan:[0,2],navig:5,ndirs_initi:0,necessari:2,need:[4,6],never:[3,6],next:6,nfev:6,nois:[3,4,6],noisi:[0,3,4],non:[3,4,6],nonconvex:4,none:[0,6],nonlinear:3,nonlinear_system:6,nonzero:0,nor:4,norm:[1,6],norm_gk:1,norm_sk:1,normal:6,note:[5,6],now:[0,6],npt:[0,1,6],nrestart:6,nrun:[1,6],nsampl:[1,6],num_extra_step:0,num_geom_step:0,num_new_dirns_each_it:0,number:[0,1,4,6],numer:3,numericalalgorithmsgroup:5,numpi:[0,2,5,6],obj:6,object:[1,2,3,4],objfun:[0,1,2,6],objfun_has_nois:[0,6],obs:4,observ:[4,6],occur:6,oct:3,octob:3,often:4,older:5,one:[0,1,4,5,6],onli:[0,1,4,6],onto:6,oper:2,oppos:0,opposit:0,opt:6,optim:[0,1,4,6],option:[0,2,3,4,5],order:6,org:5,origin:6,orthogon:0,other:0,otherwis:[0,4],our:[3,4,6],out:6,output:[1,2,3,4],outsid:6,over:[0,6],overdetermin:4,overflow:2,overflowerror:0,overlin:0,overrid:0,overridden:6,overview:[3,6],packag:[2,3,5],page:6,panda:[1,5,6],paper:[0,3,4,6],param1:6,param2:6,paramet:[2,3],part:[2,6],partial:6,particularli:[4,6],pass:[2,6],past:0,pball:6,pbox:6,per:[0,1,2,6],perform:[0,2,4,6],perturb:0,perturb_trust_region_step:0,phase:[0,6],physic:4,piec:[0,1],pip:3,pleas:3,plot:6,plt:6,point:[1,2,3,4,6],pois:1,poised:[0,1],popul:1,posisbl:0,possibl:[4,6],post:0,practic:6,precondit:0,predict:[1,4],prediction_error:6,prefer:4,preprint:[0,3,4,6],present:5,preserv:0,previou:[0,6],previous:1,print:[0,2,6],print_funct:6,print_progress:6,privat:5,privileg:5,probabl:4,problem:[0,1,2,3,4,6],proccess:4,process:[0,4],produc:[4,6],progress:[3,6],prohibit:4,project:6,provid:[4,5,6],pull:5,pure:5,purpos:6,put:4,pydata:5,pyplot:6,python:[0,5,6],quad:[3,4,6],quantit:4,quantiti:4,quit:0,quit_on_noise_level:0,r_1:[0,4,6],r_2:4,r_i:[0,4,6],r_m:[0,4,6],r_tol:0,radii:0,radiu:[0,1,2,4,6],random:[0,2,6],random_directions_make_orthogon:0,random_initial_direct:0,rang:6,rank:3,rate:0,rather:2,ratio:[0,1],reach:[0,6],recommend:0,recycl:0,reduc:[0,2],reduce_delta:0,reduct:1,refer:3,regim:2,region:[2,3,4,5,6],regress:3,regular:0,regularli:4,rel:0,rel_tol:0,relationship:4,relax:[3,6],releas:[2,3],remov:[2,5],repeat:4,replac:6,reproduc:[2,6],request:6,requir:[0,3,4,6],rerun:5,reset:0,reset_delta:0,reset_rho:0,resid:6,residu:[1,4,6],respect:0,restart:[1,2,3,6],result:[0,1,2,4,6],retriev:2,rho:[1,6],rho_:0,rho_k:[0,1],rhobeg:6,rhoend:6,rhoend_scal:0,right:[2,6],risk:0,robert:[0,3,4,6],robust:[0,3,4,6],root:5,rosenbrock:6,rosenbrock_noisi:6,rounding_error_const:0,roundoff:0,row:1,run:[0,5,6],run_in_parallel:0,runtim:2,runtimewarn:6,s_k:[0,1,4],safeti:[0,1,2],safety_step_thresh:0,sai:4,same:[0,3,4,6],sampl:6,sas:6,satisfi:[4,6],save:[0,1,6],save_diagnostic_info:[0,1],save_poised:[0,1],save_rk:[0,1],save_xk:[0,1],saw:1,scale:[0,2,6],scale_factor:0,scale_factor_for_quit:0,scaling_within_bound:6,scipi:[4,5,6],screen:0,script:6,search:[0,2,4],section:[0,1,2,6],see:[1,5,6],seed:[2,6],sens:4,sensibl:6,set:[1,2,3,4,6],set_xlabel:6,set_ylabel:6,setup:[0,5],sever:[0,4,6],shape:6,shift:[0,6],should:[0,5,6],show:6,side:[2,4],similar:4,similarli:4,simpl:[2,3,5],simpli:4,simul:4,sin:4,sinc:[0,1],singular:[0,6],site:5,situat:4,size:[4,6],slope:0,slow:[1,3,6],slow_it:1,small:[3,4,6],smaller:0,smallest:[0,1],smooth:0,soft:0,softwar:[0,3,4,5,6],soln:[1,6],solut:[2,3,4,5,6],solv:[0,1,3],solver:[0,2,3,4,6],some:[0,4,5,6],sourc:5,space:[0,2,4],specifi:[0,6],squar:[1,4],stai:4,stand:3,start:[2,4,6],statu:6,step:[0,1,2,4,6],still:0,stochast:[3,6],stop:[0,2],store:[0,2],str:6,string:6,strongli:4,structur:3,subproblem:[2,4,5,6],success:[0,1,6],successfulli:6,sudo:5,suffici:6,suitabl:4,sum:[1,4],sum_:[3,4,6],supervis:3,support:[3,5,6],suppos:[4,6],sure:[4,5],svd_max_jac_cond:0,svd_scale_factor:0,system:[0,1,3,5],t_i:6,tabl:6,take:[0,4,6],taken:6,task:4,tdata:6,technic:6,techniqu:4,termin:[1,3,6],test:[3,6],text:[0,1,3,4,6],than:[0,2,6],thei:[0,4],them:[2,4],therefor:6,thi:[0,1,2,3,4,5,6],thresh_for_slow:0,threshold:0,through:[0,2],throw_error_on_nan:0,time:[0,1,4,6],to_csv:1,toler:0,too:6,top:5,total:[1,6],tr_radiu:0,train:3,transact:[0,3,4,6],trend:0,trial:4,triangular:2,trigger:0,troubl:6,trust:[2,3,4,5,6],trustregion:[2,5],tupl:6,turn:1,two:[0,4],type:[1,4],typic:4,unabl:6,under:3,underdetermin:4,uninstal:3,unknown:4,unless:6,unpack:5,unsuccess:0,updat:2,upgrad:5,upper:[2,6],usag:[3,6],use:[0,1,2,3,5],use_full_rank_interp:0,use_old_rk:0,use_restart:0,use_soft_restart:0,used:[0,1,4,6],useful:[0,4,6],user:[0,4,5,6],user_param:[0,1,6],uses:0,using:[0,2,3,4,6],usual:0,val1:6,val2:6,valu:[1,2,3,4,6],vari:0,variabl:6,vast:4,vdot:4,vector:[0,1,4,6],veri:[0,4,6],version:[0,3,5],via:6,viewer:6,visibl:6,wai:[0,4,6],want:[5,6],warn:[2,6],well:4,what:[0,1],when:[0,2,3,6],where:[0,4,6],whether:[0,4,6],which:[0,1,2,3,4,5,6],whole:2,why:6,wish:[3,4,6],within:[0,2,6],without:[0,3],work:5,wors:6,write:6,written:5,www:5,x_0:[0,6],x_1:[4,6],x_2:[4,6],x_b:0,x_j:6,x_k:[0,4],x_n:4,xmin:6,xtol:6,y_1:4,y_i:[4,6],y_k:[0,1],y_m:4,y_t:0,ydata:6,yet:0,you:[3,4,5,6],your:[4,5,6],zenodo:2,zero:[0,2]},titles:["Advanced Usage","Diagnostic Information","Version History","DFO-LS: Derivative-Free Optimizer for Least-Squares Minimization","Overview","Installing DFO-LS","Using DFO-LS"],titleterms:{Adding:6,Using:6,acknowledg:3,advanc:0,algorithm:[0,1,4],apr:2,argument:6,bound:[],check:0,constraint:6,count:1,current:1,data:6,deriv:3,detail:4,dfo:[3,4,5,6],diagnost:1,dykstra:0,dynam:0,equat:[4,6],estim:6,evalu:6,exampl:6,feb:2,fit:[4,6],free:3,gener:0,grow:0,histori:2,how:6,inform:[0,1],initi:0,instal:5,interpol:[0,1],iter:1,jan:2,jun:2,least:[3,6],log:0,manag:0,manual:5,matrix:0,minim:[3,6],model:[0,1],more:6,multipl:0,nois:0,noisi:6,nonlinear:[4,6],object:[0,6],oct:2,optim:3,option:6,output:[0,6],overview:4,paramet:[0,4,6],pip:5,point:0,progress:[0,1],rank:0,refer:[0,4,6],region:[0,1],regress:0,requir:5,restart:0,set:0,simpl:6,slow:0,small:0,solv:[4,6],squar:[3,6],stochast:0,system:[4,6],termin:0,test:5,trust:[0,1],uninstal:5,usag:0,use:[4,6],using:5,valu:0,version:2,when:4}})
\ No newline at end of file
diff --git a/docs/build/html/userguide.html b/docs/build/html/userguide.html
index 5b40502..ab38f4f 100644
--- a/docs/build/html/userguide.html
+++ b/docs/build/html/userguide.html
@@ -8,7 +8,7 @@
- Using DFO-LS — DFO-LS v1.2.3 documentation
+ Using DFO-LS — DFO-LS v1.3.0 documentation
@@ -61,7 +61,7 @@
- 1.2.3
+ 1.3.0
@@ -94,7 +94,7 @@
How to use DFO-LS
Optional Arguments
A Simple Example
-Adding Bounds and More Output
+Adding Constraints and More Output
Example: Noisy Objective Evaluation
Example: Parameter Estimation/Data Fitting
Example: Solving a Nonlinear System of Equations
@@ -176,8 +176,8 @@ Nonlinear Least-Squares Minimization
-where the bound constraints \(a \leq x \leq b\) are optional. The upper and lower bounds on the variables are non-relaxable (i.e. DFO-LS will never ask to evaluate a point outside the bounds).
+\text{s.t.} &\quad x \in C\end{split}\]
+where the set \(C\) is an optional non-empty, closed and convex constraint set. The constraints are non-relaxable (i.e. DFO-LS will never ask to evaluate a point that is not feasible).
DFO-LS iteratively constructs an interpolation-based model for the objective, and determines a step using a trust-region framework.
For an in-depth technical description of the algorithm see the paper [CFMR2018].
@@ -239,6 +239,7 @@ Optional Arguments