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