diff --git a/ensemble/ensemble.py b/ensemble/ensemble.py
index 304b2404..d4dc46a0 100644
--- a/ensemble/ensemble.py
+++ b/ensemble/ensemble.py
@@ -80,7 +80,7 @@ def __init__(self, keys_en, sim, redund_sim=None):
 
         # If it is a restart run, we do not need to initialize anything, only load the self info. that exists in the
         # pickle save file. If it is not a restart run, we initialize everything below.
-        if 'restart' in self.keys_en and self.keys_en['restart'] == 'yes':
+        if ('restart' in self.keys_en) and (self.keys_en['restart'] == 'yes'):
             # Initiate a restart run
             self.logger.info('\033[92m--- Restart run initiated! ---\033[92m')
             # Check if the pickle save file exists in folder
@@ -117,7 +117,7 @@ def __init__(self, keys_en, sim, redund_sim=None):
                 self.disable_tqdm = False
 
             # extract information that is given for the prior model
-            self._ext_prior_info()
+            self.prior_info = self._extract_prior_info()
 
             # Calculate initial ensemble if IMPORTSTATICVAR has not been given in init. file.
             # Prior info. on state variables must be given by PRIOR_<STATICVAR-name> keyword.
@@ -172,243 +172,116 @@ def _ext_ml_info(self):
                         self.error_comp_scheme = self.keys_en['multilevel'][i][2]
                     self.ML_corr_done = False
 
-    def _ext_prior_info(self):
-        """
+    def _extract_prior_info(self) -> dict:
+        '''
         Extract prior information on STATE from keyword(s) PRIOR_<STATE entries>.
-        """
-        # Parse prior info on each state entered in STATE.
-        # Store names given in STATE
-        if not isinstance(self.keys_en['state'], list):  # Single string
-            state_names = [self.keys_en['state']]
-        else:  # List
-            state_names = self.keys_en['state']
-
-        # Check if PRIOR_<state names> exists for each entry in STATE
+        '''
+
+        # Get state names as list
+        state_names = self.keys_en['state']
+        if not isinstance(state_names, list): state_names = [state_names]
+
+        # Check if PRIOR_<state names> exists for each entry in state
         for name in state_names:
-            assert 'prior_' + name in self.keys_en, \
+            assert f'prior_{name}' in self.keys_en, \
                 'PRIOR_{0} is missing! This keyword is needed to make initial ensemble for {0} entered in ' \
                 'STATE'.format(name.upper())
+        
+        # define dict to store prior information in 
+        prior_info = {name: None for name in state_names}
 
-        # Init. prior info variable
-        self.prior_info = {keys: None for keys in state_names}
-
-        # Loop over each prior keyword and make an initial. ensemble for each state in STATE,
-        # which is subsequently stored in the state dictionary. If 3D grid dimensions are inputted, information for
-        # each layer must be inputted, else the single information will be copied to all layers.
-        grid_dim = np.array([0, 0])
+        # loop over state priors
         for name in state_names:
-            # initiallize an empty dictionary inside the dictionary.
-            self.prior_info[name] = {}
-            # List the option names inputted in prior keyword
-            # opt_list = list(zip(*self.keys_da['prior_'+name]))
-            mean = None
-            self.prior_info[name]['mean'] = mean
-            vario = [None]
-            self.prior_info[name]['vario'] = vario
-            aniso = [None]
-            self.prior_info[name]['aniso'] = aniso
-            angle = [None]
-            self.prior_info[name]['angle'] = angle
-            corr_length = [None]
-            self.prior_info[name]['corr_length'] = corr_length
-            self.prior_info[name]['nx'] = self.prior_info[name]['ny'] = self.prior_info[name]['nz'] = None
-
-            # Extract info. from the prior keyword
-            for i, opt in enumerate(list(zip(*self.keys_en['prior_' + name]))[0]):
-                if opt == 'vario':  # Variogram model
-                    if not isinstance(self.keys_en['prior_' + name][i][1], list):
-                        vario = [self.keys_en['prior_' + name][i][1]]
-                    else:
-                        vario = self.keys_en['prior_' + name][i][1]
-                elif opt == 'mean':  # Mean
-                    mean = self.keys_en['prior_' + name][i][1]
-                elif opt == 'var':  # Variance
-                    if not isinstance(self.keys_en['prior_' + name][i][1], list):
-                        variance = [self.keys_en['prior_' + name][i][1]]
-                    else:
-                        variance = self.keys_en['prior_' + name][i][1]
-                elif opt == 'aniso':  # Anisotropy factor
-                    if not isinstance(self.keys_en['prior_' + name][i][1], list):
-                        aniso = [self.keys_en['prior_' + name][i][1]]
-                    else:
-                        aniso = self.keys_en['prior_' + name][i][1]
-                elif opt == 'angle':  # Anisotropy angle
-                    if not isinstance(self.keys_en['prior_' + name][i][1], list):
-                        angle = [self.keys_en['prior_' + name][i][1]]
-                    else:
-                        angle = self.keys_en['prior_' + name][i][1]
-                elif opt == 'range':  # Correlation length
-                    if not isinstance(self.keys_en['prior_' + name][i][1], list):
-                        corr_length = [self.keys_en['prior_' + name][i][1]]
+            prior = self.keys_en[f'prior_{name}']
+            
+            # Check if is a list (old way)
+            if isinstance(prior, list):
+                # list of lists - old way of inputting prior information
+                prior_dict = {}
+                for i, opt in enumerate(list(zip(*prior))[0]):
+                    if opt == 'limits':
+                        prior_dict[opt] = prior[i][1:]
                     else:
-                        corr_length = self.keys_en['prior_' + name][i][1]
-                elif opt == 'grid':  # Grid dimensions
-                    grid_dim = self.keys_en['prior_' + name][i][1]
-                elif opt == 'limits':  # Truncation values
-                    limits = self.keys_en['prior_' + name][i][1:]
-                elif opt == 'active':  # Number of active cells (single number)
-                    active = self.keys_en['prior_' + name][i][1]
-
-            # Check if mean needs to be loaded, or if loaded
-            if type(mean) is str:
-                assert mean.endswith('.npz'), 'File name does not end with \'.npz\'!'
-                load_file = np.load(mean)
+                        prior_dict[opt] = prior[i][1]
+                prior = prior_dict
+            else:
+                assert isinstance(prior, dict), 'PRIOR_{0} must be a dictionary or list of lists!'.format(name.upper())
+
+
+            # load mean if in file
+            if isinstance(prior['mean'], str):
+                assert prior['mean'].endswith('.npz'), 'File name does not end with \'.npz\'!'
+                load_file = np.load(prior['mean'])
                 assert len(load_file.files) == 1, \
                     'More than one variable located in {0}. Only the mean vector can be stored in the .npz file!' \
-                    .format(mean)
-                mean = load_file[load_file.files[0]]
+                    .format(prior['mean'])
+                prior['mean'] = load_file[load_file.files[0]]
             else:  # Single number inputted, make it a list if not already
-                if not isinstance(mean, list):
-                    mean = [mean]
-
-            # Check if limits exists
-            try:
-                limits
-            except NameError:
-                limits = None
-
-            # check if active exists
-            try:
-                active
-            except NameError:
-                active = None
-
-            # Extract x- and y-dim
-            nx = int(grid_dim[0])
-            ny = int(grid_dim[1])
-
-            # Check if 3D grid inputted. If so, we check if info. has been given on all layers. In the case it has
-            # not been given, we just copy the info. given.
-            if len(grid_dim) == 3 and grid_dim[2] > 1:  # 3D
-                nz = int(grid_dim[2])
-
-                # Check mean when values have been inputted directly (not when mean has been loaded)
-                if isinstance(mean, list) and len(mean) < nz:
-                    # Check if it is more than one entry and give error
-                    assert len(mean) == 1, \
-                        'Information from MEAN has been given for {0} layers, whereas {1} is needed!' \
-                        .format(len(mean), nz)
-
-                    # Only 1 entry; copy this to all layers
-                    print(
-                        '\033[1;33mSingle entry for MEAN will be copied to all {0} layers\033[1;m'.format(nz))
-                    self.prior_info[name]['mean'] = mean * nz
-
-                else:
-                    self.prior_info[name]['mean'] = mean
-
-                # Check variogram model
-                if len(vario) < nz:
-                    # Check if it is more than one entry and give error
-                    assert len(vario) == 1, \
-                        'Information from VARIO has been given for {0} layers, whereas {1} is needed!' \
-                        .format(len(vario), nz)
-
-                    # Only 1 entry; copy this to all layers
-                    print(
-                        '\033[1;33mSingle entry for VARIO will be copied to all {0} layers\033[1;m'.format(nz))
-                    self.prior_info[name]['vario'] = vario * nz
-
-                else:
-                    self.prior_info[name]['vario'] = vario
-
-                # Variance
-                if len(variance) < nz:
-                    # Check if it is more than one entry and give error
-                    assert len(variance) == 1, \
-                        'Information from VAR has been given for {0} layers, whereas {1} is needed!' \
-                        .format(len(variance), nz)
-
-                    # Only 1 entry; copy this to all layers
-                    print(
-                        '\033[1;33mSingle entry for VAR will be copied to all {0} layers\033[1;m'.format(nz))
-                    self.prior_info[name]['variance'] = variance * nz
-
-                else:
-                    self.prior_info[name]['variance'] = variance
-
-                # Aniso factor
-                if len(aniso) < nz:
-                    # Check if it is more than one entry and give error
-                    assert len(aniso) == 1, \
-                        'Information from ANISO has been given for {0} layers, whereas {1} is needed!' \
-                        .format(len(aniso), nz)
-
-                    # Only 1 entry; copy this to all layers
-                    print(
-                        '\033[1;33mSingle entry for ANISO will be copied to all {0} layers\033[1;m'.format(nz))
-                    self.prior_info[name]['aniso'] = aniso * nz
-
-                else:
-                    self.prior_info[name]['aniso'] = aniso
-
-                # Aniso factor
-                if len(angle) < nz:
-                    # Check if it is more than one entry and give error
-                    assert len(angle) == 1, \
-                        'Information from ANGLE has been given for {0} layers, whereas {1} is needed!' \
-                        .format(len(angle), nz)
-
-                    # Only 1 entry; copy this to all layers
-                    print(
-                        '\033[1;33mSingle entry for ANGLE will be copied to all {0} layers\033[1;m'.format(nz))
-                    self.prior_info[name]['angle'] = angle * nz
+                if not isinstance(prior['mean'], list):
+                    prior['mean'] = [prior['mean']]
+
+            # loop over keys in prior
+            for key in prior.keys():
+                # ensure that entry is a list
+                if (not isinstance(prior[key], list)) and (key != 'mean'):
+                    prior[key] = [prior[key]]
+
+            # change the name of some keys
+            prior['variance'] = prior.pop('var', None)
+            prior['corr_length'] = prior.pop('range', None)
+
+            # process grid
+            if 'grid' in prior:
+                grid_dim = prior['grid']
+
+                # check if 3D-grid
+                if (len(grid_dim) == 3) and (grid_dim[2] > 1):
+                    nz = int(grid_dim[2])
+                    prior['nz'] = nz
+                    prior['nx'] = int(grid_dim[0])
+                    prior['ny'] = int(grid_dim[1])
+                    
+
+                    # Check mean when values have been inputted directly (not when mean has been loaded)
+                    mean = prior['mean']
+                    if isinstance(mean, list) and len(mean) < nz:
+                         # Check if it is more than one entry and give error
+                        assert len(mean) == 1, \
+                            'Information from MEAN has been given for {0} layers, whereas {1} is needed!' \
+                            .format(len(mean), nz)
+
+                        # Only 1 entry; copy this to all layers
+                        print(
+                            '\033[1;33mSingle entry for MEAN will be copied to all {0} layers\033[1;m'.format(nz))
+                        prior['mean'] = mean * nz
+
+                    #check if info. has been given on all layers. In the case it has not been given, we just copy the info. given.
+                    for key in ['vario', 'variance', 'aniso', 'angle', 'corr_length']:
+                        if key in prior.keys():
+                            val = prior[key]
+                            if len(val) < nz:
+                                # Check if it is more than one entry and give error
+                                assert len(val) == 1, \
+                                    'Information from {0} has been given for {1} layers, whereas {2} is needed!' \
+                                    .format(key.upper(), len(val), nz)
+
+                                # Only 1 entry; copy this to all layers
+                                print(
+                                    '\033[1;33mSingle entry for {0} will be copied to all {1} layers\033[1;m'.format(key.upper(), nz))
+                                prior[key] = val * nz
 
                 else:
-                    self.prior_info[name]['angle'] = angle
+                    prior['nx'] = int(grid_dim[0])
+                    prior['ny'] = int(grid_dim[1])
+                    prior['nz'] = 1
 
-                # Corr. length
-                if len(corr_length) < nz:
-                    # Check if it is more than one entry and give error
-                    assert len(corr_length) == 1, \
-                        'Information from RANGE has been given for {0} layers, whereas {1} is needed!' \
-                        .format(len(corr_length), nz)
+            prior.pop('grid', None)
 
-                    # Only 1 entry; copy this to all layers
-                    print(
-                        '\033[1;33mSingle entry for RANGE will be copied to all {0} layers\033[1;m'.format(nz))
-                    self.prior_info[name]['corr_length'] = corr_length * nz
-
-                else:
-                    self.prior_info[name]['corr_length'] = corr_length
-
-                # Limits, if exists
-                if limits is not None:
-                    self.prior_info[name]['limits'] = limits
-
-                    # if isinstance(limits[0], list) and len(limits) < nz or \
-                    #         not isinstance(limits[0], list) and len(limits) < 2 * nz:
-                    #     # Check if it is more than one entry and give error
-                    #     assert (isinstance(limits[0], list) and len(limits) == 1), \
-                    #         'Information from LIMITS has been given for {0} layers, whereas {1} is needed!' \
-                    #         .format(len(limits), nz)
-                    #     assert (not isinstance(limits[0], list) and len(limits) == 2), \
-                    #         'Information from LIMITS has been given for {0} layers, whereas {1} is needed!' \
-                    #         .format(len(limits) / 2, nz)
-                    #
-                    #     # Only 1 entry; copy this to all layers
-                    #     print(
-                    #         '\033[1;33mSingle entry for RANGE will be copied to all {0} layers\033[1;m'.format(nz))
-                    #     self.prior_info[name]['limits'] = [limits] * nz
-
-            else:  # 2D grid only, or optimization case
-                nz = 1
-                self.prior_info[name]['mean'] = mean
-                self.prior_info[name]['vario'] = vario
-                self.prior_info[name]['variance'] = variance
-                self.prior_info[name]['aniso'] = aniso
-                self.prior_info[name]['angle'] = angle
-                self.prior_info[name]['corr_length'] = corr_length
-                if limits is not None:
-                    self.prior_info[name]['limits'] = limits
-                if active is not None:
-                    self.prior_info[name]['active'] = active
-
-            self.prior_info[name]['nx'] = nx
-            self.prior_info[name]['ny'] = ny
-            self.prior_info[name]['nz'] = nz
-
-        # Loop over keys and input
+            # add prior to prior_info
+            prior_info[name] = prior
+            
+        return prior_info
+                
 
     def gen_init_ensemble(self):
         """
@@ -427,21 +300,21 @@ def gen_init_ensemble(self):
             ind_end = 0
 
             # Extract info.
-            nz = self.prior_info[name]['nz']
-            mean = self.prior_info[name]['mean']
-            nx = self.prior_info[name]['nx']
-            ny = self.prior_info[name]['ny']
+            nx = self.prior_info[name].get('nx', 0)
+            ny = self.prior_info[name].get('ny', 0)
+            nz = self.prior_info[name].get('nz', 0)
+            mean = self.prior_info[name].get('mean', None)
+
             if nx == ny == 0:  # assume ensemble will be generated elsewhere if dimensions are zero
                 break
-            variance = self.prior_info[name]['variance']
-            corr_length = self.prior_info[name]['corr_length']
-            aniso = self.prior_info[name]['aniso']
-            vario = self.prior_info[name]['vario']
-            angle = self.prior_info[name]['angle']
-            if 'limits' in self.prior_info[name]:
-                limits = self.prior_info[name]['limits']
-            else:
-                limits = None
+
+            variance = self.prior_info[name].get('variance', None)
+            corr_length = self.prior_info[name].get('corr_length', None)
+            aniso = self.prior_info[name].get('aniso', None)
+            vario = self.prior_info[name].get('vario', None)
+            angle = self.prior_info[name].get('angle', None)
+            limits= self.prior_info[name].get('limits',None)
+            
 
             # Loop over nz to make layers of 2D priors
             for i in range(self.prior_info[name]['nz']):
diff --git a/input_output/organize.py b/input_output/organize.py
index 3accf029..8e500e38 100644
--- a/input_output/organize.py
+++ b/input_output/organize.py
@@ -3,6 +3,7 @@
 from copy import deepcopy
 import csv
 import datetime as dt
+import pandas as pd
 
 
 class Organize_input():
@@ -109,6 +110,13 @@ def _org_report(self):
                         pred_prim.extend(csv_data)
                 self.keys_fwd['reportpoint'] = pred_prim
 
+            elif isinstance(self.keys_fwd['reportpoint'], dict):
+                self.keys_fwd['reportpoint'] = pd.date_range(**self.keys_fwd['reportpoint']).to_pydatetime().tolist()
+
+            else:
+                pass
+
+
         # Check if assimindex is given as a csv file. If so, we read and make a potential 2D list (if sequential).
         if 'assimindex' in self.keys_pr:
             if isinstance(self.keys_pr['assimindex'], str) and self.keys_pr['assimindex'].endswith('.csv'):
diff --git a/input_output/read_config.py b/input_output/read_config.py
index 986e7558..ba00da2b 100644
--- a/input_output/read_config.py
+++ b/input_output/read_config.py
@@ -51,6 +51,12 @@ def ndarray_constructor(loader, node):
         y = yaml.load(fid, Loader=FullLoader)
 
     # Check for dataassim and fwdsim
+    if 'ensemble' in y.keys():
+        keys_en = y['ensemble']
+        check_mand_keywords_en(keys_en)
+    else:
+        keys_en = None
+
     if 'optim' in y.keys():
         keys_pr = y['optim']
         check_mand_keywords_opt(keys_pr)
@@ -59,16 +65,17 @@ def ndarray_constructor(loader, node):
         check_mand_keywords_da(keys_pr)
     else:
         raise KeyError
+    
     if 'fwdsim' in y.keys():
         keys_fwd = y['fwdsim']
     else:
         raise KeyError
 
     # Organize keywords
-    org = Organize_input(keys_pr, keys_fwd)
+    org = Organize_input(keys_pr, keys_fwd, keys_en)
     org.organize()
 
-    return org.get_keys_pr(), org.get_keys_fwd()
+    return org.get_keys_pr(), org.get_keys_fwd(), org.get_keys_en()
 
 
 def convert_txt_to_toml(init_file):
diff --git a/pipt/loop/assimilation.py b/pipt/loop/assimilation.py
index 0ccd6dd8..f5caa72f 100644
--- a/pipt/loop/assimilation.py
+++ b/pipt/loop/assimilation.py
@@ -301,32 +301,20 @@ def _ext_max_iter(self):
         - ST 7/6-16
         """
         if 'iteration' in self.ensemble.keys_da:
-            # Make sure ITERATION is a list
-            if not isinstance(self.ensemble.keys_da['iteration'][0], list):
-                iter_opts = [self.ensemble.keys_da['iteration']]
-            else:
-                iter_opts = self.ensemble.keys_da['iteration']
-
+            iter_opts = dict(self.ensemble.keys_da['iteration'])
             # Check if 'max_iter' has been given; if not, give error (mandatory in ITERATION)
-            assert 'max_iter' in list(
-                zip(*iter_opts))[0], 'MAX_ITER has not been given in ITERATION!'
-
-            # Extract max. iter
-            max_iter = [item[1] for item in iter_opts if item[0] == 'max_iter'][0]
+            try:
+                max_iter = iter_opts['max_iter']
+            except KeyError:
+                raise AssertionError('MAX_ITER has not been given in ITERATION')
 
         elif 'mda' in self.ensemble.keys_da:
-            # Make sure ITERATION is a list
-            if not isinstance(self.ensemble.keys_da['mda'][0], list):
-                iter_opts = [self.ensemble.keys_da['mda']]
-            else:
-                iter_opts = self.ensemble.keys_da['mda']
-
-            # Check if 'max_iter' has been given; if not, give error (mandatory in ITERATION)
-            assert 'tot_assim_steps' in list(
-                zip(*iter_opts))[0], 'TOT_ASSIM_STEPS has not been given in MDA!'
-
-            # Extract max. iter
-            max_iter = [item[1] for item in iter_opts if item[0] == 'tot_assim_steps'][0]
+            iter_opts = dict(self.ensemble.keys_da['mda'])
+            # Check if 'tot_assim_steps' has been given; if not, raise error (mandatory in MDA)
+            try:
+                max_iter = iter_opts['tot_assim_steps']
+            except KeyError:
+                raise AssertionError('TOT_ASSIM_STEPS has not been given in MDA!')
 
         else:
             max_iter = 1
diff --git a/pipt/update_schemes/enrml.py b/pipt/update_schemes/enrml.py
index 632564b6..741b700e 100644
--- a/pipt/update_schemes/enrml.py
+++ b/pipt/update_schemes/enrml.py
@@ -250,36 +250,22 @@ def _ext_iter_param(self):
         file. These parameters include convergence tolerances and parameters for the damping parameter. Default
         values for these parameters have been given here, if they are not provided in ITERATION.
         """
-
-        # Predefine all the default values
-        self.data_misfit_tol = 0.01
-        self.step_tol = 0.01
-        self.lam = 100
-        self.lam_max = 1e10
-        self.lam_min = 0.01
-        self.gamma = 5
-        self.trunc_energy = 0.95
+        try:
+            options = dict(self.keys_da['iteration'])
+        except:
+            options = dict([self.keys_da['iteration']])
+
+        # unpack options
+        self.data_misfit_tol = options.get('data_misfit_tol', 0.01)
+        self.trunc_energy = options.get('energy', 0.95)
+        self.step_tol  = options.get('step_tol', 0.01)
+        self.lam       = options.get('lambda', 100)
+        self.lam_max   = options.get('lambda_max', 1e10)
+        self.lam_min   = options.get('lambda_min', 0.01)
+        self.gamma     = options.get('lambda_factor', 5)
         self.iteration = 0
 
-        # Loop over options in ITERATION and extract the parameters we want
-        for i, opt in enumerate(list(zip(*self.keys_da['iteration']))[0]):
-            if opt == 'data_misfit_tol':
-                self.data_misfit_tol = self.keys_da['iteration'][i][1]
-            if opt == 'step_tol':
-                self.step_tol = self.keys_da['iteration'][i][1]
-            if opt == 'lambda':
-                self.lam = self.keys_da['iteration'][i][1]
-            if opt == 'lambda_max':
-                self.lam_max = self.keys_da['iteration'][i][1]
-            if opt == 'lambda_min':
-                self.lam_min = self.keys_da['iteration'][i][1]
-            if opt == 'lambda_factor':
-                self.gamma = self.keys_da['iteration'][i][1]
-
-        if 'energy' in self.keys_da:
-            # initial energy (Remember to extract this)
-            self.trunc_energy = self.keys_da['energy']
-            if self.trunc_energy > 1:  # ensure that it is given as percentage
+        if self.trunc_energy > 1:  # ensure that it is given as percentage
                 self.trunc_energy /= 100.
 
 
@@ -593,33 +579,19 @@ def _ext_iter_param(self):
         file. These parameters include convergence tolerances and parameters for the damping parameter. Default
         values for these parameters have been given here, if they are not provided in ITERATION.
         """
-
-        # Predefine all the default values
-        self.data_misfit_tol = 0.01
-        self.step_tol = 0.01
-        self.gamma = 0.2
-        self.gamma_max = 0.5
-        self.gamma_factor = 2.5
-        self.trunc_energy = 0.95
-        self.iteration = 0
-
-        # Loop over options in ITERATION and extract the parameters we want
-        for i, opt in enumerate(list(zip(*self.keys_da['iteration']))[0]):
-            if opt == 'data_misfit_tol':
-                self.data_misfit_tol = self.keys_da['iteration'][i][1]
-            if opt == 'step_tol':
-                self.step_tol = self.keys_da['iteration'][i][1]
-            if opt == 'gamma':
-                self.gamma = self.keys_da['iteration'][i][1]
-            if opt == 'gamma_max':
-                self.gamma_max = self.keys_da['iteration'][i][1]
-            if opt == 'gamma_factor':
-                self.gamma_factor = self.keys_da['iteration'][i][1]
-
-        if 'energy' in self.keys_da:
-            # initial energy (Remember to extract this)
-            self.trunc_energy = self.keys_da['energy']
-            if self.trunc_energy > 1:  # ensure that it is given as percentage
+        try:
+            options = dict(self.keys_da['iteration'])
+        except:
+            options = dict([self.keys_da['iteration']])
+        
+        self.data_misfit_tol = options.get('data_misfit_tol', 0.01)
+        self.trunc_energy = options.get('energy', 0.95)
+        self.step_tol     = options.get('step_tol', 0.01)
+        self.gamma        = options.get('gamma', 0.2)
+        self.gamma_max    = options.get('gamma_max', 0.5)
+        self.gamma_factor = options.get('gamma_factor', 2.5)
+
+        if self.trunc_energy > 1:  # ensure that it is given as percentage
                 self.trunc_energy /= 100.
 
 
diff --git a/pipt/update_schemes/esmda.py b/pipt/update_schemes/esmda.py
index b848096c..7cefd164 100644
--- a/pipt/update_schemes/esmda.py
+++ b/pipt/update_schemes/esmda.py
@@ -31,7 +31,7 @@ def __init__(self, keys_da, keys_en, sim):
 
         Parameters
         ----------
-        keys_da['mda'] : list
+        keys_da['mda'] : dict
             - tot_assim_steps: total number of iterations in MDA, e.g., 3
             - inflation_param: covariance inflation factors, e.g., [2, 4, 4]
 
@@ -222,17 +222,16 @@ def _ext_inflation_param(self):
         alpha: list
             Data covariance inflation factor
         """
-        # Make sure MDA is a list
-        if not isinstance(self.keys_da['mda'][0], list):
-            mda_opts = [self.keys_da['mda']]
-        else:
-            mda_opts = self.keys_da['mda']
+        try:
+            mda_opts = dict(self.keys_da['mda'])
+        except:
+            mda_opts = dict([self.keys_da['mda']])
 
         # Check if INFLATION_PARAM has been provided, and if so, extract the value(s). If not, we set alpha to the
         # default value equal to the tot. no. assim. steps
-        if 'inflation_param' in list(zip(*mda_opts))[0]:
+        if 'inflation_param' in mda_opts:
             # Extract value
-            alpha_tmp = [item[1] for item in mda_opts if item[0] == 'inflation_param'][0]
+            alpha_tmp = mda_opts['inflation_param']
 
             # If one value is given, we copy it to all assim. steps. If multiple values are given, we check the
             # number of parameters corresponds to tot. no. assim. steps
@@ -279,22 +278,17 @@ def _ext_assim_steps(self):
         - ST 7/6-16
         - ST 1/3-17: Changed to output list of assim. steps instead of just tot. assim. steps
         """
-        # Make sure MDA is a list
-        if not isinstance(self.keys_da['mda'][0], list):
-            mda_opts = [self.keys_da['mda']]
-        else:
-            mda_opts = self.keys_da['mda']
+        try:
+            mda_opts = dict(self.keys_da['mda'])
+        except:
+            mda_opts = dict([self.keys_da['mda']])
 
+    
         # Check if 'max_iter' has been given; if not, give error (mandatory in ITERATION)
-        assert 'tot_assim_steps' in list(
-            zip(*mda_opts))[0], 'TOT_ASSIM_STEPS has not been given in MDA!'
-
-        # Extract max. iter
-        tot_no_assim = int([item[1]
-                           for item in mda_opts if item[0] == 'tot_assim_steps'][0])
-
-        # Make a list of assim. steps
-        assim_steps = list(range(tot_no_assim))
+        try:
+            assim_steps = list(range(int(mda_opts['tot_assim_steps'])))
+        except KeyError:
+            raise AssertionError('TOT_ASSIM_STEPS has not been given in MDA!')
 
         # If it is a restart run, we remove simulations already done
         if self.restart is True:
diff --git a/popt/cost_functions/ecalc_npv.py b/popt/cost_functions/ecalc_npv.py
index d13d8df6..420e6ca4 100644
--- a/popt/cost_functions/ecalc_npv.py
+++ b/popt/cost_functions/ecalc_npv.py
@@ -44,9 +44,7 @@ def ecalc_npv(pred_data, **kwargs):
     report = kwargs.get('true_order', [])
 
     # Economic values
-    npv_const = {}
-    for name, value in keys_opt['npv_const']:
-        npv_const[name] = value
+    npv_const = dict(keys_opt['npv_const'])
 
     # Collect production data
     Qop = []
diff --git a/popt/cost_functions/ecalc_pareto_npv.py b/popt/cost_functions/ecalc_pareto_npv.py
index 2847dd99..9375f9f4 100644
--- a/popt/cost_functions/ecalc_pareto_npv.py
+++ b/popt/cost_functions/ecalc_pareto_npv.py
@@ -46,9 +46,7 @@ def ecalc_pareto_npv(pred_data, kwargs):
     report = kwargs.get('true_order', [])
 
     # Economic values
-    npv_const = {}
-    for name, value in keys_opt['npv_const']:
-        npv_const[name] = value
+    npv_const = dict(keys_opt['npv_const'])
 
     # Collect production data
     Qop = []
diff --git a/popt/cost_functions/npv.py b/popt/cost_functions/npv.py
index dfb3f6f8..bb18c8cb 100644
--- a/popt/cost_functions/npv.py
+++ b/popt/cost_functions/npv.py
@@ -38,9 +38,7 @@ def npv(pred_data, **kwargs):
     report = kwargs.get('true_order', [])
 
     # Economic values
-    npv_const = {}
-    for name, value in keys_opt['npv_const']:
-        npv_const[name] = value
+    npv_const = dict(keys_opt['npv_const'])
 
     values = []
     for i in np.arange(1, len(pred_data)):
diff --git a/popt/cost_functions/ren_npv.py b/popt/cost_functions/ren_npv.py
index e76c5671..0035f4a1 100644
--- a/popt/cost_functions/ren_npv.py
+++ b/popt/cost_functions/ren_npv.py
@@ -32,9 +32,7 @@ def ren_npv(pred_data, kwargs):
     report = kwargs.get('true_order', [])
 
     # Economic values
-    npv_const = {}
-    for name, value in keys_opt['npv_const']:
-        npv_const[name] = value
+    npv_const = dict(keys_opt['npv_const'])
 
     # Loop over timesteps
     values = []
diff --git a/popt/update_schemes/enopt.py b/popt/update_schemes/enopt.py
index a397ad92..ae9ed13e 100644
--- a/popt/update_schemes/enopt.py
+++ b/popt/update_schemes/enopt.py
@@ -8,7 +8,7 @@
 # Internal imports
 from popt.misc_tools import optim_tools as ot
 from popt.loop.optimize import Optimize
-import popt.update_schemes.optimizers as opt
+import popt.update_schemes.subroutines.optimizers as opt
 
 
 class EnOpt(Optimize):
diff --git a/popt/update_schemes/genopt.py b/popt/update_schemes/genopt.py
index 2316f6d4..408baa70 100644
--- a/popt/update_schemes/genopt.py
+++ b/popt/update_schemes/genopt.py
@@ -7,8 +7,8 @@
 # Internal imports
 from popt.misc_tools import optim_tools as ot
 from popt.loop.optimize import Optimize
-import popt.update_schemes.optimizers as opt
-from popt.update_schemes.cma import CMA
+import popt.update_schemes.subroutines.optimizers as opt
+from popt.update_schemes.subroutines.cma import CMA
 
 
 class GenOpt(Optimize):
diff --git a/popt/update_schemes/linesearch.py b/popt/update_schemes/linesearch.py
index 7bf1081f..e82b305e 100644
--- a/popt/update_schemes/linesearch.py
+++ b/popt/update_schemes/linesearch.py
@@ -5,15 +5,20 @@
 import warnings
 
 from numpy import linalg as la
-from functools import cache
 from scipy.optimize import OptimizeResult
-from scipy.optimize._dcsrch import DCSRCH
-from scipy.optimize._linesearch import _zoom
 
 # Internal imports
 from popt.misc_tools import optim_tools as ot
 from popt.loop.optimize import Optimize
-from popt.update_schemes import optimizers
+from popt.update_schemes.subroutines import line_search, line_search_backtracking, bfgs_update, newton_cg
+
+# Some symbols for logger
+subk = '\u2096'
+sup2 = '\u00b2'
+jac_inf_symbol = f'‖jac(x{subk})‖\u221E'
+fun_xk_symbol  = f'fun(x{subk})'
+nabla_symbol = "\u2207"
+
 
 def LineSearch(fun, x, jac, method='GD', hess=None, args=(), bounds=None, callback=None, **options):
     '''
@@ -32,8 +37,8 @@ def LineSearch(fun, x, jac, method='GD', hess=None, args=(), bounds=None, callba
     
     method: str
         Which optimization method to use. Default is 'GD' for 'Gradient Descent'.
-        Other options are 'BFGS' for the 'Broyden–Fletcher–Goldfarb–Shanno' method,
-        and 'Newton' for the Newton method.  
+        Other options are 'BFGS' for the 'Broyden-Fletcher-Goldfarb-Shanno' method,
+        and 'Newton-CG' for the Newton-conjugate gradient method.  
     
     hess: callable, optional
         Hessian function, hess(x, *args). Default is None. 
@@ -48,72 +53,83 @@ def LineSearch(fun, x, jac, method='GD', hess=None, args=(), bounds=None, callba
         A callable called after each successful iteration. The class instance of LineSearch 
         is passed as the only argument to the callback function: callback(self) 
     
-    **options: keyword arguments, optional
+    **options: 
+        keyword arguments, optional
 
     
     LineSearch Options (**options)
     ------------------------------
-    maxiter: int
+    - maxiter: int,
         Maximum number of iterations. Default is 20.
     
-    step_size: float
-        Step-size for optimizer. Default is 0.25/inf-norm(jac(x0)).
-
-    step_size_maxiter: int
+    - lsmaxiter: int,
         Maximum number of iterations for the line search. Default is 10.
+    
+    - step_size: float,
+        Step-size for optimizer. Default is 0.25/inf-norm(jac(x0)).
 
-    step_size_max: float
-        Maximum step-size. Default is 1e5
+    - step_size_max: float,
+        Maximum step-size. Default is 1e5. If bounds are specified, 
+        the maximum step-size is set to the maximum step-size allowed by the bounds.
 
-    step_size_adapt: int
+    - step_size_adapt: int,
         Set method for choosing initial step-size for each iteration. If 0, step_size value is used.
         If 1, Equation (3.6) from "Numercal Optimization" [1] is used. If 2, the equation above Equation (3.6) is used. 
         Default is 0.
                     
-    c1: float
+    - c1: float,
         Tolerance parameter for the Armijo condition. Default is 1e-4.
 
-    c2: float
+    - c2: float,
         Tolerance parameter for the Curvature condition. Default is 0.9.
 
-    xtol: float
+    - xtol: float,
         Optimization stop whenever |dx|<xtol. Default is 1e-8. 
+    
+    - ftol: float,
+        Optimization stop whenever |f_new - f_old| < ftol * |f_old|. Default is 1e-4.
+    
+    - gtol: float,
+        Optimization stop whenever ||jac||_inf < gtol. Default is 1e-5.
 
-    line_search_method: int
+    - lsmethod: int,
         Sets method for proposing new step-size in the line search.  Default is 1.
-        If line_search_method=0: step-size is cut in half.
-        If line_search_method=1: Algorithm (3.5) from [1] with polynomial interpolation is used.
-        If line_search_method=2: the DCSRCH implementation in scipy is used.
+        If lsmethod=0: backtracking is used (step-size is cut in half).
+        If lsmethod=1: line search Algorithm (3.5) from [1] with polynomial interpolation is used.
+    
+    - 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.
 
-    saveit: bool
+    - saveit: bool,
         If True, the results from each iteration is saved. Default is True.
 
-    save_folder: str
+    - save_folder: str,
         Name of folder to save the results to. Defaul is ./ (the current directory).
 
-    fun0: float
+    - fun0: float,
         Function value of the intial control.
 
-    jac0: ndarray
+    - jac0: ndarray,
         Jacobian of the initial control.
     
-    hess0: ndarray
+    - hess0: ndarray,
         Hessian value of the initial control.
 
-    hess0_inv: ndarray
+    - hess0_inv: ndarray, 
         Initial inverse of hessian (only if method = 'BFGS').
 
-    resample: int
+    - resample: int,
         Number of jacobian re-computations allowed if a line search fails. Default is 0.
         (useful if jacobian is stochastic)
 
-    savedata: list[str]
+    - savedata: list[str],
         Further specification of which class variables to save to the result files.
 
-    restart: bool
+    - restart: bool,
         Restart optimization from a restart file. Default is False
 
-    restartsave: bool
+    - restartsave: bool,
         Save a restart file after each successful iteration. Default is False
     
                 
@@ -140,7 +156,7 @@ def LineSearch(fun, x, jac, method='GD', hess=None, args=(), bounds=None, callba
     >>> from popt.update_schemes.linesearch import LineSearch
     >>> x0 = np.random.uniform(-3, 3, 2)
     >>> kwargs = {'maxiter': 100,
-                  'line_search_maxiter': 10,
+                  'lsmaxiter': 10,
                   'step_size_adapt': 1,
                   'saveit': False}
     >>> res = LineSearch(fun=rosen, x=x0, jac=rosen_der, method='BFGS', **kwargs)
@@ -182,20 +198,28 @@ def  __init__(self, fun, x, jac, method='GD', hess=None, args=(), bounds=None, c
             self.callback = callback
         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 step-size
         self.step_size       = options.get('step_size', None)
         self.step_size_max   = options.get('step_size_max', 1e5)
         self.step_size_adapt = options.get('step_size_adapt', 0)
 
-        # Set options for line-search method
-        self.line_search_kwargs  = {
+        # Set options for line-search 
+        self.lskwargs = {
             'c1': options.get('c1', 1e-4),
             'c2': options.get('c2', 0.9),
+            'rho': options.get('rho', 0.5),
             'amax': self.step_size_max,
-            'xtol': options.get('xtol', 1e-8),
-            'maxiter': options.get('line_search_maxiter', 10),
-            'method' : options.get('line_search_method', 1)
+            'maxiter': options.get('lsmaxiter', 10),
+            'method' : options.get('lsmethod', 1),
+            'logger' : self.logger.info
         }
 
         # Set other options
@@ -203,15 +227,18 @@ def  __init__(self, fun, x, jac, method='GD', hess=None, args=(), bounds=None, c
         self.resample  = options.get('resample', 0)
         self.saveit    = options.get('saveit', True)
 
+        # set tolerance for convergence
+        self.xtol = options.get('xtol', 1e-8) # tolerance for control vector
+        self.ftol = options.get('ftol', 1e-4) # relative tolerance for function value
+        self.gtol = options.get('gtol', 1e-5) # tolerance for inf-norm of jacobian
+
         # Check method
-        valid_methods = ['GD', 'BFGS', 'Newton']
+        valid_methods = ['GD', 'BFGS', 'Newton-CG']
         if not self.method in valid_methods:
             raise ValueError(f"'{self.method}' is not a valid method. Valid methods are: {valid_methods}")
         
-        # Make sure hessian is callable if mehtod='Newton'
-        if (self.method == 'Newton') and (not callable(self.hessian)):
-            warnings.warn('Newton’s method requires a hessian, method changed to BFGS')
-            self.method = 'BFGS'
+        if (self.method == 'Newton-CG') and (self.hessian is None):
+            print(f'Warning: No hessian function provided. Finite difference approximation is used: {nabla_symbol}{sup2}f(x{subk})d ≈ ({nabla_symbol}f(x{subk}+hd)-{nabla_symbol}f(x{subk}))/h')
 
         # Calculate objective function of startpoint
         if not self.restart:
@@ -238,15 +265,15 @@ def  __init__(self, fun, x, jac, method='GD', hess=None, args=(), bounds=None, c
             self.p_old = None
         
             # Initial results
-            self.optimize_result = self.update_results()
+            self.optimize_result = self.get_intermediate_results()
             if self.saveit:
                 ot.save_optimize_results(self.optimize_result)
-
             if self.logger is not None:
                 self.logger.info(f'       ====== Running optimization - Line search ({method}) ======')
-                self.logger.info('Specified options\n'+pprint.pformat(OptimizeResult(self.options)))
-                self.logger.info(f'       {"iter.":<10} {"fun":<15} {"step-size":<15} {"|grad|":<15}')
-                self.logger.info(f'       {self.iteration:<10} {self.fk:<15.4e} {0.0:<15.4e}  {la.norm(self.jk):<15.4e}')
+                self.logger.info('\nSPECIFIED OPTIONS:\n'+pprint.pformat(OptimizeResult(self.options)))
+                self.logger.info('')
+                self.logger.info(f'       {"iter.":<10} {fun_xk_symbol:<15} {jac_inf_symbol:<15} {"step-size":<15}')
+                self.logger.info(f'       {self.iteration:<10} {self.fk:<15.4e} {la.norm(self.jk, np.inf):<15.4e} {0:<15.4e}')
                 self.logger.info('')
 
         self.run_loop() 
@@ -268,6 +295,11 @@ def _jac(self, x):
             g = self.jacobian(x)
         else:
             g = self.jacobian(x, *self.args)
+
+        # project gradient onto the feasible set
+        if self.bounds is not None:
+            g = - self._project_pk(-g, x)
+
         return g
     
     def _hess(self, x):
@@ -279,80 +311,15 @@ def _hess(self, x):
             h = self.hessian(x)
         else:
             h = self.hessian(x, *self.args)
-        return make_matrix_psd(h)
+        return h
     
-    def update_results(self):
-
-        res = {'fun': self.fk, 
-               'x': self.xk, 
-               'jac': self.jk,
-               'hess': self.Hk,
-               'hess_inv': self.Hk_inv,
-               'nfev': self.nfev,
-               'njev': self.njev,
-               'nit': self.iteration,
-               'step-size': self.step_size,
-               'method': self.method,
-               'save_folder': self.options.get('save_folder', './')}
-        
-        for a, arg in enumerate(self.args):
-            res[f'args[{a}]'] = arg
-
-        if 'savedata' in self.options:
-            # Make sure "SAVEDATA" gives a list
-            if isinstance(self.options['savedata'], list):
-                savedata = self.options['savedata']
-            else:
-                savedata = [self.options['savedata']]
-
-            # Loop over variables to store in save list
-            for save_typ in savedata:
-                if save_typ in locals():
-                    res[save_typ] = eval('{}'.format(save_typ))
-                elif hasattr(self, save_typ):
-                    res[save_typ] = eval(' self.{}'.format(save_typ))
-                else:
-                    print(f'Cannot save {save_typ}!\n\n')
-
-        return OptimizeResult(res)
     
-    def _set_step_size(self, pk):
-        ''' Sets the step-size '''
-
-        # If first iteration
-        if self.iteration == 1:
-            if self.step_size is None:
-                self.step_size = 0.25/la.norm(pk, np.inf)
-                alpha = self.step_size
-            else:
-                alpha = self.step_size
-        else:
-            if np.dot(pk, self.jk) != 0 and self.step_size_adapt != 0:
-                if self.step_size_adapt == 1:
-                    alpha = 2*(self.fk - self.f_old)/np.dot(pk, self.jk)
-                if self.step_size_adapt == 2:
-                    slope_old = np.dot(self.p_old, self.j_old)
-                    slope_new = np.dot(pk, self.jk)
-                    alpha = self.step_size*slope_old/slope_new
-            else:
-                alpha = self.step_size
-
-        if alpha < 0: 
-            alpha = abs(alpha)
-
-        #if self.method in ['BFGS', 'Newton']:
-            # From "Numerical Optimization"
-        #    alpha = min(1, 1.01*alpha)
-        
-        return min(alpha, self.step_size_max)
-
-   
     def calc_update(self, iter_resamp=0):
 
         # Initialize variables for this step
         success = False
 
-        # If in resampling mode, compute jacobian 
+        # If in resampling mode, compute jacobian
         # Else, jacobian from in __init__ or from latest line_search is used
         if self.jk is None:
             self.jk = self._jac(self.xk)
@@ -372,44 +339,53 @@ def calc_update(self, iter_resamp=0):
             pk = - self.jk
         if self.method == 'BFGS':
             pk = - np.matmul(self.Hk_inv, self.jk)
-        if self.method == 'Newton':
-            pk = - np.matmul(la.inv(self.Hk), self.jk)
+        if self.method == 'Newton-CG':
+            pk = newton_cg(self.jk, Hk=self.Hk, xk=self.xk, jac=self._jac, logger=self.logger.info)
 
-        # 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
+        # porject search direction onto the feasible set
+        if self.bounds is not None:
+            pk = self._project_pk(pk, self.xk)
         
         # Set step_size
-        step_size = self._set_step_size(pk)
-
-        # Set maximum step-size
         if self.bounds is not None:
-            mean_bound_range = np.mean([b[1]-b[0] for b in self.bounds])
-            step_size_max = mean_bound_range/np.linalg.norm(pk)
-            self.line_search_kwargs['amax'] = step_size_max
+            self.step_size_max = self._set_max_step_size(pk, self.xk)
+            self.lskwargs['amax'] = self.step_size_max
+        step_size = self._set_step_size(pk, self.step_size_max)
 
         # Perform line-search 
-        self.logger.info('Performing line search...')
-        ls_res = line_search(
-            fun=self._fun,
-            jac=self._jac,
-            xk=self.xk,
-            pk=pk,
-            ak=step_size,
-            fk=self.fk,
-            gk=self.jk,
-            logger=self.logger,
-            **self.line_search_kwargs
-        )
-        step_size, f_new, f_old, j_new, self.msg = ls_res
+        self.logger.info('Performing line search.............')
+        if self.lskwargs['method'] == 0:
+            ls_res = line_search_backtracking(
+                step_size=step_size,
+                xk=self.xk,
+                pk=pk,
+                fun=self._fun,
+                jac=self._jac,
+                fk=self.fk,
+                jk=self.jk,
+                **self.lskwargs
+            )
+        else:
+            ls_res = line_search(
+                step_size=step_size,
+                xk=self.xk,
+                pk=pk,
+                fun=self._fun,
+                jac=self._jac,
+                fk=self.fk,
+                jk=self.jk,
+                **self.lskwargs
+            )
+        step_size, f_new, j_new, _, _ = ls_res
     
         if not (step_size is None):
-
+        
+            # Save old values
             x_old = self.xk
             j_old = self.jk
+            f_old = self.fk
+
+            # Update control
             x_new = ot.clip_state(x_old + step_size*pk, self.bounds)
 
             # Update state
@@ -421,6 +397,7 @@ def calc_update(self, iter_resamp=0):
             self.j_old = j_old
             self.f_old = f_old
             self.p_old = pk
+            sk = x_new - x_old
 
             # Call the callback function
             if callable(self.callback):
@@ -428,29 +405,51 @@ def calc_update(self, iter_resamp=0):
 
             # Update BFGS
             if self.method == 'BFGS':
-                sk  = x_new - x_old
-                yk  = j_new - j_old
-                rho = 1/np.dot(yk,sk)
-                id_mat = np.eye(sk.size)
-
-                matrix1 = (id_mat - rho*np.outer(sk, yk))
-                matrix2 = (id_mat - rho*np.outer(yk, sk))
-                self.Hk_inv = matrix1@self.Hk_inv@matrix2 + rho*np.outer(sk, sk)
+                yk = j_new - j_old
+                if self.iteration == 1: self.Hk_inv = np.dot(yk,sk)/np.dot(yk,yk) * np.eye(sk.size)
+                self.Hk_inv = bfgs_update(self.Hk_inv, sk, yk)
 
             # Update status
             success = True
 
             # Save Results
-            self.optimize_result = self.update_results()
+            self.optimize_result = self.get_intermediate_results()
             if self.saveit:
                 ot.save_optimize_results(self.optimize_result)
 
             # Write logging info
             if self.logger is not None:
                 self.logger.info('')
-                self.logger.info(f'       {"iter.":<10} {"fun":<15} {"step-size":<15} {"|grad|":<15}')
-                self.logger.info(f'       {self.iteration:<10} {self.fk:<15.4e} {step_size:<15.4e}  {la.norm(self.jk):<15.4e}')
+                self.logger.info(f'       {"iter.":<10} {fun_xk_symbol:<15} {jac_inf_symbol:<15} {"step-size":<15}')
+                self.logger.info(f'       {self.iteration:<10} {self.fk:<15.4e} {la.norm(self.jk, np.inf):<15.4e} {step_size:<15.4e}')
                 self.logger.info('')
+            
+            # Check for convergence
+            if (la.norm(sk, np.inf) < self.xtol):
+                self.msg = 'Convergence criteria met: |dx| < xtol'
+                self.logger.info(self.msg)
+                success = False
+                return success
+            if (np.abs(self.fk - f_old) < self.ftol * np.abs(f_old)):
+                self.msg = 'Convergence criteria met: |f(x+dx) - f(x)| < ftol * |f(x)|'
+                self.logger.info(self.msg)
+                success = False
+                return success
+            if (la.norm(self.jk, np.inf) < self.gtol):
+                self.msg = f'Convergence criteria met: {jac_inf_symbol} < gtol'
+                self.logger.info(self.msg)
+                success = False
+                return success
+
+            # Check for custom convergence
+            if callable(self.convergence_criteria):
+                if self.convergence_criteria(self):
+                    self.logger.info('Custom convergence criteria met. Stopping optimization.')
+                    success = False
+                    return success
+
+            if self.step_size_adapt == 2:
+                self.step_size = step_size
 
             # Update iteration
             self.iteration += 1
@@ -469,225 +468,99 @@ def calc_update(self, iter_resamp=0):
                 success = False
     
         return success
+    
+    def get_intermediate_results(self):
+
+        # Define default results
+        results = {
+            'fun': self.fk, 
+            'x': self.xk, 
+            'jac': self.jk,
+            'nfev': self.nfev,
+            'njev': self.njev,
+            'nit': self.iteration,
+            'method': self.method,
+            'save_folder': self.options.get('save_folder', './')
+        }
 
+        if 'savedata' in self.options:
+            # Make sure "SAVEDATA" gives a list
+            if isinstance(self.options['savedata'], list):
+                savedata = self.options['savedata']
+            else:
+                savedata = [self.options['savedata']]
 
+            if 'args' in savedata:
+                for a, arg in enumerate(self.args):
+                    results[f'args[{a}]'] = arg
 
-def line_search(fun, jac, xk, pk, ak, fk=None, gk=None, c1=0.0001, c2=0.9, maxiter=10, **kwargs):
-    '''
-    Performs a single line search step
-    '''
-    line_search_step = LineSearchStepBase(
-        fun, 
-        jac, 
-        xk, 
-        pk, 
-        ak, 
-        fk, 
-        gk, 
-        c1, 
-        c2, 
-        maxiter, 
-        **kwargs
-    )
-    return line_search_step()
-
-class LineSearchStepBase:
-
-    def __init__(self, fun, jac, xk, pk, ak, fk=None, gk=None, c1=0.0001, c2=0.9, maxiter=10, **kwargs):
-        self.fun = fun
-        self.jac = jac
-        self.xk = xk
-        self.pk = pk
-        self.ak = ak
-        self.fk = fk
-        self.gk = gk
-        self.c1 = c1
-        self.c2 = c2
-        self.maxiter = maxiter
-        self.msg = ''
-
-        # kwargs
-        self.amax   = kwargs.get('amax', 1e5)
-        self.amin   = kwargs.get('amin', 0.0)
-        self.xtol   = kwargs.get('xtol', 1e-8)
-        self.method = kwargs.get('method', 1)
-        self.logger = kwargs.get('logger', None)
-
-        # If c2 is None, the curvature condition is not used
-        if self.c2 is None:
-            self.c2 = np.inf
-            self.method = 0
-
-        # Check for initial values
-        if self.fk is None:
-            self.phi0 = self.phi(0, eval=False)
-        else:
-            self.phi0 = self.fk
+            # Loop over variables to store in save list
+            for variable in savedata:
+                if variable in locals():
+                    results[variable] = eval('{}'.format(variable))
+                elif hasattr(self, variable):
+                    results[variable] = eval('self.{}'.format(variable))
+                else:
+                    print(f'Cannot save {variable}!\n\n')
 
-        if self.gk is None:
-            self.dphi0 = self.dphi(0, eval=False)
-        else:
-            self.dphi0 = np.dot(self.pk, self.gk)
+        return OptimizeResult(results)
+    
+    def _set_step_size(self, pk, amax):
+        ''' Sets the step-size '''
 
+        # If first iteration
+        if (self.iteration == 1):
+            if (self.step_size is None):
+                self.step_size = 0.25/la.norm(pk, np.inf)
+                alpha = self.step_size
+            else:
+                alpha = self.step_size
 
-    def __call__(self):
-        
-        if self.method == 0:
-            step_size, fnew = self._line_search_alpha_cut(step_size=self.ak)
-        
-        if self.method == 1:
-            step_size, fnew = self._line_search_alpha_interpol(step_size=self.ak)
-
-        if self.method == 2:
-            dcsrch = DCSRCH(
-                self.phi, 
-                self.dphi, 
-                self.c1, 
-                self.c2, 
-                self.xtol, 
-                self.amin, 
-                self.amax
-            )
-            dcsrch_res = dcsrch(
-                self.ak, 
-                phi0=self.phi0, 
-                derphi0=self.dphi0, 
-                maxiter=self.maxiter
-            )
-            step_size, fnew, _,  self.msg = dcsrch_res
-            self.msg = str(self.msg)
-
-        if step_size is None:
-            if self.msg is None:
-                self.msg = 'Line search did not find a solution'
-            return None, None, None, None, self.msg
-        elif la.norm(step_size*self.pk) <= self.xtol:
-            self.msg = f'|dx| < {self.xtol}'
-            return None, None, None, None, self.msg
         else:
-            step_size = min(step_size, self.amax)
-            self.msg = 'Line search was successful'
-            return step_size, fnew, self.phi0, self.jac_val, self.msg
-
-
-    def _line_search_alpha_cut(self, step_size):
-
-        ak = step_size
-        for i in range(self.maxiter):
-            phi_new = self.phi(ak)
+            if (self.step_size_adapt == 1) and (np.dot(pk, self.jk) != 0):
+                alpha = 2*(self.fk - self.f_old)/np.dot(pk, self.jk)
+            elif (self.step_size_adapt == 2) and (np.dot(pk, self.jk) == 0):
+                slope_old = np.dot(self.p_old, self.j_old)
+                slope_new = np.dot(pk, self.jk)
+                alpha = self.step_size*slope_old/slope_new
+            else:
+                alpha = self.step_size
 
-            # Check Armijo Condition
-            if phi_new < self.phi0 + self.c1*ak*self.dphi0:
-                dphi_new = self.dphi(ak)
+        if alpha < 0: 
+            alpha = abs(alpha)
 
-                # Curvature condition
-                if abs(dphi_new) <= abs(self.c2*self.dphi0):
-                    return ak, phi_new
-            
-            ak = ak/2
+        if alpha >= amax:
+            alpha = 0.75*amax
         
-        return None, None
+        return alpha
     
-    def _line_search_alpha_interpol(self, step_size):
-        ak  = step_size
-
-        # Some lists    
-        alpha = [0.0]
-        phi   = [self.phi0]
-        dphi  = [self.dphi0]
+    def _project_pk(self, pk, xk):
+        ''' Projects the jacobian onto the feasible set defined by bounds '''
+        lb = np.array(self.bounds)[:, 0]
+        ub = np.array(self.bounds)[:, 1]
+        for i, pk_val in enumerate(pk):
+            if (xk[i] <= lb[i] and pk_val < 0) or (xk[i] >= ub[i] and pk_val > 0):
+                pk[i] = 0
+        return pk
+    
+    def _set_max_step_size(self, pk, xk):
+        lb = np.array(self.bounds)[:, 0]
+        ub = np.array(self.bounds)[:, 1]
 
-        for i in range(1, self.maxiter+1):
-            
-            # Append lists
-            alpha.append(ak)
-            phi.append(self.phi(ak))
-            dphi.append(self.dphi(ak))
-
-            # Check Armijo Condition
-            if phi[i] > self.phi0 + self.c1*alpha[i]*self.dphi0 or (phi[i] >= phi[i-1] and i>1):                
-                step_size_new, phi_new = self._zoom(alpha[i-1], alpha[i], phi[i-1], phi[i], dphi[i-1])
-                return step_size_new, phi_new
-            
-            if abs(dphi[i]) < - self.c2*self.dphi0:
-                return alpha[i], phi[i]
-            
-            # Check Curvature condition
-            if dphi[i] >= 0:
-                step_size_new, phi_new = self._zoom(alpha[i], alpha[i-1], phi[i], phi[i-1], dphi[i])
-                return step_size_new, phi_new
-            
-            if alpha[i] >= self.amax:
-                return None, None
+        amax = []
+        for i, pk_val in enumerate(pk):
+            if pk_val < 0:
+                amax.append((lb[i] - xk[i])/pk_val)
+            elif pk_val > 0:
+                amax.append((ub[i] - xk[i])/pk_val)
             else:
-                ak = ak*2
-        
-        return None, None
+                continue
+
+        return max(amax)
 
 
-    def log(self, msg):
-        if self.logger is None:
-            print(msg)
-        else:
-            self.logger.info(msg)
-    
-    @cache
-    def phi(self, a, eval=True):
-        if eval:
-            self.log('  Evaluating Armijo Condition')
-        return self.fun(self.xk + a*self.pk)
-
-    @cache
-    def dphi(self, a, eval=True):
-        if eval:
-            self.log('  Evaluating Curvature Condition')
-        jval = self.jac(self.xk + a*self.pk)
-        self.jac_val = jval
-        return np.dot(self.pk, jval)
-    
-    def _zoom(self, a_lo, a_hi, phi_lo, phi_hi, dphi_low):
-        alpha_new, phi_new, _  = _zoom(a_lo=a_lo, 
-                                       a_hi=a_hi, 
-                                       phi_lo=phi_lo, 
-                                       phi_hi=phi_hi,
-                                       derphi_lo=dphi_low,
-                                       phi=self.phi,
-                                       derphi=self.dphi,
-                                       phi0=self.phi0,
-                                       derphi0=self.dphi0,
-                                       c1=self.c1,
-                                       c2=self.c2,
-                                       extra_condition=lambda *args: True)
-        return alpha_new, phi_new
-
-
-def get_near_psd(A):
-    eigval, eigvec = np.linalg.eig((A + A.T)/2)
-    eigval[eigval < 0] = 1.0
-    return eigvec.dot(np.diag(eigval)).dot(eigvec.T)
-
-def make_matrix_psd(A, maxiter=100):
-    # Set beta to Frobenius norm of A
-    beta = np.linalg.norm(A, 'fro')
-    
-    # Initialize tau
-    if np.min(np.diag(A)) > 0:
-        tau = 0
-    else:
-        tau = beta/2
-    
-    for _ in range(maxiter):
-        try:
-            M = A + tau*np.eye(A.shape[0])
-            # Attempt Cholesky 
-            np.linalg.cholesky(A + tau*np.eye(A.shape[0]))
-            return M
-        except np.linalg.LinAlgError:
-            # Set new tau
-            tau = max(2*tau, beta/2)
-    
-    return None
-                
 
+                
 
             
 
diff --git a/popt/update_schemes/smcopt.py b/popt/update_schemes/smcopt.py
index f67503b6..f0944fe8 100644
--- a/popt/update_schemes/smcopt.py
+++ b/popt/update_schemes/smcopt.py
@@ -6,7 +6,7 @@
 
 # Internal imports
 from popt.loop.optimize import Optimize
-import popt.update_schemes.optimizers as opt
+import popt.update_schemes.subroutines.optimizers as opt
 from popt.misc_tools import optim_tools as ot
 
 
diff --git a/popt/update_schemes/subroutines/__init__.py b/popt/update_schemes/subroutines/__init__.py
new file mode 100644
index 00000000..eab7b843
--- /dev/null
+++ b/popt/update_schemes/subroutines/__init__.py
@@ -0,0 +1,3 @@
+from .subroutines import *
+from .cma import *
+from .optimizers import *
\ No newline at end of file
diff --git a/popt/update_schemes/cma.py b/popt/update_schemes/subroutines/cma.py
similarity index 99%
rename from popt/update_schemes/cma.py
rename to popt/update_schemes/subroutines/cma.py
index 0a998748..bee93e5e 100644
--- a/popt/update_schemes/cma.py
+++ b/popt/update_schemes/subroutines/cma.py
@@ -2,6 +2,8 @@
 import numpy as np
 from popt.misc_tools import optim_tools as ot
 
+__all__ = ['CMA']
+
 class CMA:
 
     def __init__(self, ne, dim, alpha_mu=None, n_mu=None, alpha_1=None, alpha_c=None, corr_update=False, equal_weights=True):
diff --git a/popt/update_schemes/optimizers.py b/popt/update_schemes/subroutines/optimizers.py
similarity index 99%
rename from popt/update_schemes/optimizers.py
rename to popt/update_schemes/subroutines/optimizers.py
index 83514cdc..a2110938 100644
--- a/popt/update_schemes/optimizers.py
+++ b/popt/update_schemes/subroutines/optimizers.py
@@ -1,6 +1,8 @@
 """Gradient acceleration."""
 import numpy as np
 
+__all__ = ['GradientAscent', 'Adam', 'AdaMax', 'Steihaug', ]
+
 
 class GradientAscent:
     r"""
diff --git a/popt/update_schemes/subroutines/subroutines.py b/popt/update_schemes/subroutines/subroutines.py
new file mode 100644
index 00000000..ddd1105a
--- /dev/null
+++ b/popt/update_schemes/subroutines/subroutines.py
@@ -0,0 +1,397 @@
+import numpy as np
+import numpy.linalg as la
+from functools import lru_cache
+from scipy.optimize._linesearch import _quadmin, _cubicmin
+
+__all__ = [
+    'line_search', 
+    'zoom', 
+    'line_search_backtracking', 
+    'bfgs_update', 
+    'newton_cg'
+] 
+
+
+def line_search(step_size, xk, pk, fun, jac, fk=None, jk=None, **kwargs):
+    '''
+    Line search algorithm to find step size alpha that satisfies the Wolfe conditions.
+
+    Parameters
+    ----------
+    step_size : float
+        Initial step size to start the line search.
+
+    xk : ndarray
+        Current point in the optimization process.
+
+    pk : ndarray
+        Search direction.
+    
+    fun : callable
+        Objective function
+
+    jac : callable
+        Gradient of the objective function
+    
+    fk : float, optional
+        Function value at xk. If None, it will be computed.
+    
+    jk : ndarray, optional
+        Gradient at xk. If None, it will be computed.
+    
+    **kwargs : dict
+        Additional parameters for the line search, such as:
+        - amax : float, maximum step size (default: 1000)
+        - maxiter : int, maximum number of iterations (default: 10)
+        - c1 : float, sufficient decrease condition (default: 1e-4)
+        - c2 : float, curvature condition (default: 0.9)
+    
+    Returns
+    -------
+    alpha : float
+        Step size that satisfies the Wolfe conditions.
+    
+    fval : float
+        Function value at the new point xk + step_size*pk.
+    
+    jval : ndarray
+        Gradient at the new point xk + step_size*pk.
+    
+    nfev : int
+        Number of function evaluations.
+    
+    njev : int
+        Number of gradient evaluations.
+    '''
+
+    global ls_nfev
+    global ls_njev
+    ls_nfev = 0
+    ls_njev = 0
+
+    # Unpack some kwargs
+    amax = kwargs.get('amax', 1000)
+    maxiter = kwargs.get('maxiter', 10)
+    c1 = kwargs.get('c1', 1e-4)
+    c2 = kwargs.get('c2', 0.9)
+
+    # check for logger in kwargs
+    global logger
+    logger = kwargs.get('logger', None)
+    if logger is None:
+        logger = print
+
+    # assertions
+    assert step_size <= amax, "Initial step size must be less than or equal to amax."
+
+    # Define phi and derivative of phi
+    @lru_cache(maxsize=None)
+    def phi(alpha):
+        global ls_nfev
+        if (alpha == 0):
+            if (fk is None):
+                phi.fun_val = fun(xk)
+                ls_nfev += 1
+            else:
+                phi.fun_val = fk
+        else:
+            logger('    Evaluating Armijo condition')
+            phi.fun_val = fun(xk + alpha*pk)
+            ls_nfev += 1
+        return phi.fun_val
+    
+    @lru_cache(maxsize=None)
+    def dphi(alpha):
+        global ls_njev
+        if (alpha == 0):
+            if (jk is None):
+                dphi.jac_val = jac(xk)
+                ls_njev += 1
+            else:
+                dphi.jac_val = jk
+        else:
+            logger('    Evaluating curvature condition')
+            dphi.jac_val = jac(xk + alpha*pk)
+            ls_njev += 1
+        return np.dot(dphi.jac_val, pk)
+    
+    # Define initial values of phi and dphi
+    phi_0  = phi(0)
+    dphi_0 = dphi(0)
+
+    # Start loop
+    a = [0, step_size]
+    for i in range(1, maxiter+1):
+        logger(f'Line search iteration: {i-1}')
+
+        # Evaluate phi(ai)
+        phi_i = phi(a[i])
+
+        # Check for sufficient decrease
+        if (phi_i > phi_0 + c1*a[i]*dphi_0) or (phi_i >= phi(a[i-1]) and i>0):
+            # Call zoom function
+            step_size = zoom(a[i-1], a[i], phi, dphi, phi_0, dphi_0, maxiter+1-i, c1, c2) 
+            return step_size, phi.fun_val, dphi.jac_val, ls_nfev, ls_njev
+
+        # Evaluate dphi(ai)
+        dphi_i = dphi(a[i])
+
+        # Check curvature condition
+        if abs(dphi_i) <= -c2*dphi_0:
+            step_size = a[i]
+            return step_size, phi.fun_val, dphi.jac_val, ls_nfev, ls_njev
+        
+        # Check for posetive derivative
+        if dphi_i >= 0:
+            # Call zoom function
+            step_size = zoom(a[i], a[i-1], phi, dphi, phi_0, dphi_0, maxiter+1-i, c1, c2) 
+            return step_size, phi.fun_val, dphi.jac_val, ls_nfev, ls_njev
+        
+        # Increase ai
+        a.append(min(2*a[i], amax))
+        
+    # If we reached this point, the line search failed
+    logger('Line search failed to find a suitable step size \n')
+    return None, None, None, ls_nfev, ls_njev
+            
+
+def zoom(alo, ahi, f, df, f0, df0, maxiter, c1, c2):
+    '''Zoom function for line search algorithm. (This is the same as for scipy)'''
+
+    phi_lo = f(alo)
+    phi_hi = f(ahi)
+    dphi_lo = df(alo)
+
+    for j in range(maxiter):
+        logger(f'Line search iteration: {j+1}')
+
+        tol_cubic = 0.2*(ahi-alo)
+        tol_quad  = 0.1*(ahi-alo)
+
+        if (j > 0):
+            # cubic interpolation for alo, phi(alo), dphi(alo) and ahi, phi(ahi)
+            aj = _cubicmin(alo, phi_lo, dphi_lo, ahi, phi_hi, aold, phi_old)
+        if (j == 0) or (aj is None) or (aj < alo + tol_cubic) or (aj > ahi - tol_cubic):
+            # quadratic interpolation for alo, phi(alo), dphi(alo) and ahi, phi(ahi)
+            aj = _quadmin(alo, phi_lo, dphi_lo, ahi, phi_hi)
+
+        # Ensure aj is within bounds
+        if (aj is None) or (aj <  alo + tol_quad) or (aj > ahi - tol_quad):
+                aj = alo + 0.5*(ahi - alo)
+
+        # Evaluate phi(aj)
+        phi_j = f(aj)
+
+        # Check for sufficient decrease
+        if (phi_j > f0 + c1*aj*df0) or (phi_j >= phi_lo):
+            # store old values
+            aold = ahi
+            phi_old = phi_hi
+            # update ahi
+            ahi = aj
+            phi_hi = phi_j
+        else:
+            # check curvature condition
+            dphi_j = df(aj)
+            if abs(dphi_j) <= -c2*df0:
+                return aj
+            
+            if dphi_j*(ahi-alo) >= 0:
+                # store old values
+                aold = ahi
+                phi_old = phi_hi
+                # update alo
+                ahi = alo
+                phi_hi = phi_lo
+            else:
+                # store old values
+                aold = alo
+                phi_old = phi_lo
+
+            alo = aj
+            phi_lo = phi_j
+            dphi_lo = dphi_j
+
+    # If we reached this point, the line search failed
+    return None
+
+
+
+
+def line_search_backtracking(step_size, xk, pk, fun, jac, fk=None, jk=None, **kwargs):
+    '''
+    Backtracking line search algorithm to find step size alpha that satisfies the Wolfe conditions.
+
+    Parameters
+    ----------
+    step_size : float
+        Initial step size to start the line search.
+
+    xk : ndarray
+        Current point in the optimization process.
+
+    pk : ndarray
+        Search direction.
+    
+    fun : callable
+        Objective function
+
+    jac : callable
+        Gradient of the objective function
+    
+    fk : float, optional
+        Function value at xk. If None, it will be computed.
+    
+    jk : ndarray, optional
+        Gradient at xk. If None, it will be computed.
+    
+    **kwargs : dict
+        Additional parameters for the line search, such as:
+        - rho : float, backtracking factor (default: 0.5)
+        - maxiter : int, maximum number of iterations (default: 10)
+        - c1 : float, sufficient decrease condition (default: 1e-4)
+        - c2 : float, curvature condition (default: 0.9)
+    
+    Returns
+    -------
+    alpha : float
+        Step size that satisfies the Wolfe conditions.
+    
+    fval : float
+        Function value at the new point xk + step_size*pk.
+    
+    jval : ndarray
+        Gradient at the new point xk + step_size*pk.
+    
+    nfev : int
+        Number of function evaluations.
+    
+    njev : int
+        Number of gradient evaluations.
+    '''
+
+    global ls_nfev
+    global ls_njev
+    ls_nfev = 0
+    ls_njev = 0
+
+    # Unpack some kwargs
+    rho = kwargs.get('rho', 0.5)
+    maxiter = kwargs.get('maxiter', 10)
+    c1 = kwargs.get('c1', 1e-4)
+
+    # Define phi and derivative of phi
+    @lru_cache(maxsize=None)
+    def phi(alpha):
+        global ls_nfev
+        if (alpha == 0):
+            if (fk is None):
+                fun_val = fun(xk)
+                ls_nfev += 1
+            else:
+                fun_val = fk
+        else:
+            fun_val = fun(xk + alpha*pk)
+            ls_nfev += 1
+        return fun_val
+    
+
+    # run the backtracking line search loop
+    for i in range(maxiter):
+        # Evaluate phi(alpha)
+        phi_i = phi(step_size)
+
+        # Check for sufficient decrease
+        if (phi_i <= phi(0) + c1*step_size*np.dot(jk, pk)):
+                # Evaluate jac at new point
+                jac_new = jac(xk + step_size*pk)
+                return step_size, phi_i, jac_new, ls_nfev, ls_njev
+        
+        # Reduce step size
+        step_size *= rho  
+
+    # If we reached this point, the line search failed
+    return None, None, None, ls_nfev, ls_njev  
+
+
+def bfgs_update(Hk, sk, yk):
+    """
+    Perform the BFGS update of the inverse Hessian approximation.
+
+    Parameters:
+    - Hk: np.ndarray, current inverse Hessian approximation (n x n)
+    - sk: np.ndarray, step vector (x_{k+1} - x_k), shape (n,)
+    - yk: np.ndarray, gradient difference (grad_{k+1} - grad_k), shape (n,)
+
+    Returns:
+    - Hk_new: np.ndarray, updated inverse Hessian approximation
+    """
+    sk = sk.reshape(-1, 1)
+    yk = yk.reshape(-1, 1)
+    rho = 1.0 / (yk.T @ sk)
+
+    if rho <= 0:
+        print('Non-positive curvature detected. BFGS update skipped....')
+        return Hk
+
+    I = np.eye(Hk.shape[0])
+    Vk = I - rho * sk @ yk.T
+    Hk_new = Vk @ Hk @ Vk.T + rho * sk @ sk.T
+
+    return Hk_new
+
+def newton_cg(gk, Hk=None, maxiter=None, **kwargs):
+
+    # Check for logger
+    logger = kwargs.get('logger', None)
+    if logger is None:
+        logger = print
+    
+    logger('Running Newton-CG subroutine..........')
+
+    if Hk is None:
+        jac = kwargs.get('jac')
+        eps = kwargs.get('eps', 1e-4)
+        xk  = kwargs.get('xk')
+
+        # define a finite difference approximation of the Hessian times a vector
+        def Hessd(d):
+            return (jac(xk + eps*d) - gk)/eps
+
+    if maxiter is None:
+        maxiter = 20*gk.size # Same dfault as in scipy
+
+    tol = min(0.5, np.sqrt(la.norm(gk)))*la.norm(gk)
+    z = 0
+    r = gk
+    d = -r
+
+    for j in range(maxiter):
+        logger(f'iteration: {j}')
+        if Hk is None:
+            Hd = Hessd(d)
+        else:
+            Hd = np.matmul(Hk, d)
+
+        dTHd = np.dot(d, Hd)
+
+        if dTHd <= 0:
+            logger('Negative curvature detected, terminating subroutine')
+            logger('')
+            if j == 0:
+                return -gk
+            else:
+                return z
+            
+        rold = r
+        a = np.dot(r,r)/dTHd
+        z = z + a*d
+        r = r + a*Hd
+
+        if la.norm(r) < tol:
+            logger('Subroutine converged')
+            logger('')
+            return z
+
+        b = np.dot(r, r)/np.dot(rold, rold)
+        d = -r + b*d
\ No newline at end of file
diff --git a/popt/update_schemes/trust_region.py b/popt/update_schemes/trust_region.py
index cd616f93..b23f6b6e 100644
--- a/popt/update_schemes/trust_region.py
+++ b/popt/update_schemes/trust_region.py
@@ -239,44 +239,6 @@ def _hess(self, x):
         else:
             h = self.hessian(x, *self.args)
         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', './')
-        }
-        
-        for a, arg in enumerate(self.args):
-            res[f'args[{a}]'] = arg
-
-        if 'savedata' in self.options:
-            # Make sure "SAVEDATA" gives a list
-            if isinstance(self.options['savedata'], list):
-                savedata = self.options['savedata']
-            else:
-                savedata = [self.options['savedata']]
-
-            # Loop over variables to store in save list
-            for save_typ in savedata:
-                if save_typ in locals():
-                    res[save_typ] = eval('{}'.format(save_typ))
-                elif hasattr(self, save_typ):
-                    res[save_typ] = eval(' self.{}'.format(save_typ))
-                else:
-                    print(f'Cannot save {save_typ}!\n\n')
-
-        return OptimizeResult(res)
-
-    def _log(self, msg):
-        if self.logger is not None:
-            self.logger.info(msg)
 
     def solve_subproblem(self, g, B, delta):
         """
@@ -426,6 +388,44 @@ def calc_update(self, inner_iter=0):
                 success = False
 
         return success
+    
+    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', './')
+        }
+        
+        for a, arg in enumerate(self.args):
+            res[f'args[{a}]'] = arg
+
+        if 'savedata' in self.options:
+            # Make sure "SAVEDATA" gives a list
+            if isinstance(self.options['savedata'], list):
+                savedata = self.options['savedata']
+            else:
+                savedata = [self.options['savedata']]
+
+            # Loop over variables to store in save list
+            for save_typ in savedata:
+                if save_typ in locals():
+                    res[save_typ] = eval('{}'.format(save_typ))
+                elif hasattr(self, save_typ):
+                    res[save_typ] = eval(' self.{}'.format(save_typ))
+                else:
+                    print(f'Cannot save {save_typ}!\n\n')
+
+        return OptimizeResult(res)
+
+    def _log(self, msg):
+        if self.logger is not None:
+            self.logger.info(msg)