fix structure based constants

examples/symbolic_regression/structural_constants.py

@@ -23,7 +23,7 @@ def measure(ind):
     yhat = f(points)
     if np.isscalar(yhat):
         yhat = np.ones_like(y) * yhat
-    return nrmse(y, yhat), len(ind)
+    return nrmse(y, yhat), len(gp.individual.resolve_sc(ind))
 
 
 def update_fitness(population, map=map):
@@ -35,7 +35,7 @@ def update_fitness(population, map=map):
 
 
 def main():
-    pop_size = 100
+    pop_size = 400
 
     mate = deap.gp.cxOnePoint
     expr_mut = partial(deap.gp.genFull, min_=0, max_=2)
@@ -45,7 +45,7 @@ def main():
 
     pop = update_fitness(Individual.create_population(pop_size))
 
-    for gen in range(20):
+    for gen in range(200):
         pop = algorithm.evolve(pop)
         pop = update_fitness(pop)
         best = deap.tools.selBest(pop, 1)[0]

glyph/cli/glyph_remote.py

@@ -527,8 +527,8 @@ def send_meta_data(app):
 
 
 def main():
-    app, bc, args = make_remote_app()
     logger.info(f"Glyph-remote: Version {version}")
+    app, bc, args = make_remote_app()
     app.run(break_condition=bc)
 
 

glyph/gp/individual.py

@@ -37,25 +37,45 @@ def sc_mmqout(x, y, cmin=-1, cmax=1):
 
 
 class StructConst(deap.gp.Primitive):
-    def __init__(self, func):
+    def __init__(self, func, arity=2):
         """
         :param func: evaluate left and right subtree and assign a constant.
         """
         self.func = func
-        super().__init__("SC", [deap.gp.__type__]*2, deap.gp.__type__)
+        super().__init__("SC", [deap.gp.__type__]*arity, deap.gp.__type__)
 
     @staticmethod
     def get_len(expr, tokens=("(,")):
         regex = "|".join("\\{}".format(t) for t in tokens)
         return len(re.split(regex, expr))
 
-    def format(self, *args):
-        left, right = args
-        return str(self.func(self.get_len(left), self.get_len(right)))
+    def as_terminal(self, *args):
+        vals = list(map(len, args))
+        return deap.gp.Terminal(self.func(*vals), False, deap.gp.__type__)
+
+
+def resolve_sc(ind):
+    """Evaluate StructConst in individual top to bottom."""
+    nodes = type(ind)(ind[:])
+
+    # structure based constants need to be evaluated top to bottom first
+    is_sc = lambda n: isinstance(n, StructConst)
+    while any(n for n in nodes if is_sc(n)):
+        # get the first structure based constant
+        i, node = next(filter((lambda x: is_sc(x[1])), enumerate(nodes)))
+        slice_ = nodes.searchSubtree(i)
+        # find its subtree
+        subtree = type(ind)(nodes[slice_])
+        # replace the subtree by a Terminal representing its numeric value
+        # child_trees helper function yields a list of tree which are used as arguments by the first primitive
+        args = list(child_trees(subtree))
+        nodes[slice_] = [node.as_terminal(*args)]
+    return nodes
 
 
 def add_sc(pset, func):
     """Adds a structural constant to a given primitive set.
+
     :param func: `callable(x, y) -> float` where x and y are the expressions of the left and right subtree
     :param pset: You may want to use `sympy_primitive_set` or `numpy_primitive_set` without symbolic constants.
     :type pset: `deap.gp.PrimitiveSet`
@@ -123,7 +143,7 @@ def sympy_phenotype(individual):
     """
     pset = individual.pset
     args = _build_args_string(pset, pset.constants)
-    expr = sympy.sympify(deap.gp.compile(repr(individual), pset))
+    expr = sympy.sympify(deap.gp.compile(str(individual), pset))
     func = sympy.utilities.lambdify(args, expr, modules='numpy')
     return func
 
@@ -215,7 +235,7 @@ def numpy_phenotype(individual):
         index = []
     consts = ["c_{}".format(i) for i in range(len(index))]
     args = _build_args_string(pset, consts)
-    expr = repr(individual)
+    expr = str(individual)
     for c_ in consts:
         expr = expr.replace(c, c_, 1)
     func = sympy.utilities.lambdify(args, expr, modules=pset.context)
@@ -256,14 +276,21 @@ def __init__(self, content):
         self.fitness = Measure()
 
     def __repr__(self):
+        return self.to_polish(replace_struct=False)
+
+    def __str__(self):
+        return self.to_polish()
+
+    def to_polish(self, for_sympy=False, replace_struct=True):
         """Symbolic representation of the expression tree."""
-        repr = ''
+        nodes = resolve_sc(self) if replace_struct else self
+        repr = ""
         stack = []
-        for node in self:
+        for node in nodes:
             stack.append((node, []))
             while len(stack[-1][1]) == stack[-1][0].arity:
                 prim, args = stack.pop()
-                repr = prim.format(*args)
+                repr = prim.format(*args) if not for_sympy else convert_inverse_prim(prim, args)
                 if len(stack) == 0:
                     break
                 stack[-1][1].append(repr)
@@ -432,22 +459,6 @@ def convert_inverse_prim(prim, args):
     return prim_formatter(*args)
 
 
-def stringify_for_sympy(f):
-    """Return the expression in a human readable string.
-    """
-    string = ""
-    stack = []
-    for node in f:
-        stack.append((node, []))
-        while len(stack[-1][1]) == stack[-1][0].arity:
-            prim, args = stack.pop()
-            string = convert_inverse_prim(prim, args)
-            if len(stack) == 0:
-                break  # If stack is empty, all nodes should have been seen
-            stack[-1][1].append(string)
-    return string
-
-
 @glyph.utils.Memoize
 def simplify_this(expr, timeout=5):
     """
@@ -460,14 +471,15 @@ def simplify_this(expr, timeout=5):
     with glyph.utils.random_state(simplify_this):  # to avoid strange side effects of sympy testcases and random
         with glyph.utils.Timeout(timeout):
             try:
-                return sympy.simplify(stringify_for_sympy(expr))
+                return sympy.simplify(expr.to_polish(for_sympy=True))
             except Exception as e:
                 logger.debug(f"Exception during simplification of {expr}: {e}.")
                 return expr
         return expr
 
 
 def child_trees(ind):
+    """Yield all child tree which are used as arguments for the head node of ind."""
     start = 1
     while start < len(ind):
         slice_ = ind.searchSubtree(start)

tests/unittest/gp/test_individual.py

@@ -154,40 +154,33 @@ def test_get_len(case):
     assert StructConst.get_len(expr) == x
 
 
-def test_struct_const_format():
-    f = lambda x, y: x + y
-    sc = StructConst(f)
-
-    left = "Add(x0, x0)"
-    right = "y0"
-
-    assert sc.format(left, right) == str(4)
-
-
 @pytest.fixture
 def sc_ind():
     pset = numpy_primitive_set(1)
-    f = lambda x, y: 0
+    f = lambda x, y: x/y
     pset = add_sc(pset, f)
     pset.addPrimitive(operator.neg, 1, "Neg")
 
     return Individual(pset=pset)
 
 
-@pytest.mark.parametrize("case", get_len_case)
-def test_struct_const_repr(case, sc_ind):
-    subexpr, _ = case
-    expr = "Add(x_0, SC({subexpr}, {subexpr}))".format(subexpr=subexpr)
-    ind = sc_ind.from_string(expr)
+sc_cases = (
+    ("SC(x_0, x_0)", 1),
+    ("SC(x_0, SC(x_0, x_0))", 1/3),
+)
 
+
+@pytest.mark.parametrize("expr,res", sc_cases)
+def test_struc_const(expr, res, sc_ind):
+    ind = sc_ind.from_string(expr)
     func = numpy_phenotype(ind)
-    assert func(1) == 1
+    assert func(0) == res
 
 
 def test_simplify_this_struct_const(sc_ind):
     expr = "SC(0.0, Add(x_0, Neg(x_0)))"
     ind = sc_ind.from_string(expr)
-    assert str(ind) == str(simplify_this(ind))
+    assert ind.to_polish() in str(simplify_this(ind))  # rounding
 
 
 child_tree_cases = (