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
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.
diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py
index 9410653b..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
 
@@ -734,7 +739,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..03c21364 100644
--- a/mpisppy/utils/prox_approx.py
+++ b/mpisppy/utils/prox_approx.py
@@ -8,28 +8,32 @@
 ###############################################################################
 
 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
-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__ = ()
 
-    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,13 +41,17 @@ 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.rho = mpisppy_model.rho[ndn_i]
+        self.W = mpisppy_model.W[ndn_i]
         self.var_index = ndn_i
-        self.cuts = xsqvar_cuts
         self.cut_index = 0
+        self.cut_values = array("d")
+        self.cut_values.append(0.0)
         self._store_bounds()
 
     def _store_bounds(self):
@@ -56,61 +64,122 @@ 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 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
+
+        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):
+            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):
+            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):
+            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:
+            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=}")
+        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
-        x_bar = self.xbar.value
         y_pnt = self.xvarsqrd.value
-        f_val = x_pnt**2
+        # f_val = x_pnt**2
 
-        #print(f"y-distance: {actual_val - measured_val})")
+        # print(f"{x_pnt=}, {y_pnt=}, {f_val=}")
+        # 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:
-            '''
-            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}")
+            y_pnt = 0.0
+        # 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)
-            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
+        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
         '''
-        # handled by bound
-        if val == 0:
+        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
         # f(a) - f'(a)a = val*val - 2*val*val
@@ -145,11 +214,16 @@ 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
         '''
         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
diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py
index df9f4b26..5118bf65 100644
--- a/mpisppy/utils/sputils.py
+++ b/mpisppy/utils/sputils.py
@@ -983,8 +983,9 @@ 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. "
-                f"Variable {var} has nonlinear interactions in the objective funtion"
+                "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."
             )
     return cost_coefs