Repository: POptUS/IBCDFO (language: python)

Issue:
Fixing quantile computation when nargin==2



Discussion / hints:

Fix patch:
diff --git a/manifold_sampling/m/build_p_models.m b/manifold_sampling/m/build_p_models.m
index d815299a..e00f5d6d 100644
--- a/manifold_sampling/m/build_p_models.m
+++ b/manifold_sampling/m/build_p_models.m
@@ -7,9 +7,7 @@
     % Evaluate model-improving points if necessary
     if ~valid && ~isempty(Gres)
         [Mdir, np, valid] = formquad(X(1:nf, :), F(1:nf, :), delta, xkin, fq_pars.npmax, fq_pars.Par, 1);
-        if valid
-            error('what to do here');
-        end
+        assert(~valid, "How did this switch to being valid?");
         [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);
     end
diff --git a/manifold_sampling/m/call_user_scripts.m b/manifold_sampling/m/call_user_scripts.m
index 46095926..5f231ca3 100644
--- a/manifold_sampling/m/call_user_scripts.m
+++ b/manifold_sampling/m/call_user_scripts.m
@@ -38,7 +38,7 @@
 
         [h_dummy1, grad_dummy1, hash_dummy] = hfun(F(nf, :));
         [h_dummy2, grad_dummy2] = hfun(F(nf, :), hash_dummy);
-        assert(any(h_dummy1 == h_dummy2), "hfun values don't agree when being re-called with the same inputs");
-        assert(all(all(grad_dummy1 == grad_dummy2)), "hfun gradients don't agree when being re-called with the same inputs");
+        assert(any(h_dummy1 == h_dummy2), ['hfun values do not agree when ' func2str(hfun) ' is re-called with the same inputs']);
+        assert(all(all(grad_dummy1 == grad_dummy2)), ['hfun gradients do not agree when '  func2str(hfun) ' is re-called with the same inputs']);
     end
 end
diff --git a/manifold_sampling/m/h_examples/censored_L1_loss.m b/manifold_sampling/m/h_examples/censored_L1_loss.m
index 54f1f264..d953238c 100644
--- a/manifold_sampling/m/h_examples/censored_L1_loss.m
+++ b/manifold_sampling/m/h_examples/censored_L1_loss.m
@@ -16,7 +16,9 @@
 %   ... in the calculation of h and grads.
 
 % Get data from outside of this function
-global C D eqtol n_h hvals_mat
+global C D n_h hvals_mat
+
+eqtol = 1e-8;
 
 % Ensure column vector and collect dimensions
 z = z(:);
@@ -106,9 +108,6 @@
         end
         h(k) = sum(vals(:, k));
     end
-
-else
-    error('Too many inputs to function');
 end
 
 if ~isempty(n_h)
diff --git a/manifold_sampling/m/h_examples/censored_L1_loss_quad_MSG.m b/manifold_sampling/m/h_examples/censored_L1_loss_quad_MSG.m
index 64882d8e..519cddb3 100644
--- a/manifold_sampling/m/h_examples/censored_L1_loss_quad_MSG.m
+++ b/manifold_sampling/m/h_examples/censored_L1_loss_quad_MSG.m
@@ -16,7 +16,9 @@
 %   ... in the calculation of h and grads.
 
 % Get data from outside of this function
-global C D eqtol n_h hvals_mat
+global C D n_h hvals_mat
+
+eqtol = 1e-8;
 
 % Ensure column vector and collect dimensions
 zin = z(:);
@@ -116,9 +118,6 @@
         end
         h(k) = sum(vals(:, k));
     end
-
-else
-    error('Too many inputs to function');
 end
 
 if ~isempty(n_h)
diff --git a/manifold_sampling/m/h_examples/max_sum_beta_plus_const_viol.m b/manifold_sampling/m/h_examples/max_sum_beta_plus_const_viol.m
index 259dad80..930fc5f1 100644
--- a/manifold_sampling/m/h_examples/max_sum_beta_plus_const_viol.m
+++ b/manifold_sampling/m/h_examples/max_sum_beta_plus_const_viol.m
@@ -14,7 +14,9 @@
 %   Similarly, if H0 has a 1 in position i uses z_i in the calculation of h and grads.
 
 % Get data from outside of this function
-global h_activity_tol p1 alpha
+p1 = length(z) - 1;
+alpha = 0;
+h_activity_tol = 1e-8;
 
 p = length(z);
 
@@ -63,7 +65,4 @@
         end
         h(k) = h1 + h2;
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/piecewise_quadratic.m b/manifold_sampling/m/h_examples/piecewise_quadratic.m
index d11ba77a..723d6632 100644
--- a/manifold_sampling/m/h_examples/piecewise_quadratic.m
+++ b/manifold_sampling/m/h_examples/piecewise_quadratic.m
@@ -55,7 +55,4 @@
         h(k) = (z - zs(:, j))' * Qs(:, :, j) * (z - zs(:, j)) + cs(j);
         grads(:, k) = 2 * Qs(:, :, j) * (z - zs(:, j));
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/piecewise_quadratic_1.m b/manifold_sampling/m/h_examples/piecewise_quadratic_1.m
index 5286fd7a..c4feb867 100644
--- a/manifold_sampling/m/h_examples/piecewise_quadratic_1.m
+++ b/manifold_sampling/m/h_examples/piecewise_quadratic_1.m
@@ -1,6 +1,6 @@
 function [h, grads, Hash] = piecewise_quadratic_1(z, H0)
 % Evaluates the piecewise quadratic function
-%   max_j { 0.5*z'*Q_j*z + c_j'*z + b_j }
+%   max_j { 0.5*z'*Qs_j*z + zs_j'*z + cs_j }
 %
 % Inputs:
 %  z:              [1 x p]   point where we are evaluating h
@@ -11,14 +11,15 @@
 %  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)
 
+global Qs zs cs
+
 z = z(:);
-global Q c b
 
 if nargin == 1
-    J = length(b);
+    [n, J] = size(zs);
     manifolds = zeros(1, J);
     for j = 1:J
-        manifolds(j) = 0.5 * z' * Q(:, :, j) * z + c(:, j)' * z + b(j);
+        manifolds(j) = 0.5 * z' * Qs(:, :, j) * z + zs(:, j)' * z + cs(j);
     end
 
     h = max(manifolds);
@@ -27,11 +28,12 @@
     rtol = 1e-12;
     inds = find(abs(h - manifolds) <= atol + rtol * abs(h - manifolds));
 
-    grads = Q(:, :, inds) * z + c(:, inds);
+    grads = zeros(n, length(inds));
 
     Hash = cell(1, length(inds));
     for j = 1:length(inds)
         Hash{j} = int2str(inds(j));
+        grads(:, j) = Qs(:, :, inds(j)) * z + zs(:, inds(j));
     end
 
 elseif nargin == 2
@@ -41,10 +43,7 @@
 
     for k = 1:J
         j = str2num(H0{k});
-        h(k) = 0.5 * z' * Q(:, :, j) * z + c(:, j)' * z + b(j);
-        grads(:, k) = Q(:, :, j) * z + c(:, j);
+        h(k) = 0.5 * z' * Qs(:, :, j) * z + zs(:, j)' * z + cs(j);
+        grads(:, k) = Qs(:, :, j) * z + zs(:, j);
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/pw_maximum.m b/manifold_sampling/m/h_examples/pw_maximum.m
index 940843d4..3dd44222 100644
--- a/manifold_sampling/m/h_examples/pw_maximum.m
+++ b/manifold_sampling/m/h_examples/pw_maximum.m
@@ -40,7 +40,4 @@
         h(k) = z(j);
         grads(j, k) = 1;
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/pw_maximum_squared.m b/manifold_sampling/m/h_examples/pw_maximum_squared.m
index 0268bf3b..ba4f1bae 100644
--- a/manifold_sampling/m/h_examples/pw_maximum_squared.m
+++ b/manifold_sampling/m/h_examples/pw_maximum_squared.m
@@ -41,7 +41,4 @@
         h(k) = z(j)^2;
         grads(j, k) = 2 * z(j);
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/pw_minimum.m b/manifold_sampling/m/h_examples/pw_minimum.m
index 679b093a..1ccc9666 100644
--- a/manifold_sampling/m/h_examples/pw_minimum.m
+++ b/manifold_sampling/m/h_examples/pw_minimum.m
@@ -40,7 +40,4 @@
         h(k) = z(j);
         grads(j, k) = 1;
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/pw_minimum_squared.m b/manifold_sampling/m/h_examples/pw_minimum_squared.m
index dff1623b..790c7e1d 100644
--- a/manifold_sampling/m/h_examples/pw_minimum_squared.m
+++ b/manifold_sampling/m/h_examples/pw_minimum_squared.m
@@ -41,7 +41,4 @@
         h(k) = z(j)^2;
         grads(j, k) = 2 * z(j);
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/quantile.m b/manifold_sampling/m/h_examples/quantile.m
index 32e305f2..e9ac7d12 100644
--- a/manifold_sampling/m/h_examples/quantile.m
+++ b/manifold_sampling/m/h_examples/quantile.m
@@ -14,18 +14,17 @@
 q = 2; % hard-coded. Eventually, make global or (preferred) pass it in as an argument.
 
 z = z(:);
-n = length(z);
-z2 = z.^2;
-[sortedz2, sort_inds] = sort(z2);
 
 if nargin == 1
+    z2 = z.^2;
+    sortedz2 = sort(z2);
     h = sortedz2(q);
 
     atol = 1e-8;
     rtol = 1e-8;
     inds = find(abs(h - z2) <= atol + rtol * abs(z2));
 
-    grads = zeros(n, length(inds));
+    grads = zeros(length(z), length(inds));
 
     Hash = cell(1, length(inds));
     for j = 1:length(inds)
@@ -37,14 +36,10 @@
     J = length(H0);
     h = zeros(1, J);
     grads = zeros(length(z), J);
-    z = z(sort_inds);
 
     for k = 1:J
         j = str2num(H0{k});
         h(k) = z(j)^2;
         grads(j, k) = 2 * z(j);
     end
-
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/h_examples/sum_squared.m b/manifold_sampling/m/h_examples/sum_squared.m
index c9f14ee9..4f3d5446 100644
--- a/manifold_sampling/m/h_examples/sum_squared.m
+++ b/manifold_sampling/m/h_examples/sum_squared.m
@@ -27,6 +27,4 @@
     h = sum(z.^2);
     grads = zeros(length(z), 1);
     grads(:) = 2 * z;
-else
-    error('Too many inputs to function');
 end
diff --git a/manifold_sampling/m/manifold_sampling_primal.m b/manifold_sampling/m/manifold_sampling_primal.m
index 31c6a6b5..5b2faf0a 100644
--- a/manifold_sampling/m/manifold_sampling_primal.m
+++ b/manifold_sampling/m/manifold_sampling_primal.m
@@ -47,8 +47,8 @@
 try
     F0 = Ffun(x0);
 catch e
-    fprintf(1, 'The identifier was:\n%s', e.identifier);
-    fprintf(1, 'There was an error! The message was:\n%s', e.message);
+    fprintf(1, 'The identifier was:\n%s\n\n', e.identifier);
+    fprintf(1, 'There was an error! The message was:\n%s\n', e.message);
     warning('Problem using Ffun. Exiting');
     X = [];
     F = [];


Test patch:
diff --git a/manifold_sampling/m/tests/Testmanifoldsampling.m b/manifold_sampling/m/tests/Testmanifoldsampling.m
new file mode 100644
index 00000000..0efb8724
--- /dev/null
+++ b/manifold_sampling/m/tests/Testmanifoldsampling.m
@@ -0,0 +1,13 @@
+classdef Testmanifoldsampling < matlab.unittest.TestCase
+    methods (Test)
+
+        function shortTests(testCase)
+            test_failing_objective;
+        end
+
+        function longTest(testCase)
+            benchmark_manifold_sampling;
+        end
+
+    end
+end
diff --git a/manifold_sampling/m/tests/benchmark_manifold_sampling.m b/manifold_sampling/m/tests/benchmark_manifold_sampling.m
index 5b84828d..9c5b3d68 100644
--- a/manifold_sampling/m/tests/benchmark_manifold_sampling.m
+++ b/manifold_sampling/m/tests/benchmark_manifold_sampling.m
@@ -2,15 +2,8 @@
 % More and Wild SIOPT paper "Benchmarking derivative-free optimization algorithms"
 function [] = benchmark_manifold_sampling()
 
-addpath('../');
-addpath('../../../../BenDFO/m/');
-addpath('../../../../BenDFO/data/');
-addpath('../h_examples/');
-addpath('../../../pounders/m'); % formquad, bmpts, boxline, phi2eval
-
-mkdir('benchmark_results');
-
-nfmax = 500;
+global C D Qs zs cs
+nfmax = 100;
 factor = 10;
 
 subprob_switch = 'linprog';
@@ -20,8 +13,17 @@
 filename = ['./benchmark_results/manifold_samplingM_nfmax=' num2str(nfmax) '.mat'];
 Results = cell(1, 53);
 
+if ~exist("mpc_test_files_smaller_Q", "dir")
+    system("wget https://web.cels.anl.gov/~jmlarson/mpc_test_files_smaller_Q.zip");
+    system("unzip mpc_test_files_smaller_Q.zip");
+end
+
+C_L1_loss = load('mpc_test_files_smaller_Q/C_for_benchmark_probs.csv');
+D_L1_loss = load('mpc_test_files_smaller_Q/D_for_benchmark_probs.csv');
+Qzb = load('mpc_test_files_smaller_Q/Q_z_and_b_for_benchmark_problems_normalized_subset.mat')';
+
 % for row = find(cellfun(@length,Results)==0)
-for row = [1, 2, 7, 8, 43, 44, 45]
+for row = [2, 1, 7, 8, 43, 44, 45]
     nprob = dfo(row, 1);
     n = dfo(row, 2);
     m = dfo(row, 3);
@@ -36,18 +38,20 @@
 
     xs = dfoxs(n, nprob, factor^factor_power);
 
+    ind = find(C_L1_loss(:, 1) == row & C_L1_loss(:, 2) == 1);
+    C = C_L1_loss(ind, 4:m + 3);
+    D = D_L1_loss(ind, 4:m + 3);
+    Qs = Qzb.Q_mat{row, 1};
+    zs = Qzb.z_mat{row, 1};
+    cs = Qzb.b_mat{row, 1};
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % Manifold sampling
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-    for jj = 1:2
-        if jj == 1
-            hfun = @pw_maximum_squared;
-        end
-        if jj == 2
-            hfun = @pw_minimum_squared;
-        end
-
+    jj = 1;
+    % for hfuns = {@censored_L1_loss_quad_MSG}
+    for hfuns = {@censored_L1_loss, @max_sum_beta_plus_const_viol, @piecewise_quadratic, @piecewise_quadratic_1, @pw_maximum,  @pw_maximum_squared, @pw_minimum, @pw_minimum_squared, @quantile}
+        hfun = hfuns{1};
         Ffun = @(x)calfun_wrapper(x, BenDFO, 'smooth');
         x0 = xs';
 
@@ -58,6 +62,7 @@
         Results{jj, row}.Fvec = F;
         Results{jj, row}.H = h;
         Results{jj, row}.X = X;
+        jj = jj + 1;
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         % save(filename, 'Results');
         % save('-mat7-binary', filename, 'Results') % Octave save
diff --git a/manifold_sampling/m/tests/failing_objective.m b/manifold_sampling/m/tests/failing_objective.m
new file mode 100644
index 00000000..96a44ae5
--- /dev/null
+++ b/manifold_sampling/m/tests/failing_objective.m
@@ -0,0 +1,9 @@
+function [fvec] = failing_objective(x)
+% This function will produce a NaN in one component on average in 1 out of 10
+% queries, otherwise, it produces returns x.
+
+fvec = x;
+
+if rand < 0.1
+    fvec(1) = NaN;
+end
diff --git a/manifold_sampling/m/tests/make_coverage.m b/manifold_sampling/m/tests/make_coverage.m
new file mode 100644
index 00000000..a0599bb5
--- /dev/null
+++ b/manifold_sampling/m/tests/make_coverage.m
@@ -0,0 +1,14 @@
+clear all;
+setup_paths;
+suite = testsuite("Testmanifoldsampling.m");
+import matlab.unittest.plugins.CodeCoveragePlugin
+import matlab.unittest.plugins.codecoverage.CoverageReport
+runner = testrunner("textoutput");
+sourceCodeFolder = "../";
+% sourceCodeFolder = "../../../../minq/m/";
+reportFolder = "coverageReport";
+reportFormat = CoverageReport(reportFolder);
+% p = CodeCoveragePlugin.forFolder(sourceCodeFolder, "Producing", reportFormat, "IncludingSubfolders", true);
+p = CodeCoveragePlugin.forFolder(["../", "../h_examples/"], "Producing", reportFormat);
+runner.addPlugin(p);
+results = runner.run(suite);
diff --git a/manifold_sampling/m/tests/setup_paths.m b/manifold_sampling/m/tests/setup_paths.m
new file mode 100644
index 00000000..10469591
--- /dev/null
+++ b/manifold_sampling/m/tests/setup_paths.m
@@ -0,0 +1,7 @@
+addpath('../');
+addpath('../../../../BenDFO/m/');
+addpath('../../../../BenDFO/data/');
+addpath('../h_examples/');
+addpath('../../../pounders/m'); % formquad, bmpts, boxline, phi2eval
+
+mkdir('benchmark_results');
diff --git a/manifold_sampling/m/tests/test_failing_objective.m b/manifold_sampling/m/tests/test_failing_objective.m
new file mode 100644
index 00000000..46a6966f
--- /dev/null
+++ b/manifold_sampling/m/tests/test_failing_objective.m
@@ -0,0 +1,27 @@
+% This tests the ability of the MATLAB pounders implementation to
+% gracefully terminate when receiving a NaN as a component in the vector of
+% objective output.
+
+function [] = test_failing_objective()
+
+nfmax = 1000;
+n = 3;
+xs = [10; 20; 30];
+L = -Inf(1, n);
+U = Inf(1, n);
+X0 = xs';
+
+rand('seed', 1);
+objective = @(x)failing_objective(x);
+hfun = @pw_maximum_squared;
+
+% [X, F, h, xkin, flag] = manifold_sampling_primal(hfun, objective, X0, L, U, nfmax, 'linprog');
+% assert(flag == -3, "No NaN was encountered in this test, but (with high probability) should have been.");
+
+% Intentionally not passing a function for an objective
+[X, F, h, xkin, flag] = manifold_sampling_primal(hfun, X0, X0, L, U, nfmax, 'linprog');
+assert(flag == -1, "Should have failed");
+
+objective = @(x) x;
+[X, F, h, xkin, flag] = manifold_sampling_primal(hfun, objective, X0, [L 1], U, nfmax, 'linprog');
+assert(flag == -1, "Should have failed");