Repository: fides-dev/fides (language: python)

Issue:
Improve multithreading compatibility
Looks like having a central module specific logger instance leads to excessive lock times even when logging is not displayed.


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/flake.yaml b/.github/workflows/flake.yaml
index 8ab3e2b..44c9ec1 100644
--- a/.github/workflows/flake.yaml
+++ b/.github/workflows/flake.yaml
@@ -9,7 +9,7 @@ jobs:
       - name: Setup Python
         uses: actions/setup-python@v1
         with:
-          python-version: 3.8.6
+          python-version: 3.9.1
           architecture: x64
 
       - uses: actions/checkout@v1
diff --git a/fides/__init__.py b/fides/__init__.py
index f2ce1ce..3c498bf 100644
--- a/fides/__init__.py
+++ b/fides/__init__.py
@@ -7,6 +7,6 @@
 # flake8: noqa
 from .minimize import Optimizer
 from .hessian_approximation import SR1, BFGS, DFP, HessianApproximation
-from .logging import logger
+from .logging import create_logger
 from .version import __version__
 from .constants import Options, SubSpaceDim, StepBackStrategy, ExitFlag
diff --git a/fides/logging.py b/fides/logging.py
index 9d27c85..4d3eeaf 100644
--- a/fides/logging.py
+++ b/fides/logging.py
@@ -10,10 +10,30 @@
 
 import logging
 
-logger = logging.getLogger('fides')
-ch = logging.StreamHandler()
-formatter = logging.Formatter(
-    '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
-)
-ch.setFormatter(formatter)
-logger.addHandler(ch)
+logger_count = 0
+
+
+def create_logger(level: int) -> logging.Logger:
+    """
+    Creates a logger instance. To avoid unnecessary locks during
+    multithreading, different logger instance should be created for every
+
+    :param level:
+        logging level
+
+    :return:
+        logger instance
+    """
+    global logger_count
+    logger_count += 1
+    # add logger count to differentiate between different fides
+    # optimization instances and avoid deadlocks
+    logger = logging.getLogger(f'fides_{logger_count}')
+    ch = logging.StreamHandler()
+    formatter = logging.Formatter(
+        '%(asctime)s - fides - %(levelname)s - %(message)s'
+    )
+    ch.setFormatter(formatter)
+    logger.addHandler(ch)
+    logger.setLevel(level)
+    return logger
diff --git a/fides/minimize.py b/fides/minimize.py
index 46a425c..f1386cc 100644
--- a/fides/minimize.py
+++ b/fides/minimize.py
@@ -13,9 +13,9 @@
 from .trust_region import trust_region, Step
 from .hessian_approximation import HessianApproximation
 from .constants import Options, ExitFlag, DEFAULT_OPTIONS
-from .logging import logger
+from .logging import create_logger
 
-from typing import Callable, Dict, Optional, Tuple
+from typing import Callable, Dict, Optional, Tuple, Union
 
 
 class Optimizer:
@@ -40,7 +40,9 @@ class Optimizer:
     :ivar starttime: Time at which optimization was started
     :ivar iteration: Current iteration
     :ivar converged: Flag indicating whether optimization has converged
-    :ivar exitflag:
+    :ivar exitflag: ExitFlag to indicate reason for termination
+    :ivar verbose: Verbosity level for logging
+    :ivar logger: logger instance
     """
     def __init__(self, fun: Callable,
                  ub: np.ndarray,
@@ -112,7 +114,8 @@ def __init__(self, fun: Callable,
         self.iteration: int = 0
         self.converged: bool = False
         self.exitflag: ExitFlag = ExitFlag.DID_NOT_RUN
-        logger.setLevel(verbose)
+        self.verbose: int = verbose
+        self.logger: Union[logging.Logger, None] = None
 
     def _reset(self):
         self.starttime = time.time()
@@ -121,6 +124,7 @@ def _reset(self):
         self.delta: float = self.get_option(Options.DELTA_INIT)
         self.delta_iter: float = self.delta
         self.fval_min = np.inf
+        self.logger = create_logger(self.verbose)
 
     def minimize(self, x0: np.ndarray):
         """
@@ -232,6 +236,7 @@ def minimize(self, x0: np.ndarray):
                     subspace_dim=self.get_option(Options.SUBSPACE_DIM),
                     stepback_strategy=self.get_option(Options.STEPBACK_STRAT),
                     refine_stepback=self.get_option(Options.REFINE_STEPBACK),
+                    logger=self.logger
                 )
 
             x_new = self.x + step.s + step.s0
@@ -385,7 +390,7 @@ def check_convergence(self, fval, x, grad) -> None:
 
         if np.isclose(fval, self.fval, atol=fatol, rtol=frtol):
             self.exitflag = ExitFlag.FTOL
-            logger.warning(
+            self.logger.warning(
                 'Stopping as function difference '
                 f'{np.abs(self.fval - fval):.2E} was smaller than specified '
                 f'tolerances (atol={fatol:.2E}, rtol={frtol:.2E})'
@@ -394,7 +399,7 @@ def check_convergence(self, fval, x, grad) -> None:
 
         elif np.isclose(x, self.x, atol=xatol, rtol=xrtol).all():
             self.exitflag = ExitFlag.XTOL
-            logger.warning(
+            self.logger.warning(
                 'Stopping as step was smaller than specified tolerances ('
                 f'atol={xatol:.2E}, rtol={xrtol:.2E})'
             )
@@ -402,7 +407,7 @@ def check_convergence(self, fval, x, grad) -> None:
 
         elif gnorm <= gatol:
             self.exitflag = ExitFlag.GTOL
-            logger.warning(
+            self.logger.warning(
                 'Stopping as gradient norm satisfies absolute convergence '
                 f'criteria: {gnorm:.2E} < {gatol:.2E}'
             )
@@ -410,7 +415,7 @@ def check_convergence(self, fval, x, grad) -> None:
 
         elif gnorm <= grtol * self.fval:
             self.exitflag = ExitFlag.GTOL
-            logger.warning(
+            self.logger.warning(
                 'Stopping as gradient norm satisfies relative convergence '
                 f'criteria: {gnorm:.2E} < {grtol:.2E} * {self.fval:.2E}'
             )
@@ -433,7 +438,7 @@ def check_continue(self) -> bool:
         maxiter = self.get_option(Options.MAXITER)
         if self.iteration >= maxiter:
             self.exitflag = ExitFlag.MAXITER
-            logger.warning(
+            self.logger.warning(
                 f'Stopping as maximum number of iterations {maxiter} was '
                 f'exceeded.'
             )
@@ -445,7 +450,7 @@ def check_continue(self) -> bool:
         avg_iter_time = time_elapsed/(self.iteration + (self.iteration == 0))
         if time_remaining < avg_iter_time:
             self.exitflag = ExitFlag.MAXTIME
-            logger.warning(
+            self.logger.warning(
                 f'Stopping as maximum runtime {maxtime} is expected to be '
                 f'exceeded in the next iteration.'
             )
@@ -511,22 +516,24 @@ def log_step(self, accepted: bool, step: Step, fval: float):
             for count in [step.reflection_count, step.truncation_count]
         ]
 
-        if np.isnan(fval):
+        if not np.isfinite(fval):
             fval = self.fval
-        logger.info(f'{" " * iterspaces}{self.iteration}'
-                    f' | {fval if accepted else self.fval:.3E}'
-                    f' | {(fval - self.fval)*accepted:+.2E}'
-                    f' | {step.qpval:+.2E}'
-                    f' | {self.tr_ratio:+.2E}'
-                    f' | {self.delta_iter:.2E}'
-                    f' | {norm(self.grad):.2E}'
-                    f' | {normdx:.2E}'
-                    f' | {step.theta:.2E}'
-                    f' | {step.alpha:.2E}'
-                    f' | {step.type}{" " * steptypespaces}'
-                    f' | {" " * reflspaces}{step.reflection_count}'
-                    f' | {" " * trunspaces}{step.truncation_count}'
-                    f' | {int(accepted)}')
+        self.logger.info(
+            f'{" " * iterspaces}{self.iteration}'
+            f' | {fval if accepted else self.fval:+.3E}'
+            f' | {(fval - self.fval):+.2E}'
+            f' | {step.qpval:+.2E}'
+            f' | {self.tr_ratio:+.2E}'
+            f' | {self.delta_iter:.2E}'
+            f' | {norm(self.grad):.2E}'
+            f' | {normdx:.2E}'
+            f' | {step.theta:.2E}'
+            f' | {step.alpha:.2E}'
+            f' | {step.type}{" " * steptypespaces}'
+            f' | {" " * reflspaces}{step.reflection_count}'
+            f' | {" " * trunspaces}{step.truncation_count}'
+            f' | {int(accepted)}'
+        )
 
     def log_step_initial(self):
         """
@@ -535,20 +542,22 @@ def log_step_initial(self):
 
         iterspaces = max(len(str(self.get_option(Options.MAXITER))), 5) - \
             len(str(self.iteration))
-        logger.info(f'{" " * iterspaces}{self.iteration}'
-                    f' | {self.fval:.3E}'
-                    f' |    NaN   '
-                    f' |    NaN   '
-                    f' |    NaN   '
-                    f' | {self.delta:.2E}'
-                    f' | {norm(self.grad):.2E}'
-                    f' |   NaN   '
-                    f' |   NaN   '
-                    f' |   NaN   '
-                    f' | NaN '
-                    f' | NaN '
-                    f' | NaN '
-                    f' | {int(np.isfinite(self.fval))}')
+        self.logger.info(
+            f'{" " * iterspaces}{self.iteration}'
+            f' | {self.fval:+.3E}'
+            f' |    NaN   '
+            f' |    NaN   '
+            f' |    NaN   '
+            f' | {self.delta:.2E}'
+            f' | {norm(self.grad):.2E}'
+            f' |   NaN   '
+            f' |   NaN   '
+            f' |   NaN   '
+            f' | NaN '
+            f' | NaN '
+            f' | NaN '
+            f' | {int(np.isfinite(self.fval))}'
+        )
 
     def log_header(self):
         """
@@ -557,10 +566,12 @@ def log_header(self):
         """
         iterspaces = len(str(self.get_option(Options.MAXITER))) - 5
 
-        logger.info(f'{" " * iterspaces} iter '
-                    f'|   fval    | fval diff | pred diff | tr ratio  '
-                    f'|  delta   |  ||g||   | ||step|| |  theta   |  alpha   '
-                    f'| step | refl | trun | accept')
+        self.logger.info(
+            f'{" " * iterspaces} iter '
+            f'|    fval    | fval diff | pred diff | tr ratio  '
+            f'|  delta   |  ||g||   | ||step|| |  theta   |  alpha   '
+            f'| step | refl | trun | accept'
+        )
 
     def check_finite(self,
                      grad: Optional[np.ndarray] = None,
diff --git a/fides/stepback.py b/fides/stepback.py
new file mode 100644
index 0000000..c4f4115
--- /dev/null
+++ b/fides/stepback.py
@@ -0,0 +1,176 @@
+"""
+Trust Region StepBack
+-----------------------------
+This module provides the machinery to combine various step-back strategies
+that can be used to compute longer steps in case the initially proposed step
+had to be truncated due to non-compliance with boundary constraints.
+"""
+
+import numpy as np
+from scipy.sparse import csc_matrix
+from typing import List, Sequence
+
+from .steps import Step, TRStepReflected, TRStepTruncated, RefinedStep
+
+
+def stepback_reflect(tr_step: Step,
+                     x: np.ndarray,
+                     sg: np.ndarray,
+                     hess: np.ndarray,
+                     scaling: csc_matrix,
+                     g_dscaling: csc_matrix,
+                     delta: float,
+                     theta: float,
+                     ub: np.ndarray,
+                     lb: np.ndarray) -> List[Step]:
+    """
+    Compute new proposal steps according to a reflection strategy.
+
+    :param tr_step:
+        Reference trust region step that will be reflect
+    :param x:
+        Current values of the optimization variables
+    :param sg:
+        Rescaled objective function gradient at x
+    :param hess:
+        (Approximate) objective function Hessian at x
+    :param g_dscaling:
+        Unscaled gradient multiplied by derivative of scaling
+        transformation
+    :param scaling:
+        Scaling transformation according to distance to boundary
+    :param delta:
+        Trust region radius, note that this applies after scaling
+        transformation
+    :param theta:
+        parameter regulating stepback
+    :param lb:
+        lower optimization variable boundaries
+    :param ub:
+        upper optimization variable boundaries
+
+    :return:
+        New proposal steps
+    """
+    rtr_step = TRStepReflected(x, sg, hess, scaling, g_dscaling, delta,
+                               theta, ub, lb, tr_step)
+    rtr_step.calculate()
+    steps = [rtr_step]
+    for ireflection in range(len(x) - 1):
+        if rtr_step.alpha == 1.0:
+            break
+        # recursively add more reflections
+        rtr_old = rtr_step
+        rtr_step = TRStepReflected(x, sg, hess, scaling, g_dscaling, delta,
+                                   theta, ub, lb, rtr_old)
+        rtr_step.calculate()
+        steps.append(rtr_step)
+
+    return steps
+
+
+def stepback_truncate(tr_step: Step,
+                      x: np.ndarray,
+                      sg: np.ndarray,
+                      hess: np.ndarray,
+                      scaling: csc_matrix,
+                      g_dscaling: csc_matrix,
+                      delta: float,
+                      theta: float,
+                      ub: np.ndarray,
+                      lb: np.ndarray) -> List[Step]:
+    """
+    Compute new proposal steps according to a truncation strategy.
+
+    :param tr_step:
+        Reference trust region step that will be reflect
+    :param x:
+        Current values of the optimization variables
+    :param sg:
+        Rescaled objective function gradient at x
+    :param hess:
+        (Approximate) objective function Hessian at x
+    :param g_dscaling:
+        Unscaled gradient multiplied by derivative of scaling
+        transformation
+    :param scaling:
+        Scaling transformation according to distance to boundary
+    :param delta:
+        Trust region radius, note that this applies after scaling
+        transformation
+    :param theta:
+        parameter regulating stepback
+    :param lb:
+        lower optimization variable boundaries
+    :param ub:
+        upper optimization variable boundaries
+
+    :return:
+        New proposal steps
+    """
+    rtt_step = TRStepTruncated(x, sg, hess, scaling, g_dscaling, delta,
+                               theta, ub, lb, tr_step)
+    rtt_step.calculate()
+    steps = [rtt_step]
+    while rtt_step.subspace.shape[1] > 0:
+        if rtt_step.alpha == 1.0:
+            break
+        rtt_step = TRStepTruncated(x, sg, hess, scaling, g_dscaling, delta,
+                                   theta, ub, lb, rtt_step)
+        rtt_step.calculate()
+        steps.append(rtt_step)
+
+    return steps
+
+
+def stepback_refine(steps: Sequence[Step],
+                    x: np.ndarray,
+                    sg: np.ndarray,
+                    hess: np.ndarray,
+                    scaling: csc_matrix,
+                    g_dscaling: csc_matrix,
+                    delta: float,
+                    theta: float,
+                    ub: np.ndarray,
+                    lb: np.ndarray) -> List[Step]:
+    """
+    Refine a promising subset of the provided steps based on trust-constr
+    optimization
+
+    :param steps:
+        Reference trust region step that will be reflect
+    :param x:
+        Current values of the optimization variables
+    :param sg:
+        Rescaled objective function gradient at x
+    :param hess:
+        (Approximate) objective function Hessian at x
+    :param g_dscaling:
+        Unscaled gradient multiplied by derivative of scaling
+        transformation
+    :param scaling:
+        Scaling transformation according to distance to boundary
+    :param delta:
+        Trust region radius, note that this applies after scaling
+        transformation
+    :param theta:
+        parameter regulating stepback
+    :param lb:
+        lower optimization variable boundaries
+    :param ub:
+        upper optimization variable boundaries
+
+    :return:
+        New proposal steps
+    """
+    min_qpval = np.nanmin([step.qpval for step in steps])
+    ref_steps = [
+        RefinedStep(x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb,
+                    step)
+        for step in steps
+        if ((step.alpha == 1.0 and step.type not in ['trnd', 'tr2d', 'grad'])
+            or (step.alpha < 1.0 and step.qpval < min_qpval / 2))
+    ]
+    for step in ref_steps:
+        step.calculate()
+    return ref_steps
diff --git a/fides/steps.py b/fides/steps.py
new file mode 100644
index 0000000..541bfb8
--- /dev/null
+++ b/fides/steps.py
@@ -0,0 +1,425 @@
+"""
+Trust Region Step Calculation
+-----------------------------
+This module provides the machinery to calculate different trust-region(
+-reflective) step proposals
+"""
+
+
+import numpy as np
+import scipy.linalg as linalg
+
+from numpy.linalg import norm
+from scipy.sparse import csc_matrix
+from scipy.optimize import Bounds, NonlinearConstraint, minimize
+
+from logging import Logger
+from .subproblem import (
+    solve_1d_trust_region_subproblem, solve_nd_trust_region_subproblem
+)
+
+from typing import Union
+
+
+def quadratic_form(Q: np.ndarray, p: np.ndarray, x: np.ndarray) -> float:
+    """
+    Computes the quadratic form :math:`x^TQx + x^Tp`
+
+    :param Q: Matrix
+    :param p: Vector
+    :param x: Input
+
+    :return:
+        Value of form
+    """
+    return 0.5 * x.T.dot(Q).dot(x) + p.T.dot(x)
+
+
+def normalize(v: np.ndarray) -> None:
+    """
+    Inplace normalization of a vector
+
+    :param v:
+        vector to be normalized
+    """
+    nv = norm(v)
+    if nv > 0:
+        v[:] = v/nv  # change inplace
+
+
+class Step:
+    """
+    Base class for the computation of a proposal step
+
+    :ivar x: Current state of optimization variables
+    :ivar s: Proposed step
+    :ivar sc: Coefficients in the 1D/2D subspace that defines the affine
+        transformed step ss: `ss = subspace * sc`
+    :ivar ss: Affine transformed step: `s = scaling * ss`
+    :ivar og_s: `s` without step back
+    :ivar og_sc: `st` without step back
+    :ivar og_ss: `ss` without step back
+    :ivar sg: Rescaled gradient `scaling * g`
+    :ivar hess: Hessian of the objective function at `x`
+    :ivar g_dscaling: `diag(g) * dscaling`
+    :ivar delta: Trust region radius in the transformed space defined by
+        scaling matrix
+    :ivar theta: Controls step back, fraction of step to take if full
+        step would reach breakpoint
+    :ivar lb: Lower boundaries for x
+    :ivar ub: Upper boundaries for x
+    :ivar minbr: Maximal fraction of step s that can be taken to reach
+        first breakpoint
+    :ivar iminbr: Index of x that specifies the variable that will hit the
+        breakpoint if a step minbr * s is taken
+    :ivar qpval: Value of the quadratic subproblem for the proposed step
+    :ivar shess: Matrix of the full quadratic problem
+    :ivar cg: Projection of the g_hat to the subspace
+    :ivar chess: Projection of the B to the subspace
+    :ivar reflection_count: Number of reflections that were applied to
+        obtain this step
+    :ivar truncation_count: Number of reflections that were applied to
+        obtain this step
+
+    :cvar type: Identifier that allows identification of subclasses
+    """
+    type = 'step'
+
+    def __init__(self,
+                 x: np.ndarray,
+                 sg: np.ndarray,
+                 hess: np.ndarray,
+                 scaling: csc_matrix,
+                 g_dscaling: csc_matrix,
+                 delta: float,
+                 theta: float,
+                 ub: np.ndarray,
+                 lb: np.ndarray,
+                 logger: Logger):
+        """
+
+        :param x:
+            Reference point
+        :param sg:
+            Gradient in rescaled coordinates
+        :param hess:
+            Hessian in unscaled coordinates
+        :param scaling:
+            Matrix that defines scaling transformation
+        :param g_dscaling:
+            Unscaled gradient multiplied by derivative of scaling
+            transformation
+        :param delta:
+            Trust region Radius in scaled coordinates
+        :param theta:
+            Stepback parameter that controls how close steps are allowed to
+            get to the boundary
+        :param ub:
+            Upper boundary
+        :param lb:
+            Lower boundary
+
+        """
+        self.x: np.ndarray = x
+
+        self.s: Union[np.ndarray, None] = None
+        self.sc: Union[np.ndarray, None] = None
+        self.ss: Union[np.ndarray, None] = None
+
+        self.og_s: Union[np.ndarray, None] = None
+        self.og_sc: Union[np.ndarray, None] = None
+        self.og_ss: Union[np.ndarray, None] = None
+
+        self.sg: np.ndarray = sg
+        self.scaling: csc_matrix = scaling
+
+        self.delta: float = delta
+        self.theta: float = theta
+
+        self.lb: np.ndarray = lb
+        self.ub: np.ndarray = ub
+
+        self.br: np.ndarray = np.ones(sg.shape)
+        self.minbr: float = 1.0
+        self.alpha: float = 1.0
+        self.iminbr: np.ndarray = np.array([])
+
+        self.qpval: float = 0.0
+
+        self.shess: np.ndarray = np.asarray(scaling * hess * scaling
+                                            + g_dscaling)
+
+        self.cg: Union[np.ndarray, None] = None
+        self.chess: Union[np.ndarray, None] = None
+        self.subspace: Union[np.ndarray, None] = None
+
+        self.s0: np.ndarray = np.zeros(sg.shape)
+        self.ss0: np.ndarray = np.zeros(sg.shape)
+
+        self.reflection_count: int = 0
+        self.truncation_count: int = 0
+        self.logger: Logger = logger
+
+    def step_back(self):
+        """
+        This function truncates the step based on the distance of the
+        current point to the boundary.
+        """
+        # create copies of the calculated step
+        self.og_s = self.s.copy()
+        self.og_ss = self.ss.copy()
+        self.og_sc = self.sc.copy()
+
+        nonzero = np.abs(self.s) > 0
+        self.br = np.inf * np.ones(self.s.shape)
+        if np.any(nonzero):
+            # br quantifies the distance to the boundary normalized
+            # by the proposed step, this indicates the fraction of the step
+            # that would put the respective variable at the boundary
+            # This is defined in [Coleman-Li1996] (3.1)
+            self.br[nonzero] = np.max(np.vstack([
+                (self.ub[nonzero] - self.x[nonzero])/self.s[nonzero],
+                (self.lb[nonzero] - self.x[nonzero])/self.s[nonzero]
+            ]), axis=0)
+        self.minbr = np.min(self.br)
+        self.iminbr = np.where(self.br == self.minbr)[0]
+        # compute the minimum of the step
+        self.alpha = np.min([1, self.theta * self.minbr])
+
+        self.s *= self.alpha
+        self.sc *= self.alpha
+        self.ss *= self.alpha
+
+    def reduce_to_subspace(self) -> None:
+        """
+        This function projects the matrix shess and the vector sg to the
+        subspace
+        """
+        self.chess = self.subspace.T.dot(self.shess.dot(self.subspace))
+        self.cg = self.subspace.T.dot(self.sg)
+
+    def compute_step(self) -> None:
+        """
+        Compute the step as solution to the trust region subproblem. Special
+        code is used for the special case 1-dimensional subspace case
+        """
+        if self.subspace.shape[1] == 0:
+            self.sc = np.empty((0, 0))
+            self.ss = np.zeros(self.ss0.shape)
+            self.s = np.zeros(self.s0.shape)
+            return
+        if self.subspace.shape[1] > 1:
+            self.sc, _ = solve_nd_trust_region_subproblem(
+                self.chess, self.cg,
+                np.sqrt(max(self.delta ** 2 - norm(self.ss0) ** 2, 0.0)),
+                self.logger
+            )
+        else:
+            self.sc = solve_1d_trust_region_subproblem(
+                self.shess, self.sg, self.subspace[:, 0], self.delta, self.ss0
+            )
+        self.ss = self.subspace.dot(np.real(self.sc))
+        self.s = self.scaling.dot(self.ss)
+
+    def calculate(self):
+        """
+        Calculates step and the expected objective function value according to
+        the quadratic approximation
+        """
+        self.reduce_to_subspace()
+        self.compute_step()
+        self.step_back()
+        self.qpval = quadratic_form(self.shess, self.sg, self.ss + self.ss0)
+
+
+class TRStepFull(Step):
+    """
+    This class provides the machinery to compute an exact solution of
+    the trust region subproblem.
+    """
+
+    type = 'trnd'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, logger):
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, logger)
+        self.subspace = np.eye(hess.shape[0])
+
+
+class TRStep2D(Step):
+    """
+    This class provides the machinery to compute an approximate solution of
+    the trust region subproblem according to a 2D subproblem
+    """
+
+    type = 'tr2d'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, logger):
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, logger)
+        n = len(sg)
+
+        try:
+            s_newt = linalg.solve(hess, sg)
+        except linalg.LinAlgError:
+            s_newt = linalg.lstsq(hess, sg)[0]
+        posdef = s_newt.dot(hess.dot(s_newt)) > 0
+        normalize(s_newt)
+
+        if n > 1:
+            if not posdef:
+                # in this case we are in Case 2 of Fig 12 in
+                # [Coleman-Li1994]
+                logger.debug('Newton direction did not have negative '
+                             'curvature adding scaling * np.sign(sg) to '
+                             '2D subspace.')
+                s_grad = scaling * np.sign(sg) + (sg == 0)
+            else:
+                s_grad = sg.copy()
+
+            # orthonormalize, this ensures that S.T.dot(S) = I and we
+            # can use S/S.T for transformation
+            s_grad = s_grad - s_newt * (s_newt.dot(s_grad))
+            normalize(s_grad)
+            # if non-zero, add s_grad to subspace
+            if np.any(s_grad != 0):
+                self.subspace = np.vstack([s_newt, s_grad]).T
+                return
+            else:
+                logger.debug('Singular subspace, continuing with 1D '
+                             'subspace.')
+
+        self.subspace = np.expand_dims(s_newt, 1)
+
+
+class TRStepReflected(Step):
+    """
+    This class provides the machinery to compute a reflected step based on
+    trust region subproblem solution that hit the boundaries.
+    """
+
+    type = 'trr'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, step: Step):
+        """
+        :param step:
+            Trust-region step that is reflected
+        """
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, step.logger)
+
+        alpha = min(step.minbr, 1)
+        self.s0 = alpha * step.og_s + step.s0
+        self.ss0 = alpha * step.og_ss + step.ss0
+        # update x and at breakpoint
+        self.x = x + self.s0
+
+        # reflect the transformed step at the boundary
+        nss = step.og_ss.copy()
+        nss[step.iminbr] *= -1
+        normalize(nss)
+        self.subspace = np.expand_dims(nss, 1)
+        self.reflection_count = step.reflection_count + 1
+
+
+class TRStepTruncated(Step):
+    """
+    This class provides the machinery to compute a reduced step based on
+    trust region subproblem solution that hit the boundaries.
+    """
+
+    type = 'trt'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, step: Step):
+        """
+        :param step:
+            Trust-region step that is reduced
+        """
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, step.logger)
+
+        self.s0 = step.s0.copy()
+        self.ss0 = step.ss0.copy()
+        iminbr = step.iminbr
+        self.s0[iminbr] += step.theta * step.br[iminbr] * step.og_s[iminbr]
+        self.ss0[iminbr] += step.theta * step.br[iminbr] * step.og_ss[iminbr]
+        # update x and at breakpoint
+        self.x = x + self.s0
+
+        subspace = step.subspace.copy()
+        subspace[iminbr, :] = 0
+        # reduce subspace
+        subspace = subspace[:, (subspace != 0).any(axis=0)]
+        # normalize subspace
+        for ix in range(subspace.shape[1]):
+            normalize(subspace[:, ix])
+        self.subspace = subspace
+
+        self.truncation_count = step.truncation_count + len(iminbr)
+
+
+class GradientStep(Step):
+    """
+    This class provides the machinery to compute a gradient step.
+    """
+
+    type = 'g'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, logger):
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, logger)
+        s_grad = sg.copy()
+        normalize(s_grad)
+        self.subspace = np.expand_dims(s_grad, 1)
+
+
+class RefinedStep(Step):
+    """
+    This class provides the machinery to refine a step based on interior
+    point optimization
+    """
+
+    type = 'ref'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, step):
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, step.logger)
+        s_grad = sg.copy()
+        normalize(s_grad)
+        self.subspace = np.expand_dims(s_grad, 1)
+        self.constraints = [
+            NonlinearConstraint(
+                fun=lambda xs: (norm(xs) - delta) * np.ones((1,)),
+                jac=lambda xs: np.expand_dims(xs, 1).T / norm(xs),
+                lb=np.zeros((1,)),
+                ub=np.ones((1,)) * np.inf,
+            )
+        ]
+        self.guess = step.ss + step.ss0
+        self.bounds = Bounds(
+            step.theta * (lb - x) / scaling.diagonal(),
+            step.theta * (ub - x) / scaling.diagonal()
+        )
+        self.reflection_count = step.reflection_count
+        self.truncation_count = step.truncation_count
+
+    def calculate(self):
+        res = minimize(fun=lambda s: quadratic_form(self.shess, self.sg, s),
+                       jac=lambda s: self.shess.dot(s) + self.sg,
+                       hess=lambda s: self.shess,
+                       x0=self.guess,
+                       method='trust-constr',
+                       bounds=self.bounds,
+                       constraints=self.constraints,
+                       options={'verbose': 0, 'maxiter': 10})
+        self.ss = res.x
+        self.s = self.scaling.dot(res.x)
+        self.sc = self.ss
+        self.step_back()
+        self.qpval = quadratic_form(self.shess, self.sg, self.ss)
diff --git a/fides/subproblem.py b/fides/subproblem.py
index 143a0c0..f4d2a28 100644
--- a/fides/subproblem.py
+++ b/fides/subproblem.py
@@ -5,15 +5,15 @@
 trust-region subproblems.
 """
 
+import logging
+
 import numpy as np
 from numpy.linalg import norm
 
 from scipy import linalg
 from scipy.optimize import newton, brentq
 
-from typing import Tuple
-
-from .logging import logger
+from typing import Tuple, Union
 
 
 def solve_1d_trust_region_subproblem(B: np.ndarray,
@@ -64,9 +64,10 @@ def solve_1d_trust_region_subproblem(B: np.ndarray,
     return tau * np.ones((1,))
 
 
-def solve_nd_trust_region_subproblem(B: np.ndarray,
-                                     g: np.ndarray,
-                                     delta: float) -> Tuple[np.ndarray, str]:
+def solve_nd_trust_region_subproblem(
+        B: np.ndarray, g: np.ndarray, delta: float,
+        logger: Union[logging.Logger, None] = None
+) -> Tuple[np.ndarray, str]:
     r"""
     This function exactly solves the n-dimensional subproblem.
 
@@ -96,12 +97,17 @@ def solve_nd_trust_region_subproblem(B: np.ndarray,
         Gradient of the quadratic subproblem
     :param delta:
         Norm boundary for the solution of the quadratic subproblem
+    :param logger:
+        Logger instance to be used for logging
 
     :return:
         s: Selected step,
         step_type: Type of solution that was obtained
 
     """
+    if logger is None:
+        logger = logging.getLogger('fides')
+
     if delta == 0:
         return np.zeros(g.shape), 'zero'
 
diff --git a/fides/trust_region.py b/fides/trust_region.py
index b4bdc06..6f1a8e2 100644
--- a/fides/trust_region.py
+++ b/fides/trust_region.py
@@ -1,423 +1,19 @@
 """
-Trust Region Steps
-------------------
-This module provides the machinery to compute different trust-region(
+Trust Region Step Evaluation
+----------------------------
+This module provides the machinery to evaluate different trust-region(
 -reflective) step proposals and select among them based on to their
 performance according to the quadratic approximation of the objective function
 """
 
 import numpy as np
-import scipy.linalg as linalg
+import logging
 
-from numpy.linalg import norm
 from scipy.sparse import csc_matrix
-from scipy.optimize import Bounds, NonlinearConstraint, minimize
 
-from .logging import logger
-from .subproblem import (
-    solve_1d_trust_region_subproblem, solve_nd_trust_region_subproblem
-)
 from .constants import SubSpaceDim, StepBackStrategy
-
-from typing import List, Sequence, Union
-
-
-def normalize(v: np.ndarray) -> None:
-    """
-    Inplace normalization of a vector
-
-    :param v:
-        vector to be normalized
-    """
-    nv = norm(v)
-    if nv > 0:
-        v[:] = v/nv  # change inplace
-
-
-def quadratic_form(Q: np.ndarray, p: np.ndarray, x: np.ndarray) -> float:
-    """
-    Computes the quadratic form :math:`x^TQx + x^Tp`
-
-    :param Q: Matrix
-    :param p: Vector
-    :param x: Input
-
-    :return:
-        Value of form
-    """
-    return 0.5 * x.T.dot(Q).dot(x) + p.T.dot(x)
-
-
-class Step:
-    """
-    Base class for the computation of a proposal step
-
-    :ivar x: Current state of optimization variables
-    :ivar s: Proposed step
-    :ivar sc: Coefficients in the 1D/2D subspace that defines the affine
-        transformed step ss: `ss = subspace * sc`
-    :ivar ss: Affine transformed step: `s = scaling * ss`
-    :ivar og_s: `s` without step back
-    :ivar og_sc: `st` without step back
-    :ivar og_ss: `ss` without step back
-    :ivar sg: Rescaled gradient `scaling * g`
-    :ivar hess: Hessian of the objective function at `x`
-    :ivar g_dscaling: `diag(g) * dscaling`
-    :ivar delta: Trust region radius in the transformed space defined by
-        scaling matrix
-    :ivar theta: Controls step back, fraction of step to take if full
-        step would reach breakpoint
-    :ivar lb: Lower boundaries for x
-    :ivar ub: Upper boundaries for x
-    :ivar minbr: Maximal fraction of step s that can be taken to reach
-        first breakpoint
-    :ivar iminbr: Index of x that specifies the variable that will hit the
-        breakpoint if a step minbr * s is taken
-    :ivar qpval: Value of the quadratic subproblem for the proposed step
-    :ivar shess: Matrix of the full quadratic problem
-    :ivar cg: Projection of the g_hat to the subspace
-    :ivar chess: Projection of the B to the subspace
-    :ivar reflection_count: Number of reflections that were applied to
-        obtain this step
-    :ivar truncation_count: Number of reflections that were applied to
-        obtain this step
-
-    :cvar type: Identifier that allows identification of subclasses
-    """
-    type = 'step'
-
-    def __init__(self,
-                 x: np.ndarray,
-                 sg: np.ndarray,
-                 hess: np.ndarray,
-                 scaling: csc_matrix,
-                 g_dscaling: csc_matrix,
-                 delta: float,
-                 theta: float,
-                 ub: np.ndarray,
-                 lb: np.ndarray):
-        """
-
-        :param x:
-            Reference point
-        :param sg:
-            Gradient in rescaled coordinates
-        :param hess:
-            Hessian in unscaled coordinates
-        :param scaling:
-            Matrix that defines scaling transformation
-        :param g_dscaling:
-            Unscaled gradient multiplied by derivative of scaling
-            transformation
-        :param delta:
-            Trust region Radius in scaled coordinates
-        :param theta:
-            Stepback parameter that controls how close steps are allowed to
-            get to the boundary
-        :param ub:
-            Upper boundary
-        :param lb:
-            Lower boundary
-
-        """
-        self.x: np.ndarray = x
-
-        self.s: Union[np.ndarray, None] = None
-        self.sc: Union[np.ndarray, None] = None
-        self.ss: Union[np.ndarray, None] = None
-
-        self.og_s: Union[np.ndarray, None] = None
-        self.og_sc: Union[np.ndarray, None] = None
-        self.og_ss: Union[np.ndarray, None] = None
-
-        self.sg: np.ndarray = sg
-        self.scaling: csc_matrix = scaling
-
-        self.delta: float = delta
-        self.theta: float = theta
-
-        self.lb: np.ndarray = lb
-        self.ub: np.ndarray = ub
-
-        self.br: np.ndarray = np.ones(sg.shape)
-        self.minbr: float = 1.0
-        self.alpha: float = 1.0
-        self.iminbr: np.ndarray = np.array([])
-
-        self.qpval: float = 0.0
-
-        self.shess: np.ndarray = np.asarray(scaling * hess * scaling
-                                            + g_dscaling)
-
-        self.cg: Union[np.ndarray, None] = None
-        self.chess: Union[np.ndarray, None] = None
-        self.subspace: Union[np.ndarray, None] = None
-
-        self.s0: np.ndarray = np.zeros(sg.shape)
-        self.ss0: np.ndarray = np.zeros(sg.shape)
-
-        self.reflection_count: int = 0
-        self.truncation_count: int = 0
-
-    def step_back(self):
-        """
-        This function truncates the step based on the distance of the
-        current point to the boundary.
-        """
-        # create copies of the calculated step
-        self.og_s = self.s.copy()
-        self.og_ss = self.ss.copy()
-        self.og_sc = self.sc.copy()
-
-        nonzero = np.abs(self.s) > 0
-        self.br = np.inf * np.ones(self.s.shape)
-        if np.any(nonzero):
-            # br quantifies the distance to the boundary normalized
-            # by the proposed step, this indicates the fraction of the step
-            # that would put the respective variable at the boundary
-            # This is defined in [Coleman-Li1996] (3.1)
-            self.br[nonzero] = np.max(np.vstack([
-                (self.ub[nonzero] - self.x[nonzero])/self.s[nonzero],
-                (self.lb[nonzero] - self.x[nonzero])/self.s[nonzero]
-            ]), axis=0)
-        self.minbr = np.min(self.br)
-        self.iminbr = np.where(self.br == self.minbr)[0]
-        # compute the minimum of the step
-        self.alpha = np.min([1, self.theta * self.minbr])
-
-        self.s *= self.alpha
-        self.sc *= self.alpha
-        self.ss *= self.alpha
-
-    def reduce_to_subspace(self) -> None:
-        """
-        This function projects the matrix shess and the vector sg to the
-        subspace
-        """
-        self.chess = self.subspace.T.dot(self.shess.dot(self.subspace))
-        self.cg = self.subspace.T.dot(self.sg)
-
-    def compute_step(self) -> None:
-        """
-        Compute the step as solution to the trust region subproblem. Special
-        code is used for the special case 1-dimensional subspace case
-        """
-        if self.subspace.shape[1] == 0:
-            self.sc = np.empty((0, 0))
-            self.ss = np.zeros(self.ss0.shape)
-            self.s = np.zeros(self.s0.shape)
-            return
-        if self.subspace.shape[1] > 1:
-            self.sc, _ = solve_nd_trust_region_subproblem(
-                self.chess, self.cg,
-                np.sqrt(max(self.delta ** 2 - norm(self.ss0) ** 2, 0.0))
-            )
-        else:
-            self.sc = solve_1d_trust_region_subproblem(
-                self.shess, self.sg, self.subspace[:, 0], self.delta, self.ss0
-            )
-        self.ss = self.subspace.dot(np.real(self.sc))
-        self.s = self.scaling.dot(self.ss)
-
-    def calculate(self):
-        """
-        Calculates step and the expected objective function value according to
-        the quadratic approximation
-        """
-        self.reduce_to_subspace()
-        self.compute_step()
-        self.step_back()
-        self.qpval = quadratic_form(self.shess, self.sg, self.ss + self.ss0)
-
-
-class TRStepFull(Step):
-    """
-    This class provides the machinery to compute an exact solution of
-    the trust region subproblem.
-    """
-
-    type = 'trnd'
-
-    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
-                 ub, lb):
-        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
-                         ub, lb)
-        self.subspace = np.eye(hess.shape[0])
-
-
-class TRStep2D(Step):
-    """
-    This class provides the machinery to compute an approximate solution of
-    the trust region subproblem according to a 2D subproblem
-    """
-
-    type = 'tr2d'
-
-    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
-                 ub, lb):
-        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
-                         ub, lb)
-        n = len(sg)
-
-        s_newt = linalg.solve(hess, sg)
-        posdef = s_newt.dot(hess.dot(s_newt)) > 0
-        normalize(s_newt)
-
-        if n > 1:
-            if not posdef:
-                # in this case we are in Case 2 of Fig 12 in
-                # [Coleman-Li1994]
-                logger.debug('Newton direction did not have negative '
-                             'curvature adding scaling * np.sign(sg) to '
-                             '2D subspace.')
-                s_grad = scaling * np.sign(sg) + (sg == 0)
-            else:
-                s_grad = sg.copy()
-
-            # orthonormalize, this ensures that S.T.dot(S) = I and we
-            # can use S/S.T for transformation
-            s_grad = s_grad - s_newt * (s_newt.dot(s_grad))
-            normalize(s_grad)
-            # if non-zero, add s_grad to subspace
-            if np.any(s_grad != 0):
-                self.subspace = np.vstack([s_newt, s_grad]).T
-                return
-            else:
-                logger.debug('Singular subspace, continuing with 1D '
-                             'subspace.')
-
-        self.subspace = np.expand_dims(s_newt, 1)
-
-
-class TRStepReflected(Step):
-    """
-    This class provides the machinery to compute a reflected step based on
-    trust region subproblem solution that hit the boundaries.
-    """
-
-    type = 'trr'
-
-    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
-                 ub, lb, step: Step):
-        """
-        :param step:
-            Trust-region step that is reflected
-        """
-        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
-                         ub, lb)
-
-        alpha = min(step.minbr, 1)
-        self.s0 = alpha * step.og_s + step.s0
-        self.ss0 = alpha * step.og_ss + step.ss0
-        # update x and at breakpoint
-        self.x = x + self.s0
-
-        # reflect the transformed step at the boundary
-        nss = step.og_ss.copy()
-        nss[step.iminbr] *= -1
-        normalize(nss)
-        self.subspace = np.expand_dims(nss, 1)
-        self.reflection_count = step.reflection_count + 1
-
-
-class TRStepTruncated(Step):
-    """
-    This class provides the machinery to compute a reduced step based on
-    trust region subproblem solution that hit the boundaries.
-    """
-
-    type = 'trt'
-
-    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
-                 ub, lb, step: Step):
-        """
-        :param step:
-            Trust-region step that is reduced
-        """
-        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
-                         ub, lb)
-
-        self.s0 = step.s0.copy()
-        self.ss0 = step.ss0.copy()
-        iminbr = step.iminbr
-        self.s0[iminbr] += step.theta * step.br[iminbr] * step.og_s[iminbr]
-        self.ss0[iminbr] += step.theta * step.br[iminbr] * step.og_ss[iminbr]
-        # update x and at breakpoint
-        self.x = x + self.s0
-
-        subspace = step.subspace.copy()
-        subspace[iminbr, :] = 0
-        # reduce subspace
-        subspace = subspace[:, (subspace != 0).any(axis=0)]
-        # normalize subspace
-        for ix in range(subspace.shape[1]):
-            normalize(subspace[:, ix])
-        self.subspace = subspace
-
-        self.truncation_count = step.truncation_count + len(iminbr)
-
-
-class GradientStep(Step):
-    """
-    This class provides the machinery to compute a gradient step.
-    """
-
-    type = 'g'
-
-    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
-                 ub, lb):
-        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
-                         ub, lb)
-        s_grad = sg.copy()
-        normalize(s_grad)
-        self.subspace = np.expand_dims(s_grad, 1)
-
-
-class RefinedStep(Step):
-    """
-    This class provides the machinery to refine a step based on interior
-    point optimization
-    """
-
-    type = 'ref'
-
-    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
-                 ub, lb, step):
-        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
-                         ub, lb)
-        s_grad = sg.copy()
-        normalize(s_grad)
-        self.subspace = np.expand_dims(s_grad, 1)
-        self.constraints = [
-            NonlinearConstraint(
-                fun=lambda xs: (norm(xs) - delta) * np.ones((1,)),
-                jac=lambda xs: np.expand_dims(xs, 1).T / norm(xs),
-                lb=np.zeros((1,)),
-                ub=np.ones((1,)) * np.inf,
-            )
-        ]
-        self.guess = step.ss + step.ss0
-        self.bounds = Bounds(
-            step.theta * (lb - x) / scaling.diagonal(),
-            step.theta * (ub - x) / scaling.diagonal()
-        )
-        self.reflection_count = step.reflection_count
-        self.truncation_count = step.truncation_count
-
-    def calculate(self):
-        res = minimize(fun=lambda s: quadratic_form(self.shess, self.sg, s),
-                       jac=lambda s: self.shess.dot(s) + self.sg,
-                       hess=lambda s: self.shess,
-                       x0=self.guess,
-                       method='trust-constr',
-                       bounds=self.bounds,
-                       constraints=self.constraints,
-                       options={'verbose': 0, 'maxiter': 10})
-        self.ss = res.x
-        self.s = self.scaling.dot(res.x)
-        self.sc = self.ss
-        self.step_back()
-        self.qpval = quadratic_form(self.shess, self.sg, self.ss)
+from .steps import Step, GradientStep, TRStep2D, TRStepFull, TRStepReflected
+from .stepback import stepback_refine, stepback_reflect, stepback_truncate
 
 
 def trust_region(x: np.ndarray,
@@ -431,7 +27,8 @@ def trust_region(x: np.ndarray,
                  ub: np.ndarray,
                  subspace_dim: SubSpaceDim,
                  stepback_strategy: StepBackStrategy,
-                 refine_stepback: bool) -> Step:
+                 refine_stepback: bool,
+                 logger: logging.Logger) -> Step:
     """
     Compute a step according to the solution of the trust-region subproblem.
     If step-back is necessary, gradient and reflected trust region step are
@@ -467,6 +64,8 @@ def trust_region(x: np.ndarray,
     :param refine_stepback:
         If set to True, proposed steps that are computed via the specified
         stepback_strategy will be refined via optimization.
+    :param logger:
+        logging.Logger instance to be used for logging
 
     :return:
         s: proposed step,
@@ -481,11 +80,11 @@ def trust_region(x: np.ndarray,
 
     if subspace_dim == SubSpaceDim.TWO:
         tr_step = TRStep2D(
-            x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb
+            x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb, logger
         )
     elif subspace_dim == SubSpaceDim.FULL:
         tr_step = TRStepFull(
-            x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb,
+            x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb, logger
         )
     else:
         raise ValueError('Invalid choice of subspace dimension.')
@@ -498,7 +97,7 @@ def trust_region(x: np.ndarray,
     steps = [tr_step]
     if tr_step.alpha < 1.0 and len(g) > 1:
         g_step = GradientStep(x, sg, hess, scaling, g_dscaling, delta,
-                              theta, ub, lb)
+                              theta, ub, lb, logger)
         g_step.calculate()
 
         steps = [g_step]
@@ -540,166 +139,3 @@ def trust_region(x: np.ndarray,
 
     qpvals = [step.qpval for step in steps]
     return steps[int(np.argmin(qpvals))]
-
-
-def stepback_reflect(tr_step: Step,
-                     x: np.ndarray,
-                     sg: np.ndarray,
-                     hess: np.ndarray,
-                     scaling: csc_matrix,
-                     g_dscaling: csc_matrix,
-                     delta: float,
-                     theta: float,
-                     ub: np.ndarray,
-                     lb: np.ndarray) -> List[Step]:
-    """
-    Compute new proposal steps according to a reflection strategy.
-
-    :param tr_step:
-        Reference trust region step that will be reflect
-    :param x:
-        Current values of the optimization variables
-    :param sg:
-        Rescaled objective function gradient at x
-    :param hess:
-        (Approximate) objective function Hessian at x
-    :param g_dscaling:
-        Unscaled gradient multiplied by derivative of scaling
-        transformation
-    :param scaling:
-        Scaling transformation according to distance to boundary
-    :param delta:
-        Trust region radius, note that this applies after scaling
-        transformation
-    :param theta:
-        parameter regulating stepback
-    :param lb:
-        lower optimization variable boundaries
-    :param ub:
-        upper optimization variable boundaries
-
-    :return:
-        New proposal steps
-    """
-    rtr_step = TRStepReflected(x, sg, hess, scaling, g_dscaling, delta,
-                               theta, ub, lb, tr_step)
-    rtr_step.calculate()
-    steps = [rtr_step]
-    for ireflection in range(len(x) - 1):
-        if rtr_step.alpha == 1.0:
-            break
-        # recursively add more reflections
-        rtr_old = rtr_step
-        rtr_step = TRStepReflected(x, sg, hess, scaling, g_dscaling, delta,
-                                   theta, ub, lb, rtr_old)
-        rtr_step.calculate()
-        steps.append(rtr_step)
-
-    return steps
-
-
-def stepback_truncate(tr_step: Step,
-                      x: np.ndarray,
-                      sg: np.ndarray,
-                      hess: np.ndarray,
-                      scaling: csc_matrix,
-                      g_dscaling: csc_matrix,
-                      delta: float,
-                      theta: float,
-                      ub: np.ndarray,
-                      lb: np.ndarray) -> List[Step]:
-    """
-    Compute new proposal steps according to a truncation strategy.
-
-    :param tr_step:
-        Reference trust region step that will be reflect
-    :param x:
-        Current values of the optimization variables
-    :param sg:
-        Rescaled objective function gradient at x
-    :param hess:
-        (Approximate) objective function Hessian at x
-    :param g_dscaling:
-        Unscaled gradient multiplied by derivative of scaling
-        transformation
-    :param scaling:
-        Scaling transformation according to distance to boundary
-    :param delta:
-        Trust region radius, note that this applies after scaling
-        transformation
-    :param theta:
-        parameter regulating stepback
-    :param lb:
-        lower optimization variable boundaries
-    :param ub:
-        upper optimization variable boundaries
-
-    :return:
-        New proposal steps
-    """
-    rtt_step = TRStepTruncated(x, sg, hess, scaling, g_dscaling, delta,
-                               theta, ub, lb, tr_step)
-    rtt_step.calculate()
-    steps = [rtt_step]
-    while rtt_step.subspace.shape[1] > 0:
-        if rtt_step.alpha == 1.0:
-            break
-        rtt_step = TRStepTruncated(x, sg, hess, scaling, g_dscaling, delta,
-                                   theta, ub, lb, rtt_step)
-        rtt_step.calculate()
-        steps.append(rtt_step)
-
-    return steps
-
-
-def stepback_refine(steps: Sequence[Step],
-                    x: np.ndarray,
-                    sg: np.ndarray,
-                    hess: np.ndarray,
-                    scaling: csc_matrix,
-                    g_dscaling: csc_matrix,
-                    delta: float,
-                    theta: float,
-                    ub: np.ndarray,
-                    lb: np.ndarray) -> List[Step]:
-    """
-    Refine a promising subset of the provided steps based on trust-constr
-    optimization
-
-    :param steps:
-        Reference trust region step that will be reflect
-    :param x:
-        Current values of the optimization variables
-    :param sg:
-        Rescaled objective function gradient at x
-    :param hess:
-        (Approximate) objective function Hessian at x
-    :param g_dscaling:
-        Unscaled gradient multiplied by derivative of scaling
-        transformation
-    :param scaling:
-        Scaling transformation according to distance to boundary
-    :param delta:
-        Trust region radius, note that this applies after scaling
-        transformation
-    :param theta:
-        parameter regulating stepback
-    :param lb:
-        lower optimization variable boundaries
-    :param ub:
-        upper optimization variable boundaries
-
-    :return:
-        New proposal steps
-    """
-    min_qpval = np.nanmin([step.qpval for step in steps])
-    ref_steps = [
-        RefinedStep(x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb,
-                    step)
-        for step in steps
-        if ((step.alpha == 1.0 and step.type not in ['trnd', 'tr2d', 'grad'])
-            or (step.alpha < 1.0 and step.qpval < min_qpval / 2))
-    ]
-    for step in ref_steps:
-        step.calculate()
-    return ref_steps
diff --git a/fides/version.py b/fides/version.py
index 3ced358..b5fdc75 100644
--- a/fides/version.py
+++ b/fides/version.py
@@ -1,1 +1,1 @@
-__version__ = "0.2.1"
+__version__ = "0.2.2"


Test patch:
diff --git a/tests/test_subproblem.py b/tests/test_subproblem.py
index d4b3762..a55207a 100644
--- a/tests/test_subproblem.py
+++ b/tests/test_subproblem.py
@@ -1,7 +1,7 @@
 from fides.subproblem import (
     solve_nd_trust_region_subproblem, solve_1d_trust_region_subproblem,
 )
-from fides.trust_region import normalize
+from fides.steps import normalize
 from scipy.spatial.transform import Rotation as R
 
 import numpy as np