Lets first import dcgpy and pyaudi and set up things as to use dCGP on gduals defined over vectorized floats
In [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
In [3]:
kernels = kernel_set(["sum", "mul", "div", "diff", "log", "sin", "cos", "exp"])() # note the call operator (returns the list of kernels)
In [6]:
dCGP_example = expression(inputs=1, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = 2134)
In [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
In [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
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
In [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)
$\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)$
In [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)]
In [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)
In [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)
In [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)
In [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)
$\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)$
In [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)]
In [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
In [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))
In [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)
In [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