Add create TRSP solver function

pounders/m/create_trsp_solver.m

@@ -0,0 +1,74 @@
+function [solver] = create_trsp_solver(spsolver)
+    % Please refer to the documentation for the Python version of this
+    % function.
+
+    arguments
+        spsolver {mustBeScalarOrEmpty, mustBeNonempty, mustBeInteger}
+    end
+
+    % ----- HARDCODED VALUES
+    % Ensure that these match the analogous constants implemented for
+    % POUNDERS/Python.  These same values might be used in the code that
+    % directly tests this function.
+    %
+    % Both MATLAB and Python implementations should declare the union of all
+    % solvers available even if they don't support one or more of the solvers.
+    TRSP_SOLVER_SIMPLE = 1;
+    TRSP_SOLVER_MINQ5 = 2;
+    TRSP_SOLVER_MINQ8 = 3;
+
+    % ----- DEFINE POUNDERS-COMPATIBLE INTERFACES ON SOLVERS
+    % Stefan's crappy 10 line solver
+    function [Xsp, mdec, found_solution] = bqmin_wrapper(H, G, Low, Upp)
+        % Assume that solver error checks its arguments thoroughly and that
+        % solver always finds valid solution.
+        found_solution = true;
+        [Xsp, mdec] = bqmin(H, G, Low, Upp);
+    end
+
+    % Arnold Neumaier's minq5
+    function [Xsp, mdec, found_solution] = minq5_wrapper(H, G, Low, Upp)
+        % Assume that solver error checks its arguments thoroughly.
+        xx = zeros(size(H, 1), 1);
+        [Xsp, mdec, minq_err] = minqsw(0, G, H, Low', Upp', 0, xx);
+        % Continuous function restricted to (compact) k-cell.
+        assert(minq_err ~= 1);
+        % See comments in Python version of this function for info on handling
+        % error code 99.
+        % assert(minq_err ~= 99);
+        found_solution = (minq_err >= 0);
+    end
+
+    % Arnold Neumaier's minq8
+    function [Xsp, mdec, found_solution] = minq8_wrapper(H, G, Low, Upp)
+        % Assume that solver error checks its arguments thoroughly and that
+        % solver always finds valid solution.
+        found_solution = true;
+
+        n = size(H, 1);
+
+        data.gam = 0;
+        data.c = G;
+        data.b = zeros(n, 1);
+        [tmp1, tmp2] = ldl(H);
+        data.D = diag(tmp2);
+        data.A = tmp1';
+
+        [Xsp, mdec] = minq8(data, Low', Upp', zeros(n, 1), 10 * n);
+    end
+
+    % ----- IDENTIFY DESIRED SOLVER
+    if spsolver == TRSP_SOLVER_SIMPLE
+        warning("POUNDERS:simpleTrspSolver", ...
+                "The simple TRSP solver should only be used for testing or debugging");
+        solver = @bqmin_wrapper;
+    elseif spsolver == TRSP_SOLVER_MINQ5
+        check_minq_installation(5);
+        solver = @minq5_wrapper;
+    elseif spsolver == TRSP_SOLVER_MINQ8
+        check_minq_installation(8);
+        solver = @minq8_wrapper;
+    else
+        error('POUNDERS:badValue', sprintf("Invalid TRSP solver %d", spsolver));
+    end
+end

pounders/m/pounders.m

@@ -39,7 +39,11 @@
 %           * 1 - Debugging level of output to screen
 %           * 2 - More verbose screen output
 %
-%       * **spsolver** - Trust-region subproblem solver flag (default is 2, not recommended to change)
+%       * **spsolver** - Trust-region subproblem solver flag
+%
+%           * 2 - Arnold Neumaier's minq5 solver (default and recommended)
+%           * 3 - Arnold Neumaier's minq8 solver
+%
 %       * **hfun** - Outer function :math:`\hfun` that maps given
 %         :math:`\Ffun(\psp)` to scalars for minimization (default is
 %         sum-of-squares that yields :math:`f`.)
@@ -116,7 +120,7 @@
     Options.delta_inact = 0.75;
 end
 if ~isfield(Options, 'spsolver')
-    Options.spsolver = 2;
+    Options.spsolver = 2; % Use minq5 by default
 end
 
 if isfield(Options, 'hfun')
@@ -129,9 +133,6 @@
     hfun = @h_leastsquares;
     combinemodels = @combine_leastsquares;
 end
-if ~isfield(Options, 'spsolver')
-    Options.spsolver = 2; % Use minq5 by default
-end
 if ~isfield(Options, 'printf')
     Options.printf = 0; % Don't print by default
 end
@@ -152,7 +153,7 @@
 nfs = Prior.nfs;
 
 delta = delta_0;
-spsolver = Options.spsolver;
+solve_trsp = create_trsp_solver(Options.spsolver);
 delta_max = Options.delta_max;
 delta_min = Options.delta_min;
 gamma_dec = Options.gamma_dec;
@@ -161,12 +162,6 @@
 printf = Options.printf;
 delta_inact = Options.delta_inact;
 
-if spsolver == 2 % Arnold Neumaier's minq5
-%    check_minq_installation(5);
-elseif spsolver == 3 % Arnold Neumaier's minq8
-    check_minq_installation(8);
-end
-
 % 0. Check inputs
 [flag, X_0, np_max, F_0, Low, Upp, xk_in] = ...
     checkinputss(Ffun, X_0, n, Model.np_max, nf_max, g_tol, delta, nfs, m, Prior.F_init, Prior.xk_in, Low, Upp);
@@ -345,26 +340,12 @@
     % 3. Solve the subproblem min{G'*s+.5*s'*H*s : Lows <= s <= Upps }
     Lows = max(Low - X(xk_in, :), -delta);
     Upps = min(Upp - X(xk_in, :), delta);
-    if spsolver == 1 % Stefan's crappy 10line solver
-        [Xsp, mdec] = bqmin(H, G, Lows, Upps);
-    elseif spsolver == 2 % Arnold Neumaier's minq5
-        [Xsp, mdec, minq_err] = minqsw(0, G, H, Lows', Upps', 0, zeros(n, 1));
-        if minq_err < 0
-            [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -4);
-            return
-        end
-
-    elseif spsolver == 3 % Arnold Neumaier's minq8
-
-        data.gam = 0;
-        data.c = G;
-        data.b = zeros(n, 1);
-        [tmp1, tmp2] = ldl(H);
-        data.D = diag(tmp2);
-        data.A = tmp1';
-
-        [Xsp, mdec] = minq8(data, Lows', Upps', zeros(n, 1), 10 * n);
+    [Xsp, mdec, found_solution] = solve_trsp(H, G, Lows, Upps);
+    if ~found_solution
+        [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -4);
+        return
     end
+
     Xsp = Xsp'; % Solvers currently work with column vectors
     step_norm = norm(Xsp, inf);