Weighted dCGP for a symbolic regression task

[1]:
from pyaudi import gdual_vdouble as gdual
from dcgpy import expression_weighted_gdual_vdouble as expression
from dcgpy import kernel_set_gdual_vdouble as kernel_set
import pyaudi
import numpy as np
import math
import re
%matplotlib inline

The ES-(1+\(\lambda\)) algorithm

[2]:
def run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output):
    # The offsprings chromosome, fitness and weights
    chromosome = [1] * offsprings
    fitness = [1] * offsprings
    weights = [1] * offsprings
    # Init the best as the initial random dCGP
    best_chromosome = dCGP.get()
    best_weights = dCGP.get_weights()
    best_fitness = sum(mse(dCGP, x, yt).constant_cf)
    # Main loop over generations
    for g in range(max_gen):
        for i in range(offsprings):
            dCGP.set(best_chromosome)
            dCGP.set_weights(best_weights)
            cumsum=0
            dCGP.mutate_active(i)
            newton(dCGP, mse, x, yt, newtonParams)
            fitness[i] = sum(mse(dCGP, x, yt).constant_cf)
            chromosome[i] = dCGP.get()
            weights[i] = dCGP.get_weights()
        for i in range(offsprings):
            if fitness[i] <= best_fitness:
                dCGP.set(chromosome[i])
                dCGP.set_weights(weights[i])
                if (fitness[i] != best_fitness) and screen_output:
                    print("New best found: gen: ", g, " value: ", fitness[i], dCGP.simplify(["x"],True))
                best_chromosome = chromosome[i]
                best_fitness = fitness[i]
                best_weights = weights[i]

        if best_fitness < 1e-14:
            break
    return g, best_chromosome, best_weights

The test problems

P1: \(x^5 - \pi x^3 + x\)

P2: \(x^5 - \pi x^3 + \frac{2\pi}x\)

P3: \(\frac{e x^5 + x^3}{x + 1}\)

P4: \(\sin(\pi x) + \frac 1x\)

P5: \(e x^5 - \pi x^3 + x\)

P6: \(\frac{e x^2-1}{\pi (x + 2)}\)

P7: \(\sin(e x) + \cos(\pi x)\)

[3]:
# The following functions create the target values for a gridded input x for different test problems
def data_P1(x):
    return x**5 - np.pi*x**3 + x
def data_P2(x):
    return x**5 - np.pi*x**3 + 2*np.pi / x
def data_P3(x):
    return (np.e*x**5 + x**3)/(x + 1)
def data_P4(x):
    return pyaudi.sin(np.pi * x) + 1./x
def data_P5(x):
    return np.e * x**5 - np.pi*x**3 + np.sqrt(2) * x
def data_P5(x):
    return np.e * x**5 - np.pi*x**3 + x
def data_P6(x):
    return (np.e*x**2-1) / (np.pi*(x + 2))
def data_P7(x):
    return pyaudi.sin(np.e*x)+pyaudi.cos(np.pi*x)

The error function

[4]:
# This is used to sum over the component of a vectorized coefficient, accounting for the fact that if its dimension
# is 1, then it could represent [a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a ...] with [a]
def collapse_vectorized_coefficient(x, N):
    if len(x) == N:
        return sum(x)
    return x[0] * N

# Quadratic error of a dCGP expression. The error is computed over the input points xin (of type gdual, order 0 as
# we are not interested in expanding the program w.r.t. these). The target values are contained in yt (of type gdual,
# order 0 as we are not interested in expanding the program w.r.t. these)
def mse(dCGP, xin, yt):
    y = dCGP([xin])[0]
    return (y-yt)**2

Newton’s method

[5]:
# Newton's method for minimizing the error function f w.r.t. the weights of the dCGP expression.
# We take a specified amount of steps, each by choosing randomly 2 or 3 weights
def newton(ex, f, x, yt, p):
    n = ex.get_n()
    r = ex.get_rows()
    c = ex.get_cols()
    a = ex.get_arity()[0]
    v = np.zeros(r * c * a)

    # random initialization of weights
    w=[]
    for i in range(r*c):
        for j in range(a):
            w.append(gdual([np.random.normal(0,1)]))
    ex.set_weights(w)
    wi = ex.get_weights()

    # get active weights
    an = ex.get_active_nodes()
    is_active = [False] * (n + r * c) # bool vector of active nodes
    for k in range(len(an)):
        is_active[an[k]] = True
    aw=[] # list of active weights
    for k in range(len(an)):
        if an[k] >= n:
            for l in range(a):
                aw.append([an[k], l]) # pair node/ingoing connection
    if len(aw)<2:
        return

    for i in range(p['steps']):
        w = ex.get_weights() # initial weights

        # random choice of the weights w.r.t. which we'll minimize the error
        num_vars = np.random.randint(2, min(3, len(aw)) + 1) # number of weights (2 or 3)
        awidx = np.random.choice(len(aw), num_vars, replace = False) # indexes of chosen weights
        ss = [] # symbols
        for j in range(len(awidx)):
            ss.append("w" + str(aw[awidx[j]][0]) + "_" + str(aw[awidx[j]][1]))
            idx = (aw[awidx[j]][0] - n) * a + aw[awidx[j]][1]
            w[idx] = gdual(w[idx].constant_cf, ss[j], 2)
        ex.set_weights(w)

        # compute the error
        E = f(ex, x, yt)
        Ei = sum(E.constant_cf)

        # get gradient and Hessian
        dw = np.zeros(len(ss))
        H = np.zeros((len(ss),len(ss)))
        for k in range(len(ss)):
            dw[k] = collapse_vectorized_coefficient(E.get_derivative({"d"+ss[k]: 1}), len(x.constant_cf))
            H[k][k] = collapse_vectorized_coefficient(E.get_derivative({"d"+ss[k]: 2}), len(x.constant_cf))
            for l in range(k):
                H[k][l] = collapse_vectorized_coefficient(E.get_derivative({"d"+ss[k]: 1, "d"+ss[l]: 1}), len(x.constant_cf))
                H[l][k] = H[k][l]

        det = np.linalg.det(H)
        if det == 0: # if H is singular
            continue

        # compute the updates
        updates = - np.linalg.inv(H) @ dw

        # update the weights
        for k in range(len(updates)):
            idx = (aw[awidx[k]][0] - n) * a + aw[awidx[k]][1]
            ex.set_weight(aw[awidx[k]][0], aw[awidx[k]][1], w[idx] + updates[k])
        wfe = ex.get_weights()
        for j in range(len(awidx)):
            idx = (aw[awidx[j]][0] - n) * a + aw[awidx[j]][1]
            wfe[idx] = gdual(wfe[idx].constant_cf)
        ex.set_weights(wfe)

        # if error increased restore the initial weights
        Ef = sum(f(ex, x, yt).constant_cf)
        if not Ef < Ei:
            for j in range(len(awidx)):
                idx = (aw[awidx[j]][0] - n) * a + aw[awidx[j]][1]
                w[idx] = gdual(w[idx].constant_cf)
            ex.set_weights(w)

Problem P1: \(x^5 - \pi x^3 + x\)

[6]:
x = np.linspace(1,3,10)
x = gdual(x)
yt = data_P1(x)
[7]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(inputs = 1, outputs = 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))

res = np.array(res)
restart:         gen:    expression:
1                109     [0.999999999999997*x**5 - 3.14159265358976*x**3 + 0.999999999999932*x]
4                43      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999996*x]
7                29      [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]
8                97      [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
10               36      [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]
12               40      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999989*x]
13               79      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999998*x]
15               92      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999998*x]
16               133     [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]
18               97      [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
19               46      [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999991*x]
24               151     [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999987*x]
25               33      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999998*x]
30               131     [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]
32               23      [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]
33               81      [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
34               16      [1.0*x**5 - 3.14159265358979*x**3 + 1.00000000000001*x]
36               48      [0.999999999999999*x**5 - 3.14159265358979*x**3 + 1.0*x]
39               150     [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]
45               50      [0.999999999999995*x**5 - 3.14159265358974*x**3 + 0.999999999999874*x]
46               189     [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]
49               67      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]
55               29      [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
56               48      [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]
60               51      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999988*x]
71               90      [0.999999999999999*x**5 - 3.14159265358979*x**3 + 1004.92559245124*x**3/(1004.92559245124*x**2 - 3157.0668586492) - 3157.06685864919*x/(1004.92559245124*x**2 - 3157.0668586492)]
74               43      [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]
77               142     [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]
80               192     [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999988*x]
82               62      [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999993*x]
85               52      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]
86               32      [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999992*x]
87               83      [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999992*x]
89               149     [0.999999999999998*x**5 - 3.14159265358978*x**3 + 0.999999999999981*x]
91               141     [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999999*x]
95               96      [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
[8]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)
ERT Expected run time = avg. number of dCGP evaluations needed:  174288.8888888889

Problem P2: \(x^5 - \pi x^3 + \frac{2\pi}x\)

[6]:
x = np.linspace(0.1,5,10) # we include points close to zero here to favour learning of 1/x
x = gdual(x)
yt = data_P2(x)
[7]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))
res = np.array(res)
restart:         gen:    expression:
1                64      [1.00000000001522*x**5 - 1.27654166639185e-10*x**4 - 3.14159265332916*x**3 + 6.28318530717146/x]
3                160     [0.999999999999802*x**5 - 3.14159265358324*x**3 - 4.58664018735164e-11*x + 6.28318530718224/x]
6                43      [1.00000000000018*x**5 - 3.14159265361377*x**3 + 6.28318530717944/x]
7                35      [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
8                63      [1.0*x**5 - 3.14159265358978*x**3 + 6.28318530717959/x]
11               168     [1.0*x**5 - 3.14159265358979*x**3 - 3.01980662698043e-14*x + 6.28318530717959/x]
13               24      [0.999999999999994*x**5 - 3.14159265359058*x**3 + 2.78627727461556e-13*x + 6.28318530717957/x]
14               143     [0.99999999997825*x**5 - 3.14159265281984*x**3 - 6.06404304548391e-9*x + 6.28318530741037/x]
16               187     [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
21               82      [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
27               7       [1.0*x**5 - 3.14159265358984*x**3 + 1.84329364468004e-13*x**2 + 6.28318530717957/x]
28               67      [0.999999999996586*x**5 - 3.14159265349105*x**3 - 2.21245651961797e-9 + 6.28318530745686/x]
34               153     [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
39               115     [0.999999999999999*x**5 - 3.14159265358978*x**3 + 6.28318530717959/x]
43               84      [0.999999999999923*x**5 + 6.42597086653041e-13*x**4 - 3.14159265359109*x**3 + 6.28318530717962/x]
50               149     [0.99999999995493*x**5 - 3.14159265203297*x**3 - 1.17409352674938e-8*x + 6.28318530777391/x]
51               149     [1.00000000012198*x**5 - 3.14159265780051*x**3 + 3.16978514547372e-8*x + 6.28318530601517/x]
52               64      [1.0*x**5 - 3.1415926535898*x**3 + 6.28318530717959/x]
54               79      [-6.32238025396602e-16*x**7 + 1.00000000000002*x**5 - 3.14159265358992*x**3 + 6.28318530717958/x]
55               69      [0.999999999999999*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
57               51      [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717958/x]
58               180     [0.999999999997159*x**5 + 2.56552832733729e-11*x**4 - 3.14159265364781*x**3 + 6.28318530750492/x]
62               96      [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
64               130     [1.0*x**5 - 3.1415926535898*x**3 + 6.28318530717959/x]
67               154     [0.99999999999999*x**5 + 7.88553657836724e-14*x**4 - 3.14159265358995*x**3 - 6.63622336613371e-14*x**2 + 6.2831853071796/x]
70               26      [1.00000000000018*x**5 - 3.1415926535975*x**3 + 1.64577260147692e-11*x**2 + 6.28318530721673/x]
75               48      [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717961/x]
76               138     [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
78               82      [1.0*x**5 - 3.14159265358982*x**3 + 6.28318530719116/x]
84               165     [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]
89               136     [1.00000000000278*x**5 - 3.14159265388607*x**3 + 1.09744339194339e-9*x**2 + 6.28318530643488/x]
90               188     [1.00000000000001*x**5 - 3.91763911922868e-14*x**4 - 3.14159265358973*x**3 + 6.28318530717958/x]
91               108     [1.0*x**5 - 3.14159265358979*x**3 + 6.2831853071796/x]
99               64      [1.0*x**5 - 3.14159265358979*x**3 - 0.36315536377679/(x*(2.18996409662261e-16*x**3 - 0.0577979712554117))]
[8]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)
ERT Expected run time = avg. number of dCGP evaluations needed:  195352.94117647057

Problem P3: \(\frac{e x^5 + x^3}{x + 1}\)

[12]:
x = np.linspace(-0.9,1,10)
x = gdual(x)
yt = data_P3(x)
[13]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))
res = np.array(res)
restart:         gen:    expression:
8                186     [0.899166850493274*x**5/(0.330784998479351*x + 0.330784998479351) + 0.330784998479352*x**3/(0.330784998479351*x + 0.330784998479351)]
14               161     [1.20883912632723*x**5/(0.444707062259436*x + 0.444707062259436) + 0.444707062259438*x**3/(0.444707062259436*x + 0.444707062259436)]
73               125     [2.71828182845905*x**5/(x + 1) + 1.0*x**3/(x + 1)]
[ ]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)

Problem P4: \(\sin(\pi x) + \frac 1x\)

[14]:
x = np.linspace(-1,1,10)
x = gdual(x)
yt = data_P4(x)
[15]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div", "sin"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))
res = np.array(res)
restart:         gen:    expression:
/home/dario/.local/lib/python3.7/site-packages/numpy/linalg/linalg.py:2116: RuntimeWarning: invalid value encountered in det
  r = _umath_linalg.det(a, signature=signature)
2                1       [0.999999999999999*sin(3.14159265358979*x) + 1.0/x]
6                132     [1.0*sin(3.14159265358979*x) - 3.6500836381227e-16 + 1.0/x]
7                89      [0.999999999999999*sin(3.14159265358979*x) + 1.0/x]
8                133     [0.999999999999997*sin(3.14159265358979*x) + 1.0/x]
10               120     [0.999999999999999*sin(3.14159265358979*x) + 4.68530917602041e-17 + 1.0/x]
12               55      [-7.89209275531109e-11*x + 1.0000000000362*sin(3.14159265367275*x + 37.6991118429053) - 9.47051130607999e-10 + 1.00000000000142/x]
14               66      [1.00000000000032*sin(3.14159265358985*x + 2.53915412165946e-16) + 0.999999999999973/x]
16               12      [1.0*sin(3.14159265358979*x) + 1.0/x]
19               2       [0.999999949216629*sin(3.14159264435007*x) + 9.36443001715427e-18 + 1.00000000433431/x]
20               123     [1.0*sin(3.14159265358979*x) + (7.86221303890169e-17*sin(3.14159265358979*x) + 1.0)/x]
21               17      [9.62453199799219*x*sin(3.1415927053321*x)/(9.62453191537768*x - 1.06372808316533e-16) + 9.62453191084485/(9.62453191537768*x - 1.06372808316533e-16)]
22               13      [0.999999999998995*sin(3.14159265358961*x) + 1.00000000000009/x]
24               69      [0.999999999999953*sin(x*(9.25302298900841e-17*x + 3.14159265358979)) + 1.0/x]
26               105     [1.0*sin(x*(1.77990517211438e-16*x + 3.14159265358979)) + 5.6656141275367e-17 + 1.0/x]
27               115     [-1.24427033889047e-16*x*sin(3.14159265358979*x) + 1.0*sin(3.14159265358979*x) + 1.0/x]
/home/dario/.local/lib/python3.7/site-packages/numpy/linalg/linalg.py:2116: RuntimeWarning: overflow encountered in det
  r = _umath_linalg.det(a, signature=signature)
34               64      [0.999999999999999*sin(3.14159265358979*x) - 1.29978496930744e-16 + 1.0/x]
36               46      [1.76206042739191e-16*x*sin(3.14159265358979*x)/(x*sin(3.14159265358979*x) + 1) + 1.0*sin(3.14159265358979*x) + 1.76206042739191e-16/(x*sin(3.14159265358979*x) + 1) + 1.0/x]
38               95      [1.0*sin(x*(1.3303626169174e-16*x + 3.14159265358979)) + 1.69425620306735e-16 + 1.0/x]
40               143     [8.12378036198924e-33*x + 1.0*sin(3.14159265358981*x) - 1.80264032596514e-16 + 1.0/x]
42               16      [1.00000000004875*sin(3.14159265359866*x) + 0.542336605776001*sin(2.35272399768743e-12*x**2) + 2541.81495243938/(2541.81495245225*x - 2.15281241329688) - 2.15281241328598/(x*(2541.81495245225*x - 2.15281241329688)) + 972.092685040695*sin(2.35272399768743e-12*x**2)/(x**2*(2541.81495245225*x - 2.15281241329688))]
44               24      [1.0*sin(3.14159265358979*x) + 1.0/x]
45               162     [6.02619053589549e-16*sin(3.14159265355456*x)**2 + 0.999999999994157*sin(3.14159265355456*x) - 4.51719913513711e-16 + 1.00000000000032/x]
47               42      [-1.87473626494473e-29*x**2 - 1.8747362649445e-29*x/sin(3.14159265359048*x) + 1.00000000000011*sin(3.14159265359048*x) + 0.999999999999994/x]
48               28      [1.00000000000002*sin(3.14159265358992*x) + 0.999999999999999/x]
49               18      [1.0*sin(3.14159265358979*x) + 1.0/x]
50               41      [1.0*sin(3.14159265358979*x) - 4.08683517342902e-17 + 1.0/x]
52               61      [1.0*sin(3.14159265358979*x) + 1.0/x]
53               49      [1.0*sin(3.14159265358979*x) + 1.0/x]
55               29      [1.0*sin(3.14159265358979*x) + (1.0 - 6.66045855255763e-18*sin(3.14159265358979*x))/x]
56               15      [6.58379173825623*x + 5.64744625319335*sin(6.11418858333298*x) + 5.2879457406247*sin(16.7761307129903*x)]
57               32      [-0.974501474950861*sin(0.863186354404554*x) + 9.07750333059453*sin(1.6599296272374*x + 2.06966097571612*sin(4.75780034524026*x)) - 5.38968705261534*sin(3.48198655314495*x - 1.09415284338664*sin(4.75780034524026*x))]
59               76      [1.0*sin(3.14159265358979*x) + 4.45769125393131e-16 + 1.0/x]
60               7       [1.00000000000001*sin(3.1415926535899*x) + 4.70822994414938e-17 + 1.0/x]
63               76      [5.19135160727313*sin(3.59576183848133*x) + 8.66296848504959*sin(12.5211187903134*x - 0.987545240550494*sin(3.59576183848133*x))]
68               150     [1.0*sin(3.14159265358979*x) + 1.0/x]
72               58      [1.0*sin(3.14159265358979*x) + 1.25427228495141e-16 + 1.0/x]
73               89      [1.0*sin(3.14159265358979*x) + 1.0/x]
75               49      [1.0*sin(3.14159265358979*x) + 1.0/x]
79               76      [1.0*sin(3.14159265358979*x) + 3.94677549026297e-17 + 1.0/x]
83               76      [1.0*sin(3.14159265358979*x) + 1.0/x]
85               12      [1.0*sin(3.14159265358979*x) + 1.0/x]
87               94      [1.000000000001*sin(3.14159265359084*x) + (0.999999999999751 - 1.35525271560688e-20/sin(3.14159265359084*x))/x]
92               19      [0.999999980945392*sin(3.14159265012292*x) + 1.00000000162629/x]
93               153     [-0.999999999999999*sin(3.14159265358979*x + 3.14159265358979) + 1.0/x]
97               48      [1.0*sin(3.14159265358979*x) + 1.0/x]
[ ]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)

Problem P5: \(ex^5 - \pi x^3 + x\)

[16]:
x = np.linspace(1,3,10)
x = gdual(x)
yt = data_P5(x)
[17]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))
res = np.array(res)
restart:         gen:    expression:
2                130     [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000003*x]
5                144     [2.71828182845904*x**5 - 3.14159265358974*x**3 + 0.999999999999915*x]
6                129     [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
7                48      [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999982*x]
9                25      [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000002*x]
10               6       [2.71828182845904*x**5 - 3.14159265358976*x**3 + 0.999999999999928*x]
11               22      [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999977*x]
13               113     [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999985*x]
15               9       [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999969*x]
16               46      [2.71828182845796*x**5 - 3.1415926535816*x**3 + 0.999999999987821*x]
17               77      [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
18               189     [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
20               181     [2.71828182765082*x**5 + 4.81282166780353e-9*x**4 - 3.14159266134502*x**3 + 1.0000000050856*x]
21               163     [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000003*x]
22               16      [2.71828182844251*x**5 - 3.14159265338736*x**3 + 0.999999999457422*x]
24               37      [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]
27               119     [2.71828182845942*x**5 - 3.14159265359449*x**3 + 1.00000000000628*x]
28               4       [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
29               42      [2.71828181857779*x**5 - 3.14159237432049*x**3 - 7.49766685257696e-7*x**2 + 1.00000053196694*x]
42               33      [2.71828182845902*x**5 - 3.14159265358958*x**3 + 0.999999999999598*x]
44               8       [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]
45               189     [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000005*x]
49               5       [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]
53               21      [2.71828182845904*x**5 - 3.14159265358977*x**3 + 0.99999999999993*x]
55               27      [2.71828182845904*x**5 - 3.14159265358977*x**3 + 0.999999999999773*x]
63               83      [2.71828182845904*x**5 - 3.14159265358979*x**3 + 0.99999999999998*x]
65               29      [2.71828182845904*x**5 - 3.14159265358979*x**3 + 0.999999999999992*x]
68               149     [2.71828182845905*x**5 + x**3*(1.51490803020543 - 2.20219730313442*x**2)/(0.700981172914968*x**2 - 0.458140436777564) + x*(0.625363773202351*x**2 - 0.458140436777565)/(0.700981172914968*x**2 - 0.458140436777564)]
69               66      [2.71828182845905*x**5 - 3.14159265358978*x**3 + 0.999999999999967*x]
71               65      [2.71828182845905*x**5 - 3.14159265358988*x**3 + 1.00000000000015*x]
73               181     [2.71828182845904*x**5 - 3.14159265358973*x**3 + 0.999999999999894*x]
74               150     [2.71828182845904*x**5 - 3.14159265358979*x**3 + 0.999999999999995*x]
78               12      [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999977*x]
80               118     [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]
81               137     [2.71828182845906*x**5 - 3.14159265359007*x**3 + 1.00000000000047*x]
84               50      [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]
85               33      [2.71828182845904*x**5 - 3.1415926535897*x**3 + 0.999999999999851*x]
87               93      [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000003*x]
88               60      [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.99999999999997*x]
89               88      [2.71828182845905*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]
91               15      [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000004*x]
92               162     [2.71828182845904*x**5 - 3.14159265358977*x**3 + 0.999999999999982*x]
95               61      [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000003*x]
98               68      [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000006*x]
[18]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)
ERT Expected run time = avg. number of dCGP evaluations needed:  132245.45454545456

Problem P6: \(\frac{e x^2-1}{\pi (x + 2)}\)

[19]:
x = np.linspace(-2.1,1,10)
x = gdual(x)
yt = data_P6(x)
[20]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))
res = np.array(res)
restart:         gen:    expression:
7                121     [0.658810478707329*x**2/(0.76140528857091*x + 1.52281057714182) - 0.242362830744741/(0.76140528857091*x + 1.52281057714182)]
12               70      [x*(0.710716298622762*x - 6.5421001949062e-12)/(0.82139426426137*x + 1.64278852852271) - 0.261457914778349/(0.82139426426137*x + 1.64278852852271)]
18               184     [x*(0.000123063420861819*x**3 + 0.000387617697492767*x**2 + 0.00026373499528672*x - 3.84934325911654e-5)/((0.0332304392124186*x + 0.0664608784094327)*(0.00428004451413307*x**2 + 0.0134810245653551*x + 0.00984187107586485)) - 0.0235684911144266/(0.128798914957873*x + 0.148085179001607)]
37               56      [-0.56551064554709*x/(1.74078576306486*x + 3.48157152619754) - 1.40565254304111 + (1.50622528661617*x**2 + 3.01245057329103*x + 4.33977055108826)/(1.74078576306486*x + 3.48157152619754)]
42               154     [0.0804769112118696*x**2/(0.093009367314253*x + 0.186018734628506) - 0.0296058011238267/(0.093009367314253*x + 0.186018734628506)]
/home/dario/miniconda3/envs/dcgpy/lib/python3.7/site-packages/ipykernel_launcher.py:63: RuntimeWarning: invalid value encountered in matmul
[ ]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)

Problem P7: \(\sin(e x) + \cos(\pi x)\)

[21]:
x = np.linspace(-1,1,10)
x = gdual(x)
yt = data_P7(x)
[22]:
# We run nexp experiments and accumulate statistic for the ERT
nexp = 100
offsprings = 4
max_gen=200
res = []
kernels = kernel_set(["sum", "mul", "diff", "div", "sin", "cos"])()
newtonParams = {
    'steps': 100,
}
print("restart: \t gen: \t expression:")
for i in range(nexp):
    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))
    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):
        for k in range(dCGP.get_arity()[0]):
            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))
    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output=False)
    res.append(g)
    dCGP.set(best_chromosome)
    if g < (max_gen-1):
        print(i, "\t\t", res[i], "\t", dCGP.simplify(["x"],True))
res = np.array(res)
restart:         gen:    expression:
2                10      [1.0*sin(2.71828182845905*x) + 0.999999999999999*cos(3.14159265358979*x)]
5                177     [0.999999999990239*sin(2.71828182845446*x) + 0.99999999827576*cos(3.1415926695047*x)]
/home/dario/miniconda3/envs/dcgpy/lib/python3.7/site-packages/ipykernel_launcher.py:63: RuntimeWarning: invalid value encountered in matmul
10               52      [0.999999996606054*sin(2.71828182686597*x) + 1.0*cos(3.14159265358979*x)]
12               166     [1.00000021604097*cos(3.14159255591115*x) + 1.00000000001092*cos(2.71828182846416*x - 7.85398185936523)]
14               84      [1.00000000002564*sin(2.71828182859488*x) + 1.0*cos(3.14159265358979*x)]
15               137     [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
17               25      [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
19               145     [-0.999999999999992*sin(x*(2.59471453823685e-12*sin(2.29774027604273*x**2) - 2.718281828461)) + 1.0*cos(3.14159265358979*x)]
21               51      [1.0*sin(2.71828182845904*x) + 1.00000000000002*cos(3.14159265358978*x)]
26               8       [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
29               116     [-3.5489631288284e-16*sin(2.71828182845905*x)**2 + 1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
31               82      [0.999999999612255*sin(2.71828182743606*x + 8.09654622736853e-11*cos(3.14159265362908*x)) + 1.00000000000211*cos(3.14159265362908*x)]
32               132     [1.00000000003824*sin(2.71828182866167*x) + 1.0*cos(3.14159265358979*x)]
44               188     [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
48               122     [1.00000000001473*sin(25.5560520538422*x) + 0.999999991979214*cos(3.14159272762242*x)]
49               88      [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
50               55      [0.999999999999999*sin(2.71828182845904*x - 8.91009907439534e-10*cos(31.4159265256615*x)) - 1.00000002272797*cos(31.4159265256615*x)]
52               86      [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
54               106     [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
59               81      [1.0*sin(2.71828182845905*x) + 1.00000000226596*cos(3.14159263267476*x)]
63               110     [1.08304336385657e-7*x*sin(2.71828183307158*x)*cos(3.14159262986457*x) + 1.00000000982667*sin(2.71828183307158*x) + 0.999999952852321*cos(3.14159262986457*x)]
67               73      [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
82               53      [0.999999999999999*sin(2.71828182845904*x) - 1.0*cos(31.4159265358979*x)]
87               77      [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
93               86      [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]
98               175     [2.68964745528915e-8*x**2 + 1.0*sin(2.71828182845905*x) + 1.00000000821958*cos(3.14159270860771*x)]
[ ]:
mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']
print("ERT Expected run time = avg. number of dCGP evaluations needed: ", mean_gen * offsprings)