diff --git a/popt/loop/ensemble.py b/popt/loop/ensemble.py
index 77ff998d..507ad8ca 100644
--- a/popt/loop/ensemble.py
+++ b/popt/loop/ensemble.py
@@ -10,7 +10,6 @@
 from popt.misc_tools import optim_tools as ot
 from pipt.misc_tools import analysis_tools as at
 from ensemble.ensemble import Ensemble as PETEnsemble
-from popt.loop.extensions import GenOptExtension
 
 
 class Ensemble(PETEnsemble):
@@ -137,12 +136,6 @@ def __set__variable(var_name=None, defalut=None):
         self.bias_weights = np.ones(self.num_samples) / self.num_samples  # initialize with equal weights
         self.bias_points = None  # this is the points used to estimate the bias correction
 
-        # Setup GenOpt
-        self.genopt = GenOptExtension(self.get_state(), 
-                                      self.get_cov(), 
-                                      func=self.function, 
-                                      ne=self.num_samples)
-
     def get_state(self):
         """
         Returns
diff --git a/popt/loop/base.py b/popt/loop/ensemble_base.py
similarity index 57%
rename from popt/loop/base.py
rename to popt/loop/ensemble_base.py
index 33d4497c..2500fc2b 100644
--- a/popt/loop/base.py
+++ b/popt/loop/ensemble_base.py
@@ -9,63 +9,81 @@
 from popt.misc_tools import optim_tools as ot
 from pipt.misc_tools import analysis_tools as at
 from ensemble.ensemble import Ensemble as PETEnsemble
+from simulator.simple_models import noSimulation
 
-class EnsembleOptimizationBase(PETEnsemble):
+class EnsembleOptimizationBaseClass(PETEnsemble):
     '''
     Base class for the popt ensemble
     '''
-    def __init__(self, kwargs_ens, sim, obj_func):
+    def __init__(self, options, simulator, objective):
         '''
         Parameters
         ----------
-        kwargs_ens : dict
+        options : dict
             Options for the ensemble class
         
-        sim : callable
-            The forward simulator (e.g. flow)
+        simulator : callable
+            The forward simulator (e.g. flow). If None, no simulation is performed.
         
-        obj_func : callable
+        objective : callable
             The objective function (e.g. npv)
         '''
+        if simulator is None:
+            sim = noSimulation()
+        else:
+            sim = simulator
 
         # Initialize PETEnsemble
-        super().__init__(kwargs_ens, sim)
-
-        self.save_prediction = kwargs_ens.get('save_prediction', None)
-        self.num_models  = kwargs_ens.get('num_models', 1)
-        self.transform   = kwargs_ens.get('transform', False)
-        self.num_samples = self.ne
+        super().__init__(options, sim)
 
-        # Get bounds and varaince
-        self.upper_bound = []
-        self.lower_bound = []
+        # Unpack some options
+        self.save_prediction = options.get('save_prediction', None)
+        self.num_models      = options.get('num_models', 1)
+        self.transform       = options.get('transform', False)
+        self.num_samples     = self.ne
+        
+        # Define some variables
+        self.lb = []
+        self.ub = []
         self.bounds = []
         self.cov = np.array([])
-        for name in self.prior_info.keys():
-            self.state[name] = np.asarray(self.prior_info[name]['mean'])
-            num_state_var = len(self.state[name])
-            value_cov = self.prior_info[name]['variance'] * np.ones((num_state_var,))
-            if 'limits' in self.prior_info[name].keys():
-                lb = self.prior_info[name]['limits'][0]
-                ub = self.prior_info[name]['limits'][1]
-                self.lower_bound.append(lb)
-                self.upper_bound.append(ub)
+
+        # Get bounds and varaince, and initialize state
+        for key in self.prior_info.keys():
+            variable = self.prior_info[key]
+
+            # mean
+            self.state[key] = np.asarray(variable['mean'])
+
+            # Covariance
+            dim = self.state[key].size
+            cov = variable['variance']*np.ones(dim)
+
+            if 'limits' in variable.keys():
+                lb, ub = variable['limits']
+                self.lb(lb)
+                self.ub(ub)
+
+                # transform cov to [0, 1] if transform is True
                 if self.transform:
-                    value_cov = value_cov / (ub - lb)**2
-                    np.clip(value_cov, 0, 1, out=value_cov)
-                    self.bounds += num_state_var*[(0, 1)]
+                    cov = np.clip(cov/(ub - lb)**2, 0, 1, out=cov)
+                    self.bounds += dim*[(0, 1)]
                 else:
-                    self.bounds += num_state_var*[(lb, ub)]
-                self.cov = np.append(self.cov, value_cov)
+                    self.bounds += dim*[(lb, ub)]
             else:
-                self.bounds += num_state_var*[(None, None)]
+                self.bounds += dim*[(None, None)]
+
+            # Add to covariance
+            self.cov = np.append(self.cov, cov)
             
-        
-        self._scale_state()
+        # Make cov full covariance matrix
         self.cov = np.diag(self.cov)
 
+        # Scale the state to [0, 1] if transform is True
+        self._scale_state()
+
         # Set objective function (callable)
-        self.obj_func = obj_func
+        self.obj_func = objective
 
         # Objective function values
         self.state_func_values = None
@@ -78,8 +96,13 @@ def get_state(self):
         x : numpy.ndarray
             Control vector as ndarray, shape (number of controls, number of perturbations)
         """
-        x = ot.aug_optim_state(self.state, list(self.state.keys()))
-        return x
+        return ot.aug_optim_state(self.state, list(self.state.keys()))
+    
+    def vec_to_state(self, x):
+        """
+        Converts a control vector to the internal state representation.
+        """
+        return ot.update_optim_state(x, self.state, list(self.state.keys()))
 
     def get_bounds(self):
         """
@@ -112,7 +135,10 @@ def function(self, x, *args):
         else:
             self.ne = x.shape[1]
 
-        self.state = ot.update_optim_state(x, self.state, list(self.state.keys()))  # go from nparray to dict
+        # convert x to state
+        self.state = self.vec_to_state(x)  # go from nparray to dict
+
+        # run the simulation
         self._invert_scale_state()  # ensure that state is in [lb,ub]
         run_success = self.calc_prediction(save_prediction=self.save_prediction)  # calculate flow data
         self._scale_state()  # scale back to [0, 1]
@@ -147,17 +173,17 @@ def _scale_state(self):
         """
         Transform the internal state from [lb, ub] to [0, 1]
         """
-        if self.transform and (self.upper_bound and self.lower_bound):
+        if self.transform and (self.lb and self.ub):
             for i, key in enumerate(self.state):
-                self.state[key] = (self.state[key] - self.lower_bound[i])/(self.upper_bound[i] - self.lower_bound[i])
+                self.state[key] = (self.state[key] - self.lb[i])/(self.ub[i] - self.lb[i])
                 np.clip(self.state[key], 0, 1, out=self.state[key])
 
     def _invert_scale_state(self):
         """
         Transform the internal state from [0, 1] to [lb, ub]
         """
-        if self.transform and (self.upper_bound and self.lower_bound):
+        if self.transform and (self.lb and self.ub):
             for i, key in enumerate(self.state):
                 if self.transform:
-                    self.state[key] = self.lower_bound[i] + self.state[key]*(self.upper_bound[i] - self.lower_bound[i])
-                np.clip(self.state[key], self.lower_bound[i], self.upper_bound[i], out=self.state[key])
\ No newline at end of file
+                    self.state[key] = self.lb[i] + self.state[key]*(self.ub[i] - self.lb[i])
+                np.clip(self.state[key], self.lb[i], self.ub[i], out=self.state[key])
\ No newline at end of file
diff --git a/popt/loop/generalized_ensemble.py b/popt/loop/generalized_ensemble.py
index cacbe992..81809c26 100644
--- a/popt/loop/generalized_ensemble.py
+++ b/popt/loop/generalized_ensemble.py
@@ -10,33 +10,32 @@
 # Internal imports
 from popt.misc_tools import optim_tools as ot
 from pipt.misc_tools import analysis_tools as at
-from popt.loop.base import EnsembleOptimizationBase
+from popt.loop.ensemble_base import EnsembleOptimizationBaseClass
 
-class GeneralizedEnsemble(EnsembleOptimizationBase):
+class GeneralizedEnsemble(EnsembleOptimizationBaseClass):
 
-    def __init__(self, kwargs_ens, sim, obj_func):
+    def __init__(self, options, simulator, objective):
         '''
         Parameters
         ----------
-        kwargs_ens : dict
+        options : dict
             Options for the ensemble class
         
-        sim : callable
-            The forward simulator (e.g. flow)
+        simulator : callable
+            The forward simulator (e.g. flow). If None, no simulation is performed.
         
-        obj_func : callable
+        objective : callable
             The objective function (e.g. npv)
         '''
-        super().__init__(kwargs_ens, sim, obj_func)
-
-        self.dim   = self.get_state().size
+        super().__init__(options, simulator, objective)
 
         # construct corr matrix
         std = np.sqrt(np.diag(self.cov))
         self.corr = self.cov/np.outer(std, std)
+        self.dim  = std 
 
         # choose marginal
-        marginal = kwargs_ens.get('marginal', 'Beta')
+        marginal = options.get('marginal', 'BetaMC')
 
         if marginal in ['Beta', 'BetaMC', 'Logistic', 'TruncGaussian', 'Gaussian']:
 
@@ -45,7 +44,7 @@ def __init__(self, kwargs_ens, sim, obj_func):
 
             if marginal == 'Beta':
                 self.margs = Beta()
-                self.theta = kwargs_ens.get('theta', np.array([[20.0, 20.0] for _ in range(self.dim)]))
+                self.theta = options.get('theta', np.array([[20.0, 20.0] for _ in range(self.dim)]))
                 self.eps = self.var2eps()
                 self.grad_scale = 1/(2*self.eps)
                 self.hess_scale = 1/(4*self.eps**2)
@@ -56,20 +55,20 @@ def __init__(self, kwargs_ens, sim, obj_func):
                 var = np.diag(self.cov)
                 self.margs = BetaMC(lb, ub, 0.1*np.sqrt(var[0]))
                 default_theta = np.array([var_to_concentration(state[i], var[i], lb[i], ub[i]) for i in range(self.dim)])
-                self.theta = kwargs_ens.get('theta', default_theta)
+                self.theta = options.get('theta', default_theta)
                 
             elif marginal == 'Logistic':
                 self.margs = Logistic()
-                self.theta = kwargs_ens.get('theta', self.margs.var_to_scale(np.diag(self.cov)))
+                self.theta = options.get('theta', self.margs.var_to_scale(np.diag(self.cov)))
 
             elif marginal == 'TruncGaussian':
                 lb, ub = np.array(self.bounds).T
                 self.margs = TruncGaussian(lb,ub)
-                self.theta = kwargs_ens.get('theta', np.sqrt(np.diag(self.cov)))
+                self.theta = options.get('theta', np.sqrt(np.diag(self.cov)))
 
             elif marginal == 'Gaussian':
                 self.margs = Gaussian()
-                self.theta = kwargs_ens.get('theta', np.sqrt(np.diag(self.cov)))
+                self.theta = options.get('theta', np.sqrt(np.diag(self.cov)))
     
     def get_theta(self):
         return self.theta
diff --git a/popt/loop/optimize.py b/popt/loop/optimize.py
index 43f99c61..6fa80a44 100644
--- a/popt/loop/optimize.py
+++ b/popt/loop/optimize.py
@@ -166,48 +166,48 @@ def run_loop(self):
                     self.save()
 
             # Check if max iterations was reached
-            if self.iteration > self.max_iter:
+            if self.iteration >= self.max_iter:
                 self.optimize_result['message'] = 'Iterations stopped due to max iterations reached!'
             else:
                 if not isinstance(self.msg, str): self.msg = ''
                 self.optimize_result['message'] = self.msg
 
-                # Logging some info to screen
-                logger.info('       Optimization converged in %d iterations ', self.iteration-1)
-                logger.info('       Optimization converged with final obj_func = %.4f',
-                            np.mean(self.optimize_result['fun']))
-                logger.info('       Total number of function evaluations = %d', self.optimize_result['nfev'])
-                logger.info('       Total number of jacobi evaluations = %d', self.optimize_result['njev'])
-                if self.start_time is not None:
-                    logger.info('       Total elapsed time = %.2f minutes', (time.perf_counter()-self.start_time)/60)
-                logger.info('       ============================================')
-
-                # Test for convergence of outer epf loop
-                epf_not_converged = False
-                if self.epf:
-                    if self.epf_iteration > self.epf['max_epf_iter']:  # max epf_iterations set to 10
-                        logger.info(f'       -----> EPF-EnOpt: maximum epf iterations reached')  # print epf info
-                        break
-                    p = np.abs(previous_state-self.mean_state) / (np.abs(previous_state) + 1.0e-9)
-                    conv_crit = self.epf['conv_crit']
-                    if np.any(p > conv_crit):
-                        epf_not_converged = True
-                        previous_state = self.mean_state
-                        self.epf['r'] *= self.epf['r_factor']  # increase penalty factor
-                        self.obj_func_tol *= self.epf['tol_factor']  # decrease tolerance
-                        self.obj_func_values = self.fun(self.mean_state, **self.epf)
-                        self.iteration = 0
-                        self.epf_iteration += 1
-                        optimize_result = ot.get_optimize_result(self)
-                        ot.save_optimize_results(optimize_result)
-                        self.nfev += 1
-                        self.iteration = +1
-                        r = self.epf['r']
-                        logger.info(f'       -----> EPF-EnOpt: {self.epf_iteration}, {r} (outer iteration, penalty factor)')  # print epf info
-                    else:
-                        logger.info(f'       -----> EPF-EnOpt: converged, no variables changed more than {conv_crit*100} %')  # print epf info
-                        final_obj_no_penalty = str(round(float(self.fun(self.mean_state)),4))
-                        logger.info(f'       -----> EPF-EnOpt: objective value without penalty = {final_obj_no_penalty}') # print epf info
+            # Logging some info to screen
+            logger.info('       Optimization converged in %d iterations ', self.iteration-1)
+            logger.info('       Optimization converged with final obj_func = %.4f',
+                        np.mean(self.optimize_result['fun']))
+            logger.info('       Total number of function evaluations = %d', self.optimize_result['nfev'])
+            logger.info('       Total number of jacobi evaluations = %d', self.optimize_result['njev'])
+            if self.start_time is not None:
+                logger.info('       Total elapsed time = %.2f minutes', (time.perf_counter()-self.start_time)/60)
+            logger.info('       ============================================')
+
+            # Test for convergence of outer epf loop
+            epf_not_converged = False
+            if self.epf:
+                if self.epf_iteration > self.epf['max_epf_iter']:  # max epf_iterations set to 10
+                    logger.info(f'       -----> EPF-EnOpt: maximum epf iterations reached')  # print epf info
+                    break
+                p = np.abs(previous_state-self.mean_state) / (np.abs(previous_state) + 1.0e-9)
+                conv_crit = self.epf['conv_crit']
+                if np.any(p > conv_crit):
+                    epf_not_converged = True
+                    previous_state = self.mean_state
+                    self.epf['r'] *= self.epf['r_factor']  # increase penalty factor
+                    self.obj_func_tol *= self.epf['tol_factor']  # decrease tolerance
+                    self.obj_func_values = self.fun(self.mean_state, **self.epf)
+                    self.iteration = 0
+                    self.epf_iteration += 1
+                    optimize_result = ot.get_optimize_result(self)
+                    ot.save_optimize_results(optimize_result)
+                    self.nfev += 1
+                    self.iteration = +1
+                    r = self.epf['r']
+                    logger.info(f'       -----> EPF-EnOpt: {self.epf_iteration}, {r} (outer iteration, penalty factor)')  # print epf info
+                else:
+                    logger.info(f'       -----> EPF-EnOpt: converged, no variables changed more than {conv_crit*100} %')  # print epf info
+                    final_obj_no_penalty = str(round(float(self.fun(self.mean_state)),4))
+                    logger.info(f'       -----> EPF-EnOpt: objective value without penalty = {final_obj_no_penalty}') # print epf info
 
     def save(self):
         """
diff --git a/popt/update_schemes/linesearch.py b/popt/update_schemes/linesearch.py
index b8fe92f3..7bf1081f 100644
--- a/popt/update_schemes/linesearch.py
+++ b/popt/update_schemes/linesearch.py
@@ -13,6 +13,7 @@
 # Internal imports
 from popt.misc_tools import optim_tools as ot
 from popt.loop.optimize import Optimize
+from popt.update_schemes import optimizers
 
 def LineSearch(fun, x, jac, method='GD', hess=None, args=(), bounds=None, callback=None, **options):
     '''
@@ -373,6 +374,13 @@ def calc_update(self, iter_resamp=0):
             pk = - np.matmul(self.Hk_inv, self.jk)
         if self.method == 'Newton':
             pk = - np.matmul(la.inv(self.Hk), self.jk)
+
+        # remove components that point out of the hybercube given by [lb,ub]
+        lb = np.array(self.bounds)[:, 0]
+        ub = np.array(self.bounds)[:, 1]
+        for i in range(self.xk.size):
+            if (self.xk[i] <= lb[i] and pk[i] < 0) or (self.xk[i] >= ub[i] and pk[i] > 0):
+                pk[i] = 0
         
         # Set step_size
         step_size = self._set_step_size(pk)
diff --git a/popt/update_schemes/trust_region.py b/popt/update_schemes/trust_region.py
index 5272d750..cd616f93 100644
--- a/popt/update_schemes/trust_region.py
+++ b/popt/update_schemes/trust_region.py
@@ -11,8 +11,12 @@
 from popt.misc_tools import optim_tools as ot
 from popt.loop.optimize import Optimize
 
+# Impors from scipy
+from scipy.optimize._trustregion_ncg import CGSteihaugSubproblem
+from scipy.optimize._trustregion_exact import IterativeSubproblem
 
-def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **options):
+
+def TrustRegion(fun, x, jac, hess, method='iterative', args=(), bounds=None, callback=None, **options):
     '''
     Trust region optimization algorithm.
 
@@ -29,6 +33,10 @@ def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **option
     
     hess : callable
         Hessian of objective function. The calling signature is `hess(x, *args)`.
+
+    method : str, optional
+        Method to use for solving the trust-region subproblem. Options are 'iterative' or 'CG-Steihaug'.
+        Default is 'iterative'.
     
     args : tuple, optional
         Extra arguments passed to the objective function and its derivatives (Jacobian, Hessian).
@@ -36,9 +44,10 @@ def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **option
     bounds : sequence, optional
         Bounds for variables. Each element of the sequence must be a tuple of two scalars, 
         representing the lower and upper bounds for that variable. Use None for one of the bounds if there are no bounds.
+        Bounds are handle by clipping the state to the bounds before evaluating the objective function and its derivatives.
     
     callback: callable, optional
-        A callable called after each successful iteration. The class instance of LineSearch 
+        A callable called after each successful iteration. The class instance 
         is passed as the only argument to the callback function: callback(self) 
 
     **options : keyword arguments, optional
@@ -58,6 +67,9 @@ def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **option
         Minimum trust-region radius. Optimization is terminated if trust_radius = trust_radius_min.
         Default is trust_radius/100.
     
+    trust_radius_cuts: int
+        Number of allowed trust-region radius reductions if a step is not successful. Default is 4.
+    
     rho_tol: float
         Tolerance for rho (ratio of actual to predicted reduction). Default is 1e-6.
         
@@ -73,6 +85,13 @@ def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **option
             eta2 = 0.1 \n
             gam1 = 0.7 \n
             gam2 = 1.5 \n
+    
+    saveit: bool
+        If True, save the optimization results to a file. Default is True.
+    
+    convergence_criteria: callable
+        A callable that takes the current optimization object as an argument and returns True if the optimization should stop.
+        It can be used to implement custom convergence criteria. Default is None.
 
     save_folder: str
         Name of folder to save the results to. Defaul is ./ (the current directory).
@@ -86,9 +105,9 @@ def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **option
     hess0: ndarray
         Hessian value of the initial control.
 
-    resample: int
-        Number of jacobian re-computations allowed if a line search fails. Default is 4.
-        (useful if jacobian is stochastic)
+    resample: bool
+        If True, resample the Jacobian and Hessian if a step is not successful. Default is False.
+        (Only makes sense if the Jacobian and Hessian are stochastic).
 
     savedata: list[str]
         Further specification of which class variables to save to the result files.
@@ -110,12 +129,12 @@ def TrustRegion(fun, x, jac, hess, args=(), bounds=None, callback=None, **option
         - nfev: number of function evaluations
         - njev: number of jacobian evaluations
     '''
-    tr_obj = TrustRegionClass(fun, x, jac, hess, args, bounds, callback, **options)
+    tr_obj = TrustRegionClass(fun, x, jac, hess, method, args, bounds, callback, **options)
     return tr_obj.optimize_result
 
 class TrustRegionClass(Optimize):
 
-    def __init__(self, fun, x, jac, hess, args=(), bounds=None, callback=None, **options):
+    def __init__(self, fun, x, jac, hess, method='iterative', args=(), bounds=None, callback=None, **options):
         
         # Initialize the parent class
         super().__init__(**options)
@@ -125,6 +144,7 @@ def __init__(self, fun, x, jac, hess, args=(), bounds=None, callback=None, **opt
         self.xk       = x
         self.jacobian = jac
         self.hessian  = hess
+        self.method   = method
         self.args     = args
         self.bounds   = bounds
         self.options  = options
@@ -135,21 +155,35 @@ def __init__(self, fun, x, jac, hess, args=(), bounds=None, callback=None, **opt
         else:
             self.callback = None
 
+        # Custom convergence criteria (callable)
+        convergence_criteria = options.get('convergence_criteria', None)
+        if callable(convergence_criteria):
+            self.convergence_criteria = self.convergence_criteria
+        else:
+            self.convergence_criteria = None
+
         # Set options for trust-region radius       
-        self.trust_radius     = options.get('trust_radius', 1.0) 
-        self.trust_radius_max = options.get('trust_radius_max', 10*self.trust_radius)
-        self.trust_radius_min = options.get('trust_radius_min', self.trust_radius/100)
+        self.trust_radius      = options.get('trust_radius', 1.0) 
+        self.trust_radius_max  = options.get('trust_radius_max', 10*self.trust_radius)
+        self.trust_radius_min  = options.get('trust_radius_min', self.trust_radius/100)
+        self.trust_radius_cuts = options.get('trust_radius_cuts', 4)
 
         # Set other options
-        self.resample = options.get('resample', 3)
+        self.resample = options.get('resample', False)
         self.saveit   = options.get('saveit', True)
         self.rho_tol  = options.get('rho_tol', 1e-6)
-        self.eta1 = options.get('eta1', 0.001)
-        self.eta2 = options.get('eta2', 0.1)
-        self.gam1 = options.get('gam1', 0.7)
-        self.gam2 = options.get('gam2', 1.5)
+        self.eta1 = options.get('eta1', 0.1) # reduce raduis if rho < 10%
+        self.eta2 = options.get('eta2', 0.5)  # increase radius if rho > 50% 
+        self.gam1 = options.get('gam1', 0.5)  # reduce by 50%
+        self.gam2 = options.get('gam2', 1.5)  # increase by 50%
         self.rho  = 0.0
 
+        # Check if method is valid
+        if self.method not in ['iterative', 'CG-Steihaug']:
+            self.method = 'iterative'
+            raise ValueError(f'Method {self.method} is not valid!. Method is set to "iterative"')
+    
+            
         if not self.restart:
             self.start_time = time.perf_counter()
 
@@ -207,15 +241,17 @@ def _hess(self, x):
         return h
     
     def update_results(self):
-        res = {'fun': self.fk, 
-               'x': self.xk, 
-               'jac': self.jk,
-               'hess': self.Hk,
-               'nfev': self.nfev,
-               'njev': self.njev,
-               'nit': self.iteration,
-               'trust_radius': self.trust_radius,
-               'save_folder': self.options.get('save_folder', './')}
+        res = {
+            'fun': self.fk, 
+            'x': self.xk, 
+            'jac': self.jk,
+            'hess': self.Hk,
+            'nfev': self.nfev,
+            'njev': self.njev,
+            'nit': self.iteration,
+            'trust_radius': self.trust_radius,
+            'save_folder': self.options.get('save_folder', './')
+        }
         
         for a, arg in enumerate(self.args):
             res[f'args[{a}]'] = arg
@@ -242,22 +278,71 @@ def _log(self, msg):
         if self.logger is not None:
             self.logger.info(msg)
 
-    def calc_update(self, iter_resamp=0):
+    def solve_subproblem(self, g, B, delta):
+        """
+        Solve the trust region subproblem using the iterative method.
+        (A big thanks to copilot for the help with this implementation)
+
+        Parameters:
+        g (numpy.ndarray): Gradient vector at the current point.
+        B (numpy.ndarray): Hessian matrix at the current point.
+        delta (float): Trust region radius.
+
+        Returns:
+        pk (numpy.ndarray): Step direction.
+        pk_hits_boundary (bool): True if the step hits the boundary of the trust region.
+        """
+
+        # Define quadratic model
+        quad = lambda p: self.fk + np.dot(g,p) + np.dot(p,np.dot(B,p))/2
+
+
+        if self.method == 'iterative':
+            subproblem = IterativeSubproblem(
+                x=self.xk, 
+                fun=quad, 
+                jac=lambda _: g, 
+                hess=lambda _: B,
+            )
+            pk, pk_hits_boundary = subproblem.solve(tr_radius=delta)
+        
+        elif self.method == 'CG-Steihaug':
+            subproblem = CGSteihaugSubproblem(
+                x=self.xk, 
+                fun=quad, 
+                jac=lambda _: g, 
+                hess=lambda _: B,
+            )
+            pk, pk_hits_boundary = subproblem.solve(trust_radius=delta)
+
+        else:
+            raise ValueError(f"Method {self.method} is not valid!")
+
+        return pk, pk_hits_boundary
+
+
+    def calc_update(self, inner_iter=0):
 
         # Initialize variables for this step
         success = True
 
         # Solve subproblem
-        self._log('Solving trust region subproblem using the CG-Steihaug method')
-        sk = self.solve_sub_problem_CG_Steihaug(self.jk, self.Hk, self.trust_radius)
+        self._log('Solving trust region subproblem')
+        sk, hits_boundary = self.solve_subproblem(self.jk, self.Hk, self.trust_radius)
+
+        # truncate sk to respect bounds
+        if self.bounds is not None:
+            lb = np.array(self.bounds)[:, 0]
+            ub = np.array(self.bounds)[:, 1]
+            sk = np.clip(sk, lb - self.xk, ub - self.xk)
 
         # Calculate the actual function value
-        xk_new = ot.clip_state(self.xk + sk, self.bounds)
+        xk_new  = self.xk + sk
         fun_new = self._fun(xk_new)
 
         # Calculate rho
         actual_reduction    = self.fk - fun_new
-        predicted_reduction = - np.dot(self.jk, sk) - 0.5*np.dot(sk, np.dot(self.Hk, sk))
+        predicted_reduction = - np.dot(self.jk, sk) - np.dot(sk, np.dot(self.Hk, sk))/2
         self.rho = actual_reduction/predicted_reduction
 
         if self.rho > self.rho_tol:
@@ -283,7 +368,7 @@ def calc_update(self, iter_resamp=0):
 
             # update the trust region radius
             delta_old = self.trust_radius
-            if self.rho >= self.eta2:
+            if (self.rho >= self.eta2) and hits_boundary:
                 delta_new = min(self.gam2*delta_old, self.trust_radius_max)
             elif self.eta1 <= self.rho < self.eta2:
                 delta_new = delta_old
@@ -291,12 +376,19 @@ def calc_update(self, iter_resamp=0):
                 delta_new = self.gam1*delta_old
             
             # Log new trust-radius
-            self.trust_radius = delta_new 
+            self.trust_radius = np.clip(delta_new, self.trust_radius_min, self.trust_radius_max)
             if not (delta_old == delta_new):
                 self._log(f'Trust-radius updated: {delta_old:<10.4e} --> {delta_new:<10.4e}')
 
+            # Check for custom convergence
+            if callable(self.convergence_criteria):
+                if self.convergence_criteria(self):
+                    self._log('Custom convergence criteria met. Stopping optimization.')
+                    success = False
+                    return success
+
             # check for convergence
-            if (self.trust_radius < self.trust_radius_min) or (self.iteration==self.max_iter):
+            if self.iteration==self.max_iter:
                 success = False
             else:
                 # Calculate the jacobian and hessian
@@ -307,101 +399,34 @@ def calc_update(self, iter_resamp=0):
             self.iteration += 1
 
         else:
-            if iter_resamp < self.resample:
-
-                iter_resamp += 1
+            if inner_iter < self.trust_radius_cuts:
+                
+                # Log the failure
+                self._log(f'Step not successful: rho < {self.rho_tol:<10.4e}')
+                
+                # Reduce trust region radius to 75% of current value
+                self._log('Reducing trust-radius by 75%')
+                self.trust_radius = 0.25*self.trust_radius
 
-                # Calculate the jacobian and hessian
-                self._log('Resampling gradient and hessian')
-                self.jk = self._jac(self.xk)
-                self.Hk = self._hess(self.xk)
+                if self.trust_radius < self.trust_radius_min:
+                    self._log(f'Trust radius {self.trust_radius} is below minimum {self.trust_radius_min}. Stopping optimization.')
+                    success = False
+                    return success
 
-                # Reduce trust region radius to 50% of current value
-                self._log('Reducing trust-radius by 50%')
-                self.trust_radius = 0.5*self.trust_radius
+                # Check for resampling of Jac and Hess
+                if self.resample:
+                    self._log('Resampling gradient and hessian')
+                    self.jk = self._jac(self.xk)
+                    self.Hk = self._hess(self.xk)
 
                 # Recursivly call function
-                success = self.calc_update(iter_resamp=iter_resamp)
+                success = self.calc_update(inner_iter=inner_iter+1)
 
             else:
                 success = False
 
         return success
 
-    
-    def solve_sub_problem_CG_Steihaug(self, g, B, delta):
-        """
-        Solve the trust region subproblem using Steihaug's Conjugate Gradient method.
-        (A big thanks to copilot for the help with this implementation)
-
-        Parameters:
-        g (numpy.ndarray): Gradient vector at the current point.
-        B (numpy.ndarray): Hessian matrix at the current point.
-        delta (float): Trust region radius.
-        tol (float): Tolerance for convergence.
-        max_iter (int): Maximum number of iterations.
-
-        Returns:
-        p (numpy.ndarray): Solution vector.
-        """
-        z = np.zeros_like(g)
-        r = g
-        d = -g
-
-        # Set same default tolerance as scipy
-        tol = min(0.5, la.norm(g)**2)*la.norm(g)
-
-        if la.norm(g) <= tol:
-            return z
-
-        # make quadratic model
-        mc = lambda s: self.fk + np.dot(g,s) + np.dot(s,np.dot(B,s))/2
-
-        while True:
-            dBd = np.dot(d, np.dot(B,d))
-
-            if dBd <= 0:
-                # Solve the quadratic equation: (p + tau*d)**2 = delta**2
-                tau_lo, tau_hi = self.get_tau_at_delta(z, d, delta)
-                p_lo = z + tau_lo*d
-                p_hi = z + tau_hi*d
-
-                if mc(p_lo) < mc(p_hi): 
-                    return p_lo
-                else:
-                    return p_hi
-
-            alpha = np.dot(r,r)/dBd
-            z_new = z + alpha*d
-
-            if la.norm(z_new) >= delta:
-                # Solve the quadratic equation: (p + tau*d)**2 = delta**2, for tau > 0
-                _ , tau = self.get_tau_at_delta(z, d, delta)
-                return z + tau * d
-                
-            r_new = r + alpha*np.dot(B,d)
-
-            if la.norm(r_new) < tol:
-                return z_new
-
-            beta = np.dot(r_new,r_new)/np.dot(r,r)
-            d = -r_new + beta*d
-            r = r_new
-            z = z_new
-
-    
-    def get_tau_at_delta(self, p, d, delta):
-        """
-        Solve the quadratic equation: (p + tau*d)**2 = delta**2, for tau > 0
-        """
-        a = np.dot(d,d)
-        b = 2*np.dot(p,d)
-        c = np.dot(p,p) - delta**2
-        tau_lo = -b/(2*a) - np.sqrt(b**2 - 4*a*c)/(2*a) 
-        tau_hi = -b/(2*a) + np.sqrt(b**2 - 4*a*c)/(2*a) 
-        return tau_lo, tau_hi
-             
-