minor improvement for speed

Author: yoelcortes

thermosteam/equilibrium/lle.py

@@ -358,11 +358,12 @@ def solve_lle_liquid_mol(self, z, T, lle_chemicals, single_loop):
             stability = lle_tangential_plane_analysis(gamma, z, T, 101325, sample=sample)
             if stability.unstable:
                 y = stability.candidate
-                y[y < 1e-32] = 1e-32
-                phi = 0.999 * (z / y).min()
-                x = z - phi * y
-                x /= x.sum()
-                K = gamma(y, T) / gamma(x, T)
+                if not stability.stable_sample:
+                    y[y < 1e-32] = 1e-32
+                    phi = 0.999 * (z / y).min()
+                    x = z - phi * y
+                    x /= x.sum()
+                    K = gamma(y, T) / gamma(x, T)
             else:
                 indices = np.argsort(z * np.array([i.MW for i in lle_chemicals]))
                 x = z.copy()

thermosteam/equilibrium/tangential_plane_stability.py

@@ -19,6 +19,7 @@ class StabilityReport(NamedTuple):
     unstable: bool
     candidate: np.ndarray
     tpd: float
+    stable_sample: bool
 
 def edge_points_simplex_masked(z: np.ndarray,
                                points_per_edge: int = 12,
@@ -103,23 +104,26 @@ def objective(w, T, P, logfz, softmax=False):
     if sample is None:
         best_val = np.inf
         best_result = None
+        stable_sample = False
     else:
         best_result = sample
         best_val = objective(sample, *args)
-        result = minimize(
-            objective, 
-            best_result, 
-            method="L-BFGS-B", 
-            options=dict(maxiter=20),
-            args=(*args, True),
-        )
-        value = result.fun
-        if value < best_val:
-            w = result.x
-            w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
-            w /= w.sum()
-            best_result = w
-            best_val = value
+        stable_sample = best_val < 0 and np.abs(best_result - z).sum() > 1e-9
+        if not stable_sample:
+            result = minimize(
+                objective, 
+                best_result, 
+                method="L-BFGS-B", 
+                options=dict(maxiter=20),
+                args=(*args, True),
+            )
+            value = result.fun
+            if value < best_val:
+                w = result.x
+                w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
+                w /= w.sum()
+                best_result = w
+                best_val = value
     w = z * MW
     w /= w.sum()
     samples = edge_points_simplex_masked(w)
@@ -148,6 +152,7 @@ def objective(w, T, P, logfz, softmax=False):
         unstable=best_val < 0 and np.abs(best_result - z).sum() > 1e-9,
         candidate=best_result,
         tpd=best_val,
+        stable_sample=stable_sample,
     )
 
 class TangentPlaneStabilityAnalysis:
@@ -176,23 +181,26 @@ def __call__(self, z, T, P, reference_phase='l', potential_phase='L', sample=Non
         if sample is None:
             best_val = np.inf
             best_result = None
+            stable_sample = False
         else:
             best_result = sample
             best_val = objective(sample, *args)
-            result = minimize(
-                objective, 
-                best_result, 
-                method="L-BFGS-B", 
-                options=dict(maxiter=20),
-                args=(*args, True),
-            )
-            value = result.fun
-            if value < best_val:
-                w = result.x
-                w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
-                w /= w.sum()
-                best_result = w
-                best_val = value
+            stable_sample = best_val < 0 and np.abs(best_result - z).sum() > 1e-9
+            if not stable_sample:
+                result = minimize(
+                    objective, 
+                    best_result, 
+                    method="L-BFGS-B", 
+                    options=dict(maxiter=20),
+                    args=(*args, True),
+                )
+                value = result.fun
+                if value < best_val:
+                    w = result.x
+                    w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
+                    w /= w.sum()
+                    best_result = w
+                    best_val = value
         MW = self.MW
         w = z * MW
         w /= w.sum()
@@ -221,7 +229,8 @@ def __call__(self, z, T, P, reference_phase='l', potential_phase='L', sample=Non
         return StabilityReport(
             unstable=best_val < 0 and np.abs(best_result - z).sum() > 1e-9,
             candidate=best_result,
-            tpd=best_val
+            tpd=best_val,
+            stable_sample=stable_sample,
         )