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)