Added a try_higher_curves production rule (trying higher polynomials …

…or PER and discouraging SE if present)

Deemed the LIN with offset sufficiently better than without (it identifies the roots after all)

GPy_ABCD/KernelExpansion/grammar.py

@@ -1,5 +1,6 @@
 from GPy_ABCD.KernelExpansion.kernelExpressionOperations import *
 from GPy_ABCD.Util.genericUtil import *
+from GPy_ABCD.Kernels.baseKernels import __USE_LIN_KERNEL_HORIZONTAL_OFFSET
 
 
 ## Expansion functions
@@ -35,6 +36,8 @@ def change_point_linear(S): return [deep_apply(one_change, S, 'CP', 'LIN')] # No
 def times_shifted_base(S): return [deep_apply(multiply, S, SumKE([B, 'C'])) for B in base_kerns - {'C'}]
 def replace_with_singleton(S): return [deep_apply(replace_node, S, SumKE([B])) for B in base_kerns]
 def remove_some_term(S): return [deep_apply(remove_a_term, S)]
+def try_higher_curves(S): return [deep_apply(higher_curves, S)] # Not in original ABCD
+
 
 production_rules_by_type = {
     'basic': {
@@ -50,6 +53,7 @@ def remove_some_term(S): return [deep_apply(remove_a_term, S)]
         'times_shifted_base': times_shifted_base, # S -> S * (B + C)
         'replace_with_singleton': replace_with_singleton, # S -> B
         'remove_some_term': remove_some_term, # S + S2 -> S and S * S2 -> S
+        'try_higher_curves': try_higher_curves, # S -> S with higher polynomials or PERs, discouraging SE
     }
 }
 
@@ -68,19 +72,18 @@ def make_simple_kexs(pseudo_kexs): return [pseudo_to_real_kex(pkex)._initialise(
 
 # TODO: SumKE(['LIN', 'C']) instead of 'LIN'? Only for non-horizontal-offset-including version?
 standard_start_kernels = make_simple_kexs(list(base_kerns - {'SE'}) + # Base Kernels without SE
-                                          [SumKE(['LIN', 'C']), SumKE(['PER', 'C'])] + # More generic LIN and PER
-                                          # both_changes('LIN')) # To catch a possible changepoint or changewindow with simple enough shapes
-                                          [SumKE(['C'], [ck]) for ck in both_changes('LIN')]) # To catch a possible changepoint or changewindow with simple enough shapes
+                                          [ProductKE(['LIN', 'LIN']), ProductKE(['LIN', 'LIN', 'LIN']), SumKE(['PER', 'C'])] + # More generic LIN and PER
+                                          both_changes('LIN')) # To catch a possible changepoint or changewindow with simple enough shapes
 
-extended_start_kernels = make_simple_kexs(list(base_kerns) + # Base Kernels
-                                          [SumKE(['LIN', 'C']), SumKE(['PER', 'C'])] + # More generic LIN and PER
-                                          # both_changes(ProductKE(['LIN', 'LIN']))) # To catch a possible changepoint or changewindow with simple enough shapes
-                                          both_changes(ProductKE([], [SumKE(['LIN', 'C']), SumKE(['LIN', 'C'])]))) # To catch a possible changepoint or changewindow with simple enough shapes
+extended_start_kernels = make_simple_kexs(list(base_kerns - {'SE'}) + # Base Kernels without SE
+                                          [ProductKE(['LIN', 'LIN']), ProductKE(['LIN', 'LIN', 'LIN']), SumKE(['PER', 'C'])] + # More generic LIN and PER
+                                          [SumKE(['C'], [ck]) for ck in both_changes('LIN')]) # To catch a possible changepoint or changewindow with simple enough shapes
 
-test_start_kernels = make_simple_kexs(list(base_kerns - {'SE'}) + # Base Kernels without SE and PER
-                                          [SumKE(['PER', 'C'])] + # More generic LIN and PER
+test_start_kernels = make_simple_kexs(list(base_kerns) + # Base Kernels
+                                          [ProductKE(['LIN', 'LIN']), ProductKE(['LIN', 'LIN', 'LIN']), SumKE(['PER', 'C'])] + # More generic LIN and PER
                                           both_changes('LIN')) # To catch a possible changepoint or changewindow with simple enough shapes
-                                          # both_changes(SumKE(['LIN', 'C']))) # To catch a possible changepoint or changewindow with simple enough shapes
+
+if not __USE_LIN_KERNEL_HORIZONTAL_OFFSET: assert SumKE(['LIN', 'C']) in standard_start_kernels, 'Non-offset-LIN needs C to achieve (almost) parity with offset-LIN'
 
 
 # Production Rules

GPy_ABCD/KernelExpansion/kernelExpressionOperations.py

@@ -117,6 +117,27 @@ def remove_a_term(S):
         # elif not S.is_root(): print('THIS WAS HERE JUST TO VERIFY THAT IT NEVER HAPPENS (BY DESIGN); NOT THAT IT WOULD BREAK ANYTHING IF IT DID')
     elif isinstance(S, ChangeKE):
         res = [SumKE([branch]) if isinstance(branch, str) else branch for branch in (S.left, S.right)]
-    else: # elif isinstance(k1, str): # This never occurs in the expansion though
+    else: # elif isinstance(S, str): # This never occurs in the expansion though
+        res = [[]]
+    return res
+
+
+def higher_curves(S):
+    res = []
+    higher_kers = [ProductKE(['LIN', 'LIN']), ProductKE(['LIN', 'LIN', 'LIN'])]
+    if isinstance(S, SumKE):
+        res = [SumKE(+S.base_terms, S.composite_terms + [ct]) for ct in higher_kers]
+        # res += [SumKE((+S.base_terms) + Counter({'SE': -S.base_terms['SE'], bt: bt_c}), S.composite_terms) for bt, bt_c in [('LIN', 2), ('PER', 2)]]
+        if 'SE' in (+S.base_terms).keys(): # Discourage SE
+            for r in res: del r.base_terms['SE']
+    elif isinstance(S, ProductKE):
+        res = [ProductKE(+S.base_terms + Counter(bts), S.composite_terms) for bts in [{'LIN': 2}, {'LIN': 3}]]
+        res += [ProductKE(+S.base_terms + (Counter({'PER': 1}) if 'PER' in +S.base_terms else Counter({'PER': 2})), S.composite_terms)]
+        if 'SE' in (+S.base_terms).keys(): # Discourage SE
+            for r in res: del r.base_terms['SE']
+    elif isinstance(S, ChangeKE):
+        res += [ChangeKE(S.CP_or_CW, ct, S.right) for ct in higher_kers if isinstance(S.left, str)]
+        res += [ChangeKE(S.CP_or_CW, S.left, ct) for ct in higher_kers if isinstance(S.right, str)]
+    else: # elif isinstance(S, str): # This never occurs in the expansion though
         res = [[]]
     return res

GPy_ABCD/KernelExpansion/kernelInterpretation.py

@@ -42,13 +42,13 @@ def first_term_interpretation(bps):
     elif b == 'LIN':
         if len(ps) > 1:
             if __USE_LIN_KERNEL_HORIZONTAL_OFFSET:
-                res = 'A polynomial function of order {:d} with offsets '.format(len(ps))
+                res = 'A polynomial function of order {:d} with roots '.format(len(ps))
                 offsets = ['{:.2f}'.format(p['offset']) for p in ps]
                 if len(ps) > 1: offsets[-1] = 'and ' + offsets[-1]
                 res += ', '.join(offsets)
             else: res = 'A polynomial function of order {:d}'.format(len(ps))
         else:
-            if __USE_LIN_KERNEL_HORIZONTAL_OFFSET: res = 'A linear function with offset {:.2f}'.format(ps[0]['offset'])
+            if __USE_LIN_KERNEL_HORIZONTAL_OFFSET: res = 'A linear function with horizontal offset {:.2f}'.format(ps[0]['offset'])
             else: res = 'A linear function'
     else: raise ValueError(f'An unexpected type of first term in a pure product has arisen: {b}')
     return res
@@ -61,13 +61,13 @@ def postmodifier_interpretation(bps):
     elif b == 'LIN':
         if len(ps) > 1:
             if __USE_LIN_KERNEL_HORIZONTAL_OFFSET:
-                res = 'with polynomially varying amplitude of order {:d} with offsets '.format(len(ps))
+                res = 'with polynomially varying amplitude of order {:d} with roots '.format(len(ps))
                 offsets = ['{:.2f}'.format(p['offset']) for p in ps]
                 if len(ps) > 1: offsets[-1] = 'and ' + offsets[-1]
                 res += ', '.join(offsets)
             else: res = 'with polynomially varying amplitude of order {:d}'.format(len(ps))
         else:
-            if __USE_LIN_KERNEL_HORIZONTAL_OFFSET: res = 'with linearly varying amplitude with offset {:.2f}'.format(ps[0]['offset'])
+            if __USE_LIN_KERNEL_HORIZONTAL_OFFSET: res = 'with linearly varying amplitude with horizontal offset {:.2f}'.format(ps[0]['offset'])
             else: res = 'with linearly varying amplitude'
     elif b == 'sigmoidal_intervals':
         if not isinstance(ps, list): ps = [ps]

GPy_ABCD/Kernels/baseKernels.py

@@ -5,10 +5,10 @@
 
 
 #### CORE CONFIGURATION OF BASE KERNELS ####
-__INCLUDE_SE_KERNEL = True
-__USE_LIN_KERNEL_HORIZONTAL_OFFSET = False
-__USE_NON_PURELY_PERIODIC_PER_KERNEL = False
-__FIX_SIGMOIDAL_KERNELS_SLOPE = True
+__INCLUDE_SE_KERNEL = True # The most generic kernel; always a bargain in terms of parameters
+__USE_LIN_KERNEL_HORIZONTAL_OFFSET = True # Identifies the polynomial roots; more accurate but one extra parameter per degree
+__USE_NON_PURELY_PERIODIC_PER_KERNEL = False # Full standard periodic kernel [MacKay (1998)] instead of only its purely periodic part
+__FIX_SIGMOIDAL_KERNELS_SLOPE = True # Hence one parameter fewer for each sigmoidal and related kernel
 #### CORE CONFIGURATION OF BASE KERNELS ####
 
 
@@ -83,9 +83,11 @@ def SIOr(): return _Sk.SigmoidalIndicatorKernelOneLocation(1, True, fixed_slope
 
 CP = _Cs.ChangePointKernel
 def CP(left, right): return _Cs.ChangePointKernel(left, right, fixed_slope = __FIX_SIGMOIDAL_KERNELS_SLOPE)
-# CW = _Cs.ChangeWindowKernel
-# CW = _Cs.ChangeWindowKernelOneLocation
-CW = _Cs.ChangeWindowKernelWithWidth
+# # CW = _Cs.ChangeWindowKernel
+# def CW(left, right): return _Cs.ChangeWindowKernel(left, right, fixed_slope = __FIX_SIGMOIDAL_KERNELS_SLOPE)
+# # CW = _Cs.ChangeWindowKernelOneLocation
+# def CW(left, right): return _Cs.ChangeWindowKernelOneLocation(left, right, fixed_slope = __FIX_SIGMOIDAL_KERNELS_SLOPE)
+# CW = _Cs.ChangeWindowKernelWithWidth
 def CW(left, right): return _Cs.ChangeWindowKernelWithWidth(left, right, fixed_slope = __FIX_SIGMOIDAL_KERNELS_SLOPE)
 
 # CP = _CFs.kCP

Tests/regressLinearKernels.py

@@ -1,16 +1,29 @@
 import numpy as np
 
-from GPy_ABCD.Util.kernelUtil import doGPR
 from GPy_ABCD.Kernels.linearKernel import Linear
 from GPy_ABCD.Kernels.linearOffsetKernel import LinearWithOffset
-from GPy_ABCD.Kernels.baseKernels import C
+from GPy_ABCD.Kernels.baseKernels import C, __USE_LIN_KERNEL_HORIZONTAL_OFFSET
+from GPy_ABCD.KernelExpansion.kernelExpression import *
 from GPy_ABCD.Util.dataAndPlottingUtil import *
+from GPy_ABCD.Util.kernelUtil import doGPR, score_ps, BIC, AIC, AICc
+from GPy_ABCD.Models.modelSearch import fit_model_list_parallel
 
 
-X, Y = generate_data(lambda x: 2 * x + 20, np.linspace(-20, 20, 201), 1, 0.5)
-# kernel = Linear(1)
-kernel = Linear(1) + C()
+# X, Y = generate_data(lambda x: 2 * x + 20, np.linspace(-20, 20, 201), 1, 0.5)
+# # kernel = Linear(1)
+# # kernel = Linear(1) + C()
 # kernel = LinearWithOffset(1)
 
 
-m = doGPR(X, Y, kernel, 5)
+X, Y = generate_data(lambda x: (x - 2) * (x  + 3), np.linspace(-20, 20, 201), 1, 1)
+# kernel = Linear(1) * Linear(1)
+# kernel = Linear(1) * Linear(1) + C()
+# kernel = ProductKE(['LIN', 'LIN']).to_kernel()
+kernel = SumKE(['C'], [ProductKE(['LIN', 'LIN'])]).to_kernel()
+# kernel = LinearWithOffset(1) * LinearWithOffset(1)
+# kernel = LinearWithOffset(1) * LinearWithOffset(1) + C()
+
+
+m = doGPR(X, Y, kernel, 5, BIC)
+
+    # VERDICT: LinearWithOffset is just better: better fits (even accounting for extra params) which is interpretable (the offsets are just roots)

Tests/regressPeriodicKernel.py

@@ -1,14 +1,18 @@
+import numpy as np
+
 from GPy_ABCD.Util.kernelUtil import doGPR
+from GPy_ABCD.Util.dataAndPlottingUtil import generate_data
 from GPy_ABCD.Kernels.baseKernels import *
-import numpy as np
 
 # np.seterr(all='raise') # Raise exceptions instead of RuntimeWarnings. The exceptions can then be caught by the debugger
 
 
-X = np.linspace(-10, 10, 101)[:, None]
-
-Y = np.cos( (X - 5) / 2 )**2 * 7 + np.random.randn(101, 1) * 1 #- 100
+# X, Y = generate_data(lambda x: np.cos( (x - 5) / 2 )**2, np.linspace(-10, 10, 101), 7, 1)
+X, Y = generate_data(lambda x: np.cos( (x - 5) / 2 )**2 + 10, np.linspace(-10, 10, 101), 7, 1)
 
-doGPR(X, Y, PER + C, 10)
 
+# doGPR(X, Y, PER(), 10)
+doGPR(X, Y, PER() + C(), 10)
 
+# doGPR(X, Y, SE(), 10)
+# doGPR(X, Y, SE() + C(), 10)

Tests/test_kernelExpressionOperations.py

@@ -82,6 +82,17 @@ def test_both_changes(kex, other, res):
 def test_remove_a_term(kex, res): assert remove_a_term(kex) == res
 
 
+@pytest.mark.parametrize('kex, res', [
+    (SumKE(['PER', 'SE']), [SumKE(['PER'], [ProductKE(['LIN', 'LIN'])]), SumKE(['PER'], [ProductKE(['LIN', 'LIN', 'LIN'])])]),
+    (ProductKE(['PER', 'SE']), [ProductKE(['PER', 'LIN', 'LIN']), ProductKE(['PER', 'LIN', 'LIN', 'LIN']), ProductKE(['PER', 'PER'])]),
+    (ProductKE(['SE', 'SE']), [ProductKE(['LIN', 'LIN']), ProductKE(['LIN', 'LIN', 'LIN']), ProductKE(['PER', 'PER'])]),
+    (ChangeKE('CP', 'C', ProductKE(['LIN', 'SE'])), [ChangeKE('CP', ProductKE(['LIN', 'LIN']), ProductKE(['LIN', 'SE'])), ChangeKE('CP', ProductKE(['LIN', 'LIN', 'LIN']), ProductKE(['LIN', 'SE']))]),
+    (ChangeKE('CW', ProductKE(['LIN', 'SE']), 'PER'), [ChangeKE('CW', ProductKE(['LIN', 'SE']), ProductKE(['LIN', 'LIN'])), ChangeKE('CW', ProductKE(['LIN', 'SE']), ProductKE(['LIN', 'LIN', 'LIN']))]),
+    ('PER', [[]]) # This does not occur in expansions
+])
+def test_higher_curves(kex, res): assert higher_curves(kex) == res
+
+
 
 
 @pytest.mark.by_inspection