Updates

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

popt/loop/base.py → 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])
+                    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])

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

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):
         """

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)