{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning constants in a symbolic regression task (PART II) (deprecated)\n",
    "\n",
    "One of the long standing \"skeletons in the closet\" of GP techniques is the constant finding problem. It is widely acknoledged that the \"ephemeral random constant\", de facto the main solution proposed to this problem, is far from being satisfactory\n",
    "\n",
    "\n",
    "Using dCGP, we are here able to succesfully learn constants as well as expressions during evolution thanks to the hybridization of the evolutionary strategy with what may, essentially, be seen as a second order back-propagation algorithm\n",
    "\n",
    "NOTE: since v1.4 symbolic regression is performed via dedicated classes and not manipulating directly the dcgpy.expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets first import dcgpy and pyaudi and set up things as to compute our CGP using the type \"gdual\" and thus get for free all derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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",
    "import pyaudi\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from random import randint\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The ES-(1+$\\lambda$) algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_experiment(dCGP, offsprings, max_gen, x, yt, screen_output):\n",
    "    # The offsprings chromosome, fitness and constant\n",
    "    chromosome = [1] * offsprings\n",
    "    fitness = [1] *offsprings\n",
    "    constant = [[1.,1.1]]*offsprings\n",
    "    # Init the best as the initial random dCGP\n",
    "    best_chromosome = dCGP.get()\n",
    "    best_constants = [1,1]\n",
    "    fit, _ = err2(dCGP, x, yt, best_constants)\n",
    "    best_fitness = fit\n",
    "    # Main loop over generations\n",
    "    for g in range(max_gen):\n",
    "        for i in range(offsprings):\n",
    "            dCGP.set(best_chromosome)\n",
    "            cumsum=0\n",
    "            dCGP.mutate_active(i+1)\n",
    "            fit, constant[i] = err2(dCGP, x, yt, best_constants)\n",
    "            fitness[i] = fit\n",
    "            chromosome[i] = dCGP.get()\n",
    "        for i in range(offsprings):\n",
    "            if fitness[i] <= best_fitness:\n",
    "                if (fitness[i] != best_fitness):\n",
    "                    best_chromosome = chromosome[i]\n",
    "                    best_fitness = fitness[i]\n",
    "                    best_constants = constant[i]\n",
    "                    dCGP.set(best_chromosome)\n",
    "                    if screen_output:\n",
    "                        print(\"New best found: gen: \", g, \" value: \", fitness[i],  dCGP.simplify([\"x\",\"c1\",\"c2\"]), \"c =\", best_constants)\n",
    "\n",
    "        if best_fitness < 1e-14:\n",
    "            break\n",
    "    return g, best_chromosome, best_constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The test problems\n",
    "As target functions, we define three different problems of increasing complexity:\n",
    "\n",
    "P5: $e x^5 - \\pi x^3 + x$\n",
    "\n",
    "P6: $\\frac{e x^2-1}{\\pi (x + 2)}$\n",
    "\n",
    "P7: $\\sin(e x) + \\cos(\\pi x)$\n",
    "\n",
    "note how $\\pi$ and $e$ are present in the expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following functions create the target values for a gridded input x for different test problems\n",
    "def data_P5(x):\n",
    "    return np.e * x**5 - np.pi*x**3 + x\n",
    "def data_P6(x):\n",
    "    return (np.e*x**2-1) / (np.pi*(x + 2))\n",
    "def data_P7(x):\n",
    "    return pyaudi.sin(np.e*x)+pyaudi.cos(np.pi*x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The error functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is used to sum over the component of a vectorized coefficient, accounting for the fact that if its dimension\n",
    "# is 1, then it could represent [a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a ...] with [a]\n",
    "def collapse_vectorized_coefficient(x, N):\n",
    "    if len(x) == N:\n",
    "        return sum(x)\n",
    "    return x[0] * N\n",
    "\n",
    "# This is the quadratic error of a dCGP expression when the constant value is cin. The error is computed\n",
    "# over the input points xin (of type gdual, order 0 as we are not interested in expanding the program w.r.t. these)\n",
    "# The target values are contained in yt (of type gdual, order 0 as we are not interested in expanding the program w.r.t. these)\n",
    "def err(dCGP, xin, yt, c1, c2):\n",
    "    y = dCGP([xin,c1,c2])[0]\n",
    "    return (y-yt)**2 / len(xin.constant_cf)\n",
    "\n",
    "# This is the quadratic error of the expression when the value of the constants are learned using a, one step, \n",
    "# second order method.\n",
    "def err2(dCGP, xin, yt, constants):\n",
    "    cin1 = constants[0]\n",
    "    cin2 = constants[1]\n",
    "    # Lets comupte the Taylor expansion (order 2) of the quadratic error around the point cin1, cin2\n",
    "    c1 = gdual([cin1], \"c1\", 2)\n",
    "    c2 = gdual([cin2], \"c2\", 2)\n",
    "    error = err(dCGP, x, yt, c1, c2)\n",
    "    #initial error for the current values of constants\n",
    "    ierr = sum(error.constant_cf) \n",
    "    # Number of points used in the data\n",
    "    N = len(xin.constant_cf)\n",
    "\n",
    "    # We try one step of the Newton method\n",
    "    for i in range(1):\n",
    "        G1 = error.get_derivative({\"dc1\":1})\n",
    "        G2 = error.get_derivative({\"dc2\":1})\n",
    "        H11 = error.get_derivative({\"dc1\":2})\n",
    "        H22 = error.get_derivative({\"dc2\":2})\n",
    "        H12 = error.get_derivative({\"dc1\":1, \"dc2\":1})\n",
    "        G1 = collapse_vectorized_coefficient(G1, N)\n",
    "        G2 = collapse_vectorized_coefficient(G2, N)\n",
    "        H11 = collapse_vectorized_coefficient(H11, N)\n",
    "        H22 = collapse_vectorized_coefficient(H22, N)\n",
    "        H12 = collapse_vectorized_coefficient(H12, N)\n",
    "\n",
    "        H = [[H11,H12],[H12, H22]]\n",
    "        det = np.linalg.det(H)\n",
    "        if det != 0:\n",
    "            invH = np.linalg.inv(H)\n",
    "            # We write the Newton update formula explicitly\n",
    "            dc1 = - invH[0,0] * G1 - invH[0,1] * G2\n",
    "            dc2 = - invH[1,0] * G1 - invH[1,1] * G2\n",
    "            dc1*=1\n",
    "            dc2*=1\n",
    "            c1+=dc1\n",
    "            c2+=dc2\n",
    "            error = err(dCGP, x, yt, c1, c2)\n",
    "        else:\n",
    "            break\n",
    "            \n",
    "    # error after one Newton step\n",
    "    aerr = sum(error.constant_cf)\n",
    "\n",
    "    # if no improvement, take 5 steps of gradient descent with learing rate 0.05\n",
    "    if aerr >= ierr:\n",
    "        c1 = gdual([cin1], \"c1\", 1)\n",
    "        c2 = gdual([cin2], \"c2\", 1)\n",
    "        # We end up with a few gradient descent steps\n",
    "        for i in range(5):\n",
    "            G1 = sum(error.get_derivative({\"dc1\":1}))\n",
    "            G2 = sum(error.get_derivative({\"dc2\":1}))\n",
    "            dc1 = - G1 \n",
    "            dc2 = - G2\n",
    "            dc1*=0.05\n",
    "            dc2*=0.05\n",
    "            c1+=dc1\n",
    "            c2+=dc2\n",
    "            error = err(dCGP, x, yt, c1, c2)\n",
    "            \n",
    "    aerr = sum(error.constant_cf)      \n",
    "    return aerr, [c1.constant_cf[0], c2.constant_cf[0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Problem P5:  $ex^5 - \\pi x^3 + x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(1,3,10)\n",
    "x = gdual(x)\n",
    "yt = data_P5(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.5/site-packages/numpy/linalg/linalg.py:1776: RuntimeWarning: invalid value encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n",
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:47: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:48: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7 \t\t 1185 \t ['(x+(((x+x)*(x*x))*((c1*(x*x))-c2)))']  a.k.a  [2*c1*x**5 - 2*c2*x**3 + x] c =  [1.3591409142295225, 1.5707963267948966]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.5/site-packages/numpy/linalg/linalg.py:1776: RuntimeWarning: invalid value encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n",
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:47: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:48: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12 \t\t 1367 \t ['((x-((((x*x)*x)*(x*x))*c1))-(c2*((x*x)*x)))']  a.k.a  [-c1*x**5 - c2*x**3 + x] c =  [-2.718281828459045, 3.141592653589792]\n",
      "14 \t\t 696 \t ['(x-(((x*x)*(((x*x)*(c1-c2))-c1))*x))']  a.k.a  [-c1*x**5 + c1*x**3 + c2*x**5 + x] c =  [-3.141592653589789, -0.42331082513074425]\n",
      "18 \t\t 1865 \t ['(((c2*x)+((((c2*x)*x)-c1)*(((c2*x)*x)*x)))/c2)']  a.k.a  [-c1*x**3 + c2*x**5 + x] c =  [3.1415926535897816, 2.7182818284590433]\n",
      "21 \t\t 300 \t ['(x-((((x*x)*x)/(c2/(x*x)))*(((c2/(x*x))-c1)*c2)))']  a.k.a  [c1*x**5 - c2*x**3 + x] c =  [2.7182818284590944, 3.1415926535902035]\n",
      "33 \t\t 1349 \t ['((((x*x)*((((x*x)*c1)-c2)*(x*x)))+(x*x))/x)']  a.k.a  [c1*x**5 - c2*x**3 + x] c =  [2.7182818284590446, 3.1415926535897873]\n",
      "36 \t\t 1454 \t ['((((x*x)*x)*(c2+(c1*(x*x))))+x)']  a.k.a  [c1*x**5 + c2*x**3 + x] c =  [2.7182818284590455, -3.141592653589795]\n",
      "38 \t\t 184 \t ['(x-(((x*x)*x)*(c1-((x*x)*c2))))']  a.k.a  [-c1*x**3 + c2*x**5 + x] c =  [3.141592653589796, 2.7182818284590455]\n",
      "45 \t\t 1935 \t ['((x+((c2*x)*((x*x)*(x*x))))+(x*(c1*(x*x))))']  a.k.a  [c1*x**3 + c2*x**5 + x] c =  [-3.141592653589768, 2.7182818284590446]\n",
      "58 \t\t 544 \t ['(x+(((c1+((c2*x)*((c2*x)/c2)))*x)*(((c2*x)/c2)*((c2*x)/c2))))']  a.k.a  [c1*x**3 + c2*x**5 + x] c =  [-3.141592653589802, 2.718281828459047]\n",
      "63 \t\t 390 \t ['(((x*((((c1*x)*x)-c2)*x))*x)+((((c1*x)*x)/c1)/x))']  a.k.a  [c1*x**5 - c2*x**3 + x] c =  [2.718281828459046, 3.1415926535897984]\n",
      "64 \t\t 99 \t ['((((c2+(c1*(x*x)))*x)*(x*x))+x)']  a.k.a  [c1*x**5 + c2*x**3 + x] c =  [2.7182818284590455, -3.141592653589793]\n",
      "71 \t\t 429 \t ['(x*(((((x*x)/(x*x))*((c1*(x*x))+c2))*(x*x))+((x*x)/(x*x))))']  a.k.a  [c1*x**5 + c2*x**3 + x] c =  [2.7182818284590455, -3.1415926535897585]\n",
      "76 \t\t 1890 \t ['(((c1/(c1*x))+((((c1*x)+(c1*x))*(x*x))+(c2*x)))*((((c1*x)+(c1*x))*(x*x))/((c1*x)+(c1*x))))']  a.k.a  [2*c1*x**5 + c2*x**3 + x] c =  [1.3591409142295234, -3.1415926535898038]\n",
      "83 \t\t 1841 \t ['((((x*x)*((x*x)*x))*c1)+(((c2*(x*x))*x)+x))']  a.k.a  [c1*x**5 + c2*x**3 + x] c =  [2.7182818284590455, -3.141592653589795]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=2000\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\",\"pdiv\"])() \n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,1233456))\n",
    "    g, best_chromosome, best_constant = run_experiment(dCGP, offsprings,max_gen,x,yt, screen_output=False)\n",
    "    res.append(g)\n",
    "    dCGP.set(best_chromosome)\n",
    "    if g < (max_gen-1):\n",
    "        print(i, \"\\t\\t\", res[i], \"\\t\", dCGP([\"x\",\"c1\",\"c2\"]), \" a.k.a \", dCGP.simplify([\"x\",\"c1\",\"c2\"]), \"c = \", best_constant)\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERT Expected run time = avg. number of dCGP evaluations needed:  49451.4666667\n"
     ]
    }
   ],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1))\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Problem P6: $\\frac{e x^2-1}{\\pi (x + 2)}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-2.1,1,10)\n",
    "x = gdual(x)\n",
    "yt = data_P6(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.5/site-packages/numpy/linalg/linalg.py:1776: RuntimeWarning: invalid value encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12 \t\t 5694 \t ['((c2+(c2+c2))*(((c2/x)-(c1*x))/((c2/x)+((c2/x)+c2))))']  a.k.a  [-3*c1*c2*x/(c2 + 2*c2/x) + 3*c2**2/(c2*x + 2*c2)] c =  [-0.2884186598106795, -0.10610329539427757]\n",
      "17 \t\t 1311 \t ['(((x*x)+((c2*(x*x))+c1))*(((x*x)-(x*x))-(c1/(c1+(c1+(c1*x))))))']  a.k.a  [-c1**2/(c1*x + 2*c1) - c1*c2*x**2/(c1*x + 2*c1) - c1*x**2/(c1*x + 2*c1)] c =  [0.3183098861802687, -1.8652559794314285]\n",
      "18 \t\t 8209 \t ['((c1+(x*(x*c2)))/((c2/c2)+(x+(c2/c2))))']  a.k.a  [c1/(x + 2) + c2*x**2/(x + 2)] c =  [-0.31830988618370615, 0.8652559794322552]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:47: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:48: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/lib/python3.5/site-packages/numpy/linalg/linalg.py:1776: RuntimeWarning: overflow encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "97 \t\t 7718 \t ['(((x*(c2+c1))*((x*(c2+c1))+(c2/x)))/(((c2+c1)+(c2+c1))+(x*(c2+c1))))']  a.k.a  [c1**2*x**2/(c1*x + 2*c1 + c2*x + 2*c2) + 2*c1*c2*x**2/(c1*x + 2*c1 + c2*x + 2*c2) + c1*c2/(c1*x + 2*c1 + c2*x + 2*c2) + c2**2*x**2/(c1*x + 2*c1 + c2*x + 2*c2) + c2**2/(c1*x + 2*c1 + c2*x + 2*c2)] c =  [1.1835658656164778, -0.31830988618398554]\n",
      "99 \t\t 4569 \t ['(((c1/(c2*x))-(x+x))/((((x+x)/x)+((x+x)+((x+x)/x)))/((c2*x)+(c2*x))))']  a.k.a  [2*c1/(2*x + 4) - 4*c2*x**2/(2*x + 4)] c =  [-0.3183098861838616, -0.43262798971613803]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=10000\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\",\"pdiv\"])() \n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,1233456))\n",
    "    g, best_chromosome, best_constant = run_experiment(dCGP, offsprings,max_gen,x,yt, screen_output=False)\n",
    "    res.append(g)\n",
    "    dCGP.set(best_chromosome)\n",
    "    if g < (max_gen-1):\n",
    "        print(i, \"\\t\\t\", res[i], \"\\t\", dCGP([\"x\",\"c1\",\"c2\"]), \" a.k.a \", dCGP.simplify([\"x\",\"c1\",\"c2\"]), \"c = \", best_constant)\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERT Expected run time = avg. number of dCGP evaluations needed:  781924.8\n"
     ]
    }
   ],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1))\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Problem P7:  $\\sin(e x) + \\cos(\\pi x)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1,1,10)\n",
    "x = gdual(x)\n",
    "yt = data_P7(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:47: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/lib/python3.5/site-packages/ipykernel/__main__.py:48: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/lib/python3.5/site-packages/numpy/linalg/linalg.py:1776: RuntimeWarning: invalid value encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36 \t\t 359 \t ['(sin((c2*x))+cos(((c2*x)+(x*c1))))']  a.k.a  [sin(c2*x) + cos(c1*x + c2*x)] c =  [0.4233108251053548, 2.718281828484438]\n",
      "51 \t\t 3855 \t ['(sin((x/c1))+cos(((x/c1)-(c2*x))))']  a.k.a  [sin(x/c1) + cos(c2*x - x/c1)] c =  [0.3678794411238845, -0.4233108248124803]\n",
      "59 \t\t 1602 \t ['(sin((x*c2))-cos((c1*(((x*c2)*c2)+x))))']  a.k.a  [sin(c2*x) - cos(c1*c2**2*x + c1*x)] c =  [-3.744870239868995, 2.718281828948403]\n",
      "61 \t\t 758 \t ['(cos((x-(c1*x)))+sin((x/c2)))']  a.k.a  [sin(x/c2) + cos(c1*x - x)] c =  [-2.1415926534019243, 0.36787944117144233]\n",
      "77 \t\t 2163 \t ['(sin(((c1*(c2*x))/c2))+cos((c2*x)))']  a.k.a  [sin(c1*x) + cos(c2*x)] c =  [5508691.188571924, 251192.3238030792]\n",
      "95 \t\t 296 \t ['(cos(((c2*x)/c1))+sin((x*c1)))']  a.k.a  [sin(c1*x) + cos(c2*x/c1)] c =  [2.718281828459221, -8.539734222674024]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=5000\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\",\"div\", \"sin\",\"cos\"])() \n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(inputs=3, outputs=1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = randint(0,1233456))\n",
    "    g, best_chromosome, best_constant = run_experiment(dCGP, offsprings,max_gen,x,yt, screen_output=False)\n",
    "    res.append(g)\n",
    "    dCGP.set(best_chromosome)\n",
    "    if g < (max_gen-1):\n",
    "        print(i, \"\\t\\t\", res[i], \"\\t\", dCGP([\"x\",\"c1\",\"c2\"]), \" a.k.a \", dCGP.simplify([\"x\",\"c1\",\"c2\"]), \"c = \", best_constant)\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERT Expected run time = avg. number of dCGP evaluations needed:  319292.666667\n"
     ]
    }
   ],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1))\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}