Discovery of prime integrals with dCGP
Lets first import dcgpy and pyaudi and set up things as to use dCGP on gduals defined over vectorized floats
[1]:
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, random
np.seterr(all='ignore') # avoids numpy complaining about early on malformed expressions being evalkuated
%matplotlib inline
We consider a set of differential equations in the form:
\[\begin{split}\left\{
\begin{array}{c}
\frac{dx_1}{dt} = f_1(x_1, \cdots, x_n) \\
\vdots \\
\frac{dx_n}{dt} = f_n(x_1, \cdots, x_n)
\end{array}
\right.\end{split}\]
and we search for expressions \(P(x_1, \cdots, x_n) = 0\) which we call prime integrals of motion.
The straight forward approach to design such a search would be to represent \(P\) via a \(dCGP\) program and evolve its chromosome so that the expression, computed along points of some trajectory, evaluates to zero. This naive approach brings to the evolution of trivial programs that are identically zero and that “do not represent the intrinsic reltations between state varaibles” - Schmidt 2009.
Let us, though, differentiate \(P\) along a trajectory solution to the ODEs above. We get:
\[\frac{dP}{dt} = \sum_{i=0}^n \frac{\partial P}{\partial x_i} \frac{dx_i}{dt} = \sum_{i=0}^n \frac{\partial P}{\partial x_i} f_i = 0\]
we may try to evolve the expression \(P\) so that the above relation is satisfied on chosen points (belonging to a real trajectory or just defined on a grid). To avoid evolution to go towards trivial solutions, unlike Schmidt, we suppress all mutations that give raise to expressions for which \(\sum_{i=0}^n \left(\frac{\partial P}{\partial x_i}\right)^2 = 0\). That is, expressions that do not depend on the state.
A mass spring system
As a simple example, consider the following mass-spring system.
The ODEs are:
\[\begin{split}\left\{
\begin{array}{l}
\dot v = -kx \\
\dot x = v
\end{array}\right.\end{split}\]
We define a dCGP having three inputs (the state and the constant \(k\)) and one output (\(P\))
[2]:
kernels = kernel_set(["sum", "mul", "pdiv", "diff"])() # note the call operator (returns the list of kernels)
dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
We define 50 random control of points where we check that the prime integral holds: \(x \in [2,4]\), \(v \in [2,4]\) and \(k \in[2, 4]\)
[3]:
n_points = 50
x = []
v = []
k = []
for i in range(n_points):
x.append(random()*2 + 2)
v.append(random()*2 + 2)
k.append(random()*2 + 2)
x = gdual(x,"x",1)
v = gdual(v,"v",1)
k = gdual(k)
[4]:
def fitness_call(dCGP, x, v, k):
res = dCGP([x,v,k])[0]
dPdx = np.array(res.get_derivative({"dx": 1}))
dPdv = np.array(res.get_derivative({"dv": 1}))
xcoeff = np.array(x.constant_cf)
vcoeff = np.array(v.constant_cf)
kcoeff = np.array(k.constant_cf)
err = dPdx/dPdv - kcoeff * xcoeff / vcoeff
return sum(err * err), 3
[5]:
# We run an evolutionary strategy ES(1 + offspring)
def run_experiment(max_gen, offsprings, dCGP, x, v, k, screen_output=False):
chromosome = [1] * offsprings
fitness = [1] *offsprings
best_chromosome = dCGP.get()
best_fitness = 1e10
for g in range(max_gen):
for i in range(offsprings):
check = 0
while(check < 1e-3):
dCGP.set(best_chromosome)
dCGP.mutate_active(i+1) # we mutate a number of increasingly higher active genes
fitness[i], check = fitness_call(dCGP, x,v,k)
chromosome[i] = dCGP.get()
for i in range(offsprings):
if fitness[i] <= best_fitness:
if (fitness[i] != best_fitness) and screen_output:
dCGP.set(chromosome[i])
print("New best found: gen: ", g, " value: ", fitness[i], " ", dCGP.simplify(["x","v","k"]))
best_chromosome = chromosome[i]
best_fitness = fitness[i]
if best_fitness < 1e-12:
break
dCGP.set(best_chromosome)
return g, dCGP
[6]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
offsprings = 10
stop = 2000
res = []
print("restart: \t gen: \t expression:")
for i in range(nexp):
dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
g, dCGP = run_experiment(stop, 10, dCGP, x,v,k, False)
res.append(g)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["x","v","k"]), " a.k.a ", dCGP.simplify(["x","v","k"]))
one_sol = dCGP
res = np.array(res)
restart: gen: expression:
0 1221 ['(((k-k)+x)*(((v/k)/(x/((v/k)*k)))+((k-k)+x)))'] a.k.a [x**2 + v**2/k]
3 557 ['((v*v)+((x*k)*x))'] a.k.a [k*x**2 + v**2]
5 382 ['(((x*k)+((((v-x)/x)*v)+v))/(((v-x)/x)/(v-x)))'] a.k.a [k*x**2 + v**2]
8 1249 ['((((k/v)+v)/(k/v))+(x*x))'] a.k.a [x**2 + 1 + v**2/k]
10 9 ['((v*v)+(k*(x*x)))'] a.k.a [k*x**2 + v**2]
14 198 ['(((k-(k*x))*x)-((v*v)+(k*x)))'] a.k.a [-k*x**2 - v**2]
15 1281 ['(((k*x)/k)*((v/(x/v))+(k*x)))'] a.k.a [k*x**2 + v**2]
21 1595 ['(((k*x)*(((k/(k*x))*k)+x))+(v*v))'] a.k.a [k**2 + k*x**2 + v**2]
23 321 ['(((x*x)+v)-(((x*x)+((x*x)+v))+((v*v)/k)))'] a.k.a [-x**2 - v**2/k]
24 1804 ['((k/(x*k))/((v/(x/v))+(x*k)))'] a.k.a [1/(k*x**2 + v**2)]
31 1292 ['((v*v)+(k+(x*(k*x))))'] a.k.a [k*x**2 + k + v**2]
34 417 ['((k-v)/(((v/(k/v))+(x*x))*(k-v)))'] a.k.a [1/(x**2 + v**2/k)]
42 880 ['(((((x*k)-((v*v)/x))-(x*k))-(((x*k)*((v*v)/(v*v)))*((v*v)/(v*v))))*((x*k)*((v*v)/(v*v))))'] a.k.a [-k**2*x**2 - k*v**2]
45 849 ['((k*(x+k))*(((k-(((v/k)/(x+k))*v))-x)-((k*((v/k)/(x+k)))/(v+v))))'] a.k.a [k**3 - k**2*v**2/(k**2 + k*x) - k**2/(2*(k + x)) - k*v**2*x/(k**2 + k*x) - k*x**2 - k*x/(2*(k + x))]
46 252 ['((((x*x)+((v/k)*v))+(v/(v/k)))*((x*x)+((v/k)*v)))'] a.k.a [k*x**2 + v**2 + x**4 + 2*v**2*x**2/k + v**4/k**2]
49 1010 ['((v*v)+((v-v)-(x*(((v-v)-x)*k))))'] a.k.a [k*x**2 + v**2]
51 912 ['(((((v/x)/k)+(x+x))-((v/x)/k))*(x+(((v/x)/k)*v)))'] a.k.a [2*x**2 + 2*v**2/k]
53 1144 ['((v*(k-k))+(((v*v)/k)+(x*x)))'] a.k.a [x**2 + v**2/k]
57 700 ['(((v/k)*v)+(x*x))'] a.k.a [x**2 + v**2/k]
58 1444 ['(((x*(k*x))+(v*v))+((x*(k*x))+(v*v)))'] a.k.a [2*k*x**2 + 2*v**2]
59 842 ['(((v*v)+(k*(x*x)))+((v*v)+(k*(x*x))))'] a.k.a [2*k*x**2 + 2*v**2]
63 853 ['((k+(k+k))+((v/(k/v))+(x*x)))'] a.k.a [3*k + x**2 + v**2/k]
65 834 ['(((x*(k*x))+(v*v))*k)'] a.k.a [k**2*x**2 + k*v**2]
69 276 ['((x*(k*x))+(v*v))'] a.k.a [k*x**2 + v**2]
76 1568 ['((((x/x)/k)-(k+((x*k)*x)))-(v*v))'] a.k.a [-k*x**2 - k - v**2 + 1/k]
79 406 ['((x*(k*x))+(v*v))'] a.k.a [k*x**2 + v**2]
83 1888 ['((x*x)+((v*v)/k))'] a.k.a [x**2 + v**2/k]
87 955 ['((v/(k/v))+(v+((x*x)-v)))'] a.k.a [x**2 + v**2/k]
90 1928 ['((v*v)+((k*x)*x))'] a.k.a [k*x**2 + v**2]
93 1144 ['((((v*v)-(x*k))+((x*(k*x))+(k*x)))/k)'] a.k.a [x**2 + v**2/k]
94 690 ['((x*x)+((v*x)/((k*(x*x))/(v*x))))'] a.k.a [x**2 + v**2/k]
96 1161 ['(((((x*k)-v)-(x*k))*(k*v))+(k-((x*k)*(x*k))))'] a.k.a [-k**2*x**2 - k*v**2 + k]
99 636 ['(((v/(x/v))+(k*x))*(k*x))'] a.k.a [k**2*x**2 + k*v**2]
[7]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
ERT Expected run time - avg. number of function evaluations needed: 49888.1818182
[8]:
print(one_sol.simplify(["x","v","k"]))
plt.rcParams["figure.figsize"] = [20,20]
one_sol.visualize(["x","v","k"])
[k**2*x**2 + k*v**2]
[8]:
<matplotlib.image.AxesImage at 0x7f5acc2fc390>
Simple pendulum
Consider the simple pendulum problem. In particular its differential formulation:
The ODEs are:
\[\begin{split}\left\{
\begin{array}{l}
\dot \omega = - \frac gL\sin\theta \\
\dot \theta = \omega \\
\end{array}\right.\end{split}\]
We define a dCGP having three inputs (the state and the constant \(\frac gL\)) and one output (\(P\))
[9]:
kernels = kernel_set(["sum", "mul", "pdiv", "diff","sin","cos"])() # note the call operator (returns the list of kernels)
dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
We define 50 random control of points where we check that the prime integral holds: \(\omega \in [-1, 1]\), \(\theta \in [-1, 1]\), and \(\frac gL \in [1,2]\)
[10]:
n_points = 50
omega = []
theta = []
c = []
for i in range(n_points):
omega.append(random()*10 - 5)
theta.append(random()*10 - 5)
c.append(random()*10)
omega = gdual(omega,"omega",1)
theta = gdual(theta,"theta",1)
c = gdual(c)
[11]:
def fitness_call(dCGP, theta, omega, c):
res = dCGP([theta, omega, c])[0]
dPdtheta = np.array(res.get_derivative({"dtheta": 1}))
dPdomega = np.array(res.get_derivative({"domega": 1}))
thetacoeff = np.array(theta.constant_cf)
omegacoeff = np.array(omega.constant_cf)
ccoeff = np.array(c.constant_cf)
err = dPdtheta/dPdomega + (-ccoeff * np.sin(thetacoeff)) / omegacoeff
check = sum(dPdtheta*dPdtheta + dPdomega*dPdomega)
return sum(err * err ), check
[12]:
# We run an evolutionary strategy ES(1 + offspring)
def run_experiment(max_gen, offsprings, dCGP, theta, omega, c, screen_output=False):
chromosome = [1] * offsprings
fitness = [1] *offsprings
best_chromosome = dCGP.get()
best_fitness = 1e10
for g in range(max_gen):
for i in range(offsprings):
check = 0
while(check < 1e-3):
dCGP.set(best_chromosome)
dCGP.mutate_active(i+1) # we mutate a number of increasingly higher active genes
fitness[i], check = fitness_call(dCGP, theta, omega, c)
chromosome[i] = dCGP.get()
for i in range(offsprings):
if fitness[i] <= best_fitness:
if (fitness[i] != best_fitness) and screen_output:
dCGP.set(chromosome[i])
print("New best found: gen: ", g, " value: ", fitness[i], " ", dCGP.simplify(["theta","omega","c"]))
best_chromosome = chromosome[i]
best_fitness = fitness[i]
if best_fitness < 1e-12:
break
dCGP.set(best_chromosome)
return g, dCGP
[13]:
# We run nexp experiments to accumulate statistic for the ERT
nexp = 100
offsprings = 10
stop = 2000
res = []
print("restart: \t gen: \t expression:")
for i in range(nexp):
dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
g, dCGP = run_experiment(stop, 10, dCGP, theta, omega, c, False)
res.append(g)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["theta","omega","c"]), " a.k.a ", dCGP.simplify(["theta","omega","c"]))
one_sol = dCGP
res = np.array(res)
restart: gen: expression:
5 1087 ['(((c*cos(theta))+(c*cos(theta)))-(omega*omega))'] a.k.a [2*c*cos(theta) - omega**2]
9 1007 ['(((((omega*omega)/c)-cos(theta))-cos(theta))/c)'] a.k.a [-2*cos(theta)/c + omega**2/c**2]
11 727 ['(((cos(theta)-((omega*omega)/c))+cos(theta))-((omega*omega)/((omega*omega)/c)))'] a.k.a [-c + 2*cos(theta) - omega**2/c]
17 509 ['(cos(((theta-theta)+(theta+(theta-theta))))-(((omega*omega)-(theta-theta))/(c+c)))'] a.k.a [cos(theta) - omega**2/(2*c)]
19 1780 ['(sin(((omega*omega)/((omega*omega)/c)))*((cos(theta)-((omega*omega)/c))+cos(theta)))'] a.k.a [2*sin(c)*cos(theta) - omega**2*sin(c)/c]
31 246 ['(cos(theta)-(((omega*omega)/c)-cos(theta)))'] a.k.a [2*cos(theta) - omega**2/c]
34 154 ['(((c*cos(theta))+(cos(theta)*c))-(omega*omega))'] a.k.a [2*c*cos(theta) - omega**2]
38 506 ['((cos(theta)*c)-(((omega*omega)+c)-(cos(theta)*c)))'] a.k.a [2*c*cos(theta) - c - omega**2]
54 1063 ['((((c+c)*cos(theta))-(omega*omega))-sin((theta/theta)))'] a.k.a [2*c*cos(theta) - omega**2 - sin(1)]
62 966 ['((cos(theta)-((omega*omega)/c))+((omega*omega)/(c+c)))'] a.k.a [cos(theta) - omega**2/(2*c)]
64 1159 ['sin(cos((((omega*(omega/c))-cos(theta))-(cos(theta)+c))))'] a.k.a [sin(cos(c + 2*cos(theta) - omega**2/c))]
71 687 ['(c+((cos(theta)*c)+((((cos(theta)*c)*theta)/theta)-(omega*omega))))'] a.k.a [2*c*cos(theta) + c - omega**2]
72 1970 ['((omega*omega)-((c+c)*cos(theta)))'] a.k.a [-2*c*cos(theta) + omega**2]
86 1215 ['(((c-(cos(theta)*c))+(omega*omega))+(c-(cos(theta)*c)))'] a.k.a [-2*c*cos(theta) + 2*c + omega**2]
93 1545 ['((omega*omega)-((c+c)*cos(theta)))'] a.k.a [-2*c*cos(theta) + omega**2]
98 264 ['((omega/(c+c))*(omega-(cos(theta)/(omega/(c+c)))))'] a.k.a [-cos(theta) + omega**2/(2*c)]
[14]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
ERT Expected run time - avg. number of function evaluations needed: 114250.625
[15]:
print(one_sol.simplify(["theta","omega","c"]))
plt.rcParams["figure.figsize"] = [20,20]
one_sol.visualize(["theta","omega","c"])
[-cos(theta) + omega**2/(2*c)]
[15]:
<matplotlib.image.AxesImage at 0x7f5acc2352e8>
The two-body problem
Consider the two body problem. In particular its differential formulation in polar coordinates:
The ODEs are:
\[\begin{split}\left\{
\begin{array}{l}
\dot v = -\frac\mu{r^2} + r\omega^2 \\
\dot \omega = - 2 \frac{v\omega}{r} \\
\dot r = v \\
\dot \theta = \omega
\end{array}\right.\end{split}\]
We define a dCGP having five inputs (the state and the constant \(\mu\)) and one output (\(P\))
[16]:
kernels = kernel_set(["sum", "mul", "pdiv", "diff"])() # note the call operator (returns the list of kernels)
dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
We define 50 random control of points where we check that the prime integral holds: \(r \in [0.1,1.1]\), \(v \in [2,4]\), \(\omega \in [1,2]\) and \(\theta \in[2, 4]\) and \(\mu \in [1,2]\)
[17]:
n_points = 50
v = []
omega = []
r = []
theta = []
mu = []
for i in range(n_points):
v.append(random()*2 + 2)
omega.append(random()*1 + 1)
r.append(random() + 0.1)
theta.append(random()*2 + 2)
mu.append(random() + 1)
r = gdual(r,"r",1)
omega = gdual(omega,"omega",1)
v = gdual(v,"v",1)
theta = gdual(theta,"theta",1)
mu = gdual(mu)
[18]:
## Use this fitness if energy conservation is to be found (it basically forces the expression to depend on v)
def fitness_call(dCGP, r, v, theta, omega, mu):
res = dCGP([r, v, theta, omega, mu])[0]
dPdr = np.array(res.get_derivative({"dr": 1}))
dPdv = np.array(res.get_derivative({"dv": 1}))
dPdtheta = np.array(res.get_derivative({"dtheta": 1}))
dPdomega = np.array(res.get_derivative({"domega": 1}))
rcoeff = np.array(r.constant_cf)
vcoeff = np.array(v.constant_cf)
thetacoeff = np.array(theta.constant_cf)
omegacoeff = np.array(omega.constant_cf)
mucoeff = np.array(mu.constant_cf)
err = dPdr / dPdv + (-mucoeff/rcoeff**2 + rcoeff*omegacoeff**2) / vcoeff + dPdtheta / dPdv / vcoeff * omegacoeff + dPdomega / dPdv / vcoeff * (-2*vcoeff*omegacoeff/rcoeff)
check = sum(dPdr*dPdr + dPdv*dPdv + dPdomega*dPdomega + dPdtheta*dPdtheta)
return sum(err * err), check
## Use this fitness if any conservation is to be found (will always converge to angular momentum)
def fitness_call_free(dCGP, r, v, theta, omega, mu):
res = dCGP([r, v, theta, omega, mu])[0]
dPdr = np.array(res.get_derivative({"dr": 1}))
dPdv = np.array(res.get_derivative({"dv": 1}))
dPdtheta = np.array(res.get_derivative({"dtheta": 1}))
dPdomega = np.array(res.get_derivative({"domega": 1}))
rcoeff = np.array(r.constant_cf)
vcoeff = np.array(v.constant_cf)
thetacoeff = np.array(theta.constant_cf)
omegacoeff = np.array(omega.constant_cf)
mucoeff = np.array(mu.constant_cf)
err = dPdr * vcoeff + dPdv * (-mucoeff/rcoeff**2 + rcoeff*omegacoeff**2) + dPdtheta * omegacoeff + dPdomega * (-2*vcoeff*omegacoeff/rcoeff)
check = sum(dPdr*dPdr + dPdv*dPdv +dPdomega*dPdomega+ dPdtheta*dPdtheta)
return sum(err * err ), check
[19]:
# We run an evolutionary strategy ES(1 + offspring)
def run_experiment(max_gen, offsprings, dCGP, r, v, theta, omega, mu, obj_fun, screen_output=False):
chromosome = [1] * offsprings
fitness = [1] *offsprings
best_chromosome = dCGP.get()
best_fitness = 1e10
for g in range(max_gen):
for i in range(offsprings):
check = 0
while(check < 1e-3):
dCGP.set(best_chromosome)
dCGP.mutate_active(i+1) # we mutate a number of increasingly higher active genes
fitness[i], check = obj_fun(dCGP, r, v, theta, omega, mu)
chromosome[i] = dCGP.get()
for i in range(offsprings):
if fitness[i] <= best_fitness:
if (fitness[i] != best_fitness) and screen_output:
dCGP.set(chromosome[i])
print("New best found: gen: ", g, " value: ", fitness[i], " ", dCGP.simplify(["r","v","theta","omega","mu"]))
best_chromosome = chromosome[i]
best_fitness = fitness[i]
if best_fitness < 1e-12:
break
dCGP.set(best_chromosome)
return g, dCGP
[20]:
# We run nexp experiments to accumulate statistic for the ERT (angular momentum case)
nexp = 100
offsprings = 10
stop = 2000 #100000
res = []
print("restart: \t gen: \t expression:")
for i in range(nexp):
dCGP = expression(inputs=5, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
g, dCGP = run_experiment(stop, 10, dCGP, r, v, theta, omega, mu, fitness_call_free, False)
res.append(g)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["r","v","theta","omega","mu"]), " a.k.a ", dCGP.simplify(["r","v","theta","omega","mu"]))
one_sol = dCGP
res = np.array(res)
restart: gen: expression:
0 117 ['((r+r)/((r+r)+((mu/omega)/(r+r))))'] a.k.a [2*r/(mu/(2*omega*r) + 2*r)]
1 76 ['((omega*r)*r)'] a.k.a [omega*r**2]
2 421 ['((mu+theta)/((((mu+theta)*(theta-(mu+theta)))*r)*(r*omega)))'] a.k.a [-1/(mu*omega*r**2)]
3 22 ['((omega*((theta*(theta-theta))-r))*((theta*(theta-theta))-r))'] a.k.a [omega*r**2]
4 634 ['(((omega/omega)/(mu+((r*omega)*r)))/(omega/omega))'] a.k.a [1/(mu + omega*r**2)]
5 408 ['((r*(omega*r))/(theta/theta))'] a.k.a [omega*r**2]
6 533 ['(mu-((r*omega)*(r/mu)))'] a.k.a [mu - omega*r**2/mu]
7 73 ['((omega/((r/mu)*omega))/((r/mu)*omega))'] a.k.a [mu**2/(omega*r**2)]
8 1573 ['((r*(mu*omega))*(v-(r+v)))'] a.k.a [-mu*omega*r**2]
9 196 ['((((r/theta)*theta)/(mu*(theta/v)))*((omega*(mu*(theta/v)))*r))'] a.k.a [omega*r**2]
10 21 ['(omega*(r*r))'] a.k.a [omega*r**2]
12 316 ['((r*(omega/v))*(v*r))'] a.k.a [omega*r**2]
14 6 ['(((r*r)*omega)/mu)'] a.k.a [omega*r**2/mu]
16 68 ['((omega*r)*r)'] a.k.a [omega*r**2]
18 1660 ['((r*omega)*r)'] a.k.a [omega*r**2]
19 267 ['((theta/r)/((r*omega)*theta))'] a.k.a [1/(omega*r**2)]
20 1392 ['((((mu+theta)-theta)/omega)/((r*((mu+theta)-theta))*r))'] a.k.a [1/(omega*r**2)]
21 231 ['((((omega/v)/(((r*omega)*r)*((r*omega)*r)))/((r*omega)*r))/(omega/v))'] a.k.a [1/(omega**3*r**6)]
23 337 ['((((theta/theta)/((theta/theta)*mu))/((mu/r)/omega))/(mu/r))'] a.k.a [omega*r**2/mu**3]
24 27 ['((mu/(mu-theta))/(omega/((((mu-theta)/r)*((mu-theta)/r))*(mu/(mu-theta)))))'] a.k.a [mu**2/(omega*r**2)]
25 42 ['(((omega*(r*theta))*mu)/(((omega*(r*theta))+(theta/r))*mu))'] a.k.a [omega*r*theta/(omega*r*theta + theta/r)]
26 16 ['(theta/(((r*r)*theta)*omega))'] a.k.a [1/(omega*r**2)]
27 261 ['(r*((((v*(mu*mu))*(mu*theta))-(((v*(mu*mu))*(mu*theta))+omega))*r))'] a.k.a [-omega*r**2]
28 285 ['(mu*(omega*(r*r)))'] a.k.a [mu*omega*r**2]
29 1021 ['((mu/(((r*omega)*r)+mu))/(((r*omega)*r)+mu))'] a.k.a [mu/(mu**2 + 2*mu*omega*r**2 + omega**2*r**4)]
30 120 ['(((r*omega)*(v+r))+((r*omega)*(r-v)))'] a.k.a [2*omega*r**2]
31 528 ['(omega/((omega*r)*(omega*r)))'] a.k.a [1/(omega*r**2)]
33 931 ['(mu/(r*(omega*r)))'] a.k.a [mu/(omega*r**2)]
34 153 ['(((mu/r)/(r/mu))/omega)'] a.k.a [mu**2/(omega*r**2)]
35 113 ['((((r*r)*omega)+mu)/mu)'] a.k.a [1 + omega*r**2/mu]
36 142 ['((omega*mu)*((r*r)*mu))'] a.k.a [mu**2*omega*r**2]
37 193 ['((mu/(mu+mu))-((r*omega)*r))'] a.k.a [-omega*r**2 + 1/2]
38 5 ['(r/((mu/omega)/r))'] a.k.a [omega*r**2/mu]
39 31 ['((r*omega)*r)'] a.k.a [omega*r**2]
40 87 ['((omega/(mu/r))/(mu/r))'] a.k.a [omega*r**2/mu**2]
41 118 ['(r*(r*omega))'] a.k.a [omega*r**2]
42 301 ['((r*(omega*r))-(((theta-theta)+(theta-theta))/r))'] a.k.a [omega*r**2]
43 1147 ['(r*(omega*r))'] a.k.a [omega*r**2]
44 649 ['((r*(r*(omega*theta)))/(mu*theta))'] a.k.a [omega*r**2/mu]
46 57 ['((r*r)/(mu/omega))'] a.k.a [omega*r**2/mu]
47 309 ['((r/(((r*mu)*omega)/((r*mu)*omega)))*((r*mu)*omega))'] a.k.a [mu*omega*r**2]
48 194 ['(((mu*(omega/mu))*r)*r)'] a.k.a [omega*r**2]
49 516 ['(r*(omega*r))'] a.k.a [omega*r**2]
50 623 ['((omega*r)*r)'] a.k.a [omega*r**2]
51 291 ['(r*(r*omega))'] a.k.a [omega*r**2]
52 245 ['(omega/((mu/r)/r))'] a.k.a [omega*r**2/mu]
53 60 ['((omega/(mu*mu))*(r*(r-(v-v))))'] a.k.a [omega*r**2/mu**2]
56 202 ['(((omega*r)*(omega*r))/omega)'] a.k.a [omega*r**2]
57 85 ['((r*r)*omega)'] a.k.a [omega*r**2]
58 72 ['(r/(((theta*mu)/omega)/(r*theta)))'] a.k.a [omega*r**2/mu]
59 5 ['(((mu*mu)/(mu*r))/(omega*(mu*r)))'] a.k.a [1/(omega*r**2)]
60 7 ['(omega*((r+r)*r))'] a.k.a [2*omega*r**2]
61 545 ['(r*(r*omega))'] a.k.a [omega*r**2]
62 178 ['((((r*v)*omega)*(r/v))*(((r*v)*omega)*(r/v)))'] a.k.a [omega**2*r**4]
63 233 ['((r*((r*omega)*mu))*(r*((r*omega)*mu)))'] a.k.a [mu**2*omega**2*r**4]
64 1 ['((r/r)/(omega*(r*r)))'] a.k.a [1/(omega*r**2)]
65 100 ['(((mu/r)*(mu/r))/omega)'] a.k.a [mu**2/(omega*r**2)]
66 389 ['(((omega/(theta/((r*theta)*r)))-(mu-(r*theta)))+(mu-(r*theta)))'] a.k.a [omega*r**2]
67 17 ['(r/((mu/r)/omega))'] a.k.a [omega*r**2/mu]
68 19 ['(r*(r*omega))'] a.k.a [omega*r**2]
69 29 ['((r*r)*omega)'] a.k.a [omega*r**2]
70 882 ['(((((theta+theta)*((theta+theta)-mu))/omega)/r)/(r*((theta+theta)*((theta+theta)-mu))))'] a.k.a [1/(omega*r**2)]
71 965 ['(((mu/r)/(mu*omega))/r)'] a.k.a [1/(omega*r**2)]
72 363 ['(((r*r)/(r*r))*((r*r)*omega))'] a.k.a [omega*r**2]
73 397 ['((omega*r)*r)'] a.k.a [omega*r**2]
74 270 ['((omega/(r*(omega*mu)))/(r*(omega*mu)))'] a.k.a [1/(mu**2*omega*r**2)]
75 3 ['(r*(omega*r))'] a.k.a [omega*r**2]
76 98 ['((mu/omega)/((mu*r)*(mu*r)))'] a.k.a [1/(mu*omega*r**2)]
77 161 ['(r*(r*omega))'] a.k.a [omega*r**2]
78 947 ['((mu/r)*((mu/r)/omega))'] a.k.a [mu**2/(omega*r**2)]
79 5 ['((omega-(omega/mu))/((mu/(r*r))+(mu/(r*r))))'] a.k.a [omega*r**2/(2*mu) - omega*r**2/(2*mu**2)]
80 169 ['(omega/((((mu/r)-v)+v)*(mu/r)))'] a.k.a [omega*r**2/mu**2]
81 568 ['((r*omega)*r)'] a.k.a [omega*r**2]
82 400 ['(((r*r)*(r*r))*(omega/(r*r)))'] a.k.a [omega*r**2]
83 37 ['(omega*(r*r))'] a.k.a [omega*r**2]
84 211 ['(((v-v)-((r-omega)+omega))*(r*omega))'] a.k.a [-omega*r**2]
85 1423 ['(((omega/((mu/theta)/r))*r)/theta)'] a.k.a [omega*r**2/mu]
86 17 ['(r*((omega/(v/r))*v))'] a.k.a [omega*r**2]
87 127 ['(r*((r*mu)*omega))'] a.k.a [mu*omega*r**2]
88 1078 ['(r*(omega*r))'] a.k.a [omega*r**2]
89 7 ['((mu-((r*(r*omega))-mu))*(mu-((r*(r*omega))-mu)))'] a.k.a [4*mu**2 - 4*mu*omega*r**2 + omega**2*r**4]
90 376 ['(((mu/r)/((r+(theta+theta))-(theta+theta)))/(omega*mu))'] a.k.a [1/(omega*r**2)]
91 45 ['(((omega/(mu+r))/(omega/(mu+r)))/(((r*(r/mu))*mu)*(mu*omega)))'] a.k.a [1/(mu*omega*r**2)]
92 219 ['((r*r)*(omega*mu))'] a.k.a [mu*omega*r**2]
93 679 ['((mu/omega)/((r+r)*(r+r)))'] a.k.a [mu/(4*omega*r**2)]
94 608 ['(omega/(((theta/r)/theta)*((theta/r)/theta)))'] a.k.a [omega*r**2]
95 15 ['((omega*r)*r)'] a.k.a [omega*r**2]
96 259 ['((r*(r*omega))-mu)'] a.k.a [-mu + omega*r**2]
97 829 ['((((r/v)*v)*omega)*((r+v)-v))'] a.k.a [omega*r**2]
98 101 ['(((r+(mu*r))+(mu*r))*(omega*r))'] a.k.a [2*mu*omega*r**2 + omega*r**2]
99 18 ['((r/(mu/(omega+omega)))*r)'] a.k.a [2*omega*r**2/mu]
[21]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
ERT Expected run time - avg. number of function evaluations needed: 5270.98901099
[22]:
print(one_sol.simplify(["r","v","theta","omega","mu"]))
plt.rcParams["figure.figsize"] = [20,20]
one_sol.visualize(["r","v","theta","omega","mu"])
[2*omega*r**2/mu]
[22]:
<matplotlib.image.AxesImage at 0x7f5acc158cc0>
[23]:
# We run nexp experiments to accumulate statistic for the ERT (angular momentum case)
nexp = 100
offsprings = 10
stop = 100000
res = []
print("restart: \t gen: \t expression:")
for i in range(nexp):
dCGP = expression(inputs=5, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,234213213))
g, dCGP = run_experiment(stop, 10, dCGP, r, v, theta, omega, mu, fitness_call, False)
res.append(g)
if g < (stop-1):
print(i, "\t\t", res[i], "\t", dCGP(["r","v","theta","omega","mu"]), " a.k.a ", dCGP.simplify(["r","v","theta","omega","mu"]))
one_sol = dCGP
res = np.array(res)
restart: gen: expression:
2 25413 ['(((v-r)*(v+r))-((((mu-r)+(mu-r))/r)-((r-(r*omega))*(r-(r*omega)))))'] a.k.a [-2*mu/r + omega**2*r**2 - 2*omega*r**2 + v**2 + 2]
3 46221 ['((v*(v*(mu+mu)))-((((mu+mu)*(mu+mu))/((r*r)/r))-(((omega*omega)*(r*r))*(mu+mu))))'] a.k.a [-4*mu**2/r + 2*mu*omega**2*r**2 + 2*mu*v**2]
6 97122 ['(((((v*(v*r))-mu)+(((omega*r)*(r*(omega*r)))-mu))*(omega*r))*mu)'] a.k.a [-2*mu**2*omega*r + mu*omega**3*r**4 + mu*omega*r**2*v**2]
8 28010 ['(((mu/r)-((v*v)-(mu/r)))-((r*omega)*(r*omega)))'] a.k.a [2*mu/r - omega**2*r**2 - v**2]
16 22942 ['(((omega*r)*((mu/(mu/r))+(omega*r)))+(mu-((mu/r)+((mu/r)-(v*v)))))'] a.k.a [mu - 2*mu/r + omega**2*r**2 + omega*r**2 + v**2]
28 63318 ['(((((mu/r)+mu)-((omega*r)*(omega*r)))+(((mu/r)+mu)-(v*v)))+((((mu/r)+mu)-((omega*r)*(omega*r)))+(((mu/r)+mu)-(v*v))))'] a.k.a [4*mu + 4*mu/r - 2*omega**2*r**2 - 2*v**2]
36 12543 ['(((omega*r)*(r-(omega*r)))+((mu/r)-((v*v)-(mu/r))))'] a.k.a [2*mu/r - omega**2*r**2 + omega*r**2 - v**2]
42 64916 ['((((omega*r)*(omega*r))-mu)-(((mu/(r*v))+((mu/(r*v))-v))*v))'] a.k.a [-mu - 2*mu/r + omega**2*r**2 + v**2]
48 86173 ['(((r*(omega*omega))*r)+((((v*r)*v)+((r-mu)+(r-mu)))/r))'] a.k.a [-2*mu/r + omega**2*r**2 + v**2 + 2]
56 29174 ['(((omega*(omega+(mu/mu)))/(((mu/mu)/r)/r))-((mu/r)+((mu/r)-(v*v))))'] a.k.a [-2*mu/r + omega**2*r**2 + omega*r**2 + v**2]
59 47790 ['((((r*omega)*(r*omega))-(((mu-r)+(mu-r))/r))+((mu/mu)+(v*v)))'] a.k.a [-2*mu/r + omega**2*r**2 + v**2 + 3]
61 36982 ['((((r*omega)*omega)*r)-(((mu/r)+(mu/r))-(v*v)))'] a.k.a [-2*mu/r + omega**2*r**2 + v**2]
81 34496 ['((v*v)-((((((mu/r)+(mu/r))+(omega*(r*r)))-((omega*(r*r))*omega))+(omega*(r*r)))+mu))'] a.k.a [-mu - 2*mu/r + omega**2*r**2 - 2*omega*r**2 + v**2]
86 58369 ['(((omega*r)*(omega*r))-(((mu/r)+(mu/r))-(v*v)))'] a.k.a [-2*mu/r + omega**2*r**2 + v**2]
90 63291 ['(((omega*r)*(omega*r))+((v*v)+(((omega*(mu/(omega*r)))*r)-((mu+mu)/r))))'] a.k.a [mu - 2*mu/r + omega**2*r**2 + v**2]
[ ]:
ERT = sum(res) / sum(res<(stop-1))
print("ERT Expected run time - avg. number of function evaluations needed: ", ERT * offsprings)
[ ]: