# HG changeset patch
# User Julien Bect <julien.bect@supelec.fr>
# Date 1670531009 18000
# Node ID 972dcc46bb41db14843e974bf33d5849de12e445
# Parent  a18897e4c7b5dae4d75f31dace35d228c7dee4ff
sqp.m: Ignore missing lower/upper bounds (bug #32008)

* sqp.m: Avoid setting lower and upper bounds to +/-realmax when
unspecified in the function call. Set bounds to +/-Inf in main
function. Update constraint subfunction to correctly return empty
distance to constraint when bounds are unspecified.

diff -r a18897e4c7b5 -r 972dcc46bb41 scripts/optimization/sqp.m
--- a/scripts/optimization/sqp.m	Wed Dec 07 22:37:28 2022 -0500
+++ b/scripts/optimization/sqp.m	Thu Dec 08 15:23:29 2022 -0500
@@ -121,8 +121,7 @@
 ## inequality constraints @var{g} and @var{h}.  If the arguments are vectors
 ## then @var{x}(i) is bound by @var{lb}(i) and @var{ub}(i).  A bound can also
 ## be a scalar in which case all elements of @var{x} will share the same
-## bound.  If only one bound (lb, ub) is specified then the other will
-## default to (-@var{realmax}, +@var{realmax}).
+## bound.
 ##
 ## The seventh argument @var{maxiter} specifies the maximum number of
 ## iterations.  The default value is 100.
@@ -288,37 +287,43 @@
       ## constraint inequality function with bounds present
       lb_idx = ub_idx = true (size (x0));
       ub_grad = - (lb_grad = eye (rows (x0)));
+
+      ## if unspecified set ub and lb to +/-Inf, preserving single type
+      if (isempty (lb))
+        if (isa (x0, "single"))
+          lb = - inf (size (x0), "single");
+        else
+          lb = - inf (size (x0));
+        endif
+      endif
+
       if (isvector (lb))
-        globals.lb = tmp_lb = lb(:);
-        lb_idx(:) = tmp_idx = (lb != -Inf);
-        globals.lb = globals.lb(tmp_idx, 1);
+        lb = lb(:);
+        lb_idx(:) = (lb != -Inf);
+        globals.lb = lb(lb_idx, 1);
         lb_grad = lb_grad(lb_idx, :);
-      elseif (isempty (lb))
-        if (isa (x0, "single"))
-          globals.lb = tmp_lb = -realmax ("single");
-        else
-          globals.lb = tmp_lb = -realmax ();
-        endif
       else
         error ("sqp: invalid lower bound");
       endif
 
+      if (isempty (ub))
+        if (isa (x0, "single"))
+          ub = inf (size (x0), "single");
+        else
+          ub = inf (size (x0));
+        endif
+      endif
+
       if (isvector (ub))
-        globals.ub = tmp_ub = ub(:);
-        ub_idx(:) = tmp_idx = (ub != Inf);
-        globals.ub = globals.ub(tmp_idx, 1);
+        ub = ub(:);
+        ub_idx(:) = (ub != Inf);
+        globals.ub = ub(ub_idx, 1);
         ub_grad = ub_grad(ub_idx, :);
-      elseif (isempty (ub))
-        if (isa (x0, "single"))
-          globals.ub = tmp_ub = realmax ("single");
-        else
-          globals.ub = tmp_ub = realmax ();
-        endif
       else
         error ("sqp: invalid upper bound");
       endif
 
-      if (any (tmp_lb > tmp_ub))
+      if (any (globals.lb > globals.ub))
         error ("sqp: upper bound smaller than lower bound");
       endif
       bounds_grad = [lb_grad; ub_grad];
@@ -705,15 +710,32 @@
 
 
 function res = cf_ub_lb (x, lbidx, ubidx, globals)
+  ## function returning constraint evaluated at x, and distance between ub/lb 
+  ## and x, when they are specified (otherwise return empty)
 
-  ## combine constraint function with ub and lb
-  if (isempty (globals.cifcn))
-    res = [x(lbidx,1)-globals.lb; globals.ub-x(ubidx,1)];
+  ## inequality constraints
+  if (! isempty (globals.cifcn))
+    ci = feval (globals.cifcn, x);
+  else
+    ci = [];
+  endif
+
+  ## lower bounds
+  if (! isempty (globals.lb))
+    lb = x(lbidx, 1) - globals.lb;
   else
-    res = [feval(globals.cifcn,x); x(lbidx,1)-globals.lb;
-           globals.ub-x(ubidx,1)];
+    lb = [];
   endif
 
+  ## upper bounds
+  if (! isempty (globals.ub))
+    ub = globals.ub - x(ubidx, 1);
+  else
+    ub = [];
+  end
+
+  res = [ci; lb; ub];
+
 endfunction