From 1b0b5eda767915f6fd3dca45ee11b4b45ac603ba Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Thu, 14 Nov 2024 16:10:51 -0700
Subject: [PATCH 01/12] light refactor before changes

---
 mpisppy/phbase.py            |  2 +-
 mpisppy/utils/prox_approx.py | 14 +++++++-------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py
index 9410653b..85fc8505 100644
--- a/mpisppy/phbase.py
+++ b/mpisppy/phbase.py
@@ -734,7 +734,7 @@ def attach_PH_to_objective(self, add_duals, add_prox, add_smooth=0):
                     elif self._prox_approx:
                         xvarsqrd = scenario._mpisppy_model.xsqvar[ndn_i]
                         scenario._mpisppy_data.xsqvar_prox_approx[ndn_i] = \
-                                ProxApproxManager(xvar, xvarsqrd, xbars[ndn_i], scenario._mpisppy_model.xsqvar_cuts, ndn_i)
+                                ProxApproxManager(scenario._mpisppy_model, xvar, ndn_i)
                     else:
                         xvarsqrd = xvar**2
                     prox_expr += (scenario._mpisppy_model.rho[ndn_i] / 2.0) * \
diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index c11edb5b..1b7f0eb8 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -25,11 +25,11 @@ def _newton_step(val, x_pnt, y_pnt):
 class ProxApproxManager:
     __slots__ = ()
 
-    def __new__(cls, xvar, xvarsqrd, xbar, xsqvar_cuts, ndn_i):
+    def __new__(cls, mpisppy_model, xvar, ndn_i):
         if xvar.is_integer():
-            return ProxApproxManagerDiscrete(xvar, xvarsqrd, xbar, xsqvar_cuts, ndn_i)
+            return ProxApproxManagerDiscrete(mpisppy_model, xvar, ndn_i)
         else:
-            return ProxApproxManagerContinuous(xvar, xvarsqrd, xbar, xsqvar_cuts, ndn_i)
+            return ProxApproxManagerContinuous(mpisppy_model, xvar, ndn_i)
 
 class _ProxApproxManager:
     '''
@@ -37,12 +37,12 @@ class _ProxApproxManager:
     '''
     __slots__ = ()
 
-    def __init__(self, xvar, xvarsqrd, xbar, xsqvar_cuts, ndn_i):
+    def __init__(self, mpisppy_model, xvar, ndn_i):
         self.xvar = xvar
-        self.xvarsqrd = xvarsqrd
-        self.xbar = xbar
+        self.xvarsqrd = mpisppy_model.xsqvar[ndn_i]
+        self.cuts = mpisppy_model.xsqvar_cuts
+        self.xbar = mpisppy_model.xbars[ndn_i]
         self.var_index = ndn_i
-        self.cuts = xsqvar_cuts
         self.cut_index = 0
         self._store_bounds()
 

From 71da91b94f15dff49a8a692094aafe0ac3ab1fd3 Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Fri, 15 Nov 2024 09:16:19 -0700
Subject: [PATCH 02/12] fully inline _newton_step; add W and rho

---
 mpisppy/utils/prox_approx.py | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 1b7f0eb8..24665c89 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -11,16 +11,18 @@
 from pyomo.core.expr.numeric_expr import LinearExpression
 
 # helpers for distance from y = x**2
-def _f(val, x_pnt, y_pnt):
-    return (( val - x_pnt )**2 + ( val**2 - y_pnt )**2)/2.
-def _df(val, x_pnt, y_pnt):
-    #return 2*(val - x_pnt) + 4*(val**2 - y_pnt)*val
-    return val*(1 - 2*y_pnt + 2*val*val) - x_pnt
-def _d2f(val, x_pnt, y_pnt):
-    return 1 + 6*val*val - 2*y_pnt
+# def _f(val, x_pnt, y_pnt):
+#     return (( val - x_pnt )**2 + ( val**2 - y_pnt )**2)/2.
+# def _df(val, x_pnt, y_pnt):
+#     #return 2*(val - x_pnt) + 4*(val**2 - y_pnt)*val
+#     return val*(1 - 2*y_pnt + 2*val*val) - x_pnt
+# def _d2f(val, x_pnt, y_pnt):
+#     return 1 + 6*val*val - 2*y_pnt
+# def _newton_step(val, x_pnt, y_pnt):
+#     return val - _df(val, x_pnt, y_pnt) / _d2f(val, x_pnt, y_pnt)
 
 def _newton_step(val, x_pnt, y_pnt):
-    return val - _df(val, x_pnt, y_pnt) / _d2f(val, x_pnt, y_pnt)
+    return val - (val * (1 - 2*y_pnt + 2*val*val) - x_pnt) / (1 + 6*val*val - 2*y_pnt)
 
 class ProxApproxManager:
     __slots__ = ()
@@ -42,6 +44,8 @@ def __init__(self, mpisppy_model, xvar, ndn_i):
         self.xvarsqrd = mpisppy_model.xsqvar[ndn_i]
         self.cuts = mpisppy_model.xsqvar_cuts
         self.xbar = mpisppy_model.xbars[ndn_i]
+        self.rho = mpisppy_model.rho[ndn_i]
+        self.W = mpisppy_model.W[ndn_i]
         self.var_index = ndn_i
         self.cut_index = 0
         self._store_bounds()

From a77f3badb17a7b94916b894a9a03141309c2d0f9 Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Fri, 15 Nov 2024 14:04:33 -0700
Subject: [PATCH 03/12] setting things up for checking tolerance on the x-axis

---
 mpisppy/phbase.py            |  5 ++++
 mpisppy/utils/prox_approx.py | 56 ++++++++++++++++++++++--------------
 mpisppy/utils/sputils.py     |  3 +-
 3 files changed, 42 insertions(+), 22 deletions(-)

diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py
index 85fc8505..27751e62 100644
--- a/mpisppy/phbase.py
+++ b/mpisppy/phbase.py
@@ -8,6 +8,7 @@
 ###############################################################################
 import time
 import logging
+import math
 
 import numpy as np
 import mpisppy.MPI as MPI
@@ -689,6 +690,10 @@ def attach_PH_to_objective(self, add_duals, add_prox, add_smooth=0):
                 self.prox_approx_tol = self.options['proximal_linearization_tolerance']
             else:
                 self.prox_approx_tol = 1.e-1
+            # The proximal approximation code now checks the tolerance based on the x-coordinates
+            # as opposed to the y-coordinates. Therefore, we will use the square root of the
+            # y-coordinate tolerance.
+            self.prox_approx_tol = math.sqrt(self.prox_approx_tol)
         else:
             self._prox_approx = False
 
diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 24665c89..2aa5033a 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -60,29 +60,44 @@ def _store_bounds(self):
         else:
             self.ub = self.xvar.ub
 
-    def add_cut(self, val, persistent_solver=None):
+    def add_cut(self, val, tolerance, persistent_solver):
         '''
         create a cut at val
         '''
         pass
 
+    def add_cuts(self, x_val, tolerance, persistent_solver):
+        x_bar = self.xbar.value
+        rho = self.rho.value
+        W = self.W.value
+
+        self.add_cut(x_val, tolerance, persistent_solver)
+        # rotate x_val around x_bar, the minimizer of (\rho / 2)(x - x_bar)^2
+        # to create a vertex at this point
+        rotated_x_val_x_bar = 2*x_bar - x_val
+        if not isclose(x_val, rotated_x_val_x_bar, abs_tol=tolerance):
+            self.add_cut(rotated_x_val_x_bar, tolerance, persistent_solver)
+        # aug_lagrange_point, is the minimizer of w\cdot x + (\rho / 2)(x - x_bar)^2
+        # to create a vertex at this point
+        aug_lagrange_point = -W / rho + x_bar
+        if not isclose(x_val, aug_lagrange_point, abs_tol=tolerance):
+            self.add_cut(2*aug_lagrange_point - x_val, tolerance, persistent_solver)
+        # finally, create another vertex at the aug_lagrange_point by rotating
+        # rotated_x_val_x_bar around the aug_lagrange_point
+        if not isclose(rotated_x_val_x_bar, aug_lagrange_point, abs_tol=tolerance):
+            self.add_cut(2*aug_lagrange_point - rotated_x_val_x_bar, tolerance, persistent_solver)
+        return True
+
     def check_tol_add_cut(self, tolerance, persistent_solver=None):
         '''
         add a cut if the tolerance is not satified
         '''
         x_pnt = self.xvar.value
-        x_bar = self.xbar.value
         y_pnt = self.xvarsqrd.value
         f_val = x_pnt**2
 
         #print(f"y-distance: {actual_val - measured_val})")
-        if y_pnt is None:
-            self.add_cut(x_pnt, persistent_solver)
-            if not isclose(x_pnt, x_bar, abs_tol=1e-6):
-                self.add_cut(2*x_bar - x_pnt, persistent_solver)
-            return True
-
-        if (f_val - y_pnt) > tolerance:
+        if y_pnt is not None:
             '''
             In this case, we project the point x_pnt, y_pnt onto
             the curve y = x**2 by finding the minimum distance
@@ -96,20 +111,17 @@ def check_tol_add_cut(self, tolerance, persistent_solver=None):
             #print(f"initial distance: {_f(this_val, x_pnt, y_pnt)**(0.5)}")
             #print(f"this_val: {this_val}")
             next_val = _newton_step(this_val, x_pnt, y_pnt)
-            while not isclose(this_val, next_val, rel_tol=1e-6, abs_tol=1e-6):
-                #print(f"newton step distance: {_f(next_val, x_pnt, y_pnt)**(0.5)}")
-                #print(f"next_val: {next_val}")
-                this_val = next_val
-                next_val = _newton_step(this_val, x_pnt, y_pnt)
-            self.add_cut(next_val, persistent_solver)
-            if not isclose(next_val, x_bar, abs_tol=1e-6):
-                self.add_cut(2*x_bar - next_val, persistent_solver)
-            return True
-        return False
+            # while not isclose(this_val, next_val, rel_tol=1e-6, abs_tol=1e-6):
+            #     #print(f"newton step distance: {_f(next_val, x_pnt, y_pnt)**(0.5)}")
+            #     #print(f"next_val: {next_val}")
+            #     this_val = next_val
+            #     next_val = _newton_step(this_val, x_pnt, y_pnt)
+            x_pnt = next_val
+        return self.add_cuts(x_pnt, tolerance,  persistent_solver)
 
 class ProxApproxManagerContinuous(_ProxApproxManager):
 
-    def add_cut(self, val, persistent_solver=None):
+    def add_cut(self, val, tolerance, persistent_solver):
         '''
         create a cut at val using a taylor approximation
         '''
@@ -149,10 +161,12 @@ def _compute_mb(val):
 
 class ProxApproxManagerDiscrete(_ProxApproxManager):
 
-    def add_cut(self, val, persistent_solver=None):
+    def add_cut(self, val, tolerance, persistent_solver):
         '''
         create up to two cuts at val, exploiting integrality
         '''
+        # TODO: tolerance isn't realy used for discrete,
+        # but as long as the tolerance is <=1 it doesn't matter
         val = int(round(val))
 
         ## cuts are indexed by the x-value to the right
diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py
index df9f4b26..bb411a46 100644
--- a/mpisppy/utils/sputils.py
+++ b/mpisppy/utils/sputils.py
@@ -984,7 +984,8 @@ def nonant_cost_coeffs(s):
         if id(var) in s._mpisppy_data.varid_to_nonant_index:
             raise RuntimeError(
                 "Found nonlinear variables in the objective function. "
-                f"Variable {var} has nonlinear interactions in the objective funtion"
+                f"Variable {var} has nonlinear interactions in the objective funtion. "
+                "Consider using gradient-based rho."
             )
     return cost_coefs
 

From 6cabaa5564f80175f63efff7dd51cc7081c7e2d6 Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Fri, 15 Nov 2024 15:30:20 -0700
Subject: [PATCH 04/12] checking in new method to avoid redundant cuts

---
 mpisppy/utils/prox_approx.py | 50 +++++++++++++++++++++++++++++-------
 1 file changed, 41 insertions(+), 9 deletions(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 2aa5033a..70c0b397 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -8,6 +8,8 @@
 ###############################################################################
 
 from math import isclose
+from array import array
+from bisect import bisect
 from pyomo.core.expr.numeric_expr import LinearExpression
 
 # helpers for distance from y = x**2
@@ -48,6 +50,8 @@ def __init__(self, mpisppy_model, xvar, ndn_i):
         self.W = mpisppy_model.W[ndn_i]
         self.var_index = ndn_i
         self.cut_index = 0
+        self.cut_values = array("d")
+        self.cut_values.append(0.0)
         self._store_bounds()
 
     def _store_bounds(self):
@@ -66,27 +70,53 @@ def add_cut(self, val, tolerance, persistent_solver):
         '''
         pass
 
+    def check_and_add_value(self, val, tolerance):
+        idx = bisect(self.cut_values, val)
+        # self.cut_values is empty, has one element
+        # or we're appending to the end
+        if idx == len(self.cut_values):
+            if val - self.cut_values[idx-1] < tolerance:
+                return False
+            else:
+                self.cut_values.insert(idx, val)
+                return True
+        # here we're at the beginning
+        if idx == 0:
+            if self.cut_values[idx] - val < tolerance:
+                return False
+            else:
+                self.cut_values.insert(idx, val)
+                return True
+        # in the middle
+        if self.cut_values[idx] - val < tolerance:
+            return False
+        if val - self.cut_values[idx-1] < tolerance:
+            return False
+        self.cut_values.insert(idx, val)
+        return True
+
     def add_cuts(self, x_val, tolerance, persistent_solver):
         x_bar = self.xbar.value
         rho = self.rho.value
         W = self.W.value
 
-        self.add_cut(x_val, tolerance, persistent_solver)
+        num_cuts = self.add_cut(x_val, tolerance, persistent_solver)
         # rotate x_val around x_bar, the minimizer of (\rho / 2)(x - x_bar)^2
         # to create a vertex at this point
         rotated_x_val_x_bar = 2*x_bar - x_val
         if not isclose(x_val, rotated_x_val_x_bar, abs_tol=tolerance):
-            self.add_cut(rotated_x_val_x_bar, tolerance, persistent_solver)
+            num_cuts += self.add_cut(rotated_x_val_x_bar, tolerance, persistent_solver)
         # aug_lagrange_point, is the minimizer of w\cdot x + (\rho / 2)(x - x_bar)^2
         # to create a vertex at this point
         aug_lagrange_point = -W / rho + x_bar
         if not isclose(x_val, aug_lagrange_point, abs_tol=tolerance):
-            self.add_cut(2*aug_lagrange_point - x_val, tolerance, persistent_solver)
+            num_cuts += self.add_cut(2*aug_lagrange_point - x_val, tolerance, persistent_solver)
         # finally, create another vertex at the aug_lagrange_point by rotating
         # rotated_x_val_x_bar around the aug_lagrange_point
         if not isclose(rotated_x_val_x_bar, aug_lagrange_point, abs_tol=tolerance):
-            self.add_cut(2*aug_lagrange_point - rotated_x_val_x_bar, tolerance, persistent_solver)
-        return True
+            num_cuts += self.add_cut(2*aug_lagrange_point - rotated_x_val_x_bar, tolerance, persistent_solver)
+        # print(f"{len(self.cut_values)=}")
+        return num_cuts
 
     def check_tol_add_cut(self, tolerance, persistent_solver=None):
         '''
@@ -125,8 +155,7 @@ def add_cut(self, val, tolerance, persistent_solver):
         '''
         create a cut at val using a taylor approximation
         '''
-        # handled by bound
-        if val == 0:
+        if not self.check_and_add_value(val, tolerance):
             return 0
         # f'(a) = 2*val
         # f(a) - f'(a)a = val*val - 2*val*val
@@ -165,9 +194,12 @@ def add_cut(self, val, tolerance, persistent_solver):
         '''
         create up to two cuts at val, exploiting integrality
         '''
-        # TODO: tolerance isn't realy used for discrete,
-        # but as long as the tolerance is <=1 it doesn't matter
         val = int(round(val))
+        if tolerance > 1:
+            # TODO: We should consider how to handle this, maybe.
+            #       Tolerances less than or equal 1 won't affect
+            #       discrete cuts
+            pass
 
         ## cuts are indexed by the x-value to the right
         ## e.g., the cut for (2,3) is indexed by 3

From ad6bda4f0b1d32413f43d91cd05a359f6428057a Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Fri, 15 Nov 2024 15:34:56 -0700
Subject: [PATCH 05/12] add debugging statement

---
 mpisppy/utils/prox_approx.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 70c0b397..643a4a9a 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -115,7 +115,7 @@ def add_cuts(self, x_val, tolerance, persistent_solver):
         # rotated_x_val_x_bar around the aug_lagrange_point
         if not isclose(rotated_x_val_x_bar, aug_lagrange_point, abs_tol=tolerance):
             num_cuts += self.add_cut(2*aug_lagrange_point - rotated_x_val_x_bar, tolerance, persistent_solver)
-        # print(f"{len(self.cut_values)=}")
+        # print(f"{self.cut_index=}")
         return num_cuts
 
     def check_tol_add_cut(self, tolerance, persistent_solver=None):

From 981ce8fe88b8c5060763727d05219459da5bdeab Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Mon, 18 Nov 2024 10:07:45 -0700
Subject: [PATCH 06/12] more tweaks

---
 mpisppy/utils/prox_approx.py | 43 ++++++++++++++++++------------------
 1 file changed, 22 insertions(+), 21 deletions(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 643a4a9a..16b0fc0f 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -115,7 +115,8 @@ def add_cuts(self, x_val, tolerance, persistent_solver):
         # rotated_x_val_x_bar around the aug_lagrange_point
         if not isclose(rotated_x_val_x_bar, aug_lagrange_point, abs_tol=tolerance):
             num_cuts += self.add_cut(2*aug_lagrange_point - rotated_x_val_x_bar, tolerance, persistent_solver)
-        # print(f"{self.cut_index=}")
+        # print(f"{x_val=}, {x_bar=}, {W=}")
+        # print(f"{self.cut_values=}")
         return num_cuts
 
     def check_tol_add_cut(self, tolerance, persistent_solver=None):
@@ -126,27 +127,27 @@ def check_tol_add_cut(self, tolerance, persistent_solver=None):
         y_pnt = self.xvarsqrd.value
         f_val = x_pnt**2
 
-        #print(f"y-distance: {actual_val - measured_val})")
-        if y_pnt is not None:
-            '''
-            In this case, we project the point x_pnt, y_pnt onto
-            the curve y = x**2 by finding the minimum distance
-            between y = x**2 and x_pnt, y_pnt.
-
-            This involves solving a cubic equation, so instead
-            we start at x_pnt, y_pnt and run newtons algorithm
-            to get an approximate good-enough solution.
-            '''
-            this_val = x_pnt
-            #print(f"initial distance: {_f(this_val, x_pnt, y_pnt)**(0.5)}")
-            #print(f"this_val: {this_val}")
+        # print(f"{x_pnt=}, {y_pnt=}, {f_val=}")
+        # print(f"y-distance: {actual_val - measured_val})")
+        if y_pnt is None:
+            y_pnt = 0.0
+        # In this case, we project the point x_pnt, y_pnt onto
+        # the curve y = x**2 by finding the minimum distance
+        # between y = x**2 and x_pnt, y_pnt.
+
+        # This involves solving a cubic equation, so instead
+        # we start at x_pnt, y_pnt and run newtons algorithm
+        # to get an approximate good-enough solution.
+        this_val = x_pnt
+        # print(f"initial distance: {_f(this_val, x_pnt, y_pnt)**(0.5)}")
+        # print(f"this_val: {this_val}")
+        next_val = _newton_step(this_val, x_pnt, y_pnt)
+        while not isclose(this_val, next_val, rel_tol=1e-6, abs_tol=1e-6):
+            # print(f"newton step distance: {_f(next_val, x_pnt, y_pnt)**(0.5)}")
+            # print(f"next_val: {next_val}")
+            this_val = next_val
             next_val = _newton_step(this_val, x_pnt, y_pnt)
-            # while not isclose(this_val, next_val, rel_tol=1e-6, abs_tol=1e-6):
-            #     #print(f"newton step distance: {_f(next_val, x_pnt, y_pnt)**(0.5)}")
-            #     #print(f"next_val: {next_val}")
-            #     this_val = next_val
-            #     next_val = _newton_step(this_val, x_pnt, y_pnt)
-            x_pnt = next_val
+        x_pnt = next_val
         return self.add_cuts(x_pnt, tolerance,  persistent_solver)
 
 class ProxApproxManagerContinuous(_ProxApproxManager):

From aa8dfd8c3d26d44ae690ed346c6b7a7b1e23ca5d Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Mon, 18 Nov 2024 13:16:53 -0700
Subject: [PATCH 07/12] ignore fixed nonants; project onto bounds; add safety
 cuts

---
 mpisppy/utils/prox_approx.py | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 16b0fc0f..926157b7 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -115,26 +115,35 @@ def add_cuts(self, x_val, tolerance, persistent_solver):
         # rotated_x_val_x_bar around the aug_lagrange_point
         if not isclose(rotated_x_val_x_bar, aug_lagrange_point, abs_tol=tolerance):
             num_cuts += self.add_cut(2*aug_lagrange_point - rotated_x_val_x_bar, tolerance, persistent_solver)
+        # If we only added 0 or 1 cuts initially, add up to two more
+        # to capture something of the proximal term. This can happen
+        # when x_bar == x_val and W == 0.
+        if self.cut_index <= 1:
+            num_cuts += self.add_cut(x_val + max(1, tolerance+1e-06), tolerance, persistent_solver)
+            num_cuts += self.add_cut(x_val - max(1, tolerance-1e-06), tolerance, persistent_solver)
         # print(f"{x_val=}, {x_bar=}, {W=}")
         # print(f"{self.cut_values=}")
+        # print(f"{self.cut_index=}")
         return num_cuts
 
     def check_tol_add_cut(self, tolerance, persistent_solver=None):
         '''
         add a cut if the tolerance is not satified
         '''
+        if self.xvar.fixed:
+            # don't do anything for fixed variables
+            return 0
         x_pnt = self.xvar.value
         y_pnt = self.xvarsqrd.value
-        f_val = x_pnt**2
+        # f_val = x_pnt**2
 
         # print(f"{x_pnt=}, {y_pnt=}, {f_val=}")
         # print(f"y-distance: {actual_val - measured_val})")
         if y_pnt is None:
             y_pnt = 0.0
-        # In this case, we project the point x_pnt, y_pnt onto
+        # We project the point x_pnt, y_pnt onto
         # the curve y = x**2 by finding the minimum distance
         # between y = x**2 and x_pnt, y_pnt.
-
         # This involves solving a cubic equation, so instead
         # we start at x_pnt, y_pnt and run newtons algorithm
         # to get an approximate good-enough solution.
@@ -156,6 +165,11 @@ def add_cut(self, val, tolerance, persistent_solver):
         '''
         create a cut at val using a taylor approximation
         '''
+        lb, ub = self.xvar.bounds
+        if lb is not None and val < lb:
+            val = lb
+        if ub is not None and val > ub:
+            val = ub
         if not self.check_and_add_value(val, tolerance):
             return 0
         # f'(a) = 2*val

From 525cd2791f7021a5ce53a2f93ae3489ff84defb8 Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Mon, 18 Nov 2024 13:18:36 -0700
Subject: [PATCH 08/12] fix double negation mistake

---
 mpisppy/utils/prox_approx.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 926157b7..5ea93708 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -120,7 +120,7 @@ def add_cuts(self, x_val, tolerance, persistent_solver):
         # when x_bar == x_val and W == 0.
         if self.cut_index <= 1:
             num_cuts += self.add_cut(x_val + max(1, tolerance+1e-06), tolerance, persistent_solver)
-            num_cuts += self.add_cut(x_val - max(1, tolerance-1e-06), tolerance, persistent_solver)
+            num_cuts += self.add_cut(x_val - max(1, tolerance+1e-06), tolerance, persistent_solver)
         # print(f"{x_val=}, {x_bar=}, {W=}")
         # print(f"{self.cut_values=}")
         # print(f"{self.cut_index=}")

From ca2bf21dd58bea2300bdc010deaf694df3ace567 Mon Sep 17 00:00:00 2001
From: Bernard Knueven <Bernard.Knueven@nrel.gov>
Date: Mon, 18 Nov 2024 16:15:41 -0700
Subject: [PATCH 09/12] better prox approx with no information

---
 mpisppy/utils/prox_approx.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py
index 5ea93708..03c21364 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -119,8 +119,17 @@ def add_cuts(self, x_val, tolerance, persistent_solver):
         # to capture something of the proximal term. This can happen
         # when x_bar == x_val and W == 0.
         if self.cut_index <= 1:
-            num_cuts += self.add_cut(x_val + max(1, tolerance+1e-06), tolerance, persistent_solver)
-            num_cuts += self.add_cut(x_val - max(1, tolerance+1e-06), tolerance, persistent_solver)
+            lb, ub = self.xvar.bounds
+            upval = 1
+            dnval = 1
+            if ub is not None:
+                # after a lot of calculus, you can show that
+                # this point minimizes the error in the approximation
+                upval = 2.0 * (ub - x_val) / 3.0
+            if lb is not None:
+                dnval = 2.0 * (x_val - lb) / 3.0
+            num_cuts += self.add_cut(x_val + max(upval, tolerance+1e-06), tolerance, persistent_solver)
+            num_cuts += self.add_cut(x_val - max(dnval, tolerance+1e-06), tolerance, persistent_solver)
         # print(f"{x_val=}, {x_bar=}, {W=}")
         # print(f"{self.cut_values=}")
         # print(f"{self.cut_index=}")

From a4d6e8d46db32b931d2d6d8fbe129c3d4172bcce Mon Sep 17 00:00:00 2001
From: Maximilian Fey <38524131+maxfey@users.noreply.github.com>
Date: Mon, 2 Dec 2024 22:19:23 +0100
Subject: [PATCH 10/12] Correct `ExtensiveForm` docstring arg order and add
 missing parameter (#465)

- Update the docstring of the `ExtensiveForm` class to match the order of
  arguments in the `__init__` method signature.
- Add the missing `all_nodenames` argument to the docstring, using a
  description consistent with the parent class `SPBase` but styled to
  match `ExtensiveForm`.

Co-authored-by: Maximilian Fey <maximilian@fey-engineering.de>
Co-authored-by: David L Woodruff <DLWoodruff@users.noreply.github.com>
---
 mpisppy/opt/ef.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/mpisppy/opt/ef.py b/mpisppy/opt/ef.py
index 6c647d81..ea6dbf9d 100644
--- a/mpisppy/opt/ef.py
+++ b/mpisppy/opt/ef.py
@@ -31,10 +31,13 @@ class ExtensiveForm(mpisppy.spbase.SPBase):
         scenario_creator (callable):
             Scenario creator function, which takes as input a scenario
             name, and returns a Pyomo model of that scenario.
+        scenario_creator_kwargs (dict, optional):
+            Keyword arguments passed to `scenario_creator`.
+        all_nodenames (list, optional):
+            List of all node names, incl. leaves. Can be None for two-stage
+            problem.
         model_name (str, optional):
             Name of the resulting EF model object.
-        scenario_creator_kwargs (dict):
-            Keyword args passed to `scenario_creator`.
         suppress_warnings (bool, optional):
             Boolean to suppress warnings when building the EF. Default
             is False.

From 0a8026c9cab52282b7cfcd8eb33cfd4b5aeb3924 Mon Sep 17 00:00:00 2001
From: Maximilian Fey <38524131+maxfey@users.noreply.github.com>
Date: Mon, 2 Dec 2024 22:40:43 +0100
Subject: [PATCH 11/12] Fix typos in documentation (#466)

Co-authored-by: Maximilian Fey <maximilian@fey-engineering.de>
---
 doc/src/aph.rst                  |  6 +++---
 doc/src/confidence_intervals.rst | 12 ++++++------
 doc/src/drivers.rst              |  4 ++--
 doc/src/examples.rst             |  4 ++--
 doc/src/install_mpi.rst          |  2 +-
 doc/src/scenario_creator.rst     |  8 ++++----
 6 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/doc/src/aph.rst b/doc/src/aph.rst
index c89c25a8..eee98753 100644
--- a/doc/src/aph.rst
+++ b/doc/src/aph.rst
@@ -4,7 +4,7 @@ APH
 ===
 
 The code is based on "Algorithm 2: Asynchronous projective hedging
-(APH) -- Algorithm 1 specialize to the setup S1-S4" from "Asynchronous
+(APH) -- Algorithm 1 specialized to the setup S1-S4" from "Asynchronous
 Projective Hedging for Stochastic Programming"
 http://www.optimization-online.org/DB_HTML/2018/10/6895.html
 
@@ -100,7 +100,7 @@ farmer
 ^^^^^^
 
 The scripts for this example are currently in the paper repo in
-`AsyncPH/experiments/challange/farmer`; the driver is
+`AsyncPH/experiments/challenge/farmer`; the driver is
 `farmer_driver.py`.  The driver references the model, which is in the
 `mpi-sppy` repo.  The `aph05.bash` script is intended
 to have a dispatch fraction of 1/2 (hence the 05 for 0.5 in the name).
@@ -119,7 +119,7 @@ A Peek Under the Hood
 
 The APH implementation has a listener thread that continuously does
 MPI reductions and a worker thread that does most of the work. A wrinkle
-is that the listenter thread does a `side gig` if enough ranks have reported
+is that the listener thread does a `side gig` if enough ranks have reported
 (the "async_frac_needed" option) because after it has
 done the reductions to get u and v and needs to do some calculations, and
 then reductions to compute tau and theta.
diff --git a/doc/src/confidence_intervals.rst b/doc/src/confidence_intervals.rst
index 97e6760f..37b7b516 100644
--- a/doc/src/confidence_intervals.rst
+++ b/doc/src/confidence_intervals.rst
@@ -7,7 +7,7 @@ If we want to assess the quality of a given candidate solution ``xhat_one``
 (a first stage solution), we could try and evaluate the optimality gap, i.e. 
 the gap between the value of the objective function
 at ``xhat_one`` and the value of the solution to the problem.
-The class ``MMWConfidenceIntervals`` compute an estimator of the optimality gap
+The class ``MMWConfidenceIntervals`` computes an estimator of the optimality gap
 as described in [mmw1999]_ (Section 3.2) and an asymptotic confidence interval for
 this gap. 
 
@@ -40,7 +40,7 @@ Evaluating a candidate solution
 To evaluate a candidate solution with some scenarios, one might
 create a ``Xhat_Eval`` object and call its ``evaluate`` method 
 (resp. ``evaluate_one`` for a single scenario). It takes as
-an argument ``xhats``, a dictionary of noon-anticipative policies for all 
+an argument ``xhats``, a dictionary of non-anticipative policies for all 
 non-leaf nodes of a scenario tree. While for a 2-stage problem, ``xhats`` is
 just the candidate solution ``xhat_one``, for multistage problem the 
 dictionary can be computed using the function ``walking_tree_xhats`` 
@@ -57,16 +57,16 @@ This object has a ``run`` method that returns a gap estimator and a confidence i
 Examples
 --------
 
-There are example scrips for sequential sampling in both ``farmer`` and ``aircond``.
+There are example scripts for sequential sampling in both ``farmer`` and ``aircond``.
 
-Using stand alone ``mmw_conf.py``
+Using stand-alone ``mmw_conf.py``
 ---------------------------------
 
-(The stand-alone module is currently for use with 2-stage problem only; for multi-stage problems, instantiate an ``MMWConfidenceItervals`` object directly)
+(The stand-alone module is currently for use with 2-stage problem only; for multi-stage problems, instantiate an ``MMWConfidenceIntervals`` object directly)
 
 ``mmw_conf`` uses the ``MMWConfidenceIntervals`` class from ``mmw_ci`` in order to construct a confidence interval on the optimality gap of a particular candidate solution ``xhat`` of a model instance. 
 
-To use the stand along program a model compatible with ``Amalgamator`` and ``.npy`` file with a candidate solution to an instance of the model are required.
+To use the stand-alone program a model compatible with ``Amalgamator`` and ``.npy`` file with a candidate solution to an instance of the model are required.
 
 First, ensure that the model to be used is compatible with the
 ``Amalgamator`` class. This requires the model to have each of the
diff --git a/doc/src/drivers.rst b/doc/src/drivers.rst
index 08d1f7bc..6615fa67 100644
--- a/doc/src/drivers.rst
+++ b/doc/src/drivers.rst
@@ -40,11 +40,11 @@ Extending Examples
 Many developers
 will need to add extensions. Here are few examples:
 
-* In the ``farmer_cylinders.py`` example, there is a block of code to add a ``--crops-mult`` argument that is passed to the scenario create in the ``scenario_creator_kwargs`` dictionary.
+* In the ``farmer_cylinders.py`` example, there is a block of code to add a ``--crops-mult`` argument that is passed to the scenario creator in the ``scenario_creator_kwargs`` dictionary.
 
 * In the ``hydro_cylinders.py`` example (which has three stages). The branching factors are obtained from the command line and passed to the scenario constructor via ``scenario_creator_kwargs`` and also passed to ``sputils.create_nodenames_from_BFs`` to create a node list.
 
-* The ``uc_cylinders.py`` example adds arguments that are used to provide data or trigger the inclusion of extensions. The  extension specifications and arguments are added to the dictionaries  (e.g., ``hub_dict``) create by ``vanilla.py``.
+* The ``uc_cylinders.py`` example adds arguments that are used to provide data or trigger the inclusion of extensions. The extension specifications and arguments are added to the dictionaries (e.g., ``hub_dict``) created by ``vanilla.py``.
 
 Not Using Examples Utilities
 ----------------------------
diff --git a/doc/src/examples.rst b/doc/src/examples.rst
index cd4733f4..bdaa7619 100644
--- a/doc/src/examples.rst
+++ b/doc/src/examples.rst
@@ -3,11 +3,11 @@
 Examples
 ========
 
-If you installed directly from github, the top
+If you installed directly from GitHub, the top
 level directory  ``examples`` 
 contains some sub-directories with examples.
 
-If you did not get the code from github (e.g., if
+If you did not get the code from GitHub (e.g., if
 you installed with pip), you will need to
 get the examples from:
 https://github.com/Pyomo/mpi-sppy/tree/master/examples
diff --git a/doc/src/install_mpi.rst b/doc/src/install_mpi.rst
index 05ebb429..964b9fa0 100644
--- a/doc/src/install_mpi.rst
+++ b/doc/src/install_mpi.rst
@@ -12,7 +12,7 @@ Here are two methods that seem to work well for installation, at least when cons
 
    * ``conda install openmpi; conda install mpi4py``  (in that order)
   
-#. If you already have an existing version of MPI, it may be better compile mpi4py against it. This can be done by installing mpi4py though pip.
+#. If you already have an existing version of MPI, it may be better compile mpi4py against it. This can be done by installing mpi4py through pip.
 
    * ``pip install mpi4py``
 
diff --git a/doc/src/scenario_creator.rst b/doc/src/scenario_creator.rst
index 0198c1f2..fecbb9cc 100644
--- a/doc/src/scenario_creator.rst
+++ b/doc/src/scenario_creator.rst
@@ -5,7 +5,7 @@
 
 This function instantiates models for scenarios and usually attaches
 some information about the scenario tree. It is required, but can have
-any name, and its first argment must be the scenario name. The other
+any name, and its first argument must be the scenario name. The other
 two arguments are optional. The `scenario_creator_kwargs` option specifies data that is
 passed through from the calling program.
 `scenario_creator_kwargs` must be a dictionary, and might give, e.g., a data
@@ -51,7 +51,7 @@ illustrated in the netdes example:
 
 ::
    
-   # Add all index of model.x using wild cards
+   # Add all indexes of model.x using wild cards
    sputils.attach_root_node(model, model.FirstStageCost, [model.x[:,:], ])
 
 Scenario Probability
@@ -70,13 +70,13 @@ The function ``attach_root_node`` takes an optional argument ``nonant_ef_suppl_l
 multipliers by algorithms such as PH, but will be given non-anticipativity
 constraints when an EF is formed, either to solve the EF or when bundles are
 formed. For some problems, with the appropriate solver, adding redundant nonanticipativity constraints
-for auxilliary variables to the bundle/EF will result in a (much) smaller pre-solved model.
+for auxiliary variables to the bundle/EF will result in a (much) smaller pre-solved model.
 
 Multi-stage
 -----------
 
 When there are scenario tree nodes other than root, their names,
-although strings, must indicates their position in the tree, 
+although strings, must indicate their position in the tree, 
 like "ROOT_3_0_1". A given non-root node, which is the child number `k` of
 a node with name `parentname`, should be named `parentname_k`.
 The node constructor, ``scenario_tree.ScenarioNode`` takes as

From b3ee83a93c6af927ebdeeb54b427eb2360687034 Mon Sep 17 00:00:00 2001
From: bknueven <30801372+bknueven@users.noreply.github.com>
Date: Tue, 3 Dec 2024 09:20:44 -0700
Subject: [PATCH 12/12] More helpful error message

---
 mpisppy/utils/sputils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py
index bb411a46..5118bf65 100644
--- a/mpisppy/utils/sputils.py
+++ b/mpisppy/utils/sputils.py
@@ -983,7 +983,7 @@ def nonant_cost_coeffs(s):
     for var in repn.nonlinear_vars:
         if id(var) in s._mpisppy_data.varid_to_nonant_index:
             raise RuntimeError(
-                "Found nonlinear variables in the objective function. "
+                "A call to nonant_cost_coefficient found nonlinear variables in the objective function. "
                 f"Variable {var} has nonlinear interactions in the objective funtion. "
                 "Consider using gradient-based rho."
             )