Repository: POptUS/IBCDFO (language: python)

Issue:
New/max eig



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 327e5e47..30d30a8c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -29,8 +29,7 @@ jobs:
         - name: Install dependencies, ibcdfo
           run: |
             python -m pip install --upgrade pip
-            conda install octave
-            pip install oct2py ipdb
+            pip install ipdb
             pip install -r pounders/py/tests/test_requirements.txt
             pushd ibcdfo_pypkg
             rm -rf dist
@@ -39,7 +38,7 @@ jobs:
             popd
             git submodule update --init --recursive
             pushd ..
-            git clone https://github.com/POptUS/BenDFO.git
+            git clone -b py_calfun https://github.com/POptUS/BenDFO.git
             popd
 
         - name: Run tests, produce coverage reports
@@ -54,6 +53,9 @@ jobs:
             coverage run --source=ibcdfo --data-file=.cov_tomerge_regtest_2 test_failing_objective.py
             coverage run --source=ibcdfo --data-file=.cov_tomerge_regtest_3 benchmark_pounders4py.py
             popd
+            pushd manifold_sampling/py/tests/
+            coverage run --source=ibcdfo --data-file=.cov_tomerge_regtest_4 benchmark_manifold_sampling.py
+            popd
 
         - name: Merge coverage, run Coveralls
           env:
diff --git a/.gitignore b/.gitignore
index 0fc9778a..99ba9480 100644
--- a/.gitignore
+++ b/.gitignore
@@ -85,6 +85,9 @@ celerybeat.pid
 # SageMath parsed files
 *.sage.py
 
+# Matlab autosave files
+*.asv
+
 # Environments
 .env
 .venv
@@ -111,3 +114,5 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+*.egg-info
diff --git a/manifold_sampling/.gitignore b/manifold_sampling/.gitignore
new file mode 100644
index 00000000..e33609d2
--- /dev/null
+++ b/manifold_sampling/.gitignore
@@ -0,0 +1,1 @@
+*.png
diff --git a/manifold_sampling/m/manifold_sampling_primal.m b/manifold_sampling/m/manifold_sampling_primal.m
index a7d2735e..31c6a6b5 100644
--- a/manifold_sampling/m/manifold_sampling_primal.m
+++ b/manifold_sampling/m/manifold_sampling_primal.m
@@ -99,7 +99,7 @@
         G_k = Gres * D_k;
         beta = max(0, f_bar' - h(xkin));
 
-        % Line 6: Choose Hessions
+        % Line 6: Choose Hessians
         H_k = zeros(size(G_k, 2), n + 1, n + 1);
         for i = 1:size(G_k, 2) % would like to vectorize this tensor product ...
             for j = 1:size(Hres, 3) % p
diff --git a/manifold_sampling/py/build_p_models.py b/manifold_sampling/py/build_p_models.py
new file mode 100644
index 00000000..c61cddda
--- /dev/null
+++ b/manifold_sampling/py/build_p_models.py
@@ -0,0 +1,20 @@
+import numpy as np
+from evaluate_points_to_force_valid_model import evaluate_points_to_force_valid_model
+from ibcdfo.pounders import formquad
+from update_models import update_models
+
+
+def build_p_models(nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, tol, hfun, Ffun, Hash, L, U):
+    n = X.shape[1]
+    p = F.shape[1]
+
+    valid, Gres, Hres, X, F, h, nf, Hash = update_models(hfun, Ffun, n, p, nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, 1, Hash, tol, L, U)
+    # Evaluate model-improving points if necessary
+    if not valid and not len(Gres) == 0:
+        Mdir, np, valid, *_ = formquad(X[: nf + 1], F[: nf + 1], delta, xkin, fq_pars["npmax"], fq_pars["Par"], 1)
+        if valid:
+            raise Exception("what to do here")
+        X, F, h, nf, Hash = evaluate_points_to_force_valid_model(n, nf, xkin, delta, X, F, h, tol["gentype"], Mdir, np, hfun, Ffun, Hash, fq_pars, tol, nfmax, L, U)
+        __, Gres, Hres, X, F, h, nf, Hash = update_models(hfun, Ffun, n, p, nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, 1, Hash, tol, L, U)
+
+    return Gres, Hres, X, F, h, nf, Hash
diff --git a/manifold_sampling/py/call_user_scripts.py b/manifold_sampling/py/call_user_scripts.py
new file mode 100644
index 00000000..84d2a0bb
--- /dev/null
+++ b/manifold_sampling/py/call_user_scripts.py
@@ -0,0 +1,43 @@
+import numpy as np
+
+
+def call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, x_in, tol, L, U, allow_recalc=0):
+    if len(x_in) != len(U) or len(x_in) != len(L):
+        raise ValueError("Input vector dimensions do not match the bounds.")
+
+    if allow_recalc is None:
+        allow_recalc = 0
+
+    xnew = np.minimum(U, np.maximum(L, x_in))
+
+    for i in range(len(xnew)):
+        if U[i] - xnew[i] < np.finfo(float).eps * np.abs(U[i]) and U[i] > xnew[i]:
+            xnew[i] = U[i]
+            print("eps project!")
+        elif xnew[i] - L[i] < np.finfo(float).eps * np.abs(L[i]) and L[i] < xnew[i]:
+            xnew[i] = L[i]
+            print("eps project!")
+
+    xnew_hashable = tuple(xnew)
+
+    if not allow_recalc and xnew_hashable in [tuple(row) for row in X[:nf, :]]:
+        raise ValueError("Your method requested a point that has already been requested. FAIL!")
+
+    nf += 1
+    X[nf] = xnew
+    F[nf] = Ffun(X[nf])
+    h[nf], _, hashes_at_nf = hfun(F[nf])
+    Hash[nf] = hashes_at_nf
+    if np.any(np.isnan(F[nf, :])):
+        raise ValueError("Got a NaN. FAIL!")
+
+    if tol["hfun_test_mode"]:
+        assert isinstance(hashes_at_nf, list), "hfun must return a list of hashes"
+        # assert all(isinstance(hash_val, str) for hash_val in hashes_at_nf), "Hashes must be strings"
+
+        h_dummy1, grad_dummy1, hash_dummy = hfun(F[nf, :])
+        h_dummy2, grad_dummy2 = hfun(F[nf, :], hash_dummy)
+        assert np.any(np.abs(h_dummy1 - h_dummy2) <= 1e-16), "hfun values don't agree when being re-called with the same inputs"
+        assert np.all(grad_dummy1 == grad_dummy2), "hfun gradients don't agree when being re-called with the same inputs"
+
+    return nf, X, F, h, Hash, hashes_at_nf
diff --git a/manifold_sampling/py/check_inputs_and_initialize.py b/manifold_sampling/py/check_inputs_and_initialize.py
new file mode 100644
index 00000000..d7d62a8f
--- /dev/null
+++ b/manifold_sampling/py/check_inputs_and_initialize.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+
+def check_inputs_and_initialize(x0, F0, nfmax):
+    x0 = x0.squeeze()
+    n = int(len(x0))
+    p = int(len(F0))
+    delta = 0.1
+    printf = 1
+    h_activity_tol = min(1e-08, delta)
+
+    # Internal parameters/tolerances for formquad
+
+    fq_pars = {"Par": [np.sqrt(n), max(10, np.sqrt(n)), 0.001, 0.001], "npmax": (n + 1) * (n + 2) // 2}
+
+    # Internal parameters/tolerances for manifold sampling
+    tol = {
+        "maxdelta": 1,
+        "mindelta": 1e-13,
+        "gtol": 1e-13,
+        "norm_g_change": 0.001,
+        "kappa_d": 0.0001,
+        "eta1": 0.01,
+        "eta2": 10000.0,
+        "eta3": 0.5,
+        "gamma_dec": 0.5,
+        "gamma_inc": 2,
+        "hfun_test_mode": 1,
+        "gentype": 2,
+    }
+
+    # kappa_mh = 0;    # [dbl] > 0 that bounds the component model Hessians
+
+    assert nfmax >= n + 1, "nfmax is less than n+1, exiting"
+
+    X = np.vstack((x0, np.zeros((nfmax - 1, n))))
+    F = np.vstack((F0, np.zeros((nfmax - 1, p))))
+    h = np.zeros((nfmax, 1))
+    Hash = {}
+
+    nf = 0
+    trust_rho = 1
+    xkin = 0
+    Hres = np.zeros((n, n, p))
+
+    return n, delta, printf, fq_pars, tol, X, F, h, Hash, nf, trust_rho, xkin, Hres
diff --git a/manifold_sampling/py/choose_generator_set.py b/manifold_sampling/py/choose_generator_set.py
new file mode 100644
index 00000000..cfb6bc65
--- /dev/null
+++ b/manifold_sampling/py/choose_generator_set.py
@@ -0,0 +1,30 @@
+import numpy as np
+
+
+def choose_generator_set(X, Hash, gentype, xkin, nf, delta, F, hfun):
+    Act_Z_k = Hash[xkin]
+
+    if gentype == 2:
+        for i in [ind for ind in range(nf) if ind != xkin]:
+            if np.linalg.norm(X[xkin] - X[i]) <= delta * (1 + 1e-8):
+                Act_tmp = Hash[i]
+                Act_Z_k = np.concatenate((Act_Z_k, Act_tmp))
+
+    elif gentype == 3:
+        hxkin, _ = hfun(F[xkin, :], Act_Z_k)
+        for i in [ind for ind in range(nf) if ind != xkin]:
+            Act_tmp = Hash[i]
+            h_i, _ = hfun(F[xkin], Act_tmp)
+            if np.linalg.norm(X[xkin] - X[i], ord=np.inf) <= delta * (1 + 1e-8) and h_i[0] <= hxkin[0]:
+                Act_Z_k = np.concatenate((Act_Z_k, Act_tmp))
+            elif np.linalg.norm(X[xkin] - X[i], ord=np.inf) <= delta**2 * (1 + 1e-8) and h_i[0] > hxkin[0]:
+                Act_Z_k = np.concatenate((Act_Z_k, Act_tmp))
+
+    Act_Z_k = np.unique(Act_Z_k, axis=0)
+    f_k, D_k = hfun(F[xkin], Act_Z_k)
+    unique_indices = np.unique(D_k, axis=1, return_index=True)[1]
+    D_k = D_k[:, unique_indices]
+    Act_Z_k = Act_Z_k[unique_indices]
+    f_k = f_k[unique_indices]
+
+    return D_k, Act_Z_k, f_k
diff --git a/manifold_sampling/py/evaluate_points_to_force_valid_model.py b/manifold_sampling/py/evaluate_points_to_force_valid_model.py
new file mode 100644
index 00000000..9502e0b2
--- /dev/null
+++ b/manifold_sampling/py/evaluate_points_to_force_valid_model.py
@@ -0,0 +1,39 @@
+import numpy as np
+from call_user_scripts import call_user_scripts
+from ibcdfo.pounders import bmpts, formquad
+
+
+def evaluate_points_to_force_valid_model(n, nf, xkin, delta, X, F, h, gentype, Mdir, mp, hfun, Ffun, Hash, fq_pars, tol, nfmax, L, U):
+    # global nprob Qs zs bs p x0 nfmax h_activity_tol row_in_dfo_dat s inst
+    # A.n=n; A.nf=nf; A.xkin=xkin; A.delta=delta; A.X=X; A.F=F; A.h=h; A.gentype=gentype; A.Mdir=Mdir; A.np=np; A.hfun=hfun; A.Ffun=Ffun; A.Hash=Hash; A.fq_pars=fq_pars; A.tol=tol; A.nfmax=nfmax; A.L=L; A.U=U;
+
+    # Evaluate model-improving points to pick best one
+    # ! May eventually want to normalize Mdir first for infty norm
+    # Plus directions
+    # *** Dec 2016: THIS ASSUMES UNCONSTRAINED, proceed with caution
+    Mdir1, np1 = bmpts(X[xkin, :], Mdir[: n - mp + 1], L, U, delta, fq_pars["Par"][3])
+    # Res = zeros(n-np, 1);
+    for i in range(n - np1):
+        # if ~all(isinf(L)) || ~all(isinf(U))
+        #     D = Mdir1(i, :);
+        #     Res(i, 1) = D*(g_k+.5*H_k*D');
+        #     if Res(i, 1)> D*(-g_k+.5*H_k*D');
+        #         Mdir1(i, :)=-Mdir1(i, :); #neg dir predicted 2b better
+        #     end
+        # end
+        Xsp = Mdir1[i, :]
+        # Only do this evaluation if the point is new and nf < nfmax
+        if not np.any(np.all(X[xkin, :] + Xsp == X[: nf + 1], axis=1)) and nf + 1 < nfmax:
+            nf, X, F, h, Hash, _ = call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, X[xkin, :] + Xsp, tol, L, U, 1)
+
+    valid = formquad(X[: nf + 1], F[: nf + 1], delta, xkin, fq_pars["npmax"], fq_pars["Par"], 1)[2]
+    if not valid and nf + 1 < nfmax:
+        print(nf)
+        print(gentype)
+        print("Proceeding with nonvalid model! Report this to Stefan in Alg1")
+        # uuid = char(java.util.UUID.randomUUID);
+        # global mw_prob_num hfun
+        # save(['first_failure_for_row_in_dfo_dat=' int2str(mw_prob_num) '_hfun=' func2str(hfun{1}) ], 'A');
+        # error('a')
+
+    return X, F, h, nf, Hash
diff --git a/manifold_sampling/py/h_examples/piecewise_quadratic.py b/manifold_sampling/py/h_examples/piecewise_quadratic.py
new file mode 100644
index 00000000..ca93e567
--- /dev/null
+++ b/manifold_sampling/py/h_examples/piecewise_quadratic.py
@@ -0,0 +1,55 @@
+import numpy as np
+
+
+def piecewise_quadratic(z, H0=None, **kwargs):
+    # Evaluates the piecewise quadratic function
+    #   max_j { || z - z_j ||_{Q_j}^2 + b_j }
+
+    # Inputs:
+    #  z:              [1 x p]   point where we are evaluating h
+    #  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate
+
+    # Outputs:
+    #  h: [dbl]                       function value
+    #  grads: [p x l]                 gradients of each of the l quadratics active at z
+    #  Hash: [1 x l cell of strings]  set of hashes for each of the l quadratics active at z (in the same order as the elements of grads)
+
+    # Hashes are output (and must be input) in the following fashion:
+    #   Hash{i} = 'j' if quadratic j is active at z (or H0{i} = 'j' if the
+    #   value/gradient of quadratic j at z is desired)
+
+    Qs = kwargs["Qs"]
+    zs = kwargs["zs"]
+    cs = kwargs["cs"]
+
+    if H0 is None:
+        n, J = zs.shape
+        manifolds = np.zeros(J)
+        for j in range(J):
+            manifolds[j] = np.dot(np.dot((z - zs[:, j]), Qs[:, :, j]), (z - zs[:, j])) + cs[j]
+
+        h = np.max(manifolds)
+
+        atol = 1e-8
+        rtol = 1e-8
+        inds = np.where(np.abs(h - manifolds) <= atol + rtol * np.abs(manifolds))[0]
+
+        grads = np.zeros((n, len(inds)))
+
+        Hash = [str(ind) for ind in inds]
+        for j in range(len(inds)):
+            grads[:, j] = 2 * np.dot(Qs[:, :, inds[j]], (z - zs[:, inds[j]]))
+
+        return h, grads, Hash
+
+    else:
+        J = len(H0)
+        h = np.zeros(J)
+        grads = np.zeros((len(z), J))
+
+        for k in range(J):
+            j = int(H0[k])
+            h[k] = np.dot(np.dot((z - zs[:, j]), Qs[:, :, j]), (z - zs[:, j])) + cs[j]
+            grads[:, k] = 2 * np.dot(Qs[:, :, j], (z - zs[:, j]))
+
+        return h, grads
diff --git a/manifold_sampling/py/h_examples/pw_maximum.py b/manifold_sampling/py/h_examples/pw_maximum.py
new file mode 100644
index 00000000..54a56e6d
--- /dev/null
+++ b/manifold_sampling/py/h_examples/pw_maximum.py
@@ -0,0 +1,44 @@
+import numpy as np
+
+
+def pw_maximum(z, H0=None):
+    # Evaluates the pointwise maximum function
+    #   max_j { z_j }
+
+    # Inputs:
+    #  z:              [1 x p]   point where we are evaluating h
+    #  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate z
+
+    # Outputs:
+    #  h: [dbl]                       function value
+    #  grads: [p x l]                 gradients of each of the l manifolds active at z
+    #  Hash: [1 x l cell of strings]  set of hashes for each of the l manifolds active at z (in the same order as the elements of grads)
+
+    n = len(z)
+
+    if H0 is None:
+        h = np.max(z)
+
+        atol = 1e-8
+        rtol = 1e-8
+        inds = np.where(np.abs(h - z) <= atol + rtol * np.abs(z))[0]
+
+        grads = np.zeros((n, len(inds)))
+        Hash = [str(ind) for ind in inds]
+
+        for j in range(len(inds)):
+            grads[inds[j], j] = 1
+
+        return h, grads, Hash
+
+    else:
+        J = len(H0)
+        h = np.zeros(J)
+        grads = np.zeros((len(z), J))
+
+        for k in range(J):
+            j = int(H0[k])
+            h[k] = z[j]
+            grads[j, k] = 1
+
+        return h, grads
diff --git a/manifold_sampling/py/h_examples/pw_maximum_eig.py b/manifold_sampling/py/h_examples/pw_maximum_eig.py
new file mode 100644
index 00000000..6fe42715
--- /dev/null
+++ b/manifold_sampling/py/h_examples/pw_maximum_eig.py
@@ -0,0 +1,56 @@
+import numpy as np
+import numpy.linalg as LA
+from scipy import spatial
+
+
+def pw_maximum_eig(z, H0=None):
+    n2 = len(z)
+    n = np.int(np.sqrt(n2))
+
+    # reshape z into a matrix
+    M = np.reshape(z, (n, n))
+
+    # compute eigendecomposition
+    [eigvals, eigvecs] = LA.eig(M)
+    sorted_indices = np.argsort(eigvals)
+    eigvals = eigvals[sorted_indices]
+    eigvecs = eigvecs[sorted_indices]
+
+    if H0 is None:
+        h = np.max(eigvals)
+
+        atol = 1e-8
+        rtol = 1e-8
+        inds = np.where(np.abs(h - eigvals) <= atol + rtol * np.abs(eigvals))[0]
+
+        selection_grads = np.zeros((n, len(inds)))
+        grads = np.zeros((n2, len(inds)))
+        Hash = [[] for i in range(len(inds))]
+        for j in range(len(inds)):
+            selection_grads[inds[j], j] = 1
+            grad = eigvecs @ np.diag(selection_grads[:, j]) @ eigvecs.T
+            grads[:, j] = np.reshape(grad, (1, n2))
+            # rounding Hash to 1 digit on purpose (this will make sense, see else statement)
+            Hash[j] = [str(np.round(val, 1)) for val in eigvecs[inds[j]]]
+
+        return h, grads, Hash
+
+    else:
+        len_inds = len(H0)
+        h = np.zeros(len_inds)
+        grads = np.zeros((n2, len_inds))
+        selection_grads = np.zeros((n, len_inds))
+
+        tree = spatial.KDTree(eigvecs)
+
+        H0 = np.array(H0, dtype=float)
+
+        for k in range(len_inds):
+            _, ind = tree.query(H0[k])
+            h[k] = eigvals[ind]
+
+            selection_grads[ind, k] = 1
+            grad = eigvecs @ np.diag(selection_grads[:, k]) @ eigvecs.T
+            grads[:, k] = np.reshape(grad, (1, n2))
+
+        return h, grads
diff --git a/manifold_sampling/py/h_examples/pw_maximum_squared.py b/manifold_sampling/py/h_examples/pw_maximum_squared.py
new file mode 100644
index 00000000..603af0dc
--- /dev/null
+++ b/manifold_sampling/py/h_examples/pw_maximum_squared.py
@@ -0,0 +1,45 @@
+import numpy as np
+
+
+def pw_maximum_squared(z, H0=None):
+    # Evaluates the pointwise maximum function
+    #   max_j { z_j^2 }
+
+    # Inputs:
+    #  z:              [1 x p]   point where we are evaluating h
+    #  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate z
+
+    # Outputs:
+    #  h: [dbl]                       function value
+    #  grads: [p x l]                 gradients of each of the l manifolds active at z
+    #  Hash: [1 x l cell of strings]  set of hashes for each of the l manifolds active at z (in the same order as the elements of grads)
+
+    n = len(z)
+
+    if H0 is None:
+        z2 = z**2
+        h = np.max(z2)
+
+        atol = 1e-8
+        rtol = 1e-8
+        inds = np.where(np.abs(h - z2) <= atol + rtol * np.abs(z2))[0]
+
+        grads = np.zeros((n, len(inds)))
+        Hash = [str(ind) for ind in inds]
+
+        for j in range(len(inds)):
+            grads[inds[j], j] = 2 * z[inds[j]]
+
+        return h, grads, Hash
+
+    else:
+        J = len(H0)
+        h = np.zeros(J)
+        grads = np.zeros((len(z), J))
+
+        for k in range(J):
+            j = int(H0[k])
+            h[k] = z[j] ** 2
+            grads[j, k] = 2 * z[j]
+
+        return h, grads
diff --git a/manifold_sampling/py/h_examples/pw_minimum.py b/manifold_sampling/py/h_examples/pw_minimum.py
new file mode 100644
index 00000000..5f9737a5
--- /dev/null
+++ b/manifold_sampling/py/h_examples/pw_minimum.py
@@ -0,0 +1,44 @@
+import numpy as np
+
+
+def pw_minimum(z, H0=None):
+    # Evaluates the pointwise minimum function
+    #   min_j { z_j }
+
+    # Inputs:
+    #  z:              [1 x p]   point where we are evaluating h
+    #  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate z
+
+    # Outputs:
+    #  h: [dbl]                       function value
+    #  grads: [p x l]                 gradients of each of the l manifolds active at z
+    #  Hash: [1 x l cell of strings]  set of hashes for each of the l manifolds active at z (in the same order as the elements of grads)
+
+    n = len(z)
+
+    if H0 is None:
+        h = np.min(z)
+
+        atol = 1e-8
+        rtol = 1e-8
+        inds = np.where(np.abs(h - z) <= atol + rtol * np.abs(z))[0]
+
+        grads = np.zeros((n, len(inds)))
+        Hash = [str(ind) for ind in inds]
+
+        for j in range(len(inds)):
+            grads[inds[j], j] = 1
+
+        return h, grads, Hash
+
+    else:
+        J = len(H0)
+        h = np.zeros(J)
+        grads = np.zeros((len(z), J))
+
+        for k in range(J):
+            j = int(H0[k])
+            h[k] = z[j]
+            grads[j, k] = 1
+
+        return h, grads
diff --git a/manifold_sampling/py/h_examples/pw_minimum_squared.py b/manifold_sampling/py/h_examples/pw_minimum_squared.py
new file mode 100644
index 00000000..1125ea34
--- /dev/null
+++ b/manifold_sampling/py/h_examples/pw_minimum_squared.py
@@ -0,0 +1,45 @@
+import numpy as np
+
+
+def pw_minimum_squared(z, H0=None):
+    # Evaluates the pointwise minimum function
+    #   min_j { z_j^2 }
+
+    # Inputs:
+    #  z:              [1 x p]   point where we are evaluating h
+    #  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate z
+
+    # Outputs:
+    #  h: [dbl]                       function value
+    #  grads: [p x l]                 gradients of each of the l manifolds active at z
+    #  Hash: [1 x l cell of strings]  set of hashes for each of the l manifolds active at z (in the same order as the elements of grads)
+
+    n = len(z)
+
+    if H0 is None:
+        z2 = z**2
+        h = np.min(z2)
+
+        atol = 1e-8
+        rtol = 1e-8
+        inds = np.where(np.abs(h - z2) <= atol + rtol * np.abs(z2))[0]
+
+        grads = np.zeros((n, len(inds)))
+
+        Hash = [str(ind) for ind in inds]
+        for j in range(len(inds)):
+            grads[inds[j], j] = 2 * z[inds[j]]
+
+        return h, grads, Hash
+
+    else:
+        J = len(H0)
+        h = np.zeros(J)
+        grads = np.zeros((n, J))
+
+        for k in range(J):
+            j = int(H0[k])
+            h[k] = z[j] ** 2
+            grads[j, k] = 2 * z[j]
+
+        return h, grads
diff --git a/manifold_sampling/py/h_examples/quantile.py b/manifold_sampling/py/h_examples/quantile.py
new file mode 100644
index 00000000..c48aa215
--- /dev/null
+++ b/manifold_sampling/py/h_examples/quantile.py
@@ -0,0 +1,48 @@
+import numpy as np
+
+
+def quantile(z, H0=None):
+    # Evaluates the q^th quantile of the values
+    #  { z_j^2 }
+
+    # Inputs:
+    #  z:              [1 x p]   point where we are evaluating h
+    #  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate
+
+    # Outputs:
+    #  h: [dbl]                       function value
+    #  grads: [p x l]                 gradients of each of the l quadratics active at z
+    #  Hash: [1 x l cell of strings]  set of hashes for each of the l quadratics active at z (in the same order as the elements of grads)
+
+    q = 1
+
+    n = len(z)
+    z2 = z**2
+    sorted_inds = np.argsort(z2)
+    z2_sort = z2[sorted_inds]
+    z_sort = z[sorted_inds]
+
+    if H0 is None:
+        h = z2_sort[q]
+        atol = 1e-08
+        rtol = 1e-08
+        inds = np.where(np.abs(h - z2_sort) <= atol + rtol * np.abs(z2_sort))[0]
+
+        grads = np.zeros((n, len(inds)))
+        Hash = [str(ind) for ind in inds]
+
+        for j in range(len(inds)):
+            grads[inds[j], j] = 2 * z_sort[inds[j]]
+
+        return h, grads, Hash
+    else:
+        J = len(H0)
+        h = np.zeros(J)
+        grads = np.zeros((len(z), J))
+
+        for k in range(J):
+            j = int(H0[k])
+            h[k] = z2_sort[j]
+            grads[j, k] = 2 * z_sort[j]
+
+        return h, grads
diff --git a/manifold_sampling/py/manifold_sampling_primal.py b/manifold_sampling/py/manifold_sampling_primal.py
new file mode 100644
index 00000000..c7be913a
--- /dev/null
+++ b/manifold_sampling/py/manifold_sampling_primal.py
@@ -0,0 +1,170 @@
+# This code solves the problem
+#     minimize h(F(x))
+# where x is an [n by 1] vector, F is a blackbox function mapping from R^n to
+# R^p, and h is a nonsmooth function mapping from R^p to R.
+#
+#
+# Inputs:
+#  hfun:    [func]   Given point z, returns
+#                      - [scalar] the value h(z)
+#                      - [p x l] gradients for all l limiting gradients at z
+#                      - [1 x l list of hashes] for each manifold active at z
+#                    Given point z and l hashes H, returns
+#                      - [1 x l] the value h_i(z) for each hash in H
+#                      - [p x l] gradients of h_i(z) for each hash in H
+#  Ffun:    [func]    Evaluates F, the black box simulation, returning a [1 x p] vector.
+#  x0:      [1 x n]   Starting point
+#  nfmax:   [int]     Maximum number of function evaluations
+#
+# Outputs:
+#   X:      [nfmax x n]   Points evaluated
+#   F:      [nfmax x p]   Their simulation values
+#   h:      [nfmax x 1]   The values h(F(x))
+#   xkin:   [int]         Current trust region center
+#   flag:   [int]         Inform user why we stopped.
+#                           -1 if error
+#                            0 if nfmax function evaluations were performed
+#                            final model gradient norm otherwise
+#
+# Some other values
+#  n:       [int]     Dimension of the domain of F (deduced from x0)
+#  p:       [int]     Dimension of the domain of h (deduced from evaluating F(x0))
+#  delta:   [dbl]     Positive starting trust region radius
+# Intermediate Variables:
+#   nf    [int]         Counter for the number of function evaluations
+#   s_k   [dbl]         Step from current iterate to approx. TRSP solution
+#   norm_g [dbl]        Stationary measure ||g||
+#   Gres [n x p]        Model gradients for each of the p outputs from Ffun
+#   Hres [n x n x p]    Model Hessians for each of the p outputs from Ffun
+#   Hash [cell]         Contains the hashes active at each evaluated point in X
+#   Act_Z_k [l cell]    List of hashes for active selection functions in TR
+#   G_k  [n x l]        Matrix of model gradients composed with gradients of elements in Act_Z_k
+#   D_k  [p x l_2]      Matrix of gradients of selection functions at different points in p-space
+
+import numpy as np
+from build_p_models import build_p_models
+from call_user_scripts import call_user_scripts
+from check_inputs_and_initialize import check_inputs_and_initialize
+from choose_generator_set import choose_generator_set
+from ibcdfo.pounders import checkinputss
+from minimize_affine_envelope import minimize_affine_envelope
+
+
+def manifold_sampling_primal(hfun, Ffun, x0, L, U, nfmax, subprob_switch):
+    # Deduce p from evaluating Ffun at x0
+    try:
+        F0 = Ffun(x0)
+    finally:
+        pass
+
+    n, delta, printf, fq_pars, tol, X, F, h, Hash, nf, successful, xkin, Hres = check_inputs_and_initialize(x0, F0, nfmax)
+    flag, x0, __, F0, L, U, xkin = checkinputss(hfun, x0, n, fq_pars["npmax"], nfmax, tol["gtol"], delta, 1, len(F0), F0, xkin, L, U)
+    if flag == -1:
+        X = x0
+        F = F0
+        h = []
+        return X, F, h, xkin, flag
+
+    # Evaluate user scripts at x_0
+    h[nf], __, hashes_at_nf = hfun(F[nf])
+    Hash[nf] = hashes_at_nf
+
+    H_mm = np.zeros((n, n))
+
+    while nf + 1 < nfmax and delta > tol["mindelta"]:
+        bar_delta = delta
+
+        # Line 3: manifold sampling while loop
+        while nf + 1 < nfmax:
+            # Line 4: build models
+            Gres, Hres, X, F, h, nf, Hash = build_p_models(nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, tol, hfun, Ffun, Hash, L, U)
+            if len(Gres) == 0:
+                print(np.array(["Model building failed. Empty Gres. Delta = " + str(delta)]))
+                X = X[: nf + 1]
+                F = F[: nf + 1]
+                h = h[: nf + 1]
+                flag = -1
+                return X, F, h, xkin, flag
+            if nf + 1 >= nfmax:
+                flag = 0
+                return X, F, h, xkin, flag
+
+            # Line 5: Build set of activities Act_Z_k, gradients D_k, G_k, and beta
+            D_k, Act_Z_k, f_bar = choose_generator_set(X, Hash, 3, xkin, nf, delta, F, hfun)
+            G_k = Gres @ D_k
+            beta = np.maximum(0, f_bar - h[xkin])
+
+            # Line 6: Choose Hessians
+            H_k = np.zeros((G_k.shape[1], n + 1, n + 1))
+            for i in range(G_k.shape[1]):
+                for j in range(Hres.shape[2]):
+                    H_k[i, 1:, 1:] = np.squeeze(H_k[i, 1:, 1:]) + D_k[j, i] * Hres[:, :, j]
+
+            # Line 7: Find a candidate s_k by solving QP
+            Low = np.maximum(L - X[xkin], -delta)
+            Upp = np.minimum(U - X[xkin], delta)
+            s_k, tau_k, __, lambda_k = minimize_affine_envelope(h[xkin], f_bar, beta, G_k, H_mm, delta, Low, Upp, H_k, subprob_switch)
+
+            # Line 8: Compute stationary measure chi_k
+            Low = np.maximum(L - X[xkin], -1.0)
+            Upp = np.minimum(U - X[xkin], 1.0)
+            __, __, chi_k, __ = minimize_affine_envelope(h[xkin], f_bar, beta, G_k, np.zeros((n, n)), delta, Low, Upp, np.zeros((G_k.shape[1], n + 1, n + 1)), subprob_switch)
+
+            # Lines 9-11: Convergence test: tiny master model gradient and tiny delta
+            if chi_k <= tol["gtol"] and delta <= tol["mindelta"]:
+                print("Convergence satisfied: small stationary measure and small delta")
+                X = X[: nf + 1]
+                F = F[: nf + 1]
+                h = h[: nf + 1]
+                flag = chi_k
+                return X, F, h, xkin, flag
+
+            # Line 12: Evaluate F
+            nf, X, F, h, Hash, hashes_at_nf = call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, X[xkin] + np.transpose(s_k), tol, L, U, 1)
+            # Line 13: Compute rho_k
+            ared = h[xkin] - h[nf]
+            pred = -tau_k
+            if pred == 0:
+                rho_k = -np.inf
+            else:
+                rho_k = ared / pred
+
+            # Lines 14-16: Check for success
+            if rho_k >= tol["eta1"] and pred > 0:
+                successful = True
+                break
+            else:
+                # Line 18: Check temporary activities after adding TRSP solution to X
+                __, tmp_Act_Z_k, __ = choose_generator_set(X, Hash, 3, xkin, nf, delta, F, hfun)
+
+                # Lines 19: See if any new activities
+                if np.all(np.isin(tmp_Act_Z_k, Act_Z_k)):
+                    # Line 20: See if intersection is nonempty
+                    if np.any(np.isin(hashes_at_nf, Act_Z_k)):
+                        successful = False
+                        break
+                    else:
+                        # Line 24: Shrink delta
+                        delta = tol["gamma_dec"] * delta
+
+        if successful:
+            xkin = nf
+            if rho_k > tol["eta3"] and np.linalg.norm(s_k, ord=np.inf) > 0.8 * bar_delta:
+                # Update delta if rho is sufficiently large
+                delta = bar_delta * tol["gamma_inc"]
+                # h_activity_tol = min(1e-8, delta);
+        else:
+            # Line 21: iteration is unsuccessful; shrink Delta
+            delta = max(bar_delta * tol["gamma_dec"], tol["mindelta"])
+            # h_activity_tol = min(1e-8, delta);
+        print("nf: %8d; fval: %8e; chi: %8e; radius: %8e; \n" % (nf, h[xkin], chi_k, delta))
+
+    if nf + 1 >= nfmax:
+        flag = 0
+    else:
+        X = X[: nf + 1]
+        F = F[: nf + 1]
+        h = h[: nf + 1]
+        flag = chi_k
+
+    return X, F, h, xkin, flag
diff --git a/manifold_sampling/py/minimize_affine_envelope.py b/manifold_sampling/py/minimize_affine_envelope.py
new file mode 100644
index 00000000..5b7a3cea
--- /dev/null
+++ b/manifold_sampling/py/minimize_affine_envelope.py
@@ -0,0 +1,60 @@
+import numpy as np
+from scipy.optimize import linprog
+
+
+def minimize_affine_envelope(f, f_bar, beta, G_k, H, delta, Low, Upp, H_k, subprob_switch):
+    G_k_smaller, cols = np.unique(G_k, axis=1, return_index=True)
+
+    n, p = G_k_smaller.shape
+
+    bk = -(f_bar - f - beta)
+    bk_smaller = bk[cols]
+    H_k_smaller = H_k[cols, :, :]
+
+    A = np.hstack((-np.ones((p, 1)), G_k_smaller.T))
+    ff = np.concatenate((np.array([[1]]), np.zeros((n, 1))))
+    HH = np.block([[0, np.zeros((1, n))], [np.zeros((n, 1)), H]])
+    x0 = np.vstack((np.array([np.max(-bk_smaller)]), np.zeros((n, 1))))
+
+    if subprob_switch == "GAMS_QCP":
+        # Implement solve_matts_QCP function here
+        pass
+    elif subprob_switch == "GAMS_LP":
+        # Implement solve_matts_LP function here
+        pass
+    elif subprob_switch == "linprog":
+        options = {"disp": False}
+        try:
+            res = linprog(c=ff.flatten(), A_ub=A, b_ub=bk_smaller, bounds=list(zip([None] + list(Low), [None] + list(Upp))), options=options, x0=x0)
+            x = res.x
+            duals_g = -1.0 * res.ineqlin.marginals
+            duals_u = -1.0 * res.upper.marginals[1:]
+            duals_l = res.lower.marginals[1:]
+        except:
+            normA = np.linalg.norm(A[:, 1:])
+            rescaledA = np.zeros_like(A)
+            rescaledA[:, 0] = -np.ones(p)
+            rescaledA[:, 1:] = A[:, 1:] / normA
+            res = linprog(c=ff.flatten(), A_ub=rescaledA, b_ub=bk_smaller, bounds=list(zip([None] + list(Low), [None] + list(Upp))), options=options, x0=x0)
+            x = res.x
+            duals_g = -1.0 * res.ineqlin.marginals
+            duals_u = -1.0 * res.upper.marginals[1:] * normA
+            duals_l = res.lower.marginals[1:] * normA
+    else:
+        raise ValueError("Unrecognized subprob_switch")
+
+    lambda_star = np.zeros(G_k.shape[1])
+    lambda_star[cols] = duals_g
+
+    s = x[1:]
+    tau = max(-bk + np.dot(G_k.T, s)) + 0.5 * np.dot(s, np.dot(H, s))
+    if tau > 0:
+        tau = 0
+        s = np.zeros(n)
+
+    Low[duals_l <= 0] = 0
+    Upp[duals_u <= 0] = 0
+
+    chi = np.linalg.norm(np.dot(G_k, lambda_star) - duals_l + duals_u) + np.dot(bk, lambda_star) - np.dot(Low, duals_l) + np.dot(Upp, duals_u)
+
+    return s, tau, chi, lambda_star
diff --git a/manifold_sampling/py/update_models.py b/manifold_sampling/py/update_models.py
new file mode 100644
index 00000000..da91cd5e
--- /dev/null
+++ b/manifold_sampling/py/update_models.py
@@ -0,0 +1,36 @@
+import numpy as np
+from call_user_scripts import call_user_scripts
+from ibcdfo.pounders import bmpts, formquad
+
+
+def update_models(hfun, Ffun, n, p, nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, geometry_pt_flag, Hash, tol, L, U):
+    Cres = F[xkin, :]
+    Res = np.zeros(F.shape)  # Stores the residuals for model updates
+
+    for i in range(nf):
+        D = X[i] - X[xkin]
+        for j in range(len(Cres)):
+            Res[i, j] = (F[i, j] - Cres[j]) - 0.5 * np.dot(D, np.dot(Hres[:, :, j], D.T))
+
+    Mdir, mp, valid, Gres, Hresdel, _ = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xkin, fq_pars["npmax"], fq_pars["Par"], 0)
+
+    # Evaluate geometry points
+    if mp < n and geometry_pt_flag:  # Must obtain and evaluate bounded geometry points
+        Mdir, mp = bmpts(X[xkin], Mdir[: n - mp], L, U, delta, fq_pars["Par"][2])
+
+        for i in range(min(n - mp, nfmax - nf)):
+            nf, X, F, h, Hash, _ = call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, X[xkin, :] + Mdir[i, :], tol, L, U, 1)
+            D = Mdir[i, :]
+
+            for j in range(p):
+                Res[nf, j] = (F[nf, j] - Cres[j]) - 0.5 * np.dot(D, np.dot(Hres[:, :, j], D.T))
+
+        _, _, valid, Gres, Hresdel, _ = formquad(X[: nf + 1], Res[: nf + 1], delta, xkin, fq_pars["npmax"], fq_pars["Par"], 0)
+
+        if len(Gres) == 0:
+            return valid, Gres, Hres, X, F, h, nf, Hash
+
+    if len(Gres):  # We'll be doing evaluations; Hres will be updated after that
+        Hres = Hres + Hresdel
+
+    return valid, Gres, Hres, X, F, h, nf, Hash
diff --git a/minq b/minq
index 1da393e7..b0d325e9 160000
--- a/minq
+++ b/minq
@@ -1,1 +1,1 @@
-Subproject commit 1da393e71c105c7e14fbc0ad0893a5ca3968bc67
+Subproject commit b0d325e97e51dc85c0c2ee76c6ea5f4b7cebf1b7


Test patch:
diff --git a/manifold_sampling/py/tests/allcomb.py b/manifold_sampling/py/tests/allcomb.py
new file mode 100644
index 00000000..72cf306e
--- /dev/null
+++ b/manifold_sampling/py/tests/allcomb.py
@@ -0,0 +1,111 @@
+import warnings
+
+import numpy as np
+
+
+def allcomb(varargin):
+    # ALLCOMB - All combinations
+    #    B = ALLCOMB(A1,A2,A3,...,AN) returns all combinations of the elements
+    #    in the arrays A1, A2, ..., and AN. B is P-by-N matrix where P is the product
+    #    of the number of elements of the N inputs.
+    #    This functionality is also known as the Cartesian Product. The
+    #    arguments can be numerical and/or characters, or they can be cell arrays.
+
+    #    Examples:
+    #       allcomb([1 3 5],[-3 8],[0 1]) # numerical input:
+    #       # -> [ 1  -3   0
+    #       #      1  -3   1
+    #       #      1   8   0
+    #       #        ...
+    #       #      5  -3   1
+    #       #      5   8   1 ] ; # a 12-by-3 array
+
+    #       allcomb('abc','XY') # character arrays
+    #       # -> [ aX ; aY ; bX ; bY ; cX ; cY] # a 6-by-2 character array
+
+    #       allcomb('xy',[65 66]) # a combination -> character output
+    #       # -> ['xA' ; 'xB' ; 'yA' ; 'yB'] # a 4-by-2 character array
+
+    #       allcomb({'hello','Bye'},{'Joe', 10:12},{99999 []}) # all cell arrays
+    #       # -> {  'hello'  'Joe'        [99999]
+    #       #       'hello'  'Joe'             []
+    #       #       'hello'  [1x3 double] [99999]
+    #       #       'hello'  [1x3 double]      []
+    #       #       'Bye'    'Joe'        [99999]
+    #       #       'Bye'    'Joe'             []
+    #       #       'Bye'    [1x3 double] [99999]
+    #       #       'Bye'    [1x3 double]      [] } ; # a 8-by-3 cell array
+
+    #    ALLCOMB(..., 'matlab') causes the first column to change fastest which
+    #    is consistent with matlab indexing. Example:
+    #      allcomb(1:2,3:4,5:6,'matlab')
+    #      # -> [ 1 3 5 ; 1 4 5 ; 1 3 6 ; ... ; 2 4 6 ]
+
+    #    If one of the N arguments is empty, ALLCOMB returns a 0-by-N empty array.
+
+    #    See also NCHOOSEK, PERMS, NDGRID
+    #         and NCHOOSE, COMBN, KTHCOMBN (Matlab Central FEX)
+
+    # Tested in Matlab R2015a and up
+    # version 4.2 (apr 2018)
+    # (c) Jos van der Geest
+    # email: samelinoa@gmail.com
+
+    # History
+    # 1.1 (feb 2006), removed minor bug when entering empty cell arrays;
+    #     added option to let the first input run fastest (suggestion by JD)
+    # 1.2 (jan 2010), using ii as an index on the left-hand for the multiple
+    #     output by NDGRID. Thanks to Jan Simon, for showing this little trick
+    # 2.0 (dec 2010). Bruno Luong convinced me that an empty input should
+    # return an empty output.
+    # 2.1 (feb 2011). A cell as input argument caused the check on the last
+    #      argument (specifying the order) to crash.
+    # 2.2 (jan 2012). removed a superfluous line of code (ischar(..))
+    # 3.0 (may 2012) removed check for doubles so character arrays are accepted
+    # 4.0 (feb 2014) added support for cell arrays
+    # 4.1 (feb 2016) fixed error for cell array input with last argument being
+    #     'matlab'. Thanks to Richard for pointing this out.
+    # 4.2 (apr 2018) fixed some grammar mistakes in the help and comments
+
+    narginchk(1, Inf)
+    NC = len(varargin)
+    # check if we should flip the order
+    if ischar(varargin[end()]) and (strcmpi(varargin[end()], "matlab") or strcmpi(varargin[end()], "john")):
+        # based on a suggestion by JD on the FEX
+        NC = NC - 1
+        ii = np.arange(1, NC + 1)
+    else:
+        # default: enter arguments backwards, so last one (AN) is changing fastest
+        ii = np.arange(NC, 1 + -1, -1)
+
+    args = varargin(np.arange(1, NC + 1))
+    if np.any(cellfun("isempty", args)):
+        warnings.warn("ALLCOMB:EmptyInput", "One of more empty inputs result in an empty output.")
+        A = np.zeros((0, NC))
+    else:
+        if NC == 0:
+            A = np.zeros((0, 0))
+        else:
+            if NC == 1:
+                A = args[0]
+            else:
+                isCellInput = cellfun(iscell, args)
+                if np.any(isCellInput):
+                    if not np.all(isCellInput):
+                        raise Exception("ALLCOMB:InvalidCellInput", "For cell input, all arguments should be cell arrays.")
+                    # for cell input, we use to indices to get all combinations
+                    ix = cellfun(lambda c: np.arange(1, np.asarray(c).size + 1), args, "un", 0)
+                    # flip using ii if last column is changing fastest
+                    ix[ii] = ndgrid(ix[ii])
+                    A = cell(np.asarray(ix[0]).size, NC)
+                    for k in np.arange(1, NC + 1).reshape(-1):
+                        # combine
+                        A[:, k] = reshape(args[k](ix[k]), [], 1)
+                else:
+                    # non-cell input, assuming all numerical values or strings
+                    # flip using ii if last column is changing fastest
+                    A[ii] = ndgrid(args[ii])
+                    # concatenate
+                    A = reshape(cat(NC + 1, A[:]), [], NC)
+
+    return A
diff --git a/manifold_sampling/py/tests/benchmark_manifold_sampling.py b/manifold_sampling/py/tests/benchmark_manifold_sampling.py
new file mode 100644
index 00000000..6a09c1f9
--- /dev/null
+++ b/manifold_sampling/py/tests/benchmark_manifold_sampling.py
@@ -0,0 +1,91 @@
+# This wrapper tests various algorithms against the Benchmark functions from the
+# More and Wild SIOPT paper "Benchmarking derivative-free optimization algorithms"
+import os
+import sys
+
+import numpy as np
+import scipy as sp
+import scipy.io as sio
+
+sys.path.append("../")
+sys.path.append("../../../../BenDFO/py/")
+sys.path.append("../h_examples/")
+
+from calfun import calfun
+from dfoxs import dfoxs
+from manifold_sampling_primal import manifold_sampling_primal
+from piecewise_quadratic import piecewise_quadratic
+from pw_maximum import pw_maximum
+from pw_maximum_squared import pw_maximum_squared
+from pw_minimum import pw_minimum
+from pw_minimum_squared import pw_minimum_squared
+from quantile import quantile
+
+if not os.path.exists("benchmark_results"):
+    os.makedirs("benchmark_results")
+
+if not os.path.exists("mpc_test_files_smaller_Q"):
+    os.system("wget https://web.cels.anl.gov/~jmlarson/mpc_test_files_smaller_Q.zip")
+    os.system("unzip mpc_test_files_smaller_Q.zip")
+
+# www.mcs.anl.gov/~jlarson/mpc_test_files_smaller_Q.zip
+# Defines data for censored-L1 loss h instances
+C_L1_loss = np.loadtxt("mpc_test_files_smaller_Q/C_for_benchmark_probs.csv", delimiter=",")
+D_L1_loss = np.loadtxt("mpc_test_files_smaller_Q/D_for_benchmark_probs.csv", delimiter=",")
+# Defines data for piecewise_quadratic h instances
+Qzb = sio.loadmat("mpc_test_files_smaller_Q/Q_z_and_b_for_benchmark_problems_normalized_subset.mat")
+
+nfmax = 50
+factor = 10
+subprob_switch = "linprog"
+dfo = np.loadtxt("../../../../BenDFO/data/dfo.dat")
+filename = "./benchmark_results/manifold_sampling_py_nfmax=" + str(nfmax) + ".mat"
+
+Results = {}
+probs_to_solve = [0, 1, 6, 7, 42, 43, 44]
+
+
+subprob_switch = "linprog"
+
+
+hfuns = [pw_maximum_squared, pw_maximum, piecewise_quadratic, quantile, pw_minimum_squared, pw_minimum]
+
+for row, (nprob, n, m, factor_power) in enumerate(dfo[probs_to_solve, :]):
+    n = int(n)
+    m = int(m)
+    LB = -np.inf * np.ones((1, n))
+    UB = np.inf * np.ones((1, n))
+    x0 = dfoxs(n, nprob, factor**factor_power)
+
+    def Ffun(y):
+        out = calfun(y, m, int(nprob), "smooth", 0, vecout=True)
+        assert len(out) == m, "Incorrect output dimension"
+        return np.squeeze(out)
+
+    ind = np.where((C_L1_loss[:, 0] == probs_to_solve[row]) & (C_L1_loss[:, 1] == 0))
+    C = C_L1_loss[ind, 3 : m + 3]
+    D = D_L1_loss[ind, 3 : m + 3]
+
+    # Individual for piecewise_quadratic h instance
+    Qs = Qzb["Q_mat"][probs_to_solve[row], 0]
+    zs = Qzb["z_mat"][probs_to_solve[row], 0]
+    cs = Qzb["b_mat"][probs_to_solve[row], 0]
+
+    for i, hfun in enumerate(hfuns):
+        if hfun.__name__ == "piecewise_quadratic":
+
+            def hfun_to_pass(z, H0=None):
+                return piecewise_quadratic(z, H0, Qs=Qs, zs=zs, cs=cs)
+
+            X, F, h, xkin, flag = manifold_sampling_primal(hfun_to_pass, Ffun, x0, LB, UB, nfmax, subprob_switch)
+        else:
+            X, F, h, xkin, flag = manifold_sampling_primal(hfun, Ffun, x0, LB, UB, nfmax, subprob_switch)
+
+        Results["MSP_" + str(probs_to_solve[row] + 1) + "_" + str(i)] = {}
+        Results["MSP_" + str(probs_to_solve[row] + 1) + "_" + str(i)]["alg"] = "Manifold sampling"
+        Results["MSP_" + str(probs_to_solve[row] + 1) + "_" + str(i)]["problem"] = ["problem " + str(probs_to_solve[row] + 1) + " from More/Wild with hfun=" + str(hfun)]
+        Results["MSP_" + str(probs_to_solve[row] + 1) + "_" + str(i)]["Fvec"] = F
+        Results["MSP_" + str(probs_to_solve[row] + 1) + "_" + str(i)]["H"] = h
+        Results["MSP_" + str(probs_to_solve[row] + 1) + "_" + str(i)]["X"] = X
+
+    sp.io.savemat(filename, Results)
diff --git a/manifold_sampling/py/tests/matrix_Ffuns.py b/manifold_sampling/py/tests/matrix_Ffuns.py
new file mode 100644
index 00000000..4e18e932
--- /dev/null
+++ b/manifold_sampling/py/tests/matrix_Ffuns.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+n = 10  # Number of matrices to sum together
+m = 5  # Dimension of the square matrices
+
+np.random.seed(0)
+A = {}
+for i in range(n):
+    # Make 5 random positive-and negative-definite matricies
+    B = np.random.uniform(-1, 1, (m, m))
+    A[i] = (-1) ** i * 0.5 * (B @ B.T)
+
+
+def compute_M_and_eig(y):
+    y = y.squeeze()
+    assert len(y) == n, "Wrong input dimension"
+    M = np.zeros((m, m))
+    for i, val in enumerate(y):
+        M += val * A[i]
+
+    # eigvals, eigvecs = np.linalg.eigh(M)
+    eigvals, eigvecs = np.linalg.eig(M)
+    return eigvals, eigvecs
+
+
+def Ffun_default(y):
+    eigvals, _ = compute_M_and_eig(y)
+
+    return eigvals
+
+
+def Ffun_sort(y):
+    eigvals, _ = compute_M_and_eig(y)
+    eigvals.sort()
+
+    return eigvals
+
+
+def Ffun_mattrix(y):
+    y = y.squeeze()
+    assert len(y) == n, "Wrong input dimension"
+    M = np.zeros((m, m))
+    for i, val in enumerate(y):
+        M += val * A[i]
+
+    return np.reshape(M, (1, m**2)).squeeze()
diff --git a/manifold_sampling/py/tests/test_max_eig_hfun.py b/manifold_sampling/py/tests/test_max_eig_hfun.py
new file mode 100644
index 00000000..aefd7a24
--- /dev/null
+++ b/manifold_sampling/py/tests/test_max_eig_hfun.py
@@ -0,0 +1,128 @@
+import itertools
+import os
+import sys
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.optimize import linear_sum_assignment
+from scipy.spatial.distance import cdist
+
+sys.path.append("../")
+sys.path.append("../../../../BenDFO/py/")
+sys.path.append("../h_examples/")
+from manifold_sampling_primal import manifold_sampling_primal
+from matrix_Ffuns import Ffun_default, Ffun_sort, compute_M_and_eig, Ffun_mattrix
+from pw_maximum import pw_maximum as hfun
+from pw_maximum_eig import pw_maximum_eig as hfun2
+
+n = 10
+AllX = np.empty((0, n))
+AllEigVecs = {}
+count = 0
+
+
+def Ffun_fast_sort(y):
+    global AllX, AllEigVecs, count
+    eigvals, eigvecs = compute_M_and_eig(y)
+
+    if len(AllX):
+        ind_of_closest_past_point = np.argmin(np.linalg.norm(AllX - y, axis=1))
+        last_eig_vecs = AllEigVecs[ind_of_closest_past_point]
+        v = np.conjugate(np.transpose(eigvecs)) @ last_eig_vecs
+        b = np.argsort(v, axis=0)
+        best_perm = b[-1, :]
+        assert len(np.unique(best_perm)) == len(best_perm)
+
+        # print(best_perm)
+        # print(len(np.unique(best_perm)))
+        eigvals = eigvals[list(best_perm)]
+        eigvecs = eigvecs[:, best_perm]
+
+    AllX = np.vstack((AllX, y))
+    AllEigVecs[count] = eigvecs
+    count += 1
+
+    return eigvals
+
+
+def Ffun_slow_sort2(y):
+    global AllX, AllEigVecs, count
+    eigvals, eigvecs = compute_M_and_eig(y)
+
+    if len(AllX):
+        ind_of_closest_past_point = np.argmin(np.linalg.norm(AllX - y, axis=1))
+        last_eig_vecs = AllEigVecs[ind_of_closest_past_point]
+
+        alldists = cdist(eigvecs.T, last_eig_vecs.T)
+
+        row_ind, col_ind = linear_sum_assignment(alldists)
+
+        eigvals = eigvals[col_ind]
+        eigvecs = eigvecs[:, col_ind]
+
+    AllX = np.vstack((AllX, y))
+    AllEigVecs[count] = eigvecs
+    count += 1
+
+    return eigvals
+
+
+def Ffun_all_perms(y):
+    global AllX, AllEigVecs, count
+    eigvals, eigvecs = compute_M_and_eig(y)
+
+    if len(AllX):
+        ind_of_closest_past_point = np.argmin(np.linalg.norm(AllX - y, axis=1))
+        last_eig_vecs = AllEigVecs[ind_of_closest_past_point]
+
+        best_dist = np.inf
+        for p in itertools.permutations(np.arange(len(eigvals))):
+            this_perm_dist = np.linalg.norm(eigvecs[:, p] - last_eig_vecs)
+            if this_perm_dist < best_dist:
+                best_dist = this_perm_dist
+                best_perm = p
+
+        eigvals = eigvals[list(best_perm)]
+        eigvecs = eigvecs[:, best_perm]
+
+    AllX = np.vstack((AllX, y))
+    AllEigVecs[count] = eigvecs
+    count += 1
+
+    return eigvals
+
+
+nfmax = 80
+subprob_switch = "linprog"
+LB = -1 * np.ones((1, n))
+UB = np.ones((1, n))
+x0 = np.ones((1, n))
+# x0 = np.array([(-1)**(i) for i in range(n)])
+# x0 = np.array([(-1)**(i+1) for i in range(10)])
+
+X, F, h, xkin, flag = manifold_sampling_primal(hfun, Ffun_default, x0, LB, UB, nfmax, subprob_switch)
+plt.plot(h, linewidth=6, alpha=0.8, solid_joinstyle="miter", label="Default")
+
+X, F, h, xkin, flag = manifold_sampling_primal(hfun, Ffun_sort, x0, LB, UB, nfmax, subprob_switch)
+plt.plot(h, linewidth=6, alpha=0.8, solid_joinstyle="miter", label="Sorted")
+
+
+# plt.plot(F, linewidth=6, alpha=0.8, solid_joinstyle='miter')
+# plt.plot(F, linewidth=6, alpha=0.8)
+# plt.savefig("Initial.png", dpi=300, bbox_inches="tight")
+# plt.close()
+# X, F, h, xkin, flag = manifold_sampling_primal(hfun, Ffun_all_perms, x0, LB, UB, nfmax, subprob_switch)
+# plt.plot(h, linewidth=6, alpha=0.8, solid_joinstyle='miter', label="Complete search mapping")
+
+X, F, h, xkin, flag = manifold_sampling_primal(hfun, Ffun_slow_sort2, x0, LB, UB, nfmax, subprob_switch)
+plt.plot(h, linewidth=6, alpha=0.8, solid_joinstyle="miter", label="Hopefully smarter mapping")
+
+X, F, h, xkin, flag = manifold_sampling_primal(hfun, Ffun_all_perms, x0, LB, UB, nfmax, subprob_switch)
+plt.plot(h, linewidth=6, alpha=0.8, solid_joinstyle="miter", label="All permutations")
+
+X, F, h, xkin, flag = manifold_sampling_primal(hfun2, Ffun_mattrix, x0, LB, UB, nfmax, subprob_switch)
+plt.plot(h, linewidth=6, alpha=0.8, solid_joinstyle="miter", label="Matt")
+
+
+plt.legend()
+plt.savefig("Initial.png", dpi=300, bbox_inches="tight")
diff --git a/pounders/py/tests/regression_tests/benchmark_pounders4py.py b/pounders/py/tests/regression_tests/benchmark_pounders4py.py
index f40beb74..df6a8054 100644
--- a/pounders/py/tests/regression_tests/benchmark_pounders4py.py
+++ b/pounders/py/tests/regression_tests/benchmark_pounders4py.py
@@ -6,19 +6,21 @@
 
 import numpy as np
 import scipy as sp
-from oct2py import octave
 
 sys.path.append("../../../../minq/py/minq5/")  # Needed for spsolver=2
 import ibcdfo.pounders as pdrs
 
+BenDFO_root = "../../../../../BenDFO/"
+sys.path.append(BenDFO_root + "py/")  # Needed for spsolver=2
+from dfoxs import dfoxs
+from calfun import calfun
+
 os.makedirs("benchmark_results", exist_ok=True)
-np.seterr("raise")
+# np.seterr("raise")
 
 
 def doit():
-    bendfo_root = "../../../../../BenDFO/"
-    octave.addpath(bendfo_root + "m/")
-    dfo = np.loadtxt(bendfo_root + "data/dfo.dat")
+    dfo = np.loadtxt(BenDFO_root + "data/dfo.dat")
 
     ensure_still_solve_problems = 0
     if ensure_still_solve_problems:
@@ -36,11 +38,17 @@ def doit():
         m = int(m)
 
         def objective(y):
-            out = octave.feval("calfun_wrapper", y, m, nprob, "smooth", [], 1, 1)
+            # It is possible to have python use the same objective values via
+            # octave. This can be slow on some systems. To (for example)
+            # test difference between matlab and python, used the following
+            # line and add "from oct2py import octave" on a system with octave
+            # installed.
+            # out = octave.feval("calfun_wrapper", y, m, nprob, "smooth", [], 1, 1)
+            out = calfun(y, m, int(nprob), "smooth", 0, vecout=True)
             assert len(out) == m, "Incorrect output dimension"
             return np.squeeze(out)
 
-        X0 = octave.dfoxs(float(n), nprob, factor**factor_power).T
+        X0 = dfoxs(n, nprob, int(factor**factor_power))
         npmax = 2 * n + 1  # Maximum number of interpolation points [2*n+1]
         L = -np.inf * np.ones((1, n))  # 1-by-n Vector of lower bounds [zeros(1, n)]
         U = np.inf * np.ones((1, n))  # 1-by-n Vector of upper bounds [ones(1, n)]