Solving differential equations with dCGP
Lets first import dcgpy and pyaudi and set up things as to use dCGP on gduals defined over vectorized floats
[2]:
from dcgpy import expression_gdual_vdouble as expression
from dcgpy import kernel_set_gdual_vdouble as kernel_set
from pyaudi import gdual_vdouble as gdual
from matplotlib import pyplot as plt
import numpy as np
from numpy import sin, cos
from random import randint
np.seterr(all='ignore') # avoids numpy complaining about early on malformed expressions being evalkuated
%matplotlib inline
The kernel functions
[3]:
kernels = kernel_set(["sum", "mul", "div", "diff", "log", "sin", "cos", "exp"])() # note the call operator (returns the list of kernels)
Instantiate a (1 input, 1 output) dCGP and we inspect a randomly created program
[6]:
dCGP_example = expression(inputs=1, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = 2134)
[7]:
plt.rcParams["figure.figsize"] = [10,6]
dCGP_example.visualize() #requires graphwiz module installed (pip install graphviz --user)
print("Represented expression: ", dCGP_example(["x"])[0])
print("Simplified expression: ", dCGP_example.simplify(["x"])) #requires sympy module installed
Represented expression: ((x+x)*(x*x))
Simplified expression: [2*x**3]
Define the ES that will evolve solutions
[8]:
# We run an evolutionary strategy ES(1 + offspring)
def run_experiment(max_gen, offsprings, quadratic_error, initial_conditions_error, dCGP, screen_output=False):
chromosome = [1] * offsprings
fitness = [1] *offsprings
best_chromosome = dCGP.get()
best_fitness = quadratic_error(dCGP, grid) + initial_conditions_error(dCGP)
for g in range(max_gen):
for i in range(offsprings):
dCGP.set(best_chromosome)
dCGP.mutate_active(i+1) # we mutate a number of increasingly higher active genes
qe = quadratic_error(dCGP, grid)
ie = initial_conditions_error(dCGP)
fitness[i] = ie + qe
chromosome[i] = dCGP.get()
for i in range(offsprings):
if fitness[i] <= best_fitness:
if (fitness[i] != best_fitness) and screen_output:
print("New best found: gen: ", g, " value: ", fitness[i])
best_chromosome = chromosome[i]
best_fitness = fitness[i]
dCGP.set(best_chromosome)
if best_fitness < 1e-7:
break
return g, best_chromosome
Consider the following Ordinary Differential Equation (ODE1):
\(\frac{dy}{dx} = \frac{2x - y}{x}\), with \(y(0.1) = 20.1\) and \(x \in [0.1,1]\)
we demand its punctual validity over a grid of \(N\) equally spaced points.
The solution to the ODE is \(y = x + \frac 2x\)
[9]:
# We define the quadratic error of the dCGP in the grid points
def qe_ODE1(dCGP, grid):
retval = 0
out = dCGP([grid])[0]
y = np.array(out.constant_cf)
dydx = np.array(out.get_derivative({"dx" : 1}))
x = np.array(grid.constant_cf)
ode1 = (2. * x - y) / x
retval += (ode1 - dydx) * (ode1 - dydx)
return sum(retval)
[10]:
# We define a penalty term associated to the initial conditions violation
def ic_ODE1(dCGP):
x0 = 1
y0 = 3
out = dCGP([gdual([x0])])[0]
return (out.constant_cf[0] - y0) * (out.constant_cf[0] - y0)
[11]:
# We construct the grid of points. Since the ODE only contains first order derivatives we use truncation order 1.
# Since we are using vectorized gdual we can instantiate only one gdual
values = np.linspace(0.1,1,10)
grid = gdual(values, "x", 1)
[12]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
offsprings = 10
stop = 500
res = []
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 = randint(0,100000))
g, best_chromosome = run_experiment(max_gen = stop, offsprings = offsprings, \
quadratic_error=qe_ODE1, initial_conditions_error=ic_ODE1, dCGP = dCGP)
res.append(g)
dCGP.set(best_chromosome)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x"]), " a.k.a ", dCGP.simplify(["x"]))
res = np.array(res)
restart: gen: expression:
5 414 ['(((x/(x*x))+(x/(x*x)))+x)'] a.k.a [x + 2/x]
10 285 ['((x+(x/(x*x)))+(x/(x*x)))'] a.k.a [x + 2/x]
12 456 ['(x+(((x+x)/x)/x))'] a.k.a [x + 2/x]
13 227 ['(((x/x)/x)-((((x/x)-x)-((x/x)/x))-(x/x)))'] a.k.a [x + 2/x]
19 146 ['((((x+x)/x)/x)+x)'] a.k.a [x + 2/x]
20 117 ['log(exp((x+(((x+x)/x)/x))))'] a.k.a [x + log(exp(2/x))]
22 402 ['(x+((x+x)/(x*x)))'] a.k.a [x + 2/x]
25 38 ['(((x/(x*x))+x)+(x/(x*x)))'] a.k.a [x + 2/x]
26 69 ['((((x/x)/x)+((x/x)/x))+(x/(x/x)))'] a.k.a [x + 2/x]
28 123 ['(x+((cos((x-x))+cos((x-x)))/x))'] a.k.a [x + 2/x]
29 469 ['((((x+x)/x)/x)+x)'] a.k.a [x + 2/x]
33 31 ['(x+(((x/x)+(x/x))/x))'] a.k.a [x + 2/x]
40 204 ['(x+((x+x)/(x*x)))'] a.k.a [x + 2/x]
43 446 ['(x+((x+x)/(x*x)))'] a.k.a [x + 2/x]
44 354 ['(((x/x)/(x/(x/x)))+(((x/x)/(x/(x/x)))+(x/(x/x))))'] a.k.a [x + 2/x]
45 392 ['(((exp((x-x))/x)+x)+((x/x)/x))'] a.k.a [x + 2/x]
46 164 ['((((x*x)*(x/(x*x)))+(x/(x*x)))+(x/(x*x)))'] a.k.a [x + 2/x]
50 74 ['(((x+x)/(x*x))+x)'] a.k.a [x + 2/x]
51 248 ['((((sin(x)/sin(x))/x)+x)+((sin(x)/sin(x))/x))'] a.k.a [x + 2/x]
52 129 ['log(exp((x+((x+x)/(x*x)))))'] a.k.a [x + log(exp(2/x))]
53 346 ['((x/(x*x))+(x+(x/(x*x))))'] a.k.a [x + 2/x]
57 24 ['((x/(x*x))+(x+(x/(x*x))))'] a.k.a [x + 2/x]
59 408 ['(x+(((x/x)+(x/x))/x))'] a.k.a [x + 2/x]
66 302 ['(x+(((x/x)+(x/x))/(x*(x/x))))'] a.k.a [x + 2/x]
71 370 ['((x+(x/(x*x)))+(x/(x*x)))'] a.k.a [x + 2/x]
81 382 ['(x+((((x*x)/x)+x)/(x*x)))'] a.k.a [x + 2/x]
83 258 ['(x+((((x+x)/x)/x)/(((x+x)/x)/((x+x)/x))))'] a.k.a [x + 2/x]
84 281 ['(x+(((x+x)/x)/x))'] a.k.a [x + 2/x]
85 84 ['(x+(((x*x)+(x*x))/((x*x)*x)))'] a.k.a [x + 2/x]
86 341 ['(((x+x)/(x*x))+x)'] a.k.a [x + 2/x]
87 383 ['(x+((((x/x)/x)*(x/x))+(((x/x)/x)*((x/x)*(x/x)))))'] a.k.a [x + 2/x]
88 137 ['(x+(((x+x)/x)/x))'] a.k.a [x + 2/x]
92 433 ['((x/(x*x))+((x/(x*x))+x))'] a.k.a [x + 2/x]
93 142 ['((x/(x*x))+((x/(x*x))+x))'] a.k.a [x + 2/x]
95 426 ['(x+(((x+x)/x)/x))'] a.k.a [x + 2/x]
96 136 ['(x+(((x+x)/x)/x))'] a.k.a [x + 2/x]
98 296 ['((x+((x/x)/x))+((x/x)/x))'] a.k.a [x + 2/x]
99 230 ['(((((x+x)+(x+x))/(x+x))/x)+x)'] a.k.a [x + 2/x]
[136]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
print("Avg. number of function evaluations from Tsoulos paper: ", 653 * 200)
ERT Expected run time - avg. number of function evaluations needed: 8977.27272727
Avg. number of function evaluations from Tsoulos paper: 130600
Consider the following Ordinary Differential Equation (ODE2):
\(\frac{dy}{dx} = \frac{1 - ycos(x)}{sin(x)}\), with \(y(0.1) = \frac{2.1}{sin(0,1)}\) and \(x \in [0.1,1]\)
we demand its punctual validity over a grid of \(N\) equally spaced points.
NOTE: The solution to the ODE is \(y = \frac{x+2}{sin(x)}\)
[137]:
# We construct the grid of points. Since the ODE only contains first order derivatives we use truncation order 1.
# Since we are using vectorized gdual we can instantiate only one gdual
values = np.linspace(0.1,1,10)
grid = gdual(values, "x", 1)
[138]:
# We define the quadratic error of the dCGP in the grid points
def qe_ODE2(dCGP, grid):
retval = 0
out = dCGP([grid])[0]
y = np.array(out.constant_cf)
dydx = np.array(out.get_derivative({"dx" : 1}))
x = np.array(grid.constant_cf)
ode2 = (1. - y * cos(x)) / sin(x)
retval += (ode2 - dydx) * (ode2 - dydx)
return sum(retval)
[139]:
# We define a penalty term associated to the initial conditions violation
dummy = (2.1)/sin(0.1)
def ic_ODE2(dCGP):
x0 = 0.1
y0 = dummy
out = dCGP([gdual([x0])])[0]
return (out.constant_cf[0] - y0) * (out.constant_cf[0] - y0)
[140]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
stop = 500
offsprings = 10
res = []
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 = randint(0,100000))
g, best_chromosome = run_experiment(max_gen = stop, offsprings = offsprings, \
quadratic_error=qe_ODE2, initial_conditions_error=ic_ODE2, dCGP=dCGP)
res.append(g)
dCGP.set(best_chromosome)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x"]), " a.k.a ", dCGP.simplify(["x"]))
res = np.array(res)
restart: gen: expression:
18 269 ['((((x/sin(x))+(x/sin(x)))/(sin(x)*(x/sin(x))))+(x/sin(x)))'] a.k.a [x/sin(x) + 2/sin(x)]
30 33 ['(((x/x)+((x/x)+x))/((x/x)*sin(x)))'] a.k.a [x/sin(x) + 2/sin(x)]
33 61 ['((((x+x)/sin(x))/x)+(((x+x)/sin(x))/((x+x)/x)))'] a.k.a [x/sin(x) + 2/sin(x)]
41 456 ['((((x/sin(x))/x)+((x/sin(x))/x))+(x/sin(x)))'] a.k.a [x/sin(x) + 2/sin(x)]
48 328 ['((((x+(sin(x)/sin(x)))*((x+(sin(x)/sin(x)))/(x+(sin(x)/sin(x)))))+((x+(sin(x)/sin(x)))/(x+(sin(x)/sin(x)))))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
59 127 ['((((x/x)*(x/x))+((x/x)+x))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
60 124 ['(((x/x)+((x/x)+x))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
62 191 ['((((x/x)+x)+(x/x))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
71 8 ['(((x+(x/x))+(x/x))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
76 37 ['((((x/x)+x)+(x/x))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
90 351 ['(((x/sin(x))+((x/sin(x))/x))+((x/sin(x))/x))'] a.k.a [x/sin(x) + 2/sin(x)]
96 377 ['(((x+((x*x)+x))/x)/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
98 352 ['(((x/x)+((x+(x/x))/(x/x)))/sin(x))'] a.k.a [x/sin(x) + 2/sin(x)]
[141]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
print("Avg. number of function evaluations from Tsoulos paper: ", 742 * 200)
ERT Expected run time - avg. number of function evaluations needed: 35482.3076923
Avg. number of function evaluations from Tsoulos paper: 148400
Consider the following Ordinary Differential Equation (ODE5):
\(\frac{d^2y}{dx^2} = 6\frac{dy}{dx} - 9y\), with \(y(0) = 0\), \(\frac{dy}{dx}(0)=2\) and \(x \in [0,1]\)
we demand its punctual validity over a grid of \(N\) equally spaced points.
NOTE: The solution to the ODE is \(y = 2x \exp(3x)\)
[142]:
# We construct the grid of points. Since the ODE only contains second order derivatives we use truncation order 2.
# Since we are using vectorized gdual we can instantiate only one gdual
values = np.linspace(0,1,10)
grid = gdual(values, "x", 2)
[143]:
# We define the quadratic error of the dCGP in the grid points
def qe_ODE5(dCGP, grid):
retval = 0
out = dCGP([grid])[0]
y = np.array(out.constant_cf)
dydx = np.array(out.get_derivative({"dx" : 1}))
dydx2 = np.array(out.get_derivative({"dx" : 2}))
x = np.array(grid.constant_cf)
ode5 = 6. * dydx - 9 * y
retval += (ode5 - dydx2) * (ode5 - dydx2)
return sum(retval)
[144]:
# We define a penalty term associated to the initial conditions violation
def ic_ODE5(dCGP):
x0 = 1e-16 # avoids what seems a numerical problem with vectorized dual?
y0 = 0.
dy0 = 2.
out = dCGP([gdual([x0], "x", 1)])[0]
dCGP_y0 = out.constant_cf[0]
dCGP_dy0 = out.get_derivative({"dx" : 1})[0]
return (dCGP_y0 - y0) * (dCGP_y0 - y0) + (dCGP_dy0 - dy0) * (dCGP_dy0 - dy0)
[145]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
stop = 500
offsprings = 10
res = []
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 = randint(0,100000))
g, best_chromosome = run_experiment(max_gen = stop, offsprings = offsprings, \
quadratic_error=qe_ODE5, initial_conditions_error=ic_ODE5, dCGP=dCGP)
res.append(g)
dCGP.set(best_chromosome)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x"]), " a.k.a ", dCGP.simplify(["x"]))
res = np.array(res)
restart: gen: expression:
8 220 ['(((exp(x)*exp(x))*exp(x))*log((exp(x)*exp(x))))'] a.k.a [2*x*exp(3*x)]
12 301 ['(((x*exp(x))*exp(x))*(exp(x)+exp(x)))'] a.k.a [2*x*exp(3*x)]
16 191 ['((x+x)*exp((x+(x+x))))'] a.k.a [2*x*exp(3*x)]
31 253 ['(log((exp(x)*exp(x)))*((exp(x)*exp(x))*exp(x)))'] a.k.a [2*x*exp(3*x)]
32 359 ['((x+x)*exp((x+(x+x))))'] a.k.a [2*x*exp(3*x)]
37 173 ['(((exp(x)*x)*(exp(x)*exp(x)))+((exp(x)*x)*(exp(x)*exp(x))))'] a.k.a [2*x*exp(3*x)]
48 216 ['((exp(x)*exp(x))*(exp(x)*(x+x)))'] a.k.a [2*x*exp(3*x)]
50 444 ['(((x+x)*exp((x+x)))*exp(x))'] a.k.a [2*x*exp(3*x)]
52 473 ['((x+x)*exp(((x+x)+x)))'] a.k.a [2*x*exp(3*x)]
64 179 ['(exp(x)*((exp(x)*(exp(x)+exp(x)))*x))'] a.k.a [2*x*exp(3*x)]
65 52 ['((x+x)*exp(((x+x)+x)))'] a.k.a [2*x*exp(3*x)]
68 193 ['((x+x)*log(exp(exp((x+(x+x))))))'] a.k.a [2*x*exp(3*x)]
71 291 ['((x+x)*exp((((x+x)-x)+(x+x))))'] a.k.a [2*x*exp(3*x)]
72 468 ['((x+x)*exp(((x+x)+x)))'] a.k.a [2*x*exp(3*x)]
73 102 ['(exp((x+(x+x)))*(x+x))'] a.k.a [2*x*exp(3*x)]
78 311 ['((exp(x)+exp(x))*((x*exp(x))*exp(x)))'] a.k.a [2*x*exp(3*x)]
80 263 ['((exp(x)*((exp(x)+exp(x))*exp(x)))*x)'] a.k.a [2*x*exp(3*x)]
83 287 ['((x+x)*exp((x+(x+x))))'] a.k.a [2*x*exp(3*x)]
88 115 ['(((x+x)+(x-x))*exp((x+(x+x))))'] a.k.a [2*x*exp(3*x)]
92 410 ['(exp((x+(x+x)))*(x+x))'] a.k.a [2*x*exp(3*x)]
[146]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
print("Avg. number of function evaluations from Tsoulos paper: ", 441 * 200)
ERT Expected run time - avg. number of function evaluations needed: 22610.5
Avg. number of function evaluations from Tsoulos paper: 88200
Consider the following non linear Ordinary Differential Equation (NLODE3):
\(\frac{d^2y}{dx^2}\frac{dy}{dx} = -\frac4{x^3}\), with \(y(1) = 0\), and \(x \in [1,2]\)
we demand its punctual validity over a grid of \(N\) equally spaced points.
NOTE: The solution to the ODE is \(y = log(x^2)\)
[147]:
# We construct the grid of points. Since the ODE only contains second order derivatives we use truncation order 2.
# Since we are using vectorized gdual we can instantiate only one gdual
values = np.linspace(1,2,10)
grid = gdual(values, "x", 2)
[148]:
# We define the quadratic error of the dCGP in the grid points
def qe_NLODE3(dCGP, grid):
retval = 0
out = dCGP([grid])[0]
y = np.array(out.constant_cf)
dydx = np.array(out.get_derivative({"dx" : 1}))
dydx2 = np.array(out.get_derivative({"dx" : 2}))
x = np.array(grid.constant_cf)
nlode3 = dydx2*dydx
retval += (nlode3 + 4/x/x/x) * (nlode3 + 4/x/x/x)
return sum(retval)
[149]:
# We define a penalty term associated to the initial conditions violation
def ic_NLODE3(dCGP):
x0 = 1.
y0 = 0.
out = dCGP([gdual([x0])])[0]
dCGP_y0 = out.constant_cf[0]
dCGP_dy0 = out.get_derivative({"dx" : 1})[0]
return (dCGP_y0 - y0) * (dCGP_y0 - y0)
[150]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
stop = 500
offsprings = 10
res = []
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 = randint(0,100000))
g, best_chromosome = run_experiment(max_gen = stop, offsprings = offsprings, \
quadratic_error=qe_NLODE3, initial_conditions_error=ic_NLODE3, dCGP=dCGP)
res.append(g)
dCGP.set(best_chromosome)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x"]), " a.k.a ", dCGP.simplify(["x"]))
res = np.array(res)
restart: gen: expression:
0 13 ['log((x*x))'] a.k.a [log(x**2)]
2 4 ['log((x*x))'] a.k.a [log(x**2)]
3 8 ['log((x*x))'] a.k.a [log(x**2)]
4 21 ['log((x*x))'] a.k.a [log(x**2)]
5 39 ['(log(x)+log(x))'] a.k.a [2*log(x)]
6 138 ['log((x*x))'] a.k.a [log(x**2)]
7 0 ['log(((x*x)-log((x/x))))'] a.k.a [log(x**2)]
9 10 ['log((x*x))'] a.k.a [log(x**2)]
10 5 ['log((x*x))'] a.k.a [log(x**2)]
11 9 ['(log(x)+log(x))'] a.k.a [2*log(x)]
12 50 ['(log(x)+log(x))'] a.k.a [2*log(x)]
13 15 ['log((x*x))'] a.k.a [log(x**2)]
14 11 ['((log(x)+log(x))+(x-x))'] a.k.a [2*log(x)]
16 80 ['log((x*x))'] a.k.a [log(x**2)]
17 2 ['(log(x)+log(x))'] a.k.a [2*log(x)]
18 47 ['(log(x)+log(x))'] a.k.a [2*log(x)]
19 1 ['(log(sin(x))-log(((x*sin(x))*x)))'] a.k.a [-log(x**2*sin(x)) + log(sin(x))]
21 7 ['log((x*x))'] a.k.a [log(x**2)]
24 17 ['log((x*x))'] a.k.a [log(x**2)]
26 39 ['log((x*x))'] a.k.a [log(x**2)]
29 6 ['(log(x)+log(x))'] a.k.a [2*log(x)]
30 55 ['log((x*x))'] a.k.a [log(x**2)]
31 54 ['log((x*x))'] a.k.a [log(x**2)]
32 188 ['(log(x)+log(x))'] a.k.a [2*log(x)]
33 4 ['(log(x)+log(x))'] a.k.a [2*log(x)]
34 32 ['(log(x)+(log(x)+sin((x-x))))'] a.k.a [2*log(x)]
35 66 ['(log(x)+log(x))'] a.k.a [2*log(x)]
36 36 ['log((x*x))'] a.k.a [log(x**2)]
37 56 ['(sin(sin((x-x)))+log((x*x)))'] a.k.a [log(x**2)]
38 5 ['log((x*x))'] a.k.a [log(x**2)]
39 11 ['(log((x+(x-x)))+log((x+(x-x))))'] a.k.a [2*log(x)]
40 9 ['log((((x/x)*x)*x))'] a.k.a [log(x**2)]
41 91 ['(log((x*x))-(log((x*x))+log((x*x))))'] a.k.a [-log(x**2)]
42 44 ['log((x*x))'] a.k.a [log(x**2)]
43 0 ['log((x*x))'] a.k.a [log(x**2)]
45 83 ['(log(x)+log(x))'] a.k.a [2*log(x)]
46 23 ['log((x*x))'] a.k.a [log(x**2)]
47 6 ['(log(x)+log(x))'] a.k.a [2*log(x)]
48 13 ['(log(x)+log(x))'] a.k.a [2*log(x)]
49 14 ['(log(x)+log(x))'] a.k.a [2*log(x)]
50 28 ['(log(x)+log(x))'] a.k.a [2*log(x)]
51 6 ['(log(x)+log(x))'] a.k.a [2*log(x)]
52 10 ['(log(x)+log(x))'] a.k.a [2*log(x)]
53 5 ['(log(x)+log(x))'] a.k.a [2*log(x)]
54 13 ['log((x*x))'] a.k.a [log(x**2)]
55 44 ['log((x*x))'] a.k.a [log(x**2)]
56 35 ['(log(x)+log(x))'] a.k.a [2*log(x)]
57 25 ['(log(x)+(log(x)-(log(x)-log(x))))'] a.k.a [2*log(x)]
58 11 ['log((x*x))'] a.k.a [log(x**2)]
59 46 ['((((x/x)-log(x))-log(x))-(x/x))'] a.k.a [-2*log(x)]
60 17 ['log((x*x))'] a.k.a [log(x**2)]
61 1 ['(log(x)+log(x))'] a.k.a [2*log(x)]
62 56 ['(log(x)+log(x))'] a.k.a [2*log(x)]
63 0 ['((log(x)+log(x))*cos(((log(x)-log(x))-(exp((log(x)-log(x)))*(log(x)-log(x))))))'] a.k.a [2*log(x)]
64 21 ['(log(x)+log(x))'] a.k.a [2*log(x)]
65 52 ['log((x*x))'] a.k.a [log(x**2)]
66 178 ['log((x*(x/cos((x-x)))))'] a.k.a [log(x**2)]
67 7 ['((x-x)-(log(x)+(log(x)-(x-x))))'] a.k.a [-2*log(x)]
68 2 ['log((x*x))'] a.k.a [log(x**2)]
69 2 ['log((x*x))'] a.k.a [log(x**2)]
70 3 ['(log(x)+log(x))'] a.k.a [2*log(x)]
71 69 ['log((x*x))'] a.k.a [log(x**2)]
72 24 ['((log(x)*(x+x))/x)'] a.k.a [2*log(x)]
73 3 ['(log(x)+log(x))'] a.k.a [2*log(x)]
74 25 ['log((x*x))'] a.k.a [log(x**2)]
75 121 ['(log(x)+log(x))'] a.k.a [2*log(x)]
76 29 ['log((x*x))'] a.k.a [log(x**2)]
77 1 ['(log(x)+log(x))'] a.k.a [2*log(x)]
78 4 ['(log(x)+log(x))'] a.k.a [2*log(x)]
79 299 ['log((x*x))'] a.k.a [log(x**2)]
80 3 ['(log(x)+log(x))'] a.k.a [2*log(x)]
81 33 ['log((x*x))'] a.k.a [log(x**2)]
82 19 ['log((x*x))'] a.k.a [log(x**2)]
83 2 ['log((x*x))'] a.k.a [log(x**2)]
84 278 ['log((((x-x)-x)*((x-x)-x)))'] a.k.a [log(x**2)]
85 69 ['log((x*x))'] a.k.a [log(x**2)]
86 7 ['(log(x)+((log(x)-log(x))+log(x)))'] a.k.a [2*log(x)]
87 11 ['log((x*x))'] a.k.a [log(x**2)]
88 1 ['log((x*x))'] a.k.a [log(x**2)]
89 10 ['(log((x*x))-(((x*x)+(x*x))-((x*x)+(x*x))))'] a.k.a [log(x**2)]
90 32 ['(log(x)+log(x))'] a.k.a [2*log(x)]
91 17 ['log((x*x))'] a.k.a [log(x**2)]
92 13 ['(log(x)+log(x))'] a.k.a [2*log(x)]
93 8 ['(log(x)+log(x))'] a.k.a [2*log(x)]
94 3 ['(log(x)+log(x))'] a.k.a [2*log(x)]
95 28 ['(log(x)+log(x))'] a.k.a [2*log(x)]
96 22 ['log((x*x))'] a.k.a [log(x**2)]
97 7 ['(log(x)+log(x))'] a.k.a [2*log(x)]
98 6 ['log((x*x))'] a.k.a [log(x**2)]
99 62 ['(((log(x)+(x+x))+log(x))-(x+x))'] a.k.a [2*log(x)]
[151]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
print("Avg. number of function evaluations from Tsoulos paper: ", 191 * 200)
ERT Expected run time - avg. number of function evaluations needed: 896.666666667
Avg. number of function evaluations from Tsoulos paper: 38200
Consider the following non linear Patial Differential Equation (PDE2):
\(\nabla^2 \psi(x,y) = -\psi(x,y)\) with \(x\in[0,1]\), \(y\in[0,1]\) and boundary conditions: \(\psi(0,y) = 0\), \(\psi(1,y) = \sin(1)\cos(y)\)
we demand its punctual validity over a squared grid of \(N\) equally spaced points.
NOTE: The solution to the PDE is \(\psi(x,y) = \sin(x)\cos(y)\)
[152]:
# We construct the grid of points. Since the PDE only contains second order derivatives we use truncation order 2.
# Since we are using vectorized gdual we can instantiate only one gdual
N=10
values = np.linspace(0,1,N)
xval = np.append(values,[values]*(N-1))
yval = values.repeat(N)
grid = [gdual(xval, "x", 2), gdual(yval, "y", 2)]
[153]:
# We define the quadratic error of the dCGP in the grid points
def qe_PDE1(dCGP, grid):
retval = 0
out = dCGP([grid[0], grid[1]])[0]
psi = np.array(out.constant_cf)
dpsidx2 = np.array(out.get_derivative({"dx" : 2}))
dpsidy2 = np.array(out.get_derivative({"dy" : 2}))
x = np.array(grid[0].constant_cf)
y = np.array(grid[1].constant_cf)
pde1 = -2 * psi
retval += (pde1 - dpsidx2 - dpsidy2) * (pde1 - dpsidx2 - dpsidy2)
return sum(retval)
[154]:
# We define a penalty term associated to the initial conditions violation
sin1 = np.sin(1)
def ic_PDE1(dCGP):
x0 = gdual([0]*N)
y0 = gdual(values)
psi = dCGP([x0, y0])[0]
dCGP_psi = np.array(psi.constant_cf)
err1 = (dCGP_psi - 0.) * (dCGP_psi - 0.)
x0 = gdual([1]*N)
y0 = gdual(values)
psi = dCGP([x0, y0])[0]
dCGP_psi = psi.constant_cf
err2 = (dCGP_psi - sin1*np.cos(values)) * (dCGP_psi - sin1*np.cos(values))
return sum(err1) + sum(err2)
[155]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
stop = 500
offsprings = 10
res = []
print("restart: \t gen: \t expression:")
for i in range(nexp):
dCGP = expression(inputs=2, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,100000))
g, best_chromosome = run_experiment(max_gen = stop, offsprings = offsprings, \
quadratic_error=qe_PDE1, initial_conditions_error=ic_PDE1, screen_output=False, dCGP=dCGP)
res.append(g)
dCGP.set(best_chromosome)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x","y"]), " a.k.a ", dCGP.simplify(["x","y"]))
res = np.array(res)
restart: gen: expression:
1 291 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
7 410 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
9 487 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
16 254 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
20 203 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
30 81 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
32 409 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
35 472 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
46 270 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
49 375 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
51 245 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
55 201 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
62 53 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
70 363 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
71 481 ['(cos(y)*sin(x))'] a.k.a [sin(x)*cos(y)]
78 34 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
87 288 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
93 377 ['(cos(((y-y)-y))*sin(x))'] a.k.a [sin(x)*cos(y)]
95 252 ['(sin(x)*cos(y))'] a.k.a [sin(x)*cos(y)]
[156]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected Run Time - avg. number of function evaluations needed: ", ERT * offsprings)
print("Avg. number of function evaluations from Tsoulos paper: ", 203 * 200)
ERT Expected Run Time - avg. number of function evaluations needed: 24192.1052632
Avg. number of function evaluations from Tsoulos paper: 40600
Consider the following non linear Patial Differential Equation (PDE6):
\(\nabla^2 \psi(x,y) + \exp(\psi(x,y)) = 1 + x^2 + y^2 + \frac 4{(1 + x^2 + y^2)^2}\) with \(x\in[0.1,3]\), \(y\in[0.1,3]\) and boundary conditions: \(\psi(0,y) = \log(1+y^2)\), \(\psi(1,y) = \log(2+y^2)\)
we demand its punctual validity over a squared grid of \(N\) equally spaced points.
NOTE: The solution to the PDE is \(\psi(x,y) = \log(1 + x^2 + y^2)\)
[171]:
# We construct the grid of points. Since the PDE only contains second order derivatives we use truncation order 2.
# Since we are using vectorized gdual we can instantiate only one gdual
N=10
values = np.linspace(0.1,3,N)
xval = np.append(values,[values]*(N-1))
yval = values.repeat(N)
grid = [gdual(xval, "x", 2), gdual(yval, "y", 2)]
[172]:
# We define the quadratic error of the dCGP in the grid points
def qe_PDE6(dCGP, grid):
retval = 0
out = dCGP([grid[0], grid[1]])[0]
psi = np.array(out.constant_cf)
dpsidx2 = np.array(out.get_derivative({"dx" : 2}))
dpsidy2 = np.array(out.get_derivative({"dy" : 2}))
x = np.array(grid[0].constant_cf)
y = np.array(grid[1].constant_cf)
pde6 = 4./(1+x*x+y*y)**2
retval += (pde6 - dpsidx2 - dpsidy2) * (pde6 - dpsidx2 - dpsidy2)
return sum(retval) / N**2
[190]:
### We define a penalty term associated to the initial conditions violation
def ic_PDE6(dCGP):
x0 = gdual([0.1]*N)
y0 = gdual(values)
psi = dCGP([x0, y0])[0]
dCGP_psi = np.array(psi.constant_cf)
err1 = (dCGP_psi - np.log(1+values*values+0.01)) * (dCGP_psi - np.log(1+values*values+0.01))
x0 = gdual([3]*N)
y0 = gdual(values)
psi = dCGP([x0, y0])[0]
dCGP_psi = psi.constant_cf
err2 = (dCGP_psi - np.log(10+values*values)) * (dCGP_psi - np.log(10+values*values))
x0 = gdual(values)
y0 = gdual([0.1]*N)
psi = dCGP([x0, y0])[0]
dCGP_psi = np.array(psi.constant_cf)
err3 = (dCGP_psi - np.log(1+values*values+0.01)) * (dCGP_psi - np.log(1+values*values+0.01))
x0 = gdual(values)
y0 = gdual([3]*N)
psi = dCGP([x0, y0])[0]
dCGP_psi = psi.constant_cf
err4 = (dCGP_psi - np.log(10+values*values)) * (dCGP_psi - np.log(10+values*values))
return (sum(err1) + sum(err2) + sum(err4) + sum(err3))
[191]:
# We run nexp experiments to accumulate statistic for the ERT
kernels = kernel_set(["sum", "mul", "div", "diff", "log","exp","cos","sin"])() # note the call operator (returns the list of kernels)
nexp = 100
stop = 2000
offsprings = 10
res = []
print("restart: \t gen: \t expression:")
for i in range(nexp):
dCGP = expression(inputs=2, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,100000))
g, best_chromosome = run_experiment(max_gen = stop, offsprings = offsprings, \
quadratic_error=qe_PDE6, initial_conditions_error=ic_PDE6, screen_output=False, dCGP=dCGP)
res.append(g)
dCGP.set(best_chromosome)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x","y"]), " a.k.a ", dCGP.simplify(["x","y"]), " ", qe_PDE6(dCGP,grid)+ic_PDE6(dCGP))
res = np.array(res)
restart: gen: expression:
13 1586 ['log(((x*x)+((y*y)+((y*y)/(y*y)))))'] a.k.a [log(x**2 + y**2 + 1)] 2.05547700971e-28
36 1225 ['log((((x*x)+(y/y))+(y*y)))'] a.k.a [log(x**2 + y**2 + 1)] 7.18298744104e-32
47 511 ['log((((((x*x)/(x*x))*(y*y))+((x*x)/(x*x)))+(x*x)))'] a.k.a [log(x**2 + y**2 + 1)] 7.90590617962e-28
48 1917 ['log((((x/x)+(x*x))+(y*y)))'] a.k.a [log(x**2 + y**2 + 1)] 7.21996529597e-32
77 1230 ['log((((y/y)+(y*y))+(x*x)))'] a.k.a [log(x**2 + y**2 + 1)] 7.17836520917e-32
93 1837 ['log((((y*y)+(y/y))+(x*x)))'] a.k.a [log(x**2 + y**2 + 1)] 7.17836520917e-32
[192]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected Run Time - avg. number of function evaluations needed: ", ERT * offsprings)
print("Avg. number of function evaluations from Tsoulos paper: ", 797 * 1000) # here we assume tsoulos used the maximum population size since PDE6 is the most difficult problem
ERT Expected Run Time - avg. number of function evaluations needed: 327020.0
Avg. number of function evaluations from Tsoulos paper: 797000