Minor version bump (to 1.2.3; release 1.2.2 has been removed).

- Added a max_retries argument to explore_model_space to persist on optimiser numerical errors
- Minor refactoring and changes around the generation of the sum-of-products form of kernels (some clarifying assertions and caching of generated forms since multiple methods triggered them)
- Minor documentation changes

.github/workflows/pythonpackage.yml

@@ -10,7 +10,7 @@ jobs:
       max-parallel: 4
       matrix:
         os: [windows-latest, ubuntu-latest, macOS-latest]
-        python-version: [3.8]
+        python-version: ["3.10"]
 
     steps:
     - uses: actions/checkout@v1

GPy_ABCD/KernelExpansion/kernelOperations.py

@@ -1,3 +1,5 @@
+from copy import deepcopy
+
 from GPy_ABCD.Kernels.baseKernels import *
 from GPy_ABCD.Kernels.baseKernels import __FIX_SIGMOIDAL_KERNELS_SLOPE
 
@@ -33,7 +35,8 @@ def fit_ker_to_kex_with_params(ker, kex, verbose = False):
     return kex.match_up_fit_parameters(param_dict, 'GP_regression.')
 
 
-def init_rand_params(kex, verbose = True): # A testing function to initialise a kernel with random parameters
+def init_rand_params(kex, verbose = True):
+    '''A testing function to initialise a kernel with random parameters'''
     if verbose: print(kex)
     ker = kex._initialise().to_kernel()
     ker.randomize()

GPy_ABCD/KernelExpressions/base.py

@@ -53,10 +53,10 @@ def simplify(self):
     def extract_if_singleton(self):
         pass
 
-    # NOTE: both traverse and reduce ignore raw-string leaves (which can only happen in ChangeKEs);
-    #       care has to be taken to perform required operations on them from their parent
     @abstractmethod
     def traverse(self):
+        '''NOTE: both traverse and reduce ignore raw-string leaves (which can only happen in ChangeKEs);
+        care has to be taken to perform required operations on them from their parent'''
         pass
 
     # Same note as traverse; see Test.checkKernelExpressions for an example func
@@ -98,7 +98,8 @@ def new_tree_with_self_replaced(self, replacement_node): # NOTE: replacement_nod
             return copied_replacement.parent.reassign_child(copied_self, copied_replacement)
 
     @abstractmethod
-    def reassign_child(self, old_child, new_child): # NOTE: has to return new_child (used by new_tree_with_self_replaced)
+    def reassign_child(self, old_child, new_child):
+        '''NOTE: has to return new_child (used by new_tree_with_self_replaced)'''
         pass
 
     @abstractmethod
@@ -128,11 +129,16 @@ def _new_parameters(self, new_parameters):
         return self
 
     @abstractmethod
-    def match_up_fit_parameters(self, fit_ker, prefix): # Note: the prefix has to already contain THIS node's name followed by a dot at the end
+    def match_up_fit_parameters(self, fit_ker, prefix):
+        '''NOTE: the prefix has to already contain THIS node's name followed by a dot at the end'''
         pass
 
     @abstractmethod
-    def sum_of_prods_form(self): # Return either a ProductKE or a SumKE whose composite_terms are only ProductKEs
+    def sum_of_prods_form(self):
+        '''Return either a ProductKE or a SumKE whose composite_terms are only ProductKEs.
+
+        NOTE: this method CAN only be called when parameters are present (i.e. after .match_up_fit_parameters has been called),
+        and SHOULD only be called indirectly through GPModel.sum_of_prods_kex or GPModel.interpret()'''
         pass
 
     def get_interpretation(self, sops = None):

GPy_ABCD/KernelExpressions/change.py

@@ -115,6 +115,7 @@ def add_sum_of_prods_terms(k1, k2):
         return res._set_all_parents()
 
     def sum_of_prods_form(self):
+        assert self.parameters, 'A sum-of-products form can only be generated when parameters are present (i.e. after .match_up_fit_parameters has been triggered), and should only be called indirectly through GPModel.sum_of_prods_kex or GPModel.interpret()'
         new_children = []
         for branch, kex in (('left', self.left), ('right', self.right)):
             sigmoid_parameters = (change_k_sigmoid_names[self.CP_or_CW][branch], self.parameters[self.CP_or_CW][0])

GPy_ABCD/KernelExpressions/commutatives.py

@@ -36,6 +36,7 @@ def simplify_base_terms_params(self):
         return self
 
     def sum_of_prods_form(self):
+        assert self.parameters, 'A sum-of-products form can only be generated when parameters are present (i.e. after .match_up_fit_parameters has been triggered), and should only be called indirectly through GPModel.sum_of_prods_kex or GPModel.interpret()'
         cts = [ct.sum_of_prods_form() for ct in self.composite_terms]
         self.composite_terms.clear()
         for ct in cts: # Only SumKEs or ProductKEs now
@@ -133,6 +134,7 @@ def multiply_pure_prods_with_params(k0, ks): # k0 is meant to be the ProductKE c
         return ProductKE([]).new_bases_with_parameters([(key, p) for kex in [k0] + ks for key, ps in list(kex.parameters.items()) for p in ps])
 
     def sum_of_prods_form(self):
+        assert self.parameters, 'A sum-of-products form can only be generated when parameters are present (i.e. after .match_up_fit_parameters has been triggered), and should only be called indirectly through GPModel.sum_of_prods_kex or GPModel.interpret()'
         sops = SumKE([])
         if not self.composite_terms:
             sops.composite_terms.append(ProductKE(self.base_terms)._new_parameters(self.parameters)) # Avoid triggering simplify()

GPy_ABCD/Models/model.py

@@ -12,6 +12,7 @@ def __init__(self, X, Y, kernel_expression = SumKE(['WN'])._initialise()):
         self.X = X
         self.Y = Y
         self.kernel_expression = kernel_expression
+        self._sum_of_prods_kex = None
         self.restarts = None
         self.model = None
         self.cached_utility_function = None
@@ -30,7 +31,20 @@ def fit(self, restarts = None, optimiser = 'lbfgsb', verbose = False, robust = F
             self.model.optimize_restarts(num_restarts = self.restarts, verbose = verbose, robust = robust, optimizer = optimiser, **kwargs)
         return self
 
-    def interpret(self): return fit_ker_to_kex_with_params(self.model.kern, deepcopy(self.kernel_expression)).get_interpretation()
+    @property
+    def sum_of_prods_kex(self):
+        '''The canonical kernel form (the one described in .interpret).
+
+        NOTE: this property/method can only be called after the model has been fitted.'''
+        if self.model is None: raise ValueError('No parameters to insert into the kernel expression since the model has not yet been fitted')
+        elif self._sum_of_prods_kex is None: self._sum_of_prods_kex = fit_ker_to_kex_with_params(self.model.kern, deepcopy(self.kernel_expression)).sum_of_prods_form()
+        return self._sum_of_prods_kex
+
+    def interpret(self):
+        '''Describe the model with a few sentences (which break down the expanded kernel form, i.e. .sum_of_prods_kex).
+
+        NOTE: this method can only be called after the model has been fitted.'''
+        return self.sum_of_prods_kex.get_interpretation(sops = self._sum_of_prods_kex)
 
     def predict(self, X, quantiles = (2.5, 97.5), full_cov = False, Y_metadata = None, kern = None, likelihood = None, include_likelihood = True):
         mean, cov = self.model.predict(X, full_cov, Y_metadata, kern, likelihood, include_likelihood)

GPy_ABCD/Models/modelSearch.py

@@ -2,6 +2,7 @@
 from copy import deepcopy
 from operator import attrgetter
 from multiprocessing import Pool, cpu_count
+from warnings import warn
 from typing import Callable, List, Dict, Tuple, Union
 
 from GPy_ABCD.KernelExpressions.base import KernelExpression
@@ -13,8 +14,6 @@
 # TODO:
 #   - group input parameters into dictionaries?
 #   - make utility functions take in the kernel expression (or even full model) so that further criteria may be applied (e.g. presence of specific kernels etc)?
-#   - focus on documenting end-user and generic developer functions etc in sphinx
-#   - make the dynamic buffer configurable, or even allow inputting a list of numbers of models to keep per round
 #   - make an interactive/interruptable mode which asks whether to go further, retaining how many etc
 #   - allow the model lists in each round to be fit in batches, with interactive request to continue (timed response maybe)
 #   - show a live count of models fitted so far in each round (probably by batches)
@@ -24,24 +23,32 @@ def print_k_list(k_or_model_list):
     return ', '.join([str(m.kernel_expression) for m in k_or_model_list] if isinstance(k_or_model_list[0], GPModel) else [str(k) for k in k_or_model_list])
 
 
-def fit_one_model(X, Y, kex, restarts, optimiser = GPy_optimisers[0]): return GPModel(X, Y, kex).fit(restarts, optimiser)
+def fit_one_model(X, Y, kex, restarts, optimiser = GPy_optimisers[0], max_retries = 1):
+    retry = 0
+    while True:
+        try: return GPModel(X, Y, kex).fit(restarts, optimiser)
+        except Exception as e:
+            if retry < max_retries:
+                retry += 1
+                warn(f'Retry #{retry} (max is {max_retries}) for kernel {kex} due to error:\n\t{e}')
+            else: raise e
 
-def fit_mods_not_parallel(X, Y, k_exprs, restarts = 5, optimiser = GPy_optimisers[0]):
+def fit_mods_not_parallel(X, Y, k_exprs, restarts = 5, optimiser = GPy_optimisers[0], max_retries = 1):
     '''Fit models from a list of kernels to the same data NOT in parallel'''
-    return [fit_one_model(X, Y, kex, restarts, optimiser) for kex in k_exprs]
-def fit_mods_parallel_processes(X, Y, k_exprs, restarts = 5, optimiser = GPy_optimisers[0]):
+    return [fit_one_model(X, Y, kex, restarts, optimiser, max_retries) for kex in k_exprs]
+def fit_mods_parallel_processes(X, Y, k_exprs, restarts = 5, optimiser = GPy_optimisers[0], max_retries = 1):
     '''Concurrently fit models (from a list of kernels) on the same data in as many processes as there are available processor cores using `multiprocessing`'s `Pool`
     (this function does NOT use GPy's own `parallel` argument, which parallelises within single fits rather than across multiple ones)
 
-    NOTE: Any code calling this function (however many nested levels in) should be within a "if __name__ == '__main__':" preamble for full OS-agnostic use.
+    NOTE: Any code calling this function (however many nested levels in) should be within an "if __name__ == '__main__':" preamble for full OS-agnostic use.
     '''
-    with Pool() as pool: return pool.starmap_async(fit_one_model, [(X, Y, kex, restarts, optimiser) for kex in k_exprs], int(len(k_exprs) / cpu_count()) + 1).get()
+    with Pool() as pool: return pool.starmap_async(fit_one_model, [(X, Y, kex, restarts, optimiser, max_retries) for kex in k_exprs], int(len(k_exprs) / cpu_count()) + 1).get()
 
 
 def explore_model_space(X, Y,
                         start_kernels: List[Union[str, KernelExpression]] = start_kernels['Default'], p_rules: List[Callable] = production_rules['Default'],
                         utility_function: Callable = BIC, rounds: int = 2, beam: Union[int, List[int]] = [3, 2, 1], restarts: int = 5,
-                        model_list_fitter: Callable = fit_mods_parallel_processes, optimiser: str = GPy_optimisers[0],
+                        model_list_fitter: Callable = fit_mods_parallel_processes, optimiser: str = GPy_optimisers[0], max_retries: int = 1,
                         verbose: bool = True) -> Tuple[List[GPModel], List[List[GPModel]], List[KernelExpression], List[GPModel], List[GPModel]]:
     '''Perform `rounds` rounds of kernel expansion followed by model fit starting from the given `start_kernels` with and expanding the best `buffer` of them with `p_rules` production rules
 
@@ -57,12 +64,14 @@ def explore_model_space(X, Y,
     :type rounds: Int
     :param beam: number of best fit-models' kernels to expand each round, either an integer or a list of integers; in the latter case, if its length is less than rounds then the last value will be repeated until required
     :type beam: Union[Int, List[Int]]
-    :param restarts: number of GPy model-fitting restarts with different parameters
+    :param restarts: number of GPy model-fitting restarts with random starting parameters
     :type restarts: Int
     :param model_list_fitter: function handling the fitting of a list of kernels to the same data; this is where parallelisation implementation can be changed
     :type model_list_fitter: Callable
     :param optimiser: identifying string for the model optimiser function; GPy 1.9.9 optimiser strings (GPy > paramz > optimization > optimization.py): 'lbfgsb', 'org-bfgs', 'fmin_tnc', 'scg', 'simplex', 'adadelta', 'rprop', 'adam'
     :type optimiser: str
+    :param max_retries: number of retries in case of numerical errors; independent of number of restarts (which are always executed)
+    :type max_retries: Int
     :param verbose: produce verbose logs
     :type verbose: Boolean
 
@@ -76,7 +85,7 @@ def explore_model_space(X, Y,
     if verbose: print(f'Testing {rounds} layers of model expansion on {len(X)} datapoints, starting from: {print_k_list(start_kexs)}\nModels are fitted with {restarts} random restarts (with {optimiser} optimiser) and scored by {utility_function.__name__}')
     def score(m): return m.compute_utility(utility_function)
 
-    tested_models = [sorted(model_list_fitter(X, Y, start_kexs, restarts, optimiser), key = score)]
+    tested_models = [sorted(model_list_fitter(X, Y, start_kexs, restarts, optimiser, max_retries), key = score)]
     sorted_models = not_expanded = sorted(flatten(tested_models), key = score)
     expanded = []
     tested_k_exprs = deepcopy(start_kexs)
@@ -91,15 +100,15 @@ def score(m): return m.compute_utility(utility_function)
 
     sorted_models, tested_models, tested_k_exprs, expanded, not_expanded = model_search_rounds(X, Y,
                         sorted_models, tested_models, tested_k_exprs, expanded, not_expanded,
-                        model_list_fitter, p_rules, utility_function, rounds, beam, restarts, optimiser, verbose)
+                        model_list_fitter, p_rules, utility_function, rounds, beam, restarts, optimiser, max_retries, verbose)
 
     if verbose: print(f'\nBest models overall: {print_k_list(sorted_models[:beam[0]])}\n')
     return sorted_models, tested_models, tested_k_exprs, expanded, not_expanded
 
 
 # This function is split from the above both for tidiness and to allow the possibility of continuing a search if desired
 def model_search_rounds(X, Y, sorted_models, tested_models, tested_k_exprs, expanded, not_expanded,
-                        model_list_fitter, p_rules, utility_function, rounds, beam, restarts, optimiser, verbose):
+                        model_list_fitter, p_rules, utility_function, rounds, beam, restarts, optimiser, max_retries, verbose):
     '''
     This function can be used to continue a previous search: the 5 arguments after Y are the outputs of explore_model_space.
     See explore_model_space description and definition for argument explanation and context.
@@ -110,7 +119,7 @@ def score(m): return m.compute_utility(utility_function)
 
     for r, b in zip(range(1, rounds + 1), beam):
         new_k_exprs = [kex for kex in unique(flatten([expand(mod.kernel_expression, p_rules) for mod in not_expanded[:b]])) if kex not in tested_k_exprs]
-        tested_models.append(sorted(model_list_fitter(X, Y, new_k_exprs, restarts, optimiser), key = score))  # tested_models[d]
+        tested_models.append(sorted(model_list_fitter(X, Y, new_k_exprs, restarts, optimiser, max_retries), key = score))  # tested_models[d]
 
         sorted_models = sorted(flatten(tested_models), key = attrgetter('cached_utility_function')) # Merge-sort would be applicable
         expanded += not_expanded[:b]

GPy_ABCD/__init__.py

@@ -1,6 +1,6 @@
 """GPy-ABCD - Basic implementation with GPy of an Automatic Bayesian Covariance Discovery (ABCD) system"""
 
-__version__ = '1.2.1' # Change it in setup.py too
+__version__ = '1.2.3' # Change it in setup.py too
 __author__ = 'Thomas Fletcher <T-Fletcher@outlook.com>'
 # __all__ = []
 

README.rst

@@ -21,17 +21,25 @@ GPy-ABCD
     :target: https://github.com/T-Flet/GPy-ABCD/blob/master/LICENSE
     :alt: License
 
-Basic implementation with GPy of an Automatic Bayesian Covariance Discovery (ABCD) system
+*(Temporary note: the "failing" build badge above is due to the workflow pip not finding GPy 1.12.0 for some reason; the tests are successful)*
 
-Briefly: ABCD is a modelling system which consists in exploring a space of compositional kernels
-(i.e. covariances of Gaussian Processes) constructed by iteratively combining a small set of base ones,
-returning the best fitting models using them, and capable of generating simple text descriptions of the
+GPy-ABCD is a basic implementation with GPy of an Automatic Bayesian Covariance Discovery (ABCD) system.
+
+Briefly, ABCD is a (Gaussian Process) modelling method which consists in exploring a space of modular kernels
+(i.e. covariances) constructed by iteratively combining a small set of simple ones,
+and returning the best fitting models using them;
+due to its modularity, it is capable of generating simple text descriptions of the
 fits based on the identified functional shapes.
 
+The usefulness of ABCD is in identifying the underlying shape of data, but the process is
+computationally expensive, therefore a typical use for it is in initial data overviews
+(possibly on subsampled datasets for efficiency), then followed by more direct exploration
+of its top few results' kernels (on the full dataset if subsampled before).
+
 See the picture in `Usage` below for an example input/output and read the paper for further details:
 
 `Fletcher, T., Bundy, A., & Nuamah, K. . GPy-ABCD: A Configurable Automatic Bayesian Covariance Discovery Implementation.
-8th ICML Workshop on Automated Machine Learning (2021) <https://sites.google.com/view/automl2021/accepted-papers>`_
+8th ICML Workshop on Automated Machine Learning (2021) <https://www.research.ed.ac.uk/en/publications/gpy-abcd-a-configurable-automatic-bayesian-covariance-discovery-i>`_
 
 
 
@@ -67,8 +75,8 @@ A minimal example to showcase the various parameters follows:
         best_mods, all_mods, all_exprs, expanded, not_expanded = explore_model_space(X, Y,
             start_kernels = start_kernels['Default'], p_rules = production_rules['Default'],
             utility_function = BIC, rounds = 2, beam = [3, 2, 1], restarts = 5,
-            model_list_fitter = fit_mods_parallel_processes, optimiser = GPy_optimisers[0],
-            verbose = True)
+            model_list_fitter = fit_mods_parallel_processes,
+            optimiser = GPy_optimisers[0], max_retries = 1, verbose = True)
 
         print('\nFull lists of models by round:')
         for mod_depth in all_mods: print(', '.join([str(mod.kernel_expression) for mod in mod_depth]) + f'\n{len(mod_depth)}')
@@ -81,9 +89,12 @@ A minimal example to showcase the various parameters follows:
         from matplotlib import pyplot as plt
         plt.show() # Not required for plotly = True above
 
+        # Later fit one of the top kernels to new data (e.g. if the original data was subsampled)
+        #fit_kex(other_X, other_Y, best_mods[0].kernel_expression, 10)
 
 
-.. figure:: selected_output_example.png
+
+.. figure:: output_example.png
     :align: center
     :figclass: align-center
 
@@ -114,9 +125,7 @@ see the definitions of entries of the :code:`start_kernels` and :code:`productio
 Project Structure
 -----------------
 
-A paper on GPy-ABCD and its differences from the original ABCD is planned, but for the time being read the paper mentioned above for a full picture of what an ABCD system is.
-
-However, briefly, it consists in exploring a space of compositional kernels built from a few carefully selected base ones,
+Briefly, the ABCD process consists in exploring a space of compositional kernels built from a few carefully selected base ones,
 returning the best fitting models using them and generating simple text interpretations of the fits based
 on the functional shapes of the final composed covariance kernels and parameter values.
 
@@ -159,7 +168,10 @@ Further Notes
 Generic:
 
 - Please reach out if you have successfully used this project in your own research
-- Feel free to fork and expand this project (pull requests are welcome) since it is not the focus of my research; it was written just because I needed to use it in a broader adaptive statistical modelling context and therefore I have no need to expand its functionality in the near future
+- Feel free to fork and expand this project (pull requests are welcome) since it
+  was only written because its particular features were required in a broader
+  adaptive statistical modelling context,
+  and therefore it is unlikely that its functionality will be expanded in the near future
 
 Code-related: