Modify s2mpjlib.py

Author: OptHuang

.github/workflows/sync_s2mpj.yml

@@ -39,7 +39,7 @@ jobs:
           # 2. Clean managed files in current workspace for fresh sync
           # Only remove the specific files/folders we are syncing.
           mkdir -p src
-          rm -rf src/python_problems src/list_of_python_problems src/s2mpjlib.py
+          rm -rf src/python_problems src/list_of_python_problems
 
           # 3. Copy only the specific files/folders we want from the temporary upstream directory
           echo "Copying selected files..."
@@ -58,13 +58,6 @@ jobs:
             echo "::warning::File 'list_of_python_problems' not found in upstream"
           fi
 
-          # Copy the 's2mpjlib.py' file
-          if [ -f "temp_upstream/s2mpjlib.py" ]; then
-            cp temp_upstream/s2mpjlib.py src/
-          else
-             echo "::warning::File 's2mpjlib.py' not found in upstream"
-          fi
-
           # 4. Remove the temporary upstream directory
           rm -rf temp_upstream
 

s2mpj_tools.py

@@ -182,7 +182,6 @@ def s2mpj_load(problem_name, *args):
     beq = bx[idx_aeq] + cu[idx_aeq] if bx is not None and idx_aeq.size > 0 else None
     bub = np.concatenate([bx[idx_aub_le] + cu[idx_aub_le], -bx[idx_aub_ge] - cl[idx_aub_ge]]) if bx is not None and (idx_aub_le.size > 0 or idx_aub_ge.size > 0) else None
 
-    getidx = lambda y, idx: y[idx] if (y is not None and idx.size > 0) else None
     # Construct nonlinear constraint functions
     # Remind that in S2MPJ, the constraints are defined as:
     #  cl <= c(x) <= cu
@@ -192,32 +191,49 @@ def s2mpj_load(problem_name, *args):
     #  cub(x) = [c(x)[idx_cle] - cu[idx_cle];
     #           -c(x)[idx_cge] + cl[idx_cge]] <= 0
     def ceq(x):
+        if idx_ceq.size == 0: return None
         y = _getcx(p, x)
-        if idx_ceq.size > 0:
-            z = getidx(y, idx_ceq) - cu[idx_ceq]
-        else:
-            z = None
-        return None if z is None or (hasattr(z, "size") and z.size == 0) else z
+        if y is None: return None
+        z = y[idx_ceq] - cu[idx_ceq]
+        return None if (hasattr(z, "size") and z.size == 0) else z
     def cub(x):
+        if idx_cle.size == 0 and idx_cge.size == 0: return None
         y = _getcx(p, x)
-        le = getidx(y, idx_cle) - cu[idx_cle] if idx_cle.size > 0 else None
-        ge = getidx(y, idx_cge) - cl[idx_cge] if idx_cge.size > 0 else None
-        if le is None and ge is None:
-            return None
-        if ge is None:
-            return le
-        if le is None:
-            return -ge
+        if y is None: return None
+        le = y[idx_cle] - cu[idx_cle] if idx_cle.size > 0 else None
+        ge = y[idx_cge] - cl[idx_cge] if idx_cge.size > 0 else None
+        if le is None and ge is None: return None
+        if ge is None: return le
+        if le is None: return -ge
         return np.concatenate([le, -ge])
-
-    getidx_list = lambda y, idx: [y[i] for i in idx] if y is not None else []
-    hceq = lambda x: getidx_list(_getHx(p, x), idx_ceq)
-    hcub = lambda x: getidx_list(_getHx(p, x), idx_cle) + getidx_list(_getHx(p, x), idx_cge)
-
-    getidx_mat = lambda y, idx: y[idx, :] if y is not None else None
-    jceq = lambda x: getidx_mat(_getJx(p, x), idx_ceq)
-    jcub = lambda x: np.vstack([getidx_mat(_getJx(p, x), idx_cle),
-                                -getidx_mat(_getJx(p, x), idx_cge)]) if getidx_mat(_getJx(p, x), idx_cle) is not None else None
+    def hceq(x):
+        if idx_ceq.size == 0: return []
+        Hx = _getHx(p, x)
+        if Hx is None: raise RuntimeError("Hessian evaluation failed")
+        return [Hx[i] for i in idx_ceq]
+    def hcub(x):
+        if idx_cle.size == 0 and idx_cge.size == 0: return []
+        Hx = _getHx(p, x)
+        if Hx is None: raise RuntimeError("Hessian evaluation failed")
+        le = [Hx[i] for i in idx_cle] if idx_cle.size > 0 else []
+        ge = [Hx[i] for i in idx_cge] if idx_cge.size > 0 else []
+        return le + ge
+
+    def jceq(x):
+        if idx_ceq.size == 0: return np.empty((0, p.n))
+        Jx = _getJx(p, x)
+        if Jx is None: return None
+        return Jx[idx_ceq, :]
+    def jcub(x):
+        if idx_cle.size == 0 and idx_cge.size == 0: return np.empty((0, p.n))
+        Jx = _getJx(p, x)
+        if Jx is None: return None
+        le = Jx[idx_cle, :] if idx_cle.size > 0 else None
+        ge = Jx[idx_cge, :] if idx_cge.size > 0 else None
+        if le is None and ge is None: return np.empty((0, p.n))
+        if ge is None: return le
+        if le is None: return -ge
+        return np.vstack([le, -ge])
 
     # Construct the Problem instance.
     problem = Problem(fun, x0, name=name, xl=xl, xu=xu, aub=aub, bub=bub, aeq=aeq, beq=beq, cub=cub, ceq=ceq, grad=grad, hess=hess, jcub=jcub, jceq=jceq, hcub=hcub, hceq=hceq)
@@ -497,7 +513,8 @@ def _gethess(p, is_feasibility, x):
             try:
                 _, _, h = p.fgHx(x)
                 h = h.toarray() if hasattr(h, 'toarray') else h
-            except Exception:
+            except Exception as e:
+                print(f"Error eval Hx: {e}")
                 h = None
     return h
 
@@ -507,7 +524,8 @@ def _getcx(p, x):
         try:
             c = p.cx(x)
             c = c.toarray() if hasattr(c, 'toarray') else c
-        except Exception:
+        except Exception as e:
+            print(f"Error eval cx: {e}")
             c = None
     return c
 
@@ -517,7 +535,8 @@ def _getJx(p, x):
         try:
             _, j = p.cJx(x)[:2]
             j = j.toarray() if hasattr(j, 'toarray') else j
-        except Exception:
+        except Exception as e:
+            print(f"Error eval Jx: {e}")
             j = None
     return j
 
@@ -529,6 +548,7 @@ def _getHx(p, x):
             for i in range(len(h)):
                 if hasattr(h[i], 'toarray'):
                     h[i] = h[i].toarray()
-        except Exception:
+        except Exception as e:
+                print(f"Error eval Hx: {e}")
                 h = None
         return h

src/s2mpjlib.py

@@ -427,13 +427,13 @@ def evalgrsum( self, isobj, glist, x, nargout ):
             #  Evaluate the linear term, if any.
 
             if hasattr(self,"gconst") and ig < len(self.gconst) and not self.gconst[ig] is None:
-                fin = float(-self.gconst[ig])
+                fin = float(np.array(-self.gconst[ig]).item())
             else:
                 fin = 0
             if has_A and ig < sA1:
                 gin           = np.zeros( (n, 1) )
                 gin[:sA2, :1] = self.A[ ig, :sA2 ].T.toarray()
-                fin           = float( fin + gin.T .dot(x) )
+                fin           = float(np.array( fin + gin.T .dot(x) ).item())
             elif nargout >= 2:
                 gin =  np.zeros(( n, 1 ))
 
@@ -448,7 +448,7 @@ def evalgrsum( self, isobj, glist, x, nargout ):
                     iel    = self.grelt[ ig ][ iiel ]         #  the element's index
                     efname = self.elftype[ iel ]              #  the element's ftype
                     irange = [iv for iv in self.elvar[ iel ]] #  the elemental variable's indeces 
-                    xiel   = x[ np.array(irange) ]            #  the elemental variable's values
+                    xiel   = x[ np.array(irange) ].flatten().flatten()   #  the elemental variable's values
 
                     if  hasattr( self, 'grelw' ) and ig <= len( self.grelw ) and not self.grelw[ig] is None :
                         has_weights = True;
@@ -621,11 +621,11 @@ def evalgrsum( self, isobj, glist, x, nargout ):
                     print( "ig = ", ig, "  cx(final) = ", cx )#D
         if isobj:
             if nargout == 1:
-                return float( fx )
+                return float(np.array(fx).item())
             elif nargout == 2:
-                return float( fx ), gx.reshape(-1,1)
+                return float(np.array(fx).item()), gx.reshape(-1,1)
             elif nargout == 3:
-                return float( fx ), gx.reshape(-1,1), Hx
+                return float(np.array(fx).item()), gx.reshape(-1,1), Hx
         else:
             if nargout == 1:
                 return cx
@@ -711,7 +711,7 @@ def evalHJv( self, mode, glist, x, v, y ):
         #  Evaluate the quadratic term, if any.
 
         if mode == "Hv" and hasattr( self, "H" ):
-            HJv += self.H.dot(v);
+            HJv += self.H.dot(v).item();
 
         if debug: #D
             print( "HJv(quadratic) = ", HJv )
@@ -750,13 +750,13 @@ def evalHJv( self, mode, glist, x, v, y ):
             #  Evaluate the linear term, if any.
 
             if hasattr( self, "gconst" ) and ig < len( self.gconst ) and not self.gconst[ ig ] is None:
-                fin = float(-self.gconst[ ig ] )
+                fin = float(np.array(-self.gconst[ ig ] ).item())
             else:
                 fin = 0.0
             gin = np.zeros((n,1))
             if has_A and ig < sA1:
                 gin[:sA2, :1] = self.A[ ig, :sA2 ].T.toarray()
-                fin           = float(fin + gin.T .dot(x))
+                fin           = float(np.array(fin + gin.T .dot(x)).item())
 
             if debug:
                 print( "ig = ", ig, "  fin(linear)", fin ) #D
@@ -768,7 +768,7 @@ def evalHJv( self, mode, glist, x, v, y ):
                     iel    = self.grelt[ ig ][ iiel ]          #  the element's index
                     efname = self.elftype[ iel ];              #  the element's ftype
                     irange = [iv for iv in self.elvar[ iel ]]  #  the elemental variable's indeces
-                    xiel   = x[ np.array(irange) ]             #  the elemental variable's values
+                    xiel   = x[ np.array(irange) ].flatten().flatten()   #  the elemental variable's values
 
                     if  hasattr( self, 'grelw' ) and ig <= len( self.grelw ) and not self.grelw[ig] is None :
                         has_weights = 1;
@@ -831,7 +831,7 @@ def evalHJv( self, mode, glist, x, v, y ):
                 else:
                     fa, grada, Hessa = eval('self.'+egname+'( self, 3, fin, ig )' )
                     sgin = lil_matrix(gin);
-                    HJv  += ( ( Hessa * sgin) * (sgin.transpose().dot(v)) + grada * Hinv) / gsc
+                    HJv  += ( ( Hessa * sgin) * (sgin.transpose().dot(v).item()) + grada * Hinv) / gsc
             elif mode == "HIv":
                 if egname == "TRIVIAL":
                     ic  += 1
@@ -840,7 +840,7 @@ def evalHJv( self, mode, glist, x, v, y ):
                     fa, grada, Hessa = eval('self.'+egname+'( self, 3, fin, ig )' )
                     sgin = lil_matrix(gin);
                     ic  += 1
-                    HJv += y[ic] * ( ( Hessa * sgin) * (sgin.transpose().dot(v)) + grada * Hinv) / gsc
+                    HJv += y[ic] * ( ( Hessa * sgin) * (sgin.transpose().dot(v).item()) + grada * Hinv) / gsc
             elif  mode == "Jv":
                 ic += 1
                 if derlvl >= 1:
@@ -883,7 +883,7 @@ def evalLx( self, gobjlist, gconlist, x, y, nargout ):
             if len( gconlist ):
                 c   = self.evalgrsum( False, gconlist, x, 1 )
                 Lxy = Lxy + y.T.dot(c)
-            return float(Lxy)
+            return float(np.array(Lxy).item())
         elif nargout == 2:
             if len( gobjlist ) or hasattr( self, "H" ):
                 Lxy, Lgxy = self.evalgrsum( True, gobjlist, x, 2 )
@@ -894,7 +894,7 @@ def evalLx( self, gobjlist, gconlist, x, y, nargout ):
                 c, J = self.evalgrsum( False, gconlist, x, 2 )
                 Lxy  = Lxy + y.T.dot(c)
                 Lgxy = Lgxy  + J.T.dot(y)
-            return float(Lxy), Lgxy
+            return float(np.array(Lxy).item()), Lgxy
         elif nargout == 3:
             if len( gobjlist ) or hasattr( self, "H" ):
                 Lxy, Lgxy, LgHxy = self.evalgrsum( True, gobjlist, x, 3 )
@@ -909,7 +909,7 @@ def evalLx( self, gobjlist, gconlist, x, y, nargout ):
                 Lgxy = Lgxy  + J.T.dot(y)
                 for ig in range( len( gconlist ) ):
                     LgHxy = LgHxy + y[ig,0] * cHi[ig]
-            return float(Lxy), Lgxy.reshape(-1,1), LgHxy
+            return float(np.array(Lxy).item()), Lgxy.reshape(-1,1), LgHxy
 
     #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     #