Objective normalizations & minor tweaks

Author: Maslyaev

epde/interface/prepared_tokens.py

@@ -222,21 +222,23 @@ class TrigonometricTokens(PreparedTokens):
     """
     Class for prepared tokens, that belongs to the trigonometric family
     """
-    def __init__(self, freq: tuple = (np.pi/2., 2*np.pi), dimensionality=1):
+    def __init__(self, freq: tuple = (np.pi/2., 2*np.pi), dimensionality: int = 1, meaningful: bool = False):
         """
         Initialization of class
 
         Args:
             freq (`tuple`): optional, default - (pi/2., 2*pi)
-                interval for parameter frequency in trigonometric token
+                interval for parameter frequency in trigonometric token.
             dimensionality (`int`): optional, default - 1
-                data dimension 
+                data dimension.
+            meaningful (`bool`): optional, defaule - False
+                flag, if the token can be used as a separate term or term-forming factor in the equation.
         """
         assert freq[1] > freq[0] and len(freq) == 2, 'The tuple, defining frequncy interval, shall contain 2 elements with first - the left boundary of interval and the second - the right one. '
 
         self._token_family = TokenFamily(token_type='trigonometric')
         self._token_family.set_status(unique_specific_token=True, unique_token_type=True,
-                                      meaningful=False)
+                                      meaningful=meaningful)
 
         def latex_form(label, **params):
             '''

epde/interface/token_family.py

@@ -275,7 +275,7 @@ def chech_constancy(self, **tfkwargs):
             try:
                 constancy = np.isclose(np.min(data), np.max(data))
             except TypeError:
-                print(f"No {label} data in cache!")
+                print(f"No {label} data in cache for constancy check. Functionality of eval-d token check TBD.")
                 continue
             if constancy:
                 constant_tokens_labels.append(label)

epde/operators/common/right_part_selection.py

@@ -115,7 +115,7 @@ def simplify_equation(self, objective: Equation):
                         term.structure = [factor for factor in term.structure if factor not in factors_simplified]
                         term.reset_saved_state()
                         # If term's order became zero -- replace term
-                        if len(term.structure) == 0:
+                        if (len(term.structure) == 0 or not term.contains_meaningful()):
                             term.randomize()
                             term.reset_saved_state()
                         while objective.structure.count(term) > 1 or term == temp:

epde/operators/multiobjective/moeadd_specific.py

@@ -8,15 +8,18 @@
 import copy
 import numpy as np
 import time
-from typing import Union
+from typing import Union, Tuple
 from functools import reduce, partial
 
-from epde.optimizers.moeadd.moeadd import ParetoLevels
+from epde.optimizers.moeadd.moeadd import ParetoLevels, ObjFunNormalizer
 from epde.operators.utils.template import CompoundOperator, add_base_param_to_operator
 from epde.operators.multiobjective.mutations import get_basic_mutation
 
+from epde.structure.main_structures import SoEq
 
-def penalty_based_intersection(sol_obj, weight, ideal_obj, penalty_factor = 1.) -> float:
+
+def penalty_based_intersection(sol_obj, weight, ideal_obj, 
+                               penalty_factor = 1., obj_normalizer: ObjFunNormalizer = None) -> float:
     '''
     Calculation of the penalty pased intersection, that is minimized for the solutions inside the 
     domain, specified by **weight** vector. The calculations are held, according to the following formulas:
@@ -50,10 +53,17 @@ def penalty_based_intersection(sol_obj, weight, ideal_obj, penalty_factor = 1.)
     
     penalty_factor : float, optional, default 1.
         The penalty parameter, represents :math:`\Theta` in the equations.
+
+    obj_normalizer : ObjFunNormalizer obj., optional, defaut None.
+        Normalizer for solution objective functions.
     
     '''
-    d_1 = np.dot((sol_obj.obj_fun - ideal_obj), weight) / np.linalg.norm(weight)
-    d_2 = np.linalg.norm(sol_obj.obj_fun - (ideal_obj + d_1 * weight/np.linalg.norm(weight)))
+    print(f'Objective before normalization: {sol_obj.obj_fun} for normalizer {obj_normalizer}')
+    solution_objective = sol_obj.obj_fun if obj_normalizer is None else obj_normalizer(sol_obj.obj_fun)
+    print(f'Objective after expected normalization: {solution_objective}')
+    
+    d_1 = np.dot((solution_objective - ideal_obj), weight) / np.linalg.norm(weight)
+    d_2 = np.linalg.norm(solution_objective - (ideal_obj + d_1 * weight/np.linalg.norm(weight)))
     return d_1 + penalty_factor * d_2
 
 
@@ -87,7 +97,7 @@ def population_to_sectors(population, weights):
     return list(map(solution_selection, np.arange(len(weights))))    
 
 
-def locate_pareto_worst(levels, weights, best_obj, penalty_factor = 1.):
+def locate_pareto_worst(levels: ParetoLevels, weights: np.ndarray, best_obj: np.ndarray, penalty_factor: float = 1.):
     '''
     
     Function, dedicated to the selection of the worst solution on the Pareto levels.
@@ -114,7 +124,8 @@ def locate_pareto_worst(levels, weights, best_obj, penalty_factor = 1.):
     if len(crowded_domains) == 1:
         most_crowded_domain = crowded_domains[0]
     else:
-        PBI = lambda domain_idx: sum([penalty_based_intersection(sol_obj, weights[domain_idx], best_obj, penalty_factor) for sol_obj in domain_solutions[domain_idx]])
+        PBI = lambda domain_idx: sum([penalty_based_intersection(sol_obj, weights[domain_idx], best_obj, penalty_factor, levels.normalizer)
+                                      for sol_obj in domain_solutions[domain_idx]])
         PBIS = np.fromiter(map(PBI, crowded_domains), dtype = float)
         most_crowded_domain = crowded_domains[np.argmax(PBIS)]
 
@@ -127,21 +138,34 @@ def locate_pareto_worst(levels, weights, best_obj, penalty_factor = 1.):
     max_level = np.max(domain_solution_NDL_idxs)
     worst_NDL_section = [domain_solutions[most_crowded_domain][sol_idx] for sol_idx in np.arange(len(domain_solutions[most_crowded_domain])) 
                         if domain_solution_NDL_idxs[sol_idx] == max_level]
-    PBIS = np.fromiter(map(lambda solution: penalty_based_intersection(solution, weights[most_crowded_domain], best_obj, penalty_factor), worst_NDL_section), dtype = float)
+    PBIS = np.fromiter(map(lambda solution: penalty_based_intersection(solution, weights[most_crowded_domain], best_obj, penalty_factor, levels.normalizer),
+                           worst_NDL_section), dtype = float)
     return worst_NDL_section[np.argmax(PBIS)]
 
 
 class PopulationUpdater(CompoundOperator):
     key = 'PopulationUpdater'
 
-    def apply(self, objective : ParetoLevels, arguments : dict):
+    def apply(self, objective : Tuple[Union[SoEq, ParetoLevels]], arguments : dict):
         '''
         Update population to get the pareto-nondomiated levels with the worst element removed. 
         Here, "worst" means the solution with highest PBI value (penalty-based boundary intersection)
-        '''         
+        '''
+        assert isinstance(objective, tuple), f'Expected input of PopulationUpdater to be a Tuple of SoEq and ParetoLevels.\n'\
+                                             f'Did not get even a Tuple, instead got {type(objective)}!'
+        assert isinstance(objective[0], SoEq), f'Expected input of PopulationUpdater to be a Tuple of SoEq and ParetoLevels.\n'\
+                                               f'Did not get a SoEq obj in the first position, instead got {type(objective[0])}!'        
+        assert isinstance(objective[1], ParetoLevels), f'Expected input of PopulationUpdater to be a Tuple of SoEq and ParetoLevels.\n'\
+                                                       f'Did not get even a ParetoLevels in the second position, '\
+                                                       f'instead got {type(objective[1])}!.'
+
         self_args, subop_args = self.parse_suboperator_args(arguments = arguments)
         # print(f'PopulationUpdater.params is {self.params}')
 
+        # TODO: Init normalizer here!
+        # print('objective is ', objective)
+        # objective[1].set_normalizer()
+
         objective[1].update(objective[0])  #levels_updated = ndl_update(offspring, levels)
         if len(objective[1].levels) == 1:
             worst_solution = locate_pareto_worst(objective[1], self_args['weights'], 
@@ -168,7 +192,8 @@ def apply(self, objective : ParetoLevels, arguments : dict):
                 else:
                     PBI = lambda domain_idx: np.sum([penalty_based_intersection(sol_obj, self_args['weights'][domain_idx],
                                                                                 self_args['best_obj'], 
-                                                                                self.params['PBI_penalty'])
+                                                                                self.params['PBI_penalty'],
+                                                                                objective[1].normalizer)
                                                      for sol_obj in last_level_by_domains[domain_idx]])
                     PBIS = np.fromiter(map(PBI, crowded_domains), dtype = float)
                     most_crowded_domain = crowded_domains[np.argmax(PBIS)]
@@ -179,7 +204,8 @@ def apply(self, objective : ParetoLevels, arguments : dict):
                 else:
                     PBIS = np.fromiter(map(lambda solution: penalty_based_intersection(solution, 
                                                                                        self_args['weights'][most_crowded_domain],
-                                                                                       self_args['best_obj'], self.params['PBI_penalty']),
+                                                                                       self_args['best_obj'], self.params['PBI_penalty'],
+                                                                                       objective[1].normalizer),
                                                last_level_by_domains[most_crowded_domain]), dtype = float)
                     worst_solution = last_level_by_domains[most_crowded_domain][np.argmax(PBIS)]                    
 
@@ -242,7 +268,8 @@ def apply(self, objective : ParetoLevels, arguments : dict):
                         most_crowded_domain = crowded_domains[0]
                     else:
                         PBI = lambda domain_idx: np.sum([penalty_based_intersection(sol_obj, self_args['weights'][domain_idx], 
-                                                                                    self_args['best_obj'], self.params['PBI_penalty']) 
+                                                                                    self_args['best_obj'], self.params['PBI_penalty'],
+                                                                                    objective.normalizer) 
                                                             for sol_obj in last_level_by_domains[domain_idx]])
                         PBIS = np.fromiter(map(PBI, crowded_domains), dtype = float)
                         most_crowded_domain = crowded_domains[np.argmax(PBIS)]
@@ -417,7 +444,7 @@ def apply(self, objective : ParetoLevels, arguments : dict):
 
         '''
         self_args, subop_args = self.parse_suboperator_args(arguments = arguments)
-        
+
         if len(objective.population) == 0:
             for idx, candidate in enumerate(objective.unplaced_candidates):
                 self.suboperators['right_part_selector'].apply(objective = candidate,
@@ -432,6 +459,10 @@ def apply(self, objective : ParetoLevels, arguments : dict):
                                                                   arguments=subop_args['chromosome_fitness'])
                 objective.history.add(tuple(candidate.obj_fun))
             objective.initial_placing()
+        
+            # TODO: consider carefully, where normalizer init shall be held. If here, only the initial values are employed
+        objective.set_normalizer()
+
         return objective
 
 def get_initial_sorter(right_part_selector : CompoundOperator, 

epde/operators/multiobjective/mutations.py

@@ -115,7 +115,8 @@ def apply(self, objective : tuple, arguments : dict): #term_idx, equation):
         new_term = Term(objective[1].pool, mandatory_family = objective[1].structure[objective[0]].descr_variable_marker, 
                         create_derivs=create_derivs,
                         max_factors_in_term = objective[1].metaparameters['max_factors_in_term']['value'])
-        while not check_uniqueness(new_term, objective[1].structure[:objective[0]] + objective[1].structure[objective[0]+1:]):
+        while not (check_uniqueness(new_term, objective[1].structure[:objective[0]] + objective[1].structure[objective[0]+1:]) and
+                   new_term.contains_meaningful()):
             new_term = Term(objective[1].pool, mandatory_family = objective[1].structure[objective[0]].descr_variable_marker, 
                             create_derivs=create_derivs,
                             max_factors_in_term = objective[1].metaparameters['max_factors_in_term']['value'])

epde/optimizers/moeadd/moeadd.py

@@ -9,6 +9,7 @@
 from itertools import chain
 
 from typing import Union, List
+from functools import reduce
 
 def flatten_chain(matrix):
     return list(chain.from_iterable(matrix))
@@ -32,6 +33,15 @@ def clear_list_of_lists(inp_list) -> list:
     return [elem for elem in inp_list if len(elem) > 0]
 
 
+class ObjFunNormalizer(object):
+    def __init__(self, obj_worst_vals: np.ndarray):
+        self._worst_vals = obj_worst_vals
+
+    def __call__(self, obj_vals: np.ndarray):
+        assert obj_vals.size == self._worst_vals.size, 'Passed objective values have different length, than stored max ones.'
+        return obj_vals / self._worst_vals
+
+
 class ParetoLevels(object):
     '''
     
@@ -53,7 +63,7 @@ class ParetoLevels(object):
     
     '''
     def __init__(self, population, sorting_method = fast_non_dominated_sorting, 
-                 update_method = ndl_update, initial_sort = False):
+                 update_method = ndl_update): # , initial_sort = False
         """
         Args:
             population (`list`): List with the elements - canidate solutions of the case-specific subclass of 
@@ -67,6 +77,8 @@ def __init__(self, population, sorting_method = fast_non_dominated_sorting,
         self.population = []
         self._update_method = update_method
         self.unplaced_candidates = population
+
+        self.normalizer = None
         self.history = set()
 
     def manual_reconst(self, attribute:str, value, except_attrs:dict):
@@ -87,6 +99,15 @@ def attrs_from_dict(self, attributes, except_keys = ['obj_type']):
         self.__dict__ = {key : item for key, item in attributes.items()
                          if key not in except_keys}
 
+    def set_normalizer(self):
+        # worst_objectives = reduce(lambda x, y: x.extend(y),
+        #                           [[elem.obj_fun for elem in level] for level in self.levels], []) # : np.ndarray
+        objectives = np.stack(reduce(lambda x, y: x.extend(y) or x,
+                                     [[elem.obj_fun for elem in level] for level in self.levels]), axis = 0)
+        
+
+        self.normalizer = ObjFunNormalizer(np.max(objectives, axis = 0))
+
     @property
     def levels(self):
         return self._levels
@@ -338,8 +359,7 @@ def __init__(self, population_instruct, weights_num, pop_size, solution_params,
                               Confirmed {len(population)}/{pop_size} solutions.')
                         break
                     solution_gen_idx += 1
-            self.pareto_levels = ParetoLevels(population, sorting_method = nds_method, update_method = ndl_update,
-                                              initial_sort = False)
+            self.pareto_levels = ParetoLevels(population, sorting_method = nds_method, update_method = ndl_update) # initial_sort = False
         else:
             if not isinstance(passed_population, ParetoLevels):
                 raise TypeError(f'Incorrect type of the population passed. Expected ParetoLevels object, instead got \
@@ -449,7 +469,6 @@ def set_strategy(self, strategy_director):
         builder.assemble(True)
         self.set_sector_processer(builder.processer)
 
-        
     def optimize(self, epochs):
         """
         Method for the main unconstrained evolutionary optimization. Can be applied repeatedly to 

epde/structure/main_structures.py

@@ -327,6 +327,9 @@ def contains_deriv(self, variable=None):
     def contains_variable(self, variable):
         return any([factor.variable == variable for factor in self.structure])
 
+    def contains_meaningful(self):
+        return any([factor.status['meaningful'] for factor in self.structure])
+
     def __eq__(self, other):
         return (all([any([other_elem == self_elem for other_elem in other.structure]) for self_elem in self.structure])
                 and all([any([other_elem == self_elem for self_elem in self.structure]) for other_elem in other.structure])