Cross_Entropy_Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importance Sampling \n",
    "\n",
    "\n",
    "The expectation value of a function $H({\\bf x})$ is defined as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    " \\langle H({\\bf x})\\rangle_{f} = \\int d{\\bf x} H({\\bf x}) f({\\bf x; u}),\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $f({\\bf x})$ is the probability distribution. The unbiased estimator of this expectation value is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle H({\\bf x})\\rangle_{f} \\approx \\frac{1}{N}\\sum_{i=1} H({\\bf X_i}),\n",
    "\\end{align}\n",
    "$$\n",
    "where the sample points $X_i$ are drawn from the distribution $f$. This method may be a computationally inneficient estimate of the expectation value of $H$. We would like to modify our previous estimator to the form \n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle H({\\bf x})\\rangle_{g} \\approx \\frac{1}{N}\\sum_{i=1} w(X_i)H({\\bf X_i}),\n",
    "\\end{align}\n",
    "$$\n",
    "where $w(X_i)$ is an unknown weight factor and $X_i$ are drawn from a new probability distribution $g({\\bf x})$ that samples $H(x)$ more effectively. \n",
    "\n",
    "\n",
    "We would like to relate how the weight functions $w(x)$ relate to the new distribution $g(x)$. This relation is very simple\n",
    "$$\n",
    "\\begin{align}\n",
    " \\langle H({\\bf x})\\rangle_{f} &= \\int d{\\bf x} H({\\bf x}) f({\\bf x; u}),\\\\\n",
    "  &= \\int d{\\bf x} H({\\bf x}) \\frac{f({\\bf x; u})}{g({\\bf x})} g({\\bf x}) \\\\\n",
    "  &=  \\int d{\\bf x} H({\\bf x}) w({\\bf x; u}) g({\\bf x}) \\\\\n",
    "  &\\equiv \\langle H({\\bf x})\\rangle_{g}.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "Now we would like to choose the best distribution function $g^*({\\bf x})$ to sample $H({\\bf x})$. This can be done by choosing the $g({\\bf x})$ function that minimizes the variance of the estimate $\\langle H({\\bf x})\\rangle_g$ \n",
    "\n",
    "The variance of this estimate is given by the integral\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sigma^2 = \\int d{\\bf x} \\left( w({\\bf x})H({\\bf x})-\\mu \\right)^2 g({\\bf x}),\n",
    "\\end{align}\n",
    "$$\n",
    "where $\\mu$ is the expectation value of the function $H({\\bf x})$. Minimizing the variance of this expression with respect to $w({\\bf x})$, we get\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{d\\sigma^2}{dw} = 0 &= \\int d{\\bf x} \\left( w({\\bf x})H({\\bf x})-\\mu \\right)H(x)\\delta({\\bf x-y})  g({\\bf x}) \\\\\n",
    "&\\Rightarrow \\\\\n",
    "w({\\bf y}) &= \\frac{\\mu}{H({\\bf y})}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Re-writing $w({\\bf y})$ in terms of $f$ and $g$ gives,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "g^*({\\bf y}) = \\frac{H({\\bf y})f({\\bf y; u})}{\\mu}.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This last expression for $g^*$ represents the optimal choice of the importance sampling function $g$. If we want to find the best function $g({\\bf y})$ to carry out importance sampling, we must choose $g=g^*$ given above.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cross Entropy Estimation\n",
    "\n",
    "\n",
    "Let us return to the problem of estimating the expectation value\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle H({\\bf x})\\rangle_{f} = \\int d{\\bf x} H({\\bf x}) f({\\bf x; u}).\n",
    "\\end{align}\n",
    "$$\n",
    "The best importance sampling function to estimate this integral is\n",
    "$$\n",
    "\\begin{align}\n",
    "g^*({\\bf x}) = \\frac{H({\\bf x})f({\\bf x;u})}{\\mu}.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Instead of using this function directly with the unknown parameter $\\mu$, we approximate the function $g^*(x)$ using our parametrized function $f({\\bf x;v})$. This approximation can be carried out using the Kullback-Leibler divergence\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "D(g^*,f({\\bf x;v})) = \\int dx \\ g^*({\\bf x}) {\\rm ln}(g^*({\\bf x})) - \\int dx \\ g^*({\\bf x}) {\\rm ln}(f({\\bf x;v})).\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The above expression vanishes when the two functions $g^*$ and $f$ are equal. To minimize the above expression in terms of the parameter ${\\bf v}$, we must maximize the expression\n",
    "$$\n",
    "\\begin{align}\n",
    "D({\\bf v}) &= \\int dx \\ g^*({\\bf x}) {\\rm ln}(f({\\bf x; v})),\\\\\n",
    " &= \\frac{1}{\\mu}\\int dx \\ H({\\bf x}) f({\\bf x; u}) {\\rm ln}(f({\\bf x; v})),\\\\\n",
    "  &= \\frac{1}{\\mu}\\int dx \\ H({\\bf x}) \\frac{f({\\bf x; u})}{f({\\bf x,w})} {\\rm ln}(f({\\bf x; v})) f({\\bf x,w}),\\\\\n",
    "  &= \\frac{1}{\\mu}\\int dx \\ H({\\bf x}) w({\\bf x;u,w}) {\\rm ln}(f({\\bf x; v})) f({\\bf x,w}).\n",
    "\\end{align}\n",
    "$$\n",
    "Therefore, The optimal choice of ${\\bf v}$ is given by solving the expression\n",
    "$$\n",
    "\\begin{align}\n",
    "\\nabla_{\\bf v}D({\\bf v}) &= 0 \\\\\n",
    "\\int dx \\ H({\\bf x}) w({\\bf x;u,w}) {\\bf \\nabla}_{\\bf v}{\\rm ln}(f({\\bf x; v})) f({\\bf x,w})&=0.\n",
    "\\end{align}\n",
    "$$\n",
    "When sampling this integral, this last expression will be\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{1}{N}\\sum_i H({\\bf X_i})w({\\bf X_i;u,w})  \\nabla_{\\bf v}{\\rm ln}(f({\\bf X_i;v})) = 0\n",
    "\\end{align}\n",
    "$$\n",
    "where \n",
    "$$\n",
    "X_i \\sim  f({\\bf x,w}).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Cross Entropy for Optimization\n",
    "\n",
    "To optimize the function $H({\\bf x})$ using cross entropy optimization we will consider the integral that represents that the function $H({\\bf x})$ is greater or equal to the maximum of the function $H({\\bf x^*})=\\gamma^*$.\n",
    "\n",
    "$$\n",
    "P(H >= \\gamma^*) = \\int d{\\bf x} \\ I_{ \\lbrace H({\\bf x})  >= \\gamma^* \\rbrace }(x) f({\\bf x; u})\n",
    "$$\n",
    "\n",
    "Note that we dont know ahead of time what $x^*$ and $\\gamma^*$ are yet. These will be updated iteratively.\n",
    "\n",
    "\n",
    "\n",
    "Using the ideas of importance sampling and minimizing the KL-divergence from the previous section, \n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{1}{N}\\sum_i H({\\bf X_i})\\nabla_{\\bf v}{\\rm ln}(f({\\bf X_i;v})) = 0\n",
    "\\end{align}\n",
    "$$\n",
    "where $X_i \\sim  f({\\bf x,u})$. The algorithm to do this iteratively is below\n",
    "\n",
    "\n",
    "#### Algorithm\n",
    "\n",
    "* Choose $\\rho$  \n",
    "* Choose $N$ \n",
    "* Initilize ${\\bf u}$\n",
    "\n",
    "\n",
    "\n",
    "1. Sample $N$ points,  $X_i \\sim f({\\bf x,u})$, compute $H(X_i)$\n",
    "2. Order $(X_i,H(X_i))$ in decending order according to the values $H(X_i)$  \n",
    "3. Using the $N_\\rho = N\\rho$ elite samples determine the solution ${\\bf v}$ to \n",
    "\n",
    "\n",
    "$$  \n",
    "\\frac{1}{N_\\rho}\\sum\\limits_{i=1}^{N_\\rho} H({\\bf X_i})\\nabla_{\\bf v}{\\rm ln}(f({\\bf X_i;v})) = 0  \n",
    "$$\n",
    "\n",
    "let this new solution replace the previous result, ${\\bf v} \\rightarrow {\\bf u}$\n",
    "\n",
    "4. Go to step 1, repeat until convergence criteria is met\n",
    "\n",
    "\n",
    "#### References \n",
    "\n",
    "[ [1] ](http://web.mit.edu/6.454/www/www_fall_2003/gew/CEtutorial.pdf) **A Tutorial on the Cross-Entropy Method** Pieter-Tjerk de Boer   \n",
    "\n",
    "[[2]](https://link.springer.com/article/10.1007/s11009-006-9753-0)  D. P. Kroese, S. Porotsky and R. Y. Rubinstein, **The Cross-Entropy Method for Continuous\n",
    "Multi-Extremal Optimization**, Methodol. Comput. Appl. Probab. 8 (2006).\n",
    "\n",
    "[[3]](https://www.sciencedirect.com/science/article/pii/S0968090X1100088X?via%3Dihub) M. Maher, R. Liu, D. Ngoduy, **Signal optimization using the cross entropy method**, Transportation Research C. 27, (2013).  \n",
    "\n",
    "\n",
    "[[4]](https://www.sciencedirect.com/science/article/pii/B9780444538598000035) Zdravko I. Botev, Dirk P. Kroese, Reuven Y. Rubinstein, Pierre L’Ecuyer, **Chapter 3 - The Cross-Entropy Method for Optimization**, Handbook of Statistics, Elsevier, 31, 35-59, (2013)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let us now implement the above algorithm as a simple example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simple Implementation\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5BcZ3km8Oft2/R0z/RoZjS6WHdj2Y68ZGOY9WIgwC432yQ22SyJXZUsqVDxUuAUu8lu1hRbXgr+2AUqSVV2vUlMwQZYguNNCChExAFiAgQbNDbGWLKFZVmyJOsymlvfpm/T7/5xzulptbqnz+k+lz6t51cl3Jcz53wctZ5+5zvf+T5RVRARUfhFgm4AERG5g4FORDQkGOhEREOCgU5ENCQY6EREQyIW1IE3b96se/fuDerwRESh9OSTT15S1Zl27wUW6Hv37sXc3FxQhyciCiUROdXpPXa5EBENCQY6EdGQYKATEQ0JBjoR0ZBgoBMRDYmugS4inxWRiyLybIf3RUT+SESOi8gzIvIa95tJRETd2KnQ/wzAbRu8fzuA/eafewH8cf/NIiIip7oGuqp+B8DiBpvcBeDzangCwCYR2e5WA4locK1W1vD5x09iqVAJuikEd/rQdwA43fT8jPnaFUTkXhGZE5G5+fl5Fw5NREH60++8iAe+egQf/ZsjQTeF4PNFUVV9SFVnVXV2ZqbtnatEFCLff3EBAPDtY/PgYjnBcyPQzwLY1fR8p/kaEQ2x2lodPzmzgtF4FCurVZxaKAbdpKueG4F+EMC/M0e7vA7Aiqqec2G/RDTATlwqYLW6hnffbPSwnlpkoAfNzrDFLwF4HMANInJGRN4nIu8XkfebmxwCcALAcQCfBvABz1pLRAPjzJIR4D+/fzMA4GUGeuC6zraoqvd0eV8BfNC1FhFRKJxdLgEAbt69CYloBGcY6IHjnaJE1JNXllcRjwq2jiexc2qUFfoAYKATUU9eWV7F1kwSkYhgx6ZRvLK8GnSTrnoMdCLqybnlEq6ZGAUAzIyN4FKeNxcFjYFORD25VChjZnwEADA9lsClfJlj0QPGQCeiniwVKphMxwEAm8dGUK7VUaisBdyqqxsDnYgcW6srllermEolABiBDgCXcuUgm3XVY6ATkWPZ1SpUgcm0EejTY8Z/L+UZ6EFioBORY4tF4wLoZGuFzgujgWKgE5Fj1nS5rRX6IqfRDRQDnYgcWypWAaDRhz4xalwczZaqgbWJGOhE1AOrQt+UMoJ8NB5FLCJYWWWgB4mBTkSOWX3oU2aXi4hgYjTOQA8YA52IHFsqVpCIRZBKRBuvTYzGkWWgB4qBTkSOLRUqmEzFISKN1zKs0APHQCcix5aK1caQRUuGFXrgGOhE5Fh2tYqMObLFMjEaR7ZUC6hFBDDQiagHuVINmeTl6+NMjMbY5RIwBjoROZYv1zA2cnmgZ5JGHzpnXAwOA52IHMuVqhhPXtnlslZXFDnjYmAY6ETkiKoiV6ph/IouFyPg2e0SHAY6ETlSqtZRq2vbCh1goAeJgU5EjuTKRmCPtVTo1qgXDl0MDgOdiBzJmUMTW0e5ZJLWBF0cuhgUBjoROWIFemsfenrEmAagUGagB4WBTkSO5Mwpclv70K1hjHkGemAY6ETkSN6s0FvHoafN56zQg8NAJyJHOnW5pBJRiDDQg8RAJyJHsh26XEQEY4kYcgz0wDDQiciRXIcuF8DodmGFHhwGOhE5ki/XkE5EEY3IFe+lR6IolHnrf1AY6ETkSLt5XCxjIzGOcgmQrUAXkdtE5JiIHBeR+9u8v1tEHhORH4nIMyJyh/tNJaJBUKisNcact2KXS7C6BrqIRAE8COB2AAcA3CMiB1o2+68AHlHVmwHcDeB/u91QIhoMxXKtMUSxVZoVeqDsVOi3ADiuqidUtQLgYQB3tWyjADLm4wkAr7jXRCJn/ue3XsCbPvkYTi8Wg27KUCpW1i5bHLrZ2EgMhQoDPSh2An0HgNNNz8+YrzX7KIBfE5EzAA4B+O12OxKRe0VkTkTm5ufne2gu0caKlRp+/xs/xcuLRXz+8ZNBN2coGYHevkIfG4nxomiA3Looeg+AP1PVnQDuAPAFEbli36r6kKrOqurszMyMS4cmWvejl5cbj7//4kKALRlehUqtY4WeHok17iQl/9kJ9LMAdjU932m+1ux9AB4BAFV9HEASwGY3GkjkxOGTixAB7rllF164mEdtrR50k4ZOsbyGdMcKPYrKWh2VGs97EOwE+mEA+0Vkn4gkYFz0PNiyzcsA3goAIvIzMAKdfSrkuxcu5LFnKoXZPVOo1Oo4uVAIuklDp1CpIbXBKBeAt/8HpWugq2oNwH0AHgXwHIzRLEdE5GMicqe52e8C+C0R+TGALwH4DeVKsRSAU4sF7J5O48bt4wCA58/nAm7RcFE11gztVKGnOeNioNr/rbRQ1UMwLnY2v/ZA0+OjAN7gbtOInFFVnFoo4uZdk9g9lQIAnFlaDbhVw6Vcq2Otrh0rdGs6AI50CQbvFKWhsVysIleqYc90CuPJOMaTMZxbZqC7qVgxRrB0q9DZ5RIMBjoNDasa3zlpVOc7No3i7HIpyCYNHSuoNxqHDgB5Dl0MBAOdhsaFrBHe2yaSAIBrNo3i3AordDc1KvQOd4o2Ap1DFwPBQKehcTFXBgBsGR8BAGyfSOIsu1xcZfWNdx6HznVFg8RAp6FxMWdU6DNmoG/NJLFcrHJMtIuKZZsVOgM9EAx0GhoXsmVMpxOIR42P9VQ6AQBYKlaCbNZQ6V6h86JokBjoNDTmc6VGdQ4Am8eMQL+ULwfVpKFTNAO90yiXeDSCRCyCQoUXRYPAQKehMZ8rY0sm2Xg+lTbCfbHACt0t1sRbncahA0A6EW0EP/mLgU5DY7FYwbTZzQIA02aFvpBnoLulW4UOAKkE50QPCgOdhsZSoYpNqfWl0axwX2CF7hqrQh+Nd67Qx0ZijYun5C8GOg2FSq2OfLmGqdR6hZ5JxhGLCBYL7EN3S9GcOjfSZoFoS2okylv/A8JAp6GwbI5k2dTU5RKJCCbTCXa5uGij1Yos6USscQMS+YuBTkNhqVgFgMsqdMDodmGXi3s2Wq3IkkpEOWwxIAx0GgrWSJbJpj50wLgwusBhi64plDuvVmRJj7BCDwoDnYaC1eUymb68Qp9Kj7BCd1GxstbxLlFLeoQVelAY6DQUFq1Ab+lymUzFsbJaDaJJQ2mj9UQt6USMF0UDwkCnobBs9qFvaulymRiNI7taRb3OBbTcsNF6opZUIoZS1VgIg/zFQKehsFioIJWIItkyPjqTjKOuQJ4Voys2Wk/UYs24yLtF/cdAp6GwVKxc0d0CGBU6AGTZ7eKKjdYTtVijYHhh1H8MdBoKy8UqJtPxK17PmIHOfnR3FMr2K3Te/u8/BjoNhexqFZlku0A3qkUGev9qa3WUa/WuFbr1Pm//9x8DnYZCrlTDePLKoFnvcmG12K9i1ZxpscsoF6uC50gX/zHQaSjkSlWMt6nQ2Yfunm6rFVkaFToD3XcMdBoKnSp09qG7p9tqRZb1dUXZ5eI3BjqFXr2uyFdqbSv0sUQMEQGyJQZ6vxoVuu1RLqzQ/cZAp9DLV2pQBTJtKvRIRJAZ5d2ibmhU6F1HuVgLRbNC9xsDnUIvVzKCpl2XC2DcXMRA75+d1YqA9S6ZIoct+o6BTqFnXfBs1+UCrN/+T/0pNC6Kblyhc6Ho4DDQKfS6VegT7HJxxWrFGra4cYUOcKHooDDQKfRypY0r9MxojIHuArujXIxtYhzlEgBbgS4it4nIMRE5LiL3d9jmV0TkqIgcEZE/d7eZRJ3ZqdCzJVaL/So6qdA5J3oguv7NiEgUwIMA3g7gDIDDInJQVY82bbMfwIcBvEFVl0Rki1cNJmq1XqF3uCjKLhdXFMo1xKOCRKx7HZge4ZzoQbBTod8C4LiqnlDVCoCHAdzVss1vAXhQVZcAQFUvuttMos6s6rvdXC7W65VaHaUquwD6YWc9UQsXig6GnUDfAeB00/Mz5mvNrgdwvYj8k4g8ISK3tduRiNwrInMiMjc/P99bi4la5EpG5TjSoXKc4N2iriiUa0jb6D8HuFB0UNy6KBoDsB/AWwDcA+DTIrKpdSNVfUhVZ1V1dmZmxqVD09XOmsdFRNq+b93+n+Pdon0pVtaQ6jKPi4ULRQfDTqCfBbCr6flO87VmZwAcVNWqqr4E4KcwAp7Ic53mcbFY7/HCaH8KFVbog85OoB8GsF9E9olIAsDdAA62bPMVGNU5RGQzjC6YEy62k6gjo0LvHOjWlAA5BnpfimX7fehjvCgaiK6Brqo1APcBeBTAcwAeUdUjIvIxEbnT3OxRAAsichTAYwD+s6oueNVooma5Ug3jI+0viALrF0t5t2h/CpVa17tELVwoOhi2vm5V9RCAQy2vPdD0WAH8jvmHyFe5Ug17plMd37duOGKF3h9Ho1yaForudMMXuY93ilLodVrcwrLeh84KvR+FsrMKHeBC0X5joFPo5Uq1xtqh7aQSUUQjwlEufeqlQueFUX8x0CnUNlrcwiIiGE/G2OXSB1V1OMrFCH7O5+IvBjqF2kaLWzTLJDmFbj9K1TpUgVGnFTpHuviKgU6h1m1iLgsr9P5YwTxmsw+dC0UHg4FOodZtcQsLA70/1nqizvvQ2eXiJwY6hZrdCj2TjHOUSx+sCt35KBd+ifqJgU6h1m1xC8t4Ms4KvQ/FxuIW9mdbBLhQtN8Y6BRqTvrQWaH3zu56opbUCBeKDgIDnUKt2+IWlsxoHPlyDXXeit4TpxU6F4oOBgOdQq3b4haWTDIGVWOYIzlndZ2M2Zw+F+BC0UFgoFOodVvcwtK4/Z9j0XtSdLBAtIULRfuPgU6h1m1xC0uGE3T1Zb0P3UGFPsIK3W8MdAq1botbWDjjYn+KlRoigq6/CTVLJWLI86KorxjoFGrdFrewsMulP4XyGtKJWNffhJqNcRk63zHQKdS6LW5haawrWmag96JYqTWGItrFZej8x0CnULPf5cJl6PpRqKw1bhayiwtF+4+BTqHWbXELC7tc+mMsbuEs0FMctug7BjqFmt0KfSQWxUgswgq9R4VyzdGQRcCo0HlR1F8MdAota3GLbnOhW8aT8caNSORMsbLmuEJPc6Fo3zHQKbSsxS3sLkKcGeV8Lr0qVHqp0NcXiiZ/MNAptOxOzGXhjIu9K5adXxTlQtH+Y6BTaNld3MKSScZ4UbRHhR6GLXKhaP8x0Cm0nFbomWS8MTsj2aeqRh86K/SBx0Cn0LI7da6Fy9D1plwzLmw6vihqVugc6eIfBjqFlhXO1l2g3XCRi95YXSZ2F7ewcKFo/zHQKbScVuiZZBylah3VtbqXzRo6VpeJ3cUtLFwo2n8MdAotu4tbWHj7f28aC0Q7HLbIhaL9x0Cn0LK7uIXF6prhSBdnrAo71cONRc0/T95joFNo2V3cwsI50XtT7LVC57BF3zHQKbSyJfu3/QPNXS6s0J1oVOgO+9C5ULT/bAW6iNwmIsdE5LiI3L/Bdr8sIiois+41kag9uzMtWqy+do50cabXUS4AF4r2W9dAF5EogAcB3A7gAIB7RORAm+3GAXwIwA/cbiRRO3ZnWrQ0ptBll4sjjS4Xh33oABeK9pudCv0WAMdV9YSqVgA8DOCuNtt9HMAnAJRcbB9RR9lVe8vPWRoVOi+KOmJ1mTi9UxTgQtF+sxPoOwCcbnp+xnytQUReA2CXqv7tRjsSkXtFZE5E5ubn5x03lqhZrlSzPWQRAMY4bLEnxXINIkAy7vySWyoRYx+6j/q+KCoiEQB/AOB3u22rqg+p6qyqzs7MzPR7aLrKOe1Dj0YEYyO8/d+pQmUNqXjU0QLRlrGRGEe5+MhOoJ8FsKvp+U7zNcs4gH8G4NsichLA6wAc5IVR8lJtrY5CZc1RlwvA2/97kS/VGr/dOMWFov1lJ9APA9gvIvtEJAHgbgAHrTdVdUVVN6vqXlXdC+AJAHeq6pwnLSbC+oRPdudxsXDGRefy5RrGerggCnChaL91DXRVrQG4D8CjAJ4D8IiqHhGRj4nInV43kKgdp1PnWjjjonO5cs1R11YzLhTtL1v/GlT1EIBDLa890GHbt/TfLKKNWd0mTm4sAoxAn8+XvWjS0MqXnI0mapYe4bBFP/FOUQql7KpVoTvschnlMnRO9dPlkkpEsVpd40LRPmGgUyjlGhW6s0Bnl4tz+VLvgW79HLtd/MFAp1DqtQ89k4wju1qFKitGu3J9jXLhMnR+YqBTKDld3MIynoyjVleUqlzkwo56XZGv1DDe8ygXzrjoJwY6hVK21Fsf+vp8Lhy6aEexugZV5+fZwgrdXwx0CqVcqYpk3Jie1Qlr3DrHotuTN784e+1yseZQ50LR/mCgUygZMy06rxo546Iz+bLxxdfPjUUAL4r6hYFOoZTtcWy0NW6dMy7ak+23QudC0b5ioFMoOZ1p0ZLhMnSOWF0uvV4U5ULR/mKgUyhlHS5uYeG6os5Yfd+996FzoWg/MdAplHKlak8VOke5OJPvcTSRxepy4ReoPxjoFErZ1Royo70siRZFNCIc5WJTzqrQe+xyiUUjSCeiPN8+YaBTKDld3MIiIrz934HGsMUeAx0whoryNyJ/MNApdErVNZRrdUw4nAvdYt3+T93lStXGbzW9Ms43v0D9wECn0Fkxw7jXQGeFbl8/My1aMqNcJcovDHQKneWiEQ6bUr0HOgPGnly594m5LJlkvPElTN5ioFPoLBcrAIBNo4meft5Yho4Vuh35Hu/IbcY+dP8w0Cl0llf7rdAZ6HYZXS7RvvaRScbYh+4TBjqFzkqx/z50XhS1Z2W1t/H+zYxVoqqoc9UizzHQKXSWV80ulx4r9MxoHPlKjQFjw8pqtecvTsvEaBx1BQq8/d9zDHQKneViFdGI9Dz6IpOMQRXIM2C6ciPQrQqfM1x6j4FOobOyWsWm0ThEehsbPc4ZF20pVddQqdUbc8j3yrqjl+fbewx0Cp3l1SomeuxuATjjol0rfV58tljnm0MXvcdAp9BZKRoVeq8446I9/d7AZbEqfFbo3mOgU+gsr1awKdXbGHSAXQB2uRbo7EP3DQOdQmfZpQqdN7tsrN/hoRZ+gfqHgU6hs1Ks9lWhT5k/u1iouNWkobTsUoVujUbiF6j3GOgUKtW1OnLlWl8X6jKjMcQiwkDvwq0ul1g0grER3i3qBwY6hUrWhZEXIoLJdIKB3oUV6P3O5QKYt/+zQvccA51Cxa1ugOl0AgsM9A1lV6sYT8b6mgvdkhmNN2bJJO/YCnQRuU1EjonIcRG5v837vyMiR0XkGRH5lojscb+pRMBSwbrtv/c+dACYYoXelRt3iVomU4nGLJnkna6BLiJRAA8CuB3AAQD3iMiBls1+BGBWVX8WwF8C+KTbDSUC0Kiqp9MMdK+5GehTYwksMtA9Z6dCvwXAcVU9oaoVAA8DuKt5A1V9TFWL5tMnAOx0t5lEhoW8EQqbx0b62s90OoGFfNmNJg2thUIFU31+cVqmUonGb1fkHTuBvgPA6abnZ8zXOnkfgK+3e0NE7hWRORGZm5+ft99KIpMVwv0GzVR6BNlSDdW1uhvNGkqLhbJrgT6ZTmB5tYo1znDpKVcviorIrwGYBfCpdu+r6kOqOquqszMzM24emq4SC4UKMskYErH+PrpTY0ZQsWrsbDHvZoUehyrnc/GanX8VZwHsanq+03ztMiLyNgAfAXCnqvJ3WfLEQqGC6T67W4D1PniOdGmvVF1DobLW97UKy2SaN3P5wU6gHwawX0T2iUgCwN0ADjZvICI3A/hTGGF+0f1mEhkW8mVXQmaKAbMh67xMpfv/8gSAaXM/PN/e6hroqloDcB+ARwE8B+ARVT0iIh8TkTvNzT4FYAzA/xORp0XkYIfdEfVlIV/B9Fj/gc4KfWPrge5WhR6/bL/kDVtLvqjqIQCHWl57oOnx21xuF1FbC4UyXrt3su/9NCp0jnRpqzE81IUvT2D9fC9x6KKneKcohUa9rlgsVFzpctmUSkCEFWMniwV3RhNZJjkhmi8Y6BQaS8UK6upOyEQjgqlUApcYMG1Z4/3duiiajEeRSkQ5qshjDHQKjYs5o2rcmkm6sr+Z8RFczLLLpZ3FQgXRiDQWp3DD9FgC8+zi8hQDnULjfLYEANiacWfkxbaJJC6Y+6TLLRYqmEwlEHFhYi7L1vEkv0A9xkCn0LiwYgW6OxX61vFk40uCLncxV8aWcXe+OC1bMiO4mOP59hIDnULDCt8t4y4F+kQSl/Jl1Hj7/xXOr5SwbcKd82zZwgrdcwx0Co0L2RI2jyX6vu3fsjUzAlWwX7eNi7mSa11bli2ZEeTKNRQrXLnIKwx0Co0L2bJr1TkAbDO7bi6warxMpVbHpXzFta4ty1bz745VuncY6BQabncDWIF1foX9us2sfu5tLgf6FrPi54Vo7zDQKTQuZEuuVo3Wvnih7nLWbyyuV+iN880K3SsMdAqFcm0NC4WKq1XjdDqBeFTwyjIDvdmFrLujiSzWqBlW6N5hoFMonFlaBQDsmhp1bZ+RiGDnZAqnl4rdN76KWIHr9iiXidE4kvEIzrGLyzMMdAqFlxeN0N09lXJ1v7umUji9yEBvdn6lhEQsgsmUe3eJAoCIYNckz7eXGOgUCqc9CvTdU6ONLwsynFooYtfkKETcu0vUsmsqxfPtIQY6hcKphSKS8QhmXL57cfdUCsvFKpdGa3JqsYg902lP9r17KoUzS6tQ5dqiXmCgUyi8vFjE7qmU61WjVfGzG8Cgqnh5oeD6b0KWnZOjyJdrWC7yC9QLDHQKhdNmoLttFwP9MguFCgqVNeyZ9ibQrb9Ddrt4g4FOA69eV5xaKGL3lPvdAFbXwolLBdf3HUanFozz4FWgN75AObLIEwx0GngvLxaxWl3DDdvGXN/32EgMOzaN4vnzOdf3HUanFoyg9aoP3fqiODHPL1AvMNBp4Flhe8O2jCf7v3HbOI6dz3qy77B5cT6PWESwc9K98f7NUokYdk+lcIxfoJ5goNPAO3Y+BxHg+q3uV+gAcMO2cZyYL6BS4zS6R1/J4rotYxiJRT07xo3bxvEcv0A9wUCngff8+Sz2TKWQSsQ82f8N28ZRqytenM97sv8wOXouiwPXePObkOXGbeM4eamAUnXN0+NcjRjoNPCOnsvihm3jnu3/wHYjwI68cnVXjZfyZVzIlhvnwys3bs+grsALF/gF6jYGOg20i9kSTi0UMbtnyrNjvGpmDBOjccydXPTsGGHw3DnjC83rCt36wvjJ2RVPj3M1YqDTQPuhGbK37PMu0CMRweyeycaxrlZPnlqCCHDTNROeHmfPdApbxkfwg5cWPD3O1YiBTgPthy8tIpWI4iaPq8Z/sW8KJ+YLmL+K5+p+4sQCbromg4lRdyflaiUieN2103jixAKnAHAZA50Glqriuy9cwuzeKcSi3n5UX/+qaQDAP/503tPjDKpCuYanXl7G6/ZN+3K81107jQvZMl7iDV2uYqDTwPrphTxeulTAOw5s9fxYr94xgZ2To/jaM694fqxB9O1j86jU6ni7D+caAN50/WYAwNefPe/L8a4WDHQaWH/11BnEIoJ33rTN82OJCN71s9vxvRcuYSF/9XW7fOXps9g8lsDsXu+uVTTbOZnCLXun8OWnzrDbxUUMdBpI2VIVf3H4NN550zbXp8zt5D2v3YU1VXzmey/5crxBcXqxiG89dwG/MrsL0Yj7c6B38m9eswMvzhfw+Iu8OOoWBjoNpP/1D8exslrF+9/8Kt+Oed2WMbzr1dvxue+fvKpmX/zE3z2PeDSCX791j6/HfffNO3DNRBL/4++ex1qdVbobbAW6iNwmIsdE5LiI3N/m/RER+Qvz/R+IyF63G0pXj68+fRaf/u4J3HPLbrx6p7dD6Fr93jtvREQEH/zzp7BcrPh67CB89nsv4WvPnMMH3nIdtk94M39LJ8l4FP/l9hvxzJkVfPxrR1FnqPet673UIhIF8CCAtwM4A+CwiBxU1aNNm70PwJKqXicidwP4BIBf9aLBNJxWilUcObeCv5w7gy//6CxuvXYaD/zCAd/bsXs6hT/81Z/DB774FN71R9/Dv3/ztXjz9TPYOZnytTvCK6qK7GoNPzm7gi8dfhl/+8w5vOPAVnzwX/n3m1Czu35uB54+vYz/808n8cLFHH7zDfvw2j2TmBiNe7IE3rCTbhckRORWAB9V1Xeazz8MAKr635u2edTc5nERiQE4D2BGN9j57Oyszs3NOW7wI4dP46Hvnmg8b3eIK17RLu+32U/rNu3+n2jLVq3b2LnW0+24bffb5bjt9nPlNm3OW5fzZOdcX3kOuv9Mva4oVIx5PZLxCN576178x7dfj2Tcuwmiunny1BL+28Fn8exZ4+7JWEQwnowhGY9iJBZBpDls2j80njdtZ50LbfyP8R/V9b9R1fW/X9XLz+eG26F5W208bt1/tVZvnOvReBT3vula/Pa/vs7zYaEbUVX83ydO4Q+/+QIWC8ZvRcl4BJlkHPFoBNGIIBYRiGBoQv5Db92PX/zn1/T0syLypKrOtnvPzmxHOwCcbnp+BsC/7LSNqtZEZAXANIBLLQ25F8C9ALB7925bjW81mU7ghq0t83q0+Tve6B9Whx9B62el2z7a7ueKfbT5ma7HadO2lq2u2Efbz3mXn2n3E13a375trfvo/o/u8jwUbJ9IYt/mNG591TTSI95MwuXEa/dM4m/ueyOOX8xj7tQSTi8WkSvVUKquoVyrNwXreuJu9J2pUONcmv+/BevnyXi8fh5Fms66GOdHLvs563HT62K+0mW7aERwzcQorp1J4/Wv2ozRRHBfmhYRwa/fuhfvmd2Fp04t4cgrWVzMlZBdrWFNFWt1Ra2uQ9Ul49XNW77+y1HVhwA8BBgVei/7ePuBrb6NlaWrm4hg/9Zx7G8tIMgTyXgUr79uM15/3eagmxJadn7POgtgV9PzneZrbbcxu1wmAHAsEhGRj+wE+mEA+0Vkn4gkANwN4GDLNgcBvNd8/G8B/MNG/edEROS+rl0uZp/4fQAeBRAF8FlVPSIiHwMwp6oHAXwGwBdE5DiARaKN9/8AAAR/SURBVBihT0REPrLVh66qhwAcanntgabHJQDvcbdpRETkBO8UJSIaEgx0IqIhwUAnIhoSDHQioiHR9dZ/zw4sMg/gVI8/vhktd6EOCLbLGbbLuUFtG9vlTD/t2qOqM+3eCCzQ+yEic53mMggS2+UM2+XcoLaN7XLGq3axy4WIaEgw0ImIhkRYA/2hoBvQAdvlDNvl3KC2je1yxpN2hbIPnYiIrhTWCp2IiFow0ImIhkSoAl1EPiUiz4vIMyLy1yKyqem9D5uLVB8TkXf63K73iMgREamLyGzT63tFZFVEnjb//MkgtMt8L7Dz1dKOj4rI2aZzdEdQbTHbs+GC6EERkZMi8hPzHDlfu9G9dnxWRC6KyLNNr02JyDdE5AXzv5MD0q7AP1sisktEHhORo+a/xQ+Zr3tzzlQ1NH8AvANAzHz8CQCfMB8fAPBjACMA9gF4EUDUx3b9DIAbAHwbwGzT63sBPBvg+erUrkDPV0sbPwrgPwX92TLbEjXPxbUAEuY5OhB0u8y2nQSweQDa8SYAr2n+XAP4JID7zcf3W/8uB6BdgX+2AGwH8Brz8TiAn5r//jw5Z6Gq0FX171W1Zj59AsbqSQBwF4CHVbWsqi8BOA7gFh/b9ZyqHvPreHZt0K5Az9cAuwXAcVU9oaoVAA/DOFdkUtXvwFjzoNldAD5nPv4cgHf72ih0bFfgVPWcqj5lPs4BeA7GGsyenLNQBXqL3wTwdfNxu4Wsd/jeovb2iciPROQfReTng26MadDO131mN9png/h1vcmgnZdmCuDvReRJc7H1QbJVVc+Zj88DGKRFfwflswUR2QvgZgA/gEfnLPjl1VuIyDcBbGvz1kdU9avmNh8BUAPwxUFqVxvnAOxW1QUReS2Ar4jITaqaDbhdvtqojQD+GMDHYQTWxwH8Powva7rcG1X1rIhsAfANEXnerEoHiqqqiAzKWOiB+WyJyBiAvwLwH1Q1KyKN99w8ZwMX6Kr6to3eF5HfAPALAN6qZgcU7C1k7Wm7OvxMGUDZfPykiLwI4HoArl3U6qVd8OF8NbPbRhH5NICvedUOG3w9L06o6lnzvxdF5K9hdA8NSqBfEJHtqnpORLYDuBh0gwBAVS9Yj4P8bIlIHEaYf1FVv2y+7Mk5C1WXi4jcBuD3ANypqsWmtw4CuFtERkRkH4D9AH4YRBubiciMiETNx9fCaNeJYFsFYIDOl/lhtvwSgGc7besDOwui+05E0iIybj2GMTggyPPUqnmR+PcCGJTfDAP/bIlRin8GwHOq+gdNb3lzzoK8AtzDFePjMPo4nzb//EnTex+BMULhGIDbfW7XL8Hoby0DuADgUfP1XwZwxGzrUwB+cRDaFfT5amnjFwD8BMAz5od8e8CfsTtgjER4EUa3VWBtaWrTtTBG3PzY/DwF1i4AX4LRlVg1P1vvAzAN4FsAXgDwTQBTA9KuwD9bAN4Io8vnmabcusOrc8Zb/4mIhkSoulyIiKgzBjoR0ZBgoBMRDQkGOhHRkGCgExENCQY6EdGQYKATEQ2J/w+70FR8ZFkxJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def H(x):\n",
    "    '''\n",
    "    This function has a global maximum at x=10\n",
    "    and local maximum at x=-2\n",
    "    '''\n",
    "    return np.exp(-(x-10)**2)+0.8*np.exp(-(x+2)**2)\n",
    "\n",
    "x = np.linspace(-20,20,num=1000)\n",
    "\n",
    "plt.plot(x,H(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result from cross-entropy optimization: 10.00\n",
      "Global maximum: 10.0\n"
     ]
    }
   ],
   "source": [
    "\n",
    "N = 1000\n",
    "rho = 0.1\n",
    "mu, sigma = -20.0,1000.0\n",
    "N_rho = int(N*rho)\n",
    "\n",
    "steps =200\n",
    "\n",
    "\n",
    "for k in range(0,steps):\n",
    "    x_samples = np.random.normal(loc=mu,scale=sigma,size=N)\n",
    "    H_samples = H(x_samples) \n",
    "    data = np.array((x_samples,H_samples)).T\n",
    "\n",
    "    # Now sort these tuples according to H(x_i) values\n",
    "    sorted_data = np.array(sorted(data,key= lambda x: x[1],reverse=True))\n",
    "    x_elite_samples = sorted_data[0:N_rho,0]\n",
    "    \n",
    "    # Update the Gaussian distribution parameters with the elite samples\n",
    "    mu = np.mean(x_elite_samples)\n",
    "    sigma = np.std(x_elite_samples)\n",
    "    \n",
    "\n",
    "\n",
    "print('Result from cross-entropy optimization: {0:0=.2f}'.format(mu))\n",
    "print('Global maximum: {}'.format(10.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cirq",
   "language": "python",
   "name": "cirq"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}