{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Weighted dCGP for a symbolic regression task "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyaudi import gdual_vdouble as gdual\n",
    "from dcgpy import expression_weighted_gdual_vdouble as expression\n",
    "from dcgpy import kernel_set_gdual_vdouble as kernel_set\n",
    "import pyaudi\n",
    "import numpy as np\n",
    "import math\n",
    "import re\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The ES-(1+$\\lambda$) algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, screen_output):\n",
    "    # The offsprings chromosome, fitness and weights\n",
    "    chromosome = [1] * offsprings\n",
    "    fitness = [1] * offsprings\n",
    "    weights = [1] * offsprings\n",
    "    # Init the best as the initial random dCGP\n",
    "    best_chromosome = dCGP.get()\n",
    "    best_weights = dCGP.get_weights()\n",
    "    best_fitness = sum(mse(dCGP, x, yt).constant_cf)\n",
    "    # Main loop over generations\n",
    "    for g in range(max_gen):\n",
    "        for i in range(offsprings):\n",
    "            dCGP.set(best_chromosome)\n",
    "            dCGP.set_weights(best_weights)\n",
    "            cumsum=0\n",
    "            dCGP.mutate_active(i)\n",
    "            newton(dCGP, mse, x, yt, newtonParams)\n",
    "            fitness[i] = sum(mse(dCGP, x, yt).constant_cf)\n",
    "            chromosome[i] = dCGP.get()\n",
    "            weights[i] = dCGP.get_weights()\n",
    "        for i in range(offsprings):\n",
    "            if fitness[i] <= best_fitness:\n",
    "                dCGP.set(chromosome[i])\n",
    "                dCGP.set_weights(weights[i])\n",
    "                if (fitness[i] != best_fitness) and screen_output:\n",
    "                    print(\"New best found: gen: \", g, \" value: \", fitness[i], dCGP.simplify([\"x\"],True))\n",
    "                best_chromosome = chromosome[i]\n",
    "                best_fitness = fitness[i]\n",
    "                best_weights = weights[i]\n",
    "\n",
    "        if best_fitness < 1e-14:\n",
    "            break\n",
    "    return g, best_chromosome, best_weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The test problems\n",
    "P1: $x^5 - \\pi x^3 + x$\n",
    "\n",
    "P2: $x^5 - \\pi x^3 + \\frac{2\\pi}x$\n",
    "\n",
    "P3: $\\frac{e x^5 + x^3}{x + 1}$\n",
    "\n",
    "P4: $\\sin(\\pi x) + \\frac 1x$\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)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following functions create the target values for a gridded input x for different test problems\n",
    "def data_P1(x):\n",
    "    return x**5 - np.pi*x**3 + x\n",
    "def data_P2(x):\n",
    "    return x**5 - np.pi*x**3 + 2*np.pi / x\n",
    "def data_P3(x):\n",
    "    return (np.e*x**5 + x**3)/(x + 1)\n",
    "def data_P4(x):\n",
    "    return pyaudi.sin(np.pi * x) + 1./x\n",
    "def data_P5(x):\n",
    "    return np.e * x**5 - np.pi*x**3 + np.sqrt(2) * x\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The error function"
   ]
  },
  {
   "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",
    "# Quadratic error of a dCGP expression. The error is computed over the input points xin (of type gdual, order 0 as\n",
    "# we are not interested in expanding the program w.r.t. these). The target values are contained in yt (of type gdual,\n",
    "# order 0 as we are not interested in expanding the program w.r.t. these)\n",
    "def mse(dCGP, xin, yt):\n",
    "    y = dCGP([xin])[0]\n",
    "    return (y-yt)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newton's method for minimizing the error function f w.r.t. the weights of the dCGP expression.\n",
    "# We take a specified amount of steps, each by choosing randomly 2 or 3 weights\n",
    "def newton(ex, f, x, yt, p):\n",
    "    n = ex.get_n()\n",
    "    r = ex.get_rows()\n",
    "    c = ex.get_cols()\n",
    "    a = ex.get_arity()[0]\n",
    "    v = np.zeros(r * c * a)\n",
    "    \n",
    "    # random initialization of weights\n",
    "    w=[]\n",
    "    for i in range(r*c):\n",
    "        for j in range(a):\n",
    "            w.append(gdual([np.random.normal(0,1)]))\n",
    "    ex.set_weights(w)\n",
    "    wi = ex.get_weights()\n",
    "        \n",
    "    # get active weights\n",
    "    an = ex.get_active_nodes()\n",
    "    is_active = [False] * (n + r * c) # bool vector of active nodes\n",
    "    for k in range(len(an)):\n",
    "        is_active[an[k]] = True\n",
    "    aw=[] # list of active weights\n",
    "    for k in range(len(an)):\n",
    "        if an[k] >= n:\n",
    "            for l in range(a):\n",
    "                aw.append([an[k], l]) # pair node/ingoing connection \n",
    "    if len(aw)<2:\n",
    "        return\n",
    "    \n",
    "    for i in range(p['steps']):\n",
    "        w = ex.get_weights() # initial weights\n",
    "        \n",
    "        # random choice of the weights w.r.t. which we'll minimize the error\n",
    "        num_vars = np.random.randint(2, min(3, len(aw)) + 1) # number of weights (2 or 3)\n",
    "        awidx = np.random.choice(len(aw), num_vars, replace = False) # indexes of chosen weights\n",
    "        ss = [] # symbols\n",
    "        for j in range(len(awidx)):\n",
    "            ss.append(\"w\" + str(aw[awidx[j]][0]) + \"_\" + str(aw[awidx[j]][1]))\n",
    "            idx = (aw[awidx[j]][0] - n) * a + aw[awidx[j]][1]\n",
    "            w[idx] = gdual(w[idx].constant_cf, ss[j], 2)\n",
    "        ex.set_weights(w)\n",
    "        \n",
    "        # compute the error\n",
    "        E = f(ex, x, yt)\n",
    "        Ei = sum(E.constant_cf)\n",
    "        \n",
    "        # get gradient and Hessian\n",
    "        dw = np.zeros(len(ss))\n",
    "        H = np.zeros((len(ss),len(ss)))\n",
    "        for k in range(len(ss)):\n",
    "            dw[k] = collapse_vectorized_coefficient(E.get_derivative({\"d\"+ss[k]: 1}), len(x.constant_cf))\n",
    "            H[k][k] = collapse_vectorized_coefficient(E.get_derivative({\"d\"+ss[k]: 2}), len(x.constant_cf))\n",
    "            for l in range(k):\n",
    "                H[k][l] = collapse_vectorized_coefficient(E.get_derivative({\"d\"+ss[k]: 1, \"d\"+ss[l]: 1}), len(x.constant_cf))\n",
    "                H[l][k] = H[k][l]\n",
    "        \n",
    "        det = np.linalg.det(H)\n",
    "        if det == 0: # if H is singular\n",
    "            continue\n",
    "        \n",
    "        # compute the updates\n",
    "        updates = - np.linalg.inv(H) @ dw\n",
    "        \n",
    "        # update the weights\n",
    "        for k in range(len(updates)):\n",
    "            idx = (aw[awidx[k]][0] - n) * a + aw[awidx[k]][1]\n",
    "            ex.set_weight(aw[awidx[k]][0], aw[awidx[k]][1], w[idx] + updates[k])\n",
    "        wfe = ex.get_weights()\n",
    "        for j in range(len(awidx)):\n",
    "            idx = (aw[awidx[j]][0] - n) * a + aw[awidx[j]][1]\n",
    "            wfe[idx] = gdual(wfe[idx].constant_cf)\n",
    "        ex.set_weights(wfe)\n",
    "        \n",
    "        # if error increased restore the initial weights\n",
    "        Ef = sum(f(ex, x, yt).constant_cf)\n",
    "        if not Ef < Ei:\n",
    "            for j in range(len(awidx)):\n",
    "                idx = (aw[awidx[j]][0] - n) * a + aw[awidx[j]][1]\n",
    "                w[idx] = gdual(w[idx].constant_cf)\n",
    "            ex.set_weights(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Problem P1:  $x^5 - \\pi x^3 + x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(1,3,10)\n",
    "x = gdual(x)\n",
    "yt = data_P1(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n",
      "1 \t\t 109 \t [0.999999999999997*x**5 - 3.14159265358976*x**3 + 0.999999999999932*x]\n",
      "4 \t\t 43 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999996*x]\n",
      "7 \t\t 29 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]\n",
      "8 \t\t 97 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "10 \t\t 36 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]\n",
      "12 \t\t 40 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999989*x]\n",
      "13 \t\t 79 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999998*x]\n",
      "15 \t\t 92 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999998*x]\n",
      "16 \t\t 133 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]\n",
      "18 \t\t 97 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "19 \t\t 46 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999991*x]\n",
      "24 \t\t 151 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999987*x]\n",
      "25 \t\t 33 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999998*x]\n",
      "30 \t\t 131 \t [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]\n",
      "32 \t\t 23 \t [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]\n",
      "33 \t\t 81 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "34 \t\t 16 \t [1.0*x**5 - 3.14159265358979*x**3 + 1.00000000000001*x]\n",
      "36 \t\t 48 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 1.0*x]\n",
      "39 \t\t 150 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]\n",
      "45 \t\t 50 \t [0.999999999999995*x**5 - 3.14159265358974*x**3 + 0.999999999999874*x]\n",
      "46 \t\t 189 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.0*x]\n",
      "49 \t\t 67 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]\n",
      "55 \t\t 29 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "56 \t\t 48 \t [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]\n",
      "60 \t\t 51 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999988*x]\n",
      "71 \t\t 90 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 1004.92559245124*x**3/(1004.92559245124*x**2 - 3157.0668586492) - 3157.06685864919*x/(1004.92559245124*x**2 - 3157.0668586492)]\n",
      "74 \t\t 43 \t [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]\n",
      "77 \t\t 142 \t [1.0*x**5 - 3.14159265358979*x**3 + 1.0*x]\n",
      "80 \t\t 192 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999988*x]\n",
      "82 \t\t 62 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999993*x]\n",
      "85 \t\t 52 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]\n",
      "86 \t\t 32 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 0.999999999999992*x]\n",
      "87 \t\t 83 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999992*x]\n",
      "89 \t\t 149 \t [0.999999999999998*x**5 - 3.14159265358978*x**3 + 0.999999999999981*x]\n",
      "91 \t\t 141 \t [1.0*x**5 - 3.14159265358979*x**3 + 0.999999999999999*x]\n",
      "95 \t\t 96 \t [1.0*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\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 = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\n",
    "    \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:  174288.8888888889\n"
     ]
    }
   ],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem P2: $x^5 - \\pi x^3 + \\frac{2\\pi}x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0.1,5,10) # we include points close to zero here to favour learning of 1/x\n",
    "x = gdual(x)\n",
    "yt = data_P2(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n",
      "1 \t\t 64 \t [1.00000000001522*x**5 - 1.27654166639185e-10*x**4 - 3.14159265332916*x**3 + 6.28318530717146/x]\n",
      "3 \t\t 160 \t [0.999999999999802*x**5 - 3.14159265358324*x**3 - 4.58664018735164e-11*x + 6.28318530718224/x]\n",
      "6 \t\t 43 \t [1.00000000000018*x**5 - 3.14159265361377*x**3 + 6.28318530717944/x]\n",
      "7 \t\t 35 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "8 \t\t 63 \t [1.0*x**5 - 3.14159265358978*x**3 + 6.28318530717959/x]\n",
      "11 \t\t 168 \t [1.0*x**5 - 3.14159265358979*x**3 - 3.01980662698043e-14*x + 6.28318530717959/x]\n",
      "13 \t\t 24 \t [0.999999999999994*x**5 - 3.14159265359058*x**3 + 2.78627727461556e-13*x + 6.28318530717957/x]\n",
      "14 \t\t 143 \t [0.99999999997825*x**5 - 3.14159265281984*x**3 - 6.06404304548391e-9*x + 6.28318530741037/x]\n",
      "16 \t\t 187 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "21 \t\t 82 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "27 \t\t 7 \t [1.0*x**5 - 3.14159265358984*x**3 + 1.84329364468004e-13*x**2 + 6.28318530717957/x]\n",
      "28 \t\t 67 \t [0.999999999996586*x**5 - 3.14159265349105*x**3 - 2.21245651961797e-9 + 6.28318530745686/x]\n",
      "34 \t\t 153 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "39 \t\t 115 \t [0.999999999999999*x**5 - 3.14159265358978*x**3 + 6.28318530717959/x]\n",
      "43 \t\t 84 \t [0.999999999999923*x**5 + 6.42597086653041e-13*x**4 - 3.14159265359109*x**3 + 6.28318530717962/x]\n",
      "50 \t\t 149 \t [0.99999999995493*x**5 - 3.14159265203297*x**3 - 1.17409352674938e-8*x + 6.28318530777391/x]\n",
      "51 \t\t 149 \t [1.00000000012198*x**5 - 3.14159265780051*x**3 + 3.16978514547372e-8*x + 6.28318530601517/x]\n",
      "52 \t\t 64 \t [1.0*x**5 - 3.1415926535898*x**3 + 6.28318530717959/x]\n",
      "54 \t\t 79 \t [-6.32238025396602e-16*x**7 + 1.00000000000002*x**5 - 3.14159265358992*x**3 + 6.28318530717958/x]\n",
      "55 \t\t 69 \t [0.999999999999999*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "57 \t\t 51 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717958/x]\n",
      "58 \t\t 180 \t [0.999999999997159*x**5 + 2.56552832733729e-11*x**4 - 3.14159265364781*x**3 + 6.28318530750492/x]\n",
      "62 \t\t 96 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "64 \t\t 130 \t [1.0*x**5 - 3.1415926535898*x**3 + 6.28318530717959/x]\n",
      "67 \t\t 154 \t [0.99999999999999*x**5 + 7.88553657836724e-14*x**4 - 3.14159265358995*x**3 - 6.63622336613371e-14*x**2 + 6.2831853071796/x]\n",
      "70 \t\t 26 \t [1.00000000000018*x**5 - 3.1415926535975*x**3 + 1.64577260147692e-11*x**2 + 6.28318530721673/x]\n",
      "75 \t\t 48 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717961/x]\n",
      "76 \t\t 138 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "78 \t\t 82 \t [1.0*x**5 - 3.14159265358982*x**3 + 6.28318530719116/x]\n",
      "84 \t\t 165 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.28318530717959/x]\n",
      "89 \t\t 136 \t [1.00000000000278*x**5 - 3.14159265388607*x**3 + 1.09744339194339e-9*x**2 + 6.28318530643488/x]\n",
      "90 \t\t 188 \t [1.00000000000001*x**5 - 3.91763911922868e-14*x**4 - 3.14159265358973*x**3 + 6.28318530717958/x]\n",
      "91 \t\t 108 \t [1.0*x**5 - 3.14159265358979*x**3 + 6.2831853071796/x]\n",
      "99 \t\t 64 \t [1.0*x**5 - 3.14159265358979*x**3 - 0.36315536377679/(x*(2.18996409662261e-16*x**3 - 0.0577979712554117))]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\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:  195352.94117647057\n"
     ]
    }
   ],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem P3: $\\frac{e x^5 + x^3}{x + 1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-0.9,1,10)\n",
    "x = gdual(x)\n",
    "yt = data_P3(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n",
      "8 \t\t 186 \t [0.899166850493274*x**5/(0.330784998479351*x + 0.330784998479351) + 0.330784998479352*x**3/(0.330784998479351*x + 0.330784998479351)]\n",
      "14 \t\t 161 \t [1.20883912632723*x**5/(0.444707062259436*x + 0.444707062259436) + 0.444707062259438*x**3/(0.444707062259436*x + 0.444707062259436)]\n",
      "73 \t\t 125 \t [2.71828182845905*x**5/(x + 1) + 1.0*x**3/(x + 1)]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem P4: $\\sin(\\pi x) + \\frac 1x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1,1,10)\n",
    "x = gdual(x)\n",
    "yt = data_P4(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dario/.local/lib/python3.7/site-packages/numpy/linalg/linalg.py:2116: RuntimeWarning: invalid value encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 \t\t 1 \t [0.999999999999999*sin(3.14159265358979*x) + 1.0/x]\n",
      "6 \t\t 132 \t [1.0*sin(3.14159265358979*x) - 3.6500836381227e-16 + 1.0/x]\n",
      "7 \t\t 89 \t [0.999999999999999*sin(3.14159265358979*x) + 1.0/x]\n",
      "8 \t\t 133 \t [0.999999999999997*sin(3.14159265358979*x) + 1.0/x]\n",
      "10 \t\t 120 \t [0.999999999999999*sin(3.14159265358979*x) + 4.68530917602041e-17 + 1.0/x]\n",
      "12 \t\t 55 \t [-7.89209275531109e-11*x + 1.0000000000362*sin(3.14159265367275*x + 37.6991118429053) - 9.47051130607999e-10 + 1.00000000000142/x]\n",
      "14 \t\t 66 \t [1.00000000000032*sin(3.14159265358985*x + 2.53915412165946e-16) + 0.999999999999973/x]\n",
      "16 \t\t 12 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "19 \t\t 2 \t [0.999999949216629*sin(3.14159264435007*x) + 9.36443001715427e-18 + 1.00000000433431/x]\n",
      "20 \t\t 123 \t [1.0*sin(3.14159265358979*x) + (7.86221303890169e-17*sin(3.14159265358979*x) + 1.0)/x]\n",
      "21 \t\t 17 \t [9.62453199799219*x*sin(3.1415927053321*x)/(9.62453191537768*x - 1.06372808316533e-16) + 9.62453191084485/(9.62453191537768*x - 1.06372808316533e-16)]\n",
      "22 \t\t 13 \t [0.999999999998995*sin(3.14159265358961*x) + 1.00000000000009/x]\n",
      "24 \t\t 69 \t [0.999999999999953*sin(x*(9.25302298900841e-17*x + 3.14159265358979)) + 1.0/x]\n",
      "26 \t\t 105 \t [1.0*sin(x*(1.77990517211438e-16*x + 3.14159265358979)) + 5.6656141275367e-17 + 1.0/x]\n",
      "27 \t\t 115 \t [-1.24427033889047e-16*x*sin(3.14159265358979*x) + 1.0*sin(3.14159265358979*x) + 1.0/x]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dario/.local/lib/python3.7/site-packages/numpy/linalg/linalg.py:2116: RuntimeWarning: overflow encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "34 \t\t 64 \t [0.999999999999999*sin(3.14159265358979*x) - 1.29978496930744e-16 + 1.0/x]\n",
      "36 \t\t 46 \t [1.76206042739191e-16*x*sin(3.14159265358979*x)/(x*sin(3.14159265358979*x) + 1) + 1.0*sin(3.14159265358979*x) + 1.76206042739191e-16/(x*sin(3.14159265358979*x) + 1) + 1.0/x]\n",
      "38 \t\t 95 \t [1.0*sin(x*(1.3303626169174e-16*x + 3.14159265358979)) + 1.69425620306735e-16 + 1.0/x]\n",
      "40 \t\t 143 \t [8.12378036198924e-33*x + 1.0*sin(3.14159265358981*x) - 1.80264032596514e-16 + 1.0/x]\n",
      "42 \t\t 16 \t [1.00000000004875*sin(3.14159265359866*x) + 0.542336605776001*sin(2.35272399768743e-12*x**2) + 2541.81495243938/(2541.81495245225*x - 2.15281241329688) - 2.15281241328598/(x*(2541.81495245225*x - 2.15281241329688)) + 972.092685040695*sin(2.35272399768743e-12*x**2)/(x**2*(2541.81495245225*x - 2.15281241329688))]\n",
      "44 \t\t 24 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "45 \t\t 162 \t [6.02619053589549e-16*sin(3.14159265355456*x)**2 + 0.999999999994157*sin(3.14159265355456*x) - 4.51719913513711e-16 + 1.00000000000032/x]\n",
      "47 \t\t 42 \t [-1.87473626494473e-29*x**2 - 1.8747362649445e-29*x/sin(3.14159265359048*x) + 1.00000000000011*sin(3.14159265359048*x) + 0.999999999999994/x]\n",
      "48 \t\t 28 \t [1.00000000000002*sin(3.14159265358992*x) + 0.999999999999999/x]\n",
      "49 \t\t 18 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "50 \t\t 41 \t [1.0*sin(3.14159265358979*x) - 4.08683517342902e-17 + 1.0/x]\n",
      "52 \t\t 61 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "53 \t\t 49 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "55 \t\t 29 \t [1.0*sin(3.14159265358979*x) + (1.0 - 6.66045855255763e-18*sin(3.14159265358979*x))/x]\n",
      "56 \t\t 15 \t [6.58379173825623*x + 5.64744625319335*sin(6.11418858333298*x) + 5.2879457406247*sin(16.7761307129903*x)]\n",
      "57 \t\t 32 \t [-0.974501474950861*sin(0.863186354404554*x) + 9.07750333059453*sin(1.6599296272374*x + 2.06966097571612*sin(4.75780034524026*x)) - 5.38968705261534*sin(3.48198655314495*x - 1.09415284338664*sin(4.75780034524026*x))]\n",
      "59 \t\t 76 \t [1.0*sin(3.14159265358979*x) + 4.45769125393131e-16 + 1.0/x]\n",
      "60 \t\t 7 \t [1.00000000000001*sin(3.1415926535899*x) + 4.70822994414938e-17 + 1.0/x]\n",
      "63 \t\t 76 \t [5.19135160727313*sin(3.59576183848133*x) + 8.66296848504959*sin(12.5211187903134*x - 0.987545240550494*sin(3.59576183848133*x))]\n",
      "68 \t\t 150 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "72 \t\t 58 \t [1.0*sin(3.14159265358979*x) + 1.25427228495141e-16 + 1.0/x]\n",
      "73 \t\t 89 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "75 \t\t 49 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "79 \t\t 76 \t [1.0*sin(3.14159265358979*x) + 3.94677549026297e-17 + 1.0/x]\n",
      "83 \t\t 76 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "85 \t\t 12 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n",
      "87 \t\t 94 \t [1.000000000001*sin(3.14159265359084*x) + (0.999999999999751 - 1.35525271560688e-20/sin(3.14159265359084*x))/x]\n",
      "92 \t\t 19 \t [0.999999980945392*sin(3.14159265012292*x) + 1.00000000162629/x]\n",
      "93 \t\t 153 \t [-0.999999999999999*sin(3.14159265358979*x + 3.14159265358979) + 1.0/x]\n",
      "97 \t\t 48 \t [1.0*sin(3.14159265358979*x) + 1.0/x]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\", \"sin\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Problem P5:  $ex^5 - \\pi x^3 + x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(1,3,10)\n",
    "x = gdual(x)\n",
    "yt = data_P5(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n",
      "2 \t\t 130 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000003*x]\n",
      "5 \t\t 144 \t [2.71828182845904*x**5 - 3.14159265358974*x**3 + 0.999999999999915*x]\n",
      "6 \t\t 129 \t [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "7 \t\t 48 \t [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999982*x]\n",
      "9 \t\t 25 \t [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000002*x]\n",
      "10 \t\t 6 \t [2.71828182845904*x**5 - 3.14159265358976*x**3 + 0.999999999999928*x]\n",
      "11 \t\t 22 \t [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999977*x]\n",
      "13 \t\t 113 \t [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999985*x]\n",
      "15 \t\t 9 \t [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999969*x]\n",
      "16 \t\t 46 \t [2.71828182845796*x**5 - 3.1415926535816*x**3 + 0.999999999987821*x]\n",
      "17 \t\t 77 \t [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "18 \t\t 189 \t [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "20 \t\t 181 \t [2.71828182765082*x**5 + 4.81282166780353e-9*x**4 - 3.14159266134502*x**3 + 1.0000000050856*x]\n",
      "21 \t\t 163 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000003*x]\n",
      "22 \t\t 16 \t [2.71828182844251*x**5 - 3.14159265338736*x**3 + 0.999999999457422*x]\n",
      "24 \t\t 37 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]\n",
      "27 \t\t 119 \t [2.71828182845942*x**5 - 3.14159265359449*x**3 + 1.00000000000628*x]\n",
      "28 \t\t 4 \t [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "29 \t\t 42 \t [2.71828181857779*x**5 - 3.14159237432049*x**3 - 7.49766685257696e-7*x**2 + 1.00000053196694*x]\n",
      "42 \t\t 33 \t [2.71828182845902*x**5 - 3.14159265358958*x**3 + 0.999999999999598*x]\n",
      "44 \t\t 8 \t [2.71828182845905*x**5 - 3.1415926535898*x**3 + 1.00000000000001*x]\n",
      "45 \t\t 189 \t [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000005*x]\n",
      "49 \t\t 5 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]\n",
      "53 \t\t 21 \t [2.71828182845904*x**5 - 3.14159265358977*x**3 + 0.99999999999993*x]\n",
      "55 \t\t 27 \t [2.71828182845904*x**5 - 3.14159265358977*x**3 + 0.999999999999773*x]\n",
      "63 \t\t 83 \t [2.71828182845904*x**5 - 3.14159265358979*x**3 + 0.99999999999998*x]\n",
      "65 \t\t 29 \t [2.71828182845904*x**5 - 3.14159265358979*x**3 + 0.999999999999992*x]\n",
      "68 \t\t 149 \t [2.71828182845905*x**5 + x**3*(1.51490803020543 - 2.20219730313442*x**2)/(0.700981172914968*x**2 - 0.458140436777564) + x*(0.625363773202351*x**2 - 0.458140436777565)/(0.700981172914968*x**2 - 0.458140436777564)]\n",
      "69 \t\t 66 \t [2.71828182845905*x**5 - 3.14159265358978*x**3 + 0.999999999999967*x]\n",
      "71 \t\t 65 \t [2.71828182845905*x**5 - 3.14159265358988*x**3 + 1.00000000000015*x]\n",
      "73 \t\t 181 \t [2.71828182845904*x**5 - 3.14159265358973*x**3 + 0.999999999999894*x]\n",
      "74 \t\t 150 \t [2.71828182845904*x**5 - 3.14159265358979*x**3 + 0.999999999999995*x]\n",
      "78 \t\t 12 \t [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.999999999999977*x]\n",
      "80 \t\t 118 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]\n",
      "81 \t\t 137 \t [2.71828182845906*x**5 - 3.14159265359007*x**3 + 1.00000000000047*x]\n",
      "84 \t\t 50 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000002*x]\n",
      "85 \t\t 33 \t [2.71828182845904*x**5 - 3.1415926535897*x**3 + 0.999999999999851*x]\n",
      "87 \t\t 93 \t [2.71828182845905*x**5 - 3.14159265358981*x**3 + 1.00000000000003*x]\n",
      "88 \t\t 60 \t [2.71828182845904*x**5 - 3.14159265358978*x**3 + 0.99999999999997*x]\n",
      "89 \t\t 88 \t [2.71828182845905*x**5 - 3.14159265358979*x**3 + 0.999999999999997*x]\n",
      "91 \t\t 15 \t [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000004*x]\n",
      "92 \t\t 162 \t [2.71828182845904*x**5 - 3.14159265358977*x**3 + 0.999999999999982*x]\n",
      "95 \t\t 61 \t [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000003*x]\n",
      "98 \t\t 68 \t [2.71828182845905*x**5 - 3.14159265358982*x**3 + 1.00000000000006*x]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERT Expected run time = avg. number of dCGP evaluations needed:  132245.45454545456\n"
     ]
    }
   ],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\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)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-2.1,1,10)\n",
    "x = gdual(x)\n",
    "yt = data_P6(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n",
      "7 \t\t 121 \t [0.658810478707329*x**2/(0.76140528857091*x + 1.52281057714182) - 0.242362830744741/(0.76140528857091*x + 1.52281057714182)]\n",
      "12 \t\t 70 \t [x*(0.710716298622762*x - 6.5421001949062e-12)/(0.82139426426137*x + 1.64278852852271) - 0.261457914778349/(0.82139426426137*x + 1.64278852852271)]\n",
      "18 \t\t 184 \t [x*(0.000123063420861819*x**3 + 0.000387617697492767*x**2 + 0.00026373499528672*x - 3.84934325911654e-5)/((0.0332304392124186*x + 0.0664608784094327)*(0.00428004451413307*x**2 + 0.0134810245653551*x + 0.00984187107586485)) - 0.0235684911144266/(0.128798914957873*x + 0.148085179001607)]\n",
      "37 \t\t 56 \t [-0.56551064554709*x/(1.74078576306486*x + 3.48157152619754) - 1.40565254304111 + (1.50622528661617*x**2 + 3.01245057329103*x + 4.33977055108826)/(1.74078576306486*x + 3.48157152619754)]\n",
      "42 \t\t 154 \t [0.0804769112118696*x**2/(0.093009367314253*x + 0.186018734628506) - 0.0296058011238267/(0.093009367314253*x + 0.186018734628506)]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dario/miniconda3/envs/dcgpy/lib/python3.7/site-packages/ipykernel_launcher.py:63: RuntimeWarning: invalid value encountered in matmul\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\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)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1,1,10)\n",
    "x = gdual(x)\n",
    "yt = data_P7(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restart: \t gen: \t expression:\n",
      "2 \t\t 10 \t [1.0*sin(2.71828182845905*x) + 0.999999999999999*cos(3.14159265358979*x)]\n",
      "5 \t\t 177 \t [0.999999999990239*sin(2.71828182845446*x) + 0.99999999827576*cos(3.1415926695047*x)]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dario/miniconda3/envs/dcgpy/lib/python3.7/site-packages/ipykernel_launcher.py:63: RuntimeWarning: invalid value encountered in matmul\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 \t\t 52 \t [0.999999996606054*sin(2.71828182686597*x) + 1.0*cos(3.14159265358979*x)]\n",
      "12 \t\t 166 \t [1.00000021604097*cos(3.14159255591115*x) + 1.00000000001092*cos(2.71828182846416*x - 7.85398185936523)]\n",
      "14 \t\t 84 \t [1.00000000002564*sin(2.71828182859488*x) + 1.0*cos(3.14159265358979*x)]\n",
      "15 \t\t 137 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "17 \t\t 25 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "19 \t\t 145 \t [-0.999999999999992*sin(x*(2.59471453823685e-12*sin(2.29774027604273*x**2) - 2.718281828461)) + 1.0*cos(3.14159265358979*x)]\n",
      "21 \t\t 51 \t [1.0*sin(2.71828182845904*x) + 1.00000000000002*cos(3.14159265358978*x)]\n",
      "26 \t\t 8 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "29 \t\t 116 \t [-3.5489631288284e-16*sin(2.71828182845905*x)**2 + 1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "31 \t\t 82 \t [0.999999999612255*sin(2.71828182743606*x + 8.09654622736853e-11*cos(3.14159265362908*x)) + 1.00000000000211*cos(3.14159265362908*x)]\n",
      "32 \t\t 132 \t [1.00000000003824*sin(2.71828182866167*x) + 1.0*cos(3.14159265358979*x)]\n",
      "44 \t\t 188 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "48 \t\t 122 \t [1.00000000001473*sin(25.5560520538422*x) + 0.999999991979214*cos(3.14159272762242*x)]\n",
      "49 \t\t 88 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "50 \t\t 55 \t [0.999999999999999*sin(2.71828182845904*x - 8.91009907439534e-10*cos(31.4159265256615*x)) - 1.00000002272797*cos(31.4159265256615*x)]\n",
      "52 \t\t 86 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "54 \t\t 106 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "59 \t\t 81 \t [1.0*sin(2.71828182845905*x) + 1.00000000226596*cos(3.14159263267476*x)]\n",
      "63 \t\t 110 \t [1.08304336385657e-7*x*sin(2.71828183307158*x)*cos(3.14159262986457*x) + 1.00000000982667*sin(2.71828183307158*x) + 0.999999952852321*cos(3.14159262986457*x)]\n",
      "67 \t\t 73 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "82 \t\t 53 \t [0.999999999999999*sin(2.71828182845904*x) - 1.0*cos(31.4159265358979*x)]\n",
      "87 \t\t 77 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "93 \t\t 86 \t [1.0*sin(2.71828182845905*x) + 1.0*cos(3.14159265358979*x)]\n",
      "98 \t\t 175 \t [2.68964745528915e-8*x**2 + 1.0*sin(2.71828182845905*x) + 1.00000000821958*cos(3.14159270860771*x)]\n"
     ]
    }
   ],
   "source": [
    "# We run nexp experiments and accumulate statistic for the ERT\n",
    "nexp = 100\n",
    "offsprings = 4\n",
    "max_gen=200\n",
    "res = []\n",
    "kernels = kernel_set([\"sum\", \"mul\", \"diff\", \"div\", \"sin\", \"cos\"])()\n",
    "newtonParams = {\n",
    "    'steps': 100,\n",
    "}\n",
    "print(\"restart: \\t gen: \\t expression:\")\n",
    "for i in range(nexp):\n",
    "    dCGP = expression(1, 1, rows=1, cols=15, levels_back=16, arity=2, kernels=kernels, seed = np.random.randint(1233456))\n",
    "    for j in range(dCGP.get_n(), dCGP.get_n() + dCGP.get_rows() * dCGP.get_cols()):\n",
    "        for k in range(dCGP.get_arity()[0]):\n",
    "            dCGP.set_weight(j, k, gdual([np.random.normal(0,1)]))\n",
    "    g, best_chromosome, best_weights = run_experiment(dCGP, offsprings, max_gen, x, yt, newtonParams, 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.simplify([\"x\"],True))\n",
    "res = np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_gen = sum(res) / sum(res<(max_gen-1)) * newtonParams['steps']\n",
    "print(\"ERT Expected run time = avg. number of dCGP evaluations needed: \", mean_gen * offsprings)"
   ]
  }
 ],
 "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
}