# Weighted dCGP for a symbolic regression task



In [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



In [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)$



In [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



In [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



In [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)

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

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

# update the weights
idx = (aw[awidx[k]][0] - n) * a + aw[awidx[k]][1]
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$



In [6]:

x = np.linspace(1,3,10)
x = gdual(x)
yt = data_P1(x)




In [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]




In [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$



In [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)




In [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))]




In [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}$



In [12]:

x = np.linspace(-0.9,1,10)
x = gdual(x)
yt = data_P3(x)




In [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)]




In [ ]:

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$



In [14]:

x = np.linspace(-1,1,10)
x = gdual(x)
yt = data_P4(x)




In [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]




In [ ]:

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$



In [16]:

x = np.linspace(1,3,10)
x = gdual(x)
yt = data_P5(x)




In [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]




In [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)}$



In [19]:

x = np.linspace(-2.1,1,10)
x = gdual(x)
yt = data_P6(x)




In [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




In [ ]:

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)$



In [21]:

x = np.linspace(-1,1,10)
x = gdual(x)
yt = data_P7(x)




In [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)]




In [ ]:

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)