System discovery

epde/cache/cache.py

@@ -20,6 +20,9 @@
 
 import numpy as np
 import psutil
+from functools import partial
+
+from typing import Union, Callable
 
 from copy import deepcopy
 from collections import OrderedDict
@@ -31,7 +34,6 @@ def upload_simple_tokens(labels, cache, tensors, grid_setting = False):
         else:
             label_completed = (label, (1.0,))
         cache.add(label_completed, tensors[idx])
-        # cache.add(label_completed, tensors[idx])        
         cache.add_base_matrix(label_completed)
 
 
@@ -51,61 +53,55 @@ def upload_grids(grids, cache):
 
 def np_ndarray_section(matrix, boundary = None, except_idx : list = []):
     #Добавить проверку на случай, когда boundary = 0, чтобы тогда возвращало матрицу без изменения
-    if isinstance(boundary, int):
-        if boundary != 0:
-            for idx in except_idx:
-                if idx < 0:
-                    except_idx.append(matrix.ndim - 1)
-            for dim_idx in np.arange(matrix.ndim):
-                if not dim_idx in except_idx:
-                    matrix = np.moveaxis(matrix, source = dim_idx, destination=0)
-                    matrix = matrix[boundary:-boundary, ...]
-                    matrix = np.moveaxis(matrix, source = 0, destination=dim_idx)
-    elif isinstance(boundary, (list, tuple)):
-        if len(boundary) != matrix.ndim:
-            raise IndexError('Sizes of boundary do not match the dimensionality of data')
-        for idx in except_idx:
-            if idx < 0:
-                except_idx.append(matrix.ndim - 1)
-        for dim_idx in np.arange(matrix.ndim):
-            if not dim_idx in except_idx and boundary[dim_idx] > 0:
-                matrix = np.moveaxis(matrix, source = dim_idx, destination=0)
-                matrix = matrix[boundary[dim_idx]:-boundary[dim_idx], ...]
-                matrix = np.moveaxis(matrix, source = 0, destination=dim_idx)    
-    return matrix
+    raise DeprecationWarning('No data section shall be used during equation discovery!')
+    # if isinstance(boundary, int):
+    #     if boundary != 0:
+    #         for idx in except_idx:
+    #             if idx < 0:
+    #                 except_idx.append(matrix.ndim - 1)
+    #         for dim_idx in np.arange(matrix.ndim):
+    #             if not dim_idx in except_idx:
+    #                 matrix = np.moveaxis(matrix, source = dim_idx, destination=0)
+    #                 matrix = matrix[boundary:-boundary, ...]
+    #                 matrix = np.moveaxis(matrix, source = 0, destination=dim_idx)
+    # elif isinstance(boundary, (list, tuple)):
+    #     if len(boundary) != matrix.ndim:
+    #         raise IndexError('Sizes of boundary do not match the dimensionality of data')
+    #     for idx in except_idx:
+    #         if idx < 0:
+    #             except_idx.append(matrix.ndim - 1)
+    #     for dim_idx in np.arange(matrix.ndim):
+    #         if not dim_idx in except_idx and boundary[dim_idx] > 0:
+    #             matrix = np.moveaxis(matrix, source = dim_idx, destination=0)
+    #             matrix = matrix[boundary[dim_idx]:-boundary[dim_idx], ...]
+    #             matrix = np.moveaxis(matrix, source = 0, destination=dim_idx)    
+    # return matrix
 
 
-def prepare_var_tensor(var_tensor, derivs_tensor, time_axis, boundary, cut_except = []):
+def prepare_var_tensor(var_tensor, derivs_tensor, time_axis, cut_except = []):
     initial_shape = var_tensor.shape
     print('initial_shape', initial_shape, 'derivs_tensor.shape', derivs_tensor.shape)
     var_tensor = np.moveaxis(var_tensor, time_axis, 0)
-    if isinstance(boundary, int):
-        result = np.ones((1 + derivs_tensor.shape[-1], ) + tuple([shape - 2*boundary for shape in var_tensor.shape]))
-    elif isinstance(boundary, (tuple, list)):
-        result = np.ones((1 + derivs_tensor.shape[-1], ) + tuple([shape - 2*boundary[dim_idx] for dim_idx, 
-                                                                  shape in enumerate(var_tensor.shape)]))
-    else:
-        raise TypeError('')
+    result = np.ones((1 + derivs_tensor.shape[-1], ) + tuple([shape for shape in var_tensor.shape])) #  - 2*boundary
+
     increment = 1
-    print('var_tensor.shape', var_tensor.shape, np_ndarray_section(var_tensor, boundary, cut_except).shape, boundary)
-    result[increment - 1, :] = np_ndarray_section(var_tensor, boundary, cut_except)
+    result[increment - 1, :] = var_tensor#, boundary, cut_except)
     if derivs_tensor.ndim == 2:
         for i_outer in range(0, derivs_tensor.shape[1]):
-            result[i_outer+increment, ...] = np.moveaxis(np_ndarray_section(derivs_tensor[:, i_outer].reshape(initial_shape), boundary, cut_except),
+            result[i_outer+increment, ...] = np.moveaxis(derivs_tensor[:, i_outer].reshape(initial_shape), # np_ndarray_section( , boundary, cut_except)
                          source=time_axis, destination=0)
     else:
         for i_outer in range(0, derivs_tensor.shape[-1]):
-#            print(np_ndarray_section(derivs_tensor[..., i_outer], boundary, []).shape, result[i_outer+1, ...].shape)
             assert derivs_tensor.ndim == var_tensor.ndim + increment, 'The shape of tensor of derivatives does not correspond '
-            result[i_outer+increment, ...] = np.moveaxis(np_ndarray_section(derivs_tensor[..., i_outer], boundary, []), #-1,
+            result[i_outer+increment, ...] = np.moveaxis(derivs_tensor[..., i_outer], # np_ndarray_section( -1, , boundary, [])
                              source=time_axis, destination=0)    
     return result
 
 
-def download_variable(var_filename, deriv_filename, boundary, time_axis):
+def download_variable(var_filename, deriv_filename, time_axis):
     var = np.load(var_filename)
     derivs = np.load(deriv_filename)
-    tokens_tensor = prepare_var_tensor(var, derivs, time_axis, boundary)
+    tokens_tensor = prepare_var_tensor(var, derivs, time_axis)
     return tokens_tensor
 
 
@@ -117,6 +113,7 @@ def __init__(self):
         self.memory_structural = dict()
         self.mem_prop_set = False
         self.base_tensors = [] #storage of non-normalized tensors, that will not be affected by change of variables
+
 
     def use_structural(self, use_base_data = True, label = None, replacing_data = None):
         assert use_base_data or not type(replacing_data) == type(None), 'Structural data must be declared with base data or by additional tensors.'
@@ -160,11 +157,22 @@ def use_structural(self, use_base_data = True, label = None, replacing_data = No
                         raise ValueError('Shapes of tensors in new structural data do not match their counterparts in the base data')
                 self.memory_structural[label] = replacing_data
 
+    @property
+    def g_func(self): # , g_func : Union[Callable, type(None)] = None
+        assert '0' in self.memory_default.keys()    # Check if we are working with the grid cache
+        return self._g_func(self.get_all()[1])
+
+    @g_func.setter
+    def g_func(self, function : Callable):
+        self._g_func = function
+
     def add_base_matrix(self, label):
         assert label in self.memory_default.keys()
-        self.base_tensors.append(label)#$] = deepcopy(self.memory[label])
-#        if self.scale_used:
-#            self.base_tensors_scaled[label] = deepcopy(self.memory_scaled[label])
+        self.base_tensors.append(label)
+
+    def set_boundaries(self, boundary_width : Union[int, list]):
+        assert '0' in self.memory_default.keys(), 'Boundaries should be specified for grid cache.'
+        self.boundary_width = boundary_width
 
     def memory_usage_properties(self, obj_test_case = None, mem_for_cache_frac = None, mem_for_cache_abs = None):
         '''
@@ -216,14 +224,12 @@ def change_variables(self, increment, increment_structral = None):
         '''
         assert not (increment_structral is None and not all(self.structural_and_base_merged)), 'Not all structural data taken from the default, and the increment for structural was not sent' 
         random_key = list(self.memory_default.keys())[0]
-#        print(random_key)
         increment = np.reshape(increment, newshape=self.memory_default[random_key].shape)
         del self.memory_normalized
         self.memory_default = {key : self.memory_default[key] for key in self.base_tensors} #deepcopy(self.base_tensors); 
         self.memory_structural = {key : self.memory_structural[key] for key in self.base_tensors} #deepcopy(self.base_tensors_structural) 
         self.memory_normalized = dict()
         for key in self.memory_default.keys():
-#            print(self.memory[key].shape, increment.shape)
             assert np.all(self.memory_default[key].shape == increment.shape) 
             self.memory_default[key] = self.memory_default[key] - increment
             if not self.structural_and_base_merged[key]:
@@ -235,8 +241,6 @@ def add(self, label, tensor, normalized : bool = False, structural : bool = Fals
         Method for addition of a new tensor into the cache. Returns True if there was enough memory and the tensor was save, and False otherwise. 
         '''
         assert not (normalized and structural), 'The added matrix can not be simultaneously normalized and structural. Possibly, bug in token/term saving'
-#        assert not scaled or self.scale_used, 'Trying to add scaled data, while the cache it not allowed to get it'
-#        assert not structural, 'the structural data must be added with cache.use_structural method'
         if normalized:
             if self.max_allowed_tensors is None:
                 self.memory_usage_properties(obj_test_case = tensor, mem_for_cache_frac = 5)
@@ -293,6 +297,7 @@ def get(self, label, normalized = False, structural = False, saved_as = None):
         assert not (normalized and structural), 'The added matrix can not be simultaneously normalized and scaled'
 #        assert not scaled or self.scale_used, 'Trying to add scaled data, while the cache it not allowed to get it'
         if label is None:
+            print(self.memory_default.keys())
             return np.random.choice(list(self.memory_default.values()))
         if normalized:
             try:

epde/decorators.py

@@ -6,12 +6,17 @@
 @author: mike_ubuntu
 """
 
+import numpy as np
+
 from functools import wraps
+from typing import Union
+
+import epde.globals as global_var
 
 changelog_entry_templates = {}
 
-class Reset_equation_status():
-    def __init__(self, reset_input : bool = True, reset_output : bool = False, 
+class ResetEquationStatus():
+    def __init__(self, reset_input : bool = False, reset_output : bool = False, 
                  reset_right_part : bool = True):
         self.reset_input = reset_input; self.reset_output = reset_output
         self.reset_right_part = reset_right_part
@@ -35,17 +40,32 @@ def wrapper(obj, *args, **kwargs):
                         except AttributeError:
                             pass
             if self.reset_output:
-                try:
+                if isinstance(result, (list, tuple, set)):
                     for equation in result:
                         try: 
-                            equation.reset_state(self.reset_right_part )
+                            equation.reset_state(self.reset_right_part)
                         except AttributeError:
                             pass
-                except TypeError:
-                    try: 
-                        result.reset_state(self.reset_right_part )
+                else:
+                    try:
+                        result.reset_state(self.reset_right_part)
                     except AttributeError:
-                        pass                    
+                        pass              
+
+                # print(f'result in resetter: {result}')
+                # try:
+                #     for equation in result:
+                #         try: 
+                #             print(f'Here I must reset state of {equation}')
+                #             equation.reset_state(self.reset_right_part)
+                #         except AttributeError:
+                #             pass
+                # except TypeError:
+                #     try:
+                #         print(f'Here I must reset state of {result}')
+                #         result.reset_state(self.reset_right_part)
+                #     except AttributeError:
+                #         pass                    
             return result
         return wrapper
 
@@ -89,17 +109,26 @@ def historized(h_obj):
             return result
         return wrapper
 
-# class Parallelize_method():
-#     def __init__(self, num_cpus = 1):
-#         self.num_cpus = num_cpus
-
-#     def __call__(self, method):
-#         @wraps(method)
-#         @ray.remote(num_cpus = self.num_cpus)
-#         def wrapper(*args, **kwargs):
-#             return method(*args, **kwargs)
-
+class Boundary_exclusion():
+    def __init__(self, boundary_width = 0):
+        self.boundary_width = boundary_width
 
-#class Ray_parallelizer():
-#    def __init__(self):
-#        
+    def __call__(self, func):
+        @wraps(func)
+        def wrapper(grids, boundary_width : Union[int, list] = 0):
+            assert len(grids) == grids[0].ndim
+            if isinstance(self.boundary_width, int): self.boundary_width = len(grids)*[self.boundary_width,]
+            indexes_shape = grids[0].shape
+            indexes = np.indices(indexes_shape)
+
+            mask_partial = np.array([np.where((indexes[idx, ...] >= self.boundary_width[idx]) &
+                                              (indexes[idx, ...] < indexes_shape[idx] - self.boundary_width[idx]),
+                                              1, 0)
+                                     for idx in np.arange(indexes.shape[0])])
+
+            mask = np.multiply.reduce(mask_partial, axis = 0)
+            g_function_res = func(grids)
+            assert np.shape(g_function_res) == np.shape(mask)
+            return func(grids) * mask
+
+        return wrapper

epde/eq_mo_objectives.py

@@ -9,46 +9,48 @@
 import numpy as np
 from functools import partial
 
-def generate_partial(obj_function, equation_idx):
-    return partial(obj_function, equation_idx = equation_idx)
+def generate_partial(obj_function, equation_key):
+    return partial(obj_function, equation_key = equation_key)
 
-def equation_discrepancy(system, equation_idx):
+def equation_fitness(system, equation_key):
     '''
-    Evaluate the discrepancy of the system of PDEs, using sum of the L2 norm of the discrepancy 
-    of each equation in the system from zero.
+    Evaluate the quality of the system of PDEs, using the individual values of fitness function for equations.
 
     Parameters:
     -----------
-        system - ``epde.structure.SoEq`` object
+        system - ``epde.structure.main_structures.SoEq`` object
         The system, that is to be evaluated.
 
     Returns:
     ----------
-        discrepancy : float.
+        error : float.
         The value of the error metric.
     '''
-    res = np.sum(np.abs(system.structure[equation_idx].evaluate(normalize = False, return_val = True)[0]))
+    # print(f'System, for which we evaluate fitness: {system.text_form}')
+    # print(f'For equation key {equation_key}, {system.vals[equation_key].fitness_calculated}')
+    assert system.vals[equation_key].fitness_calculated, 'Trying to call fitness before its evaluation.'
+    res = system.vals[equation_key].fitness_value
     return res
 
-def equation_complexity_by_terms(system, equation_idx):
+def equation_complexity_by_terms(system, equation_key):
     '''
     Evaluate the complexity of the system of PDEs, evaluating a number of terms for each equation. 
     In the evaluation, we consider only terms with non-zero weights, and the target term with the free 
     coefficient are not included in the final metric due to their ubiquty in the equations.
 
     Parameters:
     -----------
-        system - ``epde.structure.SoEq`` object
+        system - ``epde.structure.main_structures.SoEq`` object
         The system, that is to be evaluated.
 
     Returns:
     ----------
         discrepancy : list of integers.
         The values of the error metric: list entry for each of the equations.
     '''    
-    return np.count_nonzero(system.structure[equation_idx].weights_internal)
+    return np.count_nonzero(system.vals[equation_key].weights_internal)
 
-def equation_complexity_by_factors(system, equation_idx):
+def equation_complexity_by_factors(system, equation_key):
     '''
     Evaluate the complexity of the system of PDEs, evaluating a number of factors in terms for each 
     equation. In the evaluation, we consider only terms with non-zero weights and target, while
@@ -57,7 +59,7 @@ def equation_complexity_by_factors(system, equation_idx):
 
     Parameters:
     -----------
-        system - ``epde.structure.SoEq`` object
+        system - ``epde.structure.main_structures.SoEq`` object
         The system, that is to be evaluated.
 
     Returns:
@@ -66,15 +68,13 @@ def equation_complexity_by_factors(system, equation_idx):
         The values of the error metric: list entry for each of the equations.
     '''    
     eq_compl = 0
-    # print(equation)
-    # print(equation.text_form)
 
-    for idx, term in enumerate(system.structure[equation_idx].structure):
-        if idx < system.structure[equation_idx].target_idx:
-            if not system.structure[equation_idx].weights_final[idx] == 0: 
+    for idx, term in enumerate(system.vals[equation_key].structure):
+        if idx < system.vals[equation_key].target_idx:
+            if not system.vals[equation_key].weights_final[idx] == 0: 
                 eq_compl += len(term.structure)
-        elif idx > system.structure[equation_idx].target_idx:
-            if not system.structure[equation_idx].weights_final[idx-1] == 0: eq_compl += len(term.structure)
+        elif idx > system.vals[equation_key].target_idx:
+            if not system.vals[equation_key].weights_final[idx-1] == 0: eq_compl += len(term.structure)
         else:
             eq_compl += len(term.structure)
     return eq_compl