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