From 89e52daa177ac03bd2aab3cead8276b2147a03bb Mon Sep 17 00:00:00 2001 From: didier amyot Date: Tue, 15 Dec 2015 21:32:08 -0500 Subject: [PATCH] First commit --- .gitignore | 2 + .../accep_reject-checkpoint.ipynb | 144 ++++++++++ .../brute_force-checkpoint.ipynb | 189 +++++++++++++ accep_reject.ipynb | 160 +++++++++++ brute_force.ipynb | 187 +++++++++++++ fuck_me.ipynb | 112 ++++++++ gibbs_sampler.ipynb | 258 ++++++++++++++++++ importance_sampling.ipynb | 179 ++++++++++++ makefile | 2 + network/__init__.py | 0 network/acceptance_rejection.py | 22 ++ network/brute_force.py | 27 ++ network/gibbs.py | 36 +++ network/importance_sampling.py | 92 +++++++ network/network.py | 97 +++++++ play_around.ipynb | 220 +++++++++++++++ requirements.txt | 35 +++ tests/test_brute_force.py | 26 ++ tests/test_gibbs.py | 35 +++ tests/test_importance_sampling.py | 121 ++++++++ tests/test_network.py | 89 ++++++ 21 files changed, 2033 insertions(+) create mode 100644 .gitignore create mode 100644 .ipynb_checkpoints/accep_reject-checkpoint.ipynb create mode 100644 .ipynb_checkpoints/brute_force-checkpoint.ipynb create mode 100644 accep_reject.ipynb create mode 100644 brute_force.ipynb create mode 100644 fuck_me.ipynb create mode 100644 gibbs_sampler.ipynb create mode 100644 importance_sampling.ipynb create mode 100644 makefile create mode 100644 network/__init__.py create mode 100644 network/acceptance_rejection.py create mode 100644 network/brute_force.py create mode 100644 network/gibbs.py create mode 100644 network/importance_sampling.py create mode 100644 network/network.py create mode 100644 play_around.ipynb create mode 100644 requirements.txt create mode 100644 tests/test_brute_force.py create mode 100644 tests/test_gibbs.py create mode 100644 tests/test_importance_sampling.py create mode 100644 tests/test_network.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5f1c35b --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/env +__pycache__ diff --git a/.ipynb_checkpoints/accep_reject-checkpoint.ipynb b/.ipynb_checkpoints/accep_reject-checkpoint.ipynb new file mode 100644 index 0000000..9be6c7e --- /dev/null +++ b/.ipynb_checkpoints/accep_reject-checkpoint.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import is_valid_network, get_number_of_activated_antenna, pp_network\n", + "from numpy import mean\n", + "from random import Random" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "random = Random(666)\n", + "def get_some_possible_network(n, box_size):\n", + " number_of_square = box_size * box_size\n", + " for i in range(n):\n", + " yield random.randint(0, 2 ** number_of_square - 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "box_size = 5\n", + "n = 1000000\n", + "\n", + "valid_networks = []\n", + "for network in get_some_possible_network(n, box_size):\n", + " if is_valid_network(network, box_size):\n", + " valid_networks.append(network)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFfBJREFUeJzt3X/sXXd93/HnK3jJAmmiFJp8WQxxolAwTCxhYLp5XS+j\nmAS2OENTauhG0kCHlGVEZULYkSp/VVWUICWo2hZpbVjkRMkyw4bibFViQrhUIM1ZiV3i2A3WIhvH\nw9+WEshYELOT9/64x+ba/sbf+/1xfa6Pnw/pq5zvueee+/I3X7987uecez6pKiRJ3XVW2wEkSeNl\n0UtSx1n0ktRxFr0kdZxFL0kdZ9FLUseNXPRJzkqyPcmW5vuNSZ5L8mTzdfXQthuS7EmyO8macQSX\nJI1m2Ty2vRV4Gjh/aN2dVXXn8EZJVgLXAyuB5cBjSd5UXrAvSa0Y6Yg+yXLgA8Ddxz80y+ZrgQer\n6nBV7QX2AKsWE1KStHCjDt18Afg0cPxR+S1JdiS5O8kFzbpLgP1D2xxo1kmSWjBn0Sf5IDBTVTs4\n9gj+LuDyqroSOAjcMZ6IkqTFGGWMfjVwbZIPAOcCv5Dk3qr66NA2fww83CwfAN4w9NjyZt0xkjhm\nL0kLUFWzDZu/ojmP6Kvqtqp6Y1VdDqwDHq+qjyaZGtrsQ8DOZnkLsC7J2UkuA64AnniFfU/c18aN\nG1vPYCYznYm5zDTa10LM56qb430+yZXAy8Be4BNNee9KshnYBRwCbq6FppMkLdq8ir6qvgF8o1n+\n6Em2+wPgDxYXTZK0FPxk7HF6vV7bEU5gptGYaXSTmMtM45O2RlWSOKIjSfOUhFrqk7GSpNObRS9J\nHWfRS1LHWfSS1HEWvSR1nEUvSR1n0UtSx1n0ktRxFr0kdZxFL0kdZ9GrE6amVpCEJExNrWg7jjRR\nvNeNOiEJP5/pMgu+b7c06bzXjSTpBBa9JHXcyEWf5KwkTybZ0nx/YZKtSZ5J8miSC4a23ZBkT5Ld\nSdaMI7gkaTTzOaK/lcH0gEesBx6rqjcDjwMbAJK8FbgeWAlcA9yVwQCqJKkFIxV9kuXAB4C7h1av\nBTY1y5uA65rla4EHq+pwVe0F9gCrliStJGneRj2i/wLwaX5+WQPAxVU1A1BVB4GLmvWXAPuHtjvQ\nrJMktWDOycGTfBCYqaodSXon2XTe17NNT08fXe71ep2Zn1GSlkq/36ff7y9qH3NeR5/ks8A/Bw4D\n5wK/AHwFeCfQq6qZJFPA16tqZZL1QFXV7c3zHwE2VtW24/brdfRaMl5HrzPFWK6jr6rbquqNVXU5\nsA54vKr+BfAwcGOz2Q3AQ83yFmBdkrOTXAZcATwxn1CSpKUz59DNSXwO2JzkJmAfgyttqKpdSTYz\nuELnEHCzh+6S1B5vgaBOcOhGZwpvgSBJOoFFL0kdZ9FLUsdZ9JLUcRa9JHWcRS9JHWfRS1LHWfTS\nGAzPYes8tmqbH5hSJ0zaB6aOzQOTkEnd4AemJEknsOglqeMseknqOItekjrOopekjrPoJanjLHpJ\n6rg5iz7JOUm2Jdme5OlmDlmSbEzyXJInm6+rh56zIcmeJLuTrBnnH0CSdHJzTiVYVT9L8p6qejHJ\nq4BvJVndPHxnVd05vH2SlQymFVwJLAceS/ImPx0lSe0Yaeimql5sFs9pnvN88/1sn85aCzxYVYer\nai+wB1i1yJySpAUaqeiTnJVkO3AQ6FfVruahW5LsSHJ3kguadZcA+4eefqBZJ0lqwZxDNwBV9TJw\nVZLzga1Jfg24C/i9qqokvw/cAXx8Pi8+PT19dLnX69Hr9ebzdEnqvH6/T7/fX9Q+5n1TsyS/C7xY\nVXcMrbsUeLiq3p5kPVBVdXvz2CPAxqradtx+HLbXkvGmZjpTjOWmZkled2RYJsm5wPuAHUmmhjb7\nELCzWd4CrEtydpLLgCuAJ+YTSpK0dEYZunk9sCmDQ5SzgPuq6mtJ7k1yJfAysBf4BEBV7UqyGdgF\nHAJu9tBdktrj/ejVCQ7d6Ezh/eglSSew6CWp4yx6Seo4i16SOs6il6SOs+glqeMseknqOItekjrO\nopekjrPoJanjLHpJ6jiLXpI6zqKXpI6z6CWp4yx6Seo4i16SOm6UqQTPSbItyfYkTyf5bLP+wiRb\nkzyT5NEj0w02j21IsifJ7iRrxvkHkCSd3EgzTCV5dVW9mORVwLeAfwNcC/x1VX0+yWeAC6tqfZK3\nAvcD7wKWA48Bbzp+OilnmNJScoYpnSnGNsNUVb3YLJ7TPOd5YC2wqVm/CbiuWb4WeLCqDlfVXmAP\nsGo+oSRJS2ekok9yVpLtwEGgX1W7gIuragagqg4CFzWbXwLsH3r6gWadJKkFy0bZqKpeBq5Kcj7w\naJIex74vZZbv5zQ9PX10udfr0ev15rsLSeq0fr9Pv99f1D5GGqM/5gnJ7wI/BT4G9KpqJskU8PWq\nWplkPVBVdXuz/SPAxqradtx+HKPXknGMXmeKsYzRJ3ndkStqkpwLvA/YDmwBbmw2uwF4qFneAqxL\ncnaSy4ArgCfmE0qStHRGGbp5PbApg0OUs4D7quprzZj95iQ3AfuA6wGqaleSzcAu4BBws4fuktSe\neQ/dLNkLO3Rz2pqaWsHMzD4ALr74Ug4e3NtuIBy60ZljIUM3Fr3mbdJKFSYvk0WvcRnbdfSSpNOX\nRS9JHWfRS1LHWfSS1HEWvSR1nEUvSR1n0UtSx1n0ktRxFr0kdZxFL0kdZ9FLUsdZ9JLUcRa9JHWc\nRS9JHTfKDFPLkzye5OkkTyX51836jUmeS/Jk83X10HM2JNmTZHeSNeP8A0iSTm7O+9E388FOVdWO\nJOcB3wbWAr8B/J+quvO47VcCDwDvApYDjwFvOv7m896P/vQ1afd+h8nL5P3oNS5juR99VR2sqh3N\n8k+A3cAlR15zlqesBR6sqsNVtRfYA6yaTyhJ0tKZ1xh9khXAlcC2ZtUtSXYkufvIBOIM/hHYP/S0\nA/z8HwZJ0ik2ctE3wzZfBm5tjuzvAi6vqiuBg8Ad44koSVqMZaNslGQZg5K/r6oeAqiqvxra5I+B\nh5vlA8Abhh5b3qw7wfT09NHlXq9Hr9cbMbYknRn6/T79fn9R+xhpcvAk9wI/qKpPDa2bqqqDzfLv\nAO+qqo8keStwP/BuBkM2X8WTsZ0yaSc+YfIyeTJW47KQk7FzHtEnWQ38JvBUku0MfntvAz6S5Erg\nZWAv8AmAqtqVZDOwCzgE3GyjS1J7RjqiH8sLe0R/2pq0o2eYvEwe0WtcxnJ5pSTp9GbRS1LHWfSS\n1HEWvSR1nEUvnSGmplaQ5OjX1NSKtiPpFPGqG83bpF3hApOXaRKvupnETJo/r7qRJJ3AopekjrPo\nJanjLHpJ6jiLXpI6zqKXpI6z6CWp4yx6Seo4i16SOs6il6SOm7PokyxP8niSp5M8leSTzfoLk2xN\n8kySR5NcMPScDUn2JNmdZM04/wCSpJOb8143SaaAqarakeQ84NvAWuC3gL+uqs8n+QxwYVWtH5oz\n9l0MJgZ/DOeM7ZRJu68MTF6mSbyvzCRm0vyN5V43VXWwqnY0yz8BdjMo8LXApmazTcB1zfK1wINV\ndbiq9gJ7gFXzCSVJWjrzGqNPsgK4EvgfwMVVNQODfwyAi5rNLgH2Dz3tQLNOktSCkYu+Gbb5MnBr\nc2R//Hs+3wNK0gRaNspGSZYxKPn7quqhZvVMkouraqYZx//LZv0B4A1DT1/erDvB9PT00eVer0ev\n15tXeEnqun6/T7/fX9Q+Rpp4JMm9wA+q6lND624HflhVt7/Cydh3Mxiy+SqejO2USTvxCZOXaRJP\nfE5iJs3fQk7GjnLVzWrgT4GnGPyWFHAb8ASwmcHR+z7g+qr6UfOcDcDHgEMMhnq2zrJfi/40NWml\nCpOXaRJLdRIzaf7GUvTjYtGfviatVGHyMk1iqU5iJs2fUwlKkk5g0UtSx1n0ktRxFr0kdZxFL0kd\nZ9FLUsdZ9JLUcRa9JHWcRS9JHWfRS1LHWfSS1HEWvSR1nEUvSR1n0UtSx1n0ktRxFr0kddycRZ/k\ni0lmknxnaN3GJM8lebL5unrosQ1J9iTZnWTNuIJLkkYzyhH9PcD7Z1l/Z1W9o/l6BCDJSuB6YCVw\nDXBXBtPaSJJaMmfRV9U3gedneWi2Al8LPFhVh6tqL7AHWLWohJKkRVnMGP0tSXYkuTvJBc26S4D9\nQ9scaNZpgaamVpCEJExNrWg7jqTT0LIFPu8u4PeqqpL8PnAH8PH57mR6evrocq/Xo9frLTBOd83M\n7OPIhM4zM46CSWeafr9Pv99f1D4yyizwSS4FHq6qt5/ssSTrgaqq25vHHgE2VtW2WZ5XzkA/t8Ep\njiM/pzAJPzMzze3YPGAmLZUkVNW8jvpGHboJQ2PySaaGHvsQsLNZ3gKsS3J2ksuAK4An5hNIkrS0\n5hy6SfIA0ANem+R7wEbgPUmuBF4G9gKfAKiqXUk2A7uAQ8DNHrZLUrtGGroZyws7dDOSSRuSADON\nYhKHSSYxk+ZvnEM3kqTTlEUvSR1n0UtSx1n0ktRxFr0kdZxFL0kdZ9FLUsdZ9JLUcRa9JHWcRS9J\nHWfRS1LHWfSS1HEWvSR1nEUvSR1n0UtSx1n0ktRxcxZ9ki8mmUnynaF1FybZmuSZJI8muWDosQ1J\n9iTZnWTNuIJLkkYzyhH9PcD7j1u3Hnisqt4MPA5sAEjyVuB6YCVwDXBXBtPaSJJaMmfRV9U3geeP\nW70W2NQsbwKua5avBR6sqsNVtRfYA6xamqiSpIVY6Bj9RVU1A1BVB4GLmvWXAPuHtjvQrJMktWTZ\nEu1nQTMMT09PH13u9Xr0er0liiNJ3dDv9+n3+4vaR0aZBT7JpcDDVfX25vvdQK+qZpJMAV+vqpVJ\n1gNVVbc32z0CbKyqbbPss5yBfm6DUxxHfk5hEn5mZprbsXnATFoqSaiqeZ37HHXoJs3XEVuAG5vl\nG4CHhtavS3J2ksuAK4An5hNIkrS05hy6SfIA0ANem+R7wEbgc8CXktwE7GNwpQ1VtSvJZmAXcAi4\n2cN2SWrXSEM3Y3lhh25GMmlDEmCmUUziMMkkZtL8jXPoRpJ0mrLoJanjLHpJ6jiLXpI6zqKXpI6z\n6CWp4yx6Sa2ZmlpBEpIwNbWi7Tid5XX0E27Srg8HM41iEq9Zn/xM7ec5HXgdvSTpBBa9JHWcRS9J\nHWfRS1LHWfSS1HEWvSR1nEUvSR1n0UtSxy1qcvAke4EfAy8Dh6pqVZILgf8MXArsBa6vqh8vMqck\naYEWe0T/MoNJwq+qqlXNuvXAY1X1ZuBxYMMiX0OStAiLLfrMso+1wKZmeRNw3SJf45QZvu+G996Q\n1BWLutdNkmeBHwEvAf+hqu5O8nxVXTi0zQ+r6hdnee7E3evGe4GMxkxzm/zfJZi8TO3nOR0s5F43\nixqjB1ZX1feT/BKwNckzHPubxCzfHzU9PX10udfr0ev1FhlHkrql3+/T7/cXtY8lu3tlko3AT4CP\nMxi3n0kyBXy9qlbOsr1H9COYxCMeM81t8n+XYPIytZ/ndHBK716Z5NVJzmuWXwOsAZ4CtgA3Npvd\nADy00NeQJC3eYoZuLga+kqSa/dxfVVuT/BmwOclNwD7g+iXIKUlaICceGeJb29GYaW6T/7sEk5ep\n/TynAycekSSdwKKXpI6z6CWp4yx6Seo4i16SOs6il6SOs+glqeMseknqOItekjrOopekjrPoJanj\nLHpJGjI801xXZpnzpmZDvOnTaMw0t8n/XYLJy9R+HpjMTMO8qZkk6QQWvSR13NiKPsnVSf4iyXeT\nfGa2bXbu3MnOnTt59tlnxxVDks54Yyn6JGcB/w54P/A24MNJ3nL8dqtXr2P16nW85S1/e4LKvt92\ngFn02w4wi37bAWbRbzvALPptB3gF/bYD6BQa1xH9KmBPVe2rqkPAg8Da4zd64YWdvPDCTs4993Je\nfPHFMUWZr37bAWbRbzvALPptB5hFv+0As+i3HeAV9NsOoFNoXEV/CbB/6PvnmnWSpHkavuRzIRYz\nOfiinX/+PwHgpz/dx7JlrUaRpIk1M7OP4Us+52ss19En+RVguqqubr5fD1RV3T60zWRdnCpJp4n5\nXkc/rqJ/FfAM8F7g+8ATwIeraveSv5gk6aTGMl5SVS8luQXYyuA8wBcteUlqR2u3QJAknRqn/JOx\nSZYneTzJ00meSvLJU51hlkznJNmWZHuT67NtZzoiyVlJnkyype0sRyTZm+TPm5/XE23nAUhyQZIv\nJdnd/D98d8t5frn5+TzZ/PfHE/K7vqH5+Xwnyf1Jzp6ATLc2XdBqHyT5YpKZJN8ZWndhkq1Jnkny\naJILJiDTP0uyM8lLSd4xyn7auAXCYeBTVfU24O8B/2q2D1OdSlX1M+A9VXUV8HbgHyVZ3WamIbcC\nu9oOcZyXgV5VXVVVq9oO0/hD4E+qaiXwd4BWhwqr6rvNz+cdwN8F/i/wlTYzJbkU+G3gqqp6O4Oh\n23UtZ3ob8DHgncCVwD9OcnlLce5h8CHPYeuBx6rqzcDjwIYJyPQU8E+Bb4y6k1Ne9FV1sKp2NMs/\nYfAXsvVr7KvqyCe2zmHwc3m+xTjA4N0P8AHg7razHCdM0H2SkpwP/GpV3QNQVYer6oWWYw37deB/\nVdX+ObccrxeA/we8Jsky4NXA/243EiuBbVX1s6p6CfhT4ENtBKmqb3Li3/u1wKZmeRNwXduZquqZ\nqtrDPK6zbPUva5IVDP4V39ZmDjg6RLIdOAj0q2oSjqK/AHyaY+8tOwkK+GqS/5nkt9sOA1wG/CDJ\nPc1QyR8lObftUEN+A/hPbYeoqueBO4DvAQeAH1XVY+2mYifwq80QyasZHNi8oeVMwy6qqhkYHKQC\nF7WcZ0FaK/ok5wFfBm5tjuxbVVUvN0M3y4F/mOTX2syT5IPATPPuJyzkUxLjs7oZkvgAg6G3f9By\nnmXAO4B/3+R6kcFb7tYl+RvAtcCXJiDL5cDvAJcCfws4L8lH2sxUVX8B3A58FfgTYDvwUpuZ5jBp\nB10jaaXom7eNXwbuq6qH2sjwSpq3/P+dwZhhm1YD1yZ5lsHR4HuS3NtyJgCq6vvNf/+Kwbhz2+P0\nzwH7q+rPmu+/zKD4J8E1wLebn1Xb3gl8q6p+2AyT/Ffg77eciaq6p6reWVU94EfAd1uONGwmycUA\nSaaAv2w5z4K0dUT/H4FdVfWHLb3+MZK87sjZ9OYt//uAHW1mqqrbquqNVXU5gxNmj1fVR9vMBJDk\n1c27MZK8BljD4O13a5q31vuT/HKz6r1MzgnsDzMBwzaNZ4BfSfI3M7hpyntp+aQ1QJJfav77RgYn\nGR9oMw7HvnveAtzYLN8AtHFgerJ39CO90z/lN5hprmb5TeCpZky8gNuq6pFTnWXI64FNzS//WQze\naXytxTyT7GLgK80tLJYB91fV1pYzAXwSuL8ZKnkW+K2W89CMOf868C/bzgJQVX/evCv8NoPhke3A\nH7WbCoD/kuQXgUPAzW2dSE/yANADXpvke8BG4HPAl5LcBOwDrp+ATM8D/xZ4HfDfkuyoqmtOuh8/\nMCVJ3TYxl8hJksbDopekjrPoJanjLHpJ6jiLXpI6zqKXpI6z6CWp4yx6Seq4/w9KwjDpHZH3KQAA\nAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(list(map(get_number_of_activated_antenna, valid_networks)), bins=100)\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "6.2328431372549016" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean(list(map(get_number_of_activated_antenna, valid_networks)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "valid_networks[0]" + ] + }, + { + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/.ipynb_checkpoints/brute_force-checkpoint.ipynb b/.ipynb_checkpoints/brute_force-checkpoint.ipynb new file mode 100644 index 0000000..9415974 --- /dev/null +++ b/.ipynb_checkpoints/brute_force-checkpoint.ipynb @@ -0,0 +1,189 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import is_valid_network, get_number_of_activated_antenna, pp_network\n", + "from numpy import mean\n", + "\n", + "def get_all_possible_network(box_size):\n", + " number_of_square = box_size * box_size\n", + " return iter(range(2**number_of_square))\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "box_size = 5\n", + "\n", + "valid_networks = []\n", + "for network in get_all_possible_network(box_size):\n", + " if is_valid_network(network, box_size):\n", + " valid_networks.append(network)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADJxJREFUeJzt3V+spPVdx/H3h64SKYGQKhzTbXtCGv/EtFlb7c0anWoE\n4p/S2KT2TyL2ouFCLGmMAbnZc2VSLzAkpjcClRqxUZJawLSllU4MmAICK7BsKYlCRdljoyDd9KaB\nrxdnllnZsztzzplznv3uvF/JhOc8Z57z/PgxvPPMM/PMpKqQJPV03tADkCRtnxGXpMaMuCQ1ZsQl\nqTEjLkmNGXFJamxmxJPsT3J/kiNJnkzy+5P1h5K8kOSxye2q3R+uJOlkmfU+8SQrwEpVHU5yIfAo\ncDXw28D3qurm3R+mJGkz+2bdoaqOAccmy8eTHAXeOvl1dnFskqQZtnROPMkqcAB4aLLquiSHk9ya\n5OIFj02SNMPcEZ+cSrkLuL6qjgOfBS6vqgNsHKl7WkWS9tjMc+IASfYB9wJfrqpbNvn9O4B7qurd\nm/zOD2eRpG2oqpmnrOc9Er8dePrkgE9e8Dzht4CnzjAQb1UcOnRo8DGcLTfnwrlwLs58m9fMFzaT\nHAQ+DjyZ5HGggJuAjyU5ALwGPAdcO/deJUkLMc+7Ux4E3rTJr76y+OFIkrbCKzb30Gg0GnoIZw3n\nYsq5mHIutm6uFzZ3tIOkdnsfknSuSUIt8IVNSdJZyIhLUmNGXJIaM+IaxMrKKklIwsrK6tDDkdry\nhU0NIgkblxwAZEsXN0jLwBc2JWkJGHFJasyIS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0Z\ncUlqzIhLUmNGXJIaM+KS1JgRl6TGjLgkNWbEpYH5BRnaCb8UQoPwSyGmnAttxi+FkKQlYMQlqTEj\nLkmNGXFJasyIS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqbGbEk+xPcn+SI0meTPKp\nyfpLktyX5JkkX01y8e4PV5J0spmfYphkBVipqsNJLgQeBa4GPgH8d1X9SZIbgEuq6sZNtvdTDHUK\nP7lvyrnQZhb2KYZVdayqDk+WjwNHgf1shPyOyd3uAD64/eFKkrZjS+fEk6wCB4BvApdV1TpshB64\ndNGDkySd2b557zg5lXIXcH1VHU/yxud8p30OuLa29vryaDRiNBptbZSSdI4bj8eMx+MtbzfXN/sk\n2QfcC3y5qm6ZrDsKjKpqfXLe/BtV9dObbOs5cZ3C88BTzoU2s+hv9rkdePpEwCfuBn53snwN8KUt\njVCStGPzvDvlIPCPwJNsHC4UcBPwMPA3wNuA54EPV9XLm2zvkbhO4dHnlHOhzcx7JO4XJWsQhmvK\nudBm/KJkSVoCRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1JgRl6TG\njLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEjLkmNGXFJasyIS1Jj\nRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1NjMiCe5Lcl6kidOWnco\nyQtJHpvcrtrdYUqSNjPPkfjngCs3WX9zVb1ncvvKgsclSZrDzIhX1QPAS5v8KosfjiRpK3ZyTvy6\nJIeT3Jrk4oWNSJI0t+1G/LPA5VV1ADgG3Ly4IUmS5rVvOxtV1XdP+vHPgXvOdP+1tbXXl0ejEaPR\naDu7laRz1ng8Zjweb3m7VNXsOyWrwD1V9a7JzytVdWyy/Gng56vqY6fZtubZh5ZLEuDE4yIs82PE\nudBmklBVM197nHkknuROYAS8Jcl3gEPA+5McAF4DngOu3dFoJUnbMteR+I524JG4NuHR55Rzoc3M\neyTuFZuS1JgRl6TGjLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEj\nLkmNGXFJasyIS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1JgR\nl6TGjLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqbGZEU9yW5L1JE+c\ntO6SJPcleSbJV5NcvLvDlCRtZp4j8c8BV75h3Y3A16vqJ4H7gT9a9MAkSbPNjHhVPQC89IbVVwN3\nTJbvAD644HFJkuaw3XPil1bVOkBVHQMuXdyQJEnzWtQLm7WgvyNJ2oJ929xuPcllVbWeZAX4rzPd\neW1t7fXl0WjEaDTa5m4l6dw0Ho8Zj8db3i5Vsw+ik6wC91TVuyY/fwb4n6r6TJIbgEuq6sbTbFvz\n7EPLJQnTJ3BhmR8jzoU2k4Sqysz7zXrAJLkTGAFvAdaBQ8DfAX8LvA14HvhwVb18mu2NuE5huKac\nC21mYRFfwECMuE5huKacC21m3oh7xaYkNWbEJakxIy5JjRlxSWrMiO+hlZVVkpCElZXVoYcj6Rzg\nu1P2kO9CmHIuppwLbcZ3p0jSEjDiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS\n1JgRl6TGjLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEjLkmNGXFJ\nasyIS1JjRlySGjPiks4aKyurJCEJKyurQw+nhVTV7u4gqd3eRxdJgBNzEZZ5XpyLKediyrmYSkJV\nZdb9PBKXpMaMuCQ1ZsQlqbF9O9k4yXPA/wKvAT+oqvctYlCSpPnsKOJsxHtUVS8tYjCSpK3Z6emU\nLOBvSJK2aacBLuBrSR5J8slFDEiSNL+dnk45WFUvJvkxNmJ+tKoeeOOd1tbWXl8ejUaMRqMd7laS\nzi3j8ZjxeLzl7RZ2sU+SQ8D3qurmN6z3Yp8JL2SYci6mnIsp52Jq1y/2SXJBkgsny28GrgCe2u7f\nkyRt3U5Op1wGfDFJTf7OX1XVfYsZliRpHn52yh7yqeKUczHlXEw5F1N+dookLQEjLkmNGXFJasyI\nS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1JgRl6TGjLgkNWbE\nJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEjLkmNGXFJasyIS1JjRlySGjPi\nktTYvr3YyYc+dA0AV1wx4tprP7EXu5SkpZCq2t0dJAV/Afwb73znvTz77D/v6v7OZkmAE/Mddnvu\nz2bOxZRzMbXMc7Gyssr6+vP/b11VZdZ2e3Q65RrgN/ZmV5LU0EbA66TbfDwnLkmNGXFJamxHEU9y\nVZJvJfl2khsWNShJ0ny2HfEk5wF/BlwJ/Azw0SQ/taiBSZJm28mR+PuAZ6vq+ar6AfAF4OrFDEuS\nNI+dRPytwL+f9PMLk3WSpD2yJxf7XHTRb/Lqqy9z/vk/tBe7k6SlsZOI/wfw9pN+3j9Zd4pXXrkX\ngCNHTryZf5lN//2dC+diyrmYWua52Pq/77av2EzyJuAZ4FeAF4GHgY9W1dFt/UFJ0pZt+0i8ql5N\nch1wHxvn1m8z4JK0t3b9s1MkSbtn167Y9EKgqSS3JVlP8sTQYxlSkv1J7k9yJMmTST419JiGkuT8\nJA8leXwyH3889JiGluS8JI8luXvosQwpyXNJ/mXy2Hh45v1340h8ciHQt9k4X/6fwCPAR6rqWwvf\nWQNJfgE4Dny+qt499HiGkmQFWKmqw0kuBB4Frl7ix8UFVfX9yetLDwJ/UFUPDj2uoST5NPBe4KKq\n+sDQ4xlKkn8F3ltVL81z/906EvdCoJNU1QPAXP9BzmVVdayqDk+WjwNHWeJrC6rq+5PF89n4f3Fp\nHyNJ9gO/Btw69FjOAmELbd6tiHshkM4oySpwAHho2JEMZ3L64HHgGDCuqqeHHtOA/hT4Q7byGazn\nrgK+luSRJJ+cdWc/xVB7bnIq5S7g+skR+VKqqteq6mfZuMbiF5P80tBjGkKSXwfWJ8/SwnbeLH1u\nOVhV72HjmcnvTU7HntZuRXzuC4G0XJLsYyPgf1lVXxp6PGeDqnoF+Hvg54Yey0AOAh+YnAv+a+D9\nST4/8JgGU1UvTv75XeCLbJyePq3divgjwDuTvCPJDwMfAZb6FWc8wjjhduDpqrpl6IEMKcmPJrl4\nsvwjwK8Ch4cd1TCq6qaqentVXc5GK+6vqt8ZelxDSHLB5JkqSd4MXAE8daZtdiXiVfUqcOJCoCPA\nF5b5QqAkdwL/BPxEku8kWcpvi05yEPg48MuTt089luSqocc1kB8HvjE5J/5N4O6q+oeBx6ThXQY8\ncNLj4p6quu9MG3ixjyQ15gubktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIa+z8zS46e\nC5aY7wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(list(map(get_number_of_activated_antenna, valid_networks)), bins=100)\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2.4126984126984126" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean(list(map(get_number_of_activated_antenna, valid_networks)))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "■-□-■-□-■\n", + "□-■-□-■-□\n", + "■-□-■-□-■\n", + "□-■-□-■-□\n", + "■-□-■-□-■\n", + "\n" + ] + } + ], + "source": [ + "print(pp_network(valid_networks[-1], box_size))" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "63" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(valid_networks)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "18446744073709551616" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "2**64" + ] + } + ], + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/accep_reject.ipynb b/accep_reject.ipynb new file mode 100644 index 0000000..051a66c --- /dev/null +++ b/accep_reject.ipynb @@ -0,0 +1,160 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import is_valid_network, get_number_of_activated_antenna, pp_network\n", + "from numpy import mean\n", + "from random import Random" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "random = Random(666)\n", + "def get_some_possible_network(n, box_size):\n", + " number_of_square = box_size * box_size\n", + " for i in range(n):\n", + " yield random.randint(0, 2 ** number_of_square - 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "box_size = 8\n", + "n = 10**7\n", + "\n", + "valid_networks = []\n", + "for network in get_some_possible_network(n, box_size):\n", + " if is_valid_network(network, box_size):\n", + " valid_networks.append(network)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADrJJREFUeJzt3X+M5Hddx/Hn6zgqkUINbVLlaqsUmiZEqAjHRSEMVmGL\niWf4Q65NAJsY7w+Lxj+wxWi6/xgl0aiISk7PVlQ4YiF6mhKq0IkBoT3SXwh7vSuYcnctNVh+BBQ9\nLm//mPGcLLs7czvf3Zmdz/ORbDrfnc9+5/PJDM/98pmdNlWFJGnx7Zr1BCRJ28PgS1IjDL4kNcLg\nS1IjDL4kNcLgS1IjxgY/yeEkTyV5ZIMx70pyMslDSa7rdoqSpC5McoV/B/D69e5McgNwdVW9CDgI\nvKejuUmSOjQ2+FX1ceArGwzZD7x3OPY+4JIkl3czPUlSV7rYw98DnBo5PjP8niRpjvimrSQ1YncH\n5zgDfP/I8RXD732HJP6LeyRpE6oq055j0iv8DL/WchR4C0CSfcBXq+qp9U5UVQv7dfvtt898Dq5v\n81/DV+hCvk4X/blb9PV1ZewVfpL3AT3g0iRfBG4HLhr8b6IOVdXdSd6Q5DHgm8DNnc1OktSZscGv\nqpsmGHNLN9ORJG0V37TtUK/Xm/UUttSir2+RLfpzt+jr60q63B8a+2BJbefjSRciCYM9/HS6bypN\nKwm1jW/aSpJ2OIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMv\nSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w\n+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY0w+JLUCIMvSY2YKPhJlpIcT3Iiya1r\n3H9pkg8neSjJZ5L8XOczlSRNJVW18YBkF3ACuB54AjgGHKiq4yNjbgeeVVXvSHIZ8ChweVV9e9W5\natzjSbOSBCgg+DrVPElCVWXa80xyhb8XOFlVj1fVWeAIsH/VmC8Bzxnefg7wH6tjL0mard0TjNkD\nnBo5Ps3gl8CoPwU+muQJ4GLgTd1MT5LUlUmCP4l3AA9X1WuTXA38Y5KXVNU3Vg9cXl4+f7vX69Hr\n9TqagiQthn6/T7/f7/y8k+zh7wOWq2ppeHwbUFX1zpExdwO/WVWfGB5/FLi1qj696lzu4WtuuYev\nebWde/jHgBcmuSrJRcAB4OiqMSvATwwndjlwDfCFaScnSerO2C2dqjqX5BbgHga/IA5X1UqSg4O7\n6xDwW8AdSR4GAvxqVT29lROXJF2YsVs6nT6YWzqaY27paF5t55aOJGkBGHxJaoTBl6RGGHxJaoTB\nl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RG\nGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJ\naoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGTBT8JEtJjic5keTWdcb0kjyY5F+T3NvtNCVJ00pVbTwg\n2QWcAK4HngCOAQeq6vjImEuAfwFeV1VnklxWVV9e41w17vGkWUkCFBB8nWqeJKGqMu15JrnC3wuc\nrKrHq+oscATYv2rMTcAHq+oMwFqxlyTN1iTB3wOcGjk+PfzeqGuA5yW5N8mxJG/uaoKSpG7s7vA8\nLwN+HHg28Mkkn6yqxzo6vyRpSpME/wxw5cjxFcPvjToNfLmqvgV8K8k/Ay8FviP4y8vL52/3ej16\nvd6FzViSFly/36ff73d+3knetH0G8CiDN22fBO4HbqyqlZEx1wJ/CCwB3wXcB7ypqj636ly+aau5\n5Zu2mlddvWk79gq/qs4luQW4h8Ge/+GqWklycHB3Haqq40k+AjwCnAMOrY69JGm2xl7hd/pgXuFr\njnmFr3m1nX+WKUlaAAZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWp\nEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZf\nkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhox\nUfCTLCU5nuREkls3GPeKJGeTvLG7KUqSujA2+El2Ae8GXg+8GLgxybXrjPtt4CNdT1KSNL1JrvD3\nAier6vGqOgscAfavMe5twF3Av3c4P0lSRyYJ/h7g1Mjx6eH3zkvyfOBnqupPgHQ3PUlSV7p60/b3\ngdG9faMvSXNm9wRjzgBXjhxfMfzeqJcDR5IEuAy4IcnZqjq6+mTLy8vnb/d6PXq93gVOWZIWW7/f\np9/vd37eVNXGA5JnAI8C1wNPAvcDN1bVyjrj7wD+vqo+tMZ9Ne7xpFkZXK8UEHydap4koaqm3jkZ\ne4VfVeeS3ALcw2AL6HBVrSQ5OLi7Dq3+kWknJUnq3tgr/E4fzCt8zTGv8DWvurrC95O2ktQIgy9J\njTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4\nktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQI\ngy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjTD4ktQIgy9JjZgo+EmWkhxPciLJrWvcf1OS\nh4dfH0/yQ91PVZI0jVTVxgOSXcAJ4HrgCeAYcKCqjo+M2QesVNXXkiwBy1W1b41z1bjHk2YlCVBA\n8HWqeZKEqsq055nkCn8vcLKqHq+qs8ARYP/ogKr6VFV9bXj4KWDPtBOTJHVrkuDvAU6NHJ9m46D/\nPPDhaSYlSere7i5PluS1wM3Aq9Ybs7y8fP52r9ej1+t1OQVJ2vH6/T79fr/z806yh7+PwZ780vD4\nNqCq6p2rxr0E+CCwVFWfX+dc7uFrbrmHr3m1nXv4x4AXJrkqyUXAAeDoqslcySD2b14v9pKk2Rq7\npVNV55LcAtzD4BfE4apaSXJwcHcdAn4DeB7wxxlcJp2tqr1bOXFJ0oUZu6XT6YO5paM55paO5tV2\nbulIkhaAwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8\nSWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqE\nwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWqEwZekRhh8SWrERMFPspTkeJITSW5d\nZ8y7kpxM8lCS67qdpiRpWmODn2QX8G7g9cCLgRuTXLtqzA3A1VX1IuAg8J4tmOvc6/f7s57Cllr0\n9S2yRX/uFn19XZnkCn8vcLKqHq+qs8ARYP+qMfuB9wJU1X3AJUku73SmO8Civ+gWfX2LbNGfu0Vf\nX1cmCf4e4NTI8enh9zYac2aNMZKkGdo96wlspbe+9a08/fTTXHrppdx5552zno4kzVSqauMByT5g\nuaqWhse3AVVV7xwZ8x7g3qr6wPD4OPCaqnpq1bk2fjBJ0pqqKtOeY5Ir/GPAC5NcBTwJHABuXDXm\nKPCLwAeGvyC+ujr20M2EJUmbMzb4VXUuyS3APQz2/A9X1UqSg4O761BV3Z3kDUkeA74J3Ly105Yk\nXaixWzqSpMXQySdtkxxO8lSSR0a+dyTJA8Ovf0vywAY/v2s47mgX8+naNOtLckmSv0mykuSzSV65\nfTOfzJTre8dwXY8k+eskF23fzCezzvpekeT+JA8O//nydX527IcOZ2mza0tyRZKPDZ+7zyT5pe2d\n+WSmee6GY3diWyZ9bV54W6pq6i/gVcB1wCPr3P87wK9v8PO/AvwVcLSL+XT9Nc36gDuBm4e3dwPP\nnfV6ulofcBXwBeCi4fEHgLfMej2TrA+4F3jd8PYNDP7oYPXP7QIeG67zmcBDwLWzXk9Ha/te4Lrh\n7YuBR+dtbdOsb2TsjmvLpOvbTFs6ucKvqo8DX9lgyM8C71/rjiRXAG8A/qyLuWyFza4vyXOBV1fV\nHcPzfLuqvr41s9y8KZ6/rwP/Azw7yW7gu4Enup/hdNZZ35PAJcPb38PgsyOrTfKhw5na7Nqq6ktV\n9dDw9jeAFebwszNTPHc7uS1j17fZtmz53+EneTXwpar6/DpDfg94O/+/wB1lzPp+EPhykjuAlwKf\nBn65qv5rO+c4jY3WV1VfSfK7wBeB/wTuqap/2u45btJtwCeG8w/wo2uMWetDh3u3YW7TmmRt5yX5\nAQZXmfdt+cy6Men6dmpbJlnfptqyHf+2zBtZ/+r+p4CnhlcaGX7tNOuuj8Ev1JcBf1RVL2MQxdu2\na2Id2ej5ewGD/8t8FfB84OIkN23j3KZxGHhbVV3JYA1/PuP5dGnitSW5GLiLQSy+sU3zm9bY9e3w\ntkzy/G2qLVsa/CTPAN7IYG93LT8G/HSSLzCIymuTvHcr59SlCdZ3GjhVVZ8eHt/F4EnaESZY38uB\nT1TV01V1DvgQY64m58grq+pvAarqLta+cj8DXDlyfAXrbB/MmUnWxnAb7i7gL6vq77ZxftOaZH07\nuS2TrG9Tbeky+Gv9Fv1JYKWq1tzXrapfq6orq+oFDD7Q9bGqekuHc+rSZtb3FHAqyTXDb10PfG7r\npjiVC14fgzf69iV5VpIwWN/KFs5xGqvXdzLJawCSXA+cWONnzn/ocPjXRwcYfMhw3mxmbTC4cvxc\nVf3BFs9vWhe8vh3elknWt7m2dPRO8/sYvFn33wz2c//vneM7gF9YNfb7gH9Y4xyvYX7fSd/0+hjs\nrx1j8BceHwIumfV6Ol7f24HPAo8AfwE8c9brmWR9wI8w2LN+EPgk8MPrrG+JwS+2k8Bts15LV2tj\ncAV8bvi6fBB4AFia9Xq6fO5GzrGj2nIBr80LbosfvJKkRvifOJSkRhh8SWqEwZekRhh8SWqEwZek\nRhh8SWqEwZekRhh8SWrE/wK3PGNUV5ra5AAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(list(map(get_number_of_activated_antenna, valid_networks)), bins=100)\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "18.0" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean(list(map(get_number_of_activated_antenna, valid_networks)))" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "□-□-□-□-■-□-□-□\n", + "■-□-□-□-□-■-□-■\n", + "□-■-□-■-□-□-■-□\n", + "□-□-■-□-□-■-□-■\n", + "■-□-□-□-■-□-■-□\n", + "□-■-□-□-□-□-□-□\n", + "□-□-□-□-□-□-■-□\n", + "■-□-□-□-□-■-□-■\n", + "\n" + ] + } + ], + "source": [ + "print(pp_network(valid_networks[0], 8))" + ] + }, + { + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/brute_force.ipynb b/brute_force.ipynb new file mode 100644 index 0000000..37f80f8 --- /dev/null +++ b/brute_force.ipynb @@ -0,0 +1,187 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import is_valid_network, get_number_of_activated_antenna, pp_network\n", + "from numpy import mean\n", + "\n", + "def get_all_possible_network(box_size):\n", + " number_of_square = box_size * box_size\n", + " return range(2**number_of_square)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "box_size = 3\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADJxJREFUeJzt3V+spPVdx/H3h64SKYGQKhzTbXtCGv/EtFlb7c0anWoE\n4p/S2KT2TyL2ouFCLGmMAbnZc2VSLzAkpjcClRqxUZJawLSllU4MmAICK7BsKYlCRdljoyDd9KaB\nrxdnllnZsztzzplznv3uvF/JhOc8Z57z/PgxvPPMM/PMpKqQJPV03tADkCRtnxGXpMaMuCQ1ZsQl\nqTEjLkmNGXFJamxmxJPsT3J/kiNJnkzy+5P1h5K8kOSxye2q3R+uJOlkmfU+8SQrwEpVHU5yIfAo\ncDXw28D3qurm3R+mJGkz+2bdoaqOAccmy8eTHAXeOvl1dnFskqQZtnROPMkqcAB4aLLquiSHk9ya\n5OIFj02SNMPcEZ+cSrkLuL6qjgOfBS6vqgNsHKl7WkWS9tjMc+IASfYB9wJfrqpbNvn9O4B7qurd\nm/zOD2eRpG2oqpmnrOc9Er8dePrkgE9e8Dzht4CnzjAQb1UcOnRo8DGcLTfnwrlwLs58m9fMFzaT\nHAQ+DjyZ5HGggJuAjyU5ALwGPAdcO/deJUkLMc+7Ux4E3rTJr76y+OFIkrbCKzb30Gg0GnoIZw3n\nYsq5mHIutm6uFzZ3tIOkdnsfknSuSUIt8IVNSdJZyIhLUmNGXJIaM+IaxMrKKklIwsrK6tDDkdry\nhU0NIgkblxwAZEsXN0jLwBc2JWkJGHFJasyIS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0Z\ncUlqzIhLUmNGXJIaM+KS1JgRl6TGjLgkNWbEpYH5BRnaCb8UQoPwSyGmnAttxi+FkKQlYMQlqTEj\nLkmNGXFJasyIS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqbGbEk+xPcn+SI0meTPKp\nyfpLktyX5JkkX01y8e4PV5J0spmfYphkBVipqsNJLgQeBa4GPgH8d1X9SZIbgEuq6sZNtvdTDHUK\nP7lvyrnQZhb2KYZVdayqDk+WjwNHgf1shPyOyd3uAD64/eFKkrZjS+fEk6wCB4BvApdV1TpshB64\ndNGDkySd2b557zg5lXIXcH1VHU/yxud8p30OuLa29vryaDRiNBptbZSSdI4bj8eMx+MtbzfXN/sk\n2QfcC3y5qm6ZrDsKjKpqfXLe/BtV9dObbOs5cZ3C88BTzoU2s+hv9rkdePpEwCfuBn53snwN8KUt\njVCStGPzvDvlIPCPwJNsHC4UcBPwMPA3wNuA54EPV9XLm2zvkbhO4dHnlHOhzcx7JO4XJWsQhmvK\nudBm/KJkSVoCRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1JgRl6TG\njLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEjLkmNGXFJasyIS1Jj\nRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1NjMiCe5Lcl6kidOWnco\nyQtJHpvcrtrdYUqSNjPPkfjngCs3WX9zVb1ncvvKgsclSZrDzIhX1QPAS5v8KosfjiRpK3ZyTvy6\nJIeT3Jrk4oWNSJI0t+1G/LPA5VV1ADgG3Ly4IUmS5rVvOxtV1XdP+vHPgXvOdP+1tbXXl0ejEaPR\naDu7laRz1ng8Zjweb3m7VNXsOyWrwD1V9a7JzytVdWyy/Gng56vqY6fZtubZh5ZLEuDE4yIs82PE\nudBmklBVM197nHkknuROYAS8Jcl3gEPA+5McAF4DngOu3dFoJUnbMteR+I524JG4NuHR55Rzoc3M\neyTuFZuS1JgRl6TGjLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEj\nLkmNGXFJasyIS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1JgR\nl6TGjLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqbGZEU9yW5L1JE+c\ntO6SJPcleSbJV5NcvLvDlCRtZp4j8c8BV75h3Y3A16vqJ4H7gT9a9MAkSbPNjHhVPQC89IbVVwN3\nTJbvAD644HFJkuaw3XPil1bVOkBVHQMuXdyQJEnzWtQLm7WgvyNJ2oJ929xuPcllVbWeZAX4rzPd\neW1t7fXl0WjEaDTa5m4l6dw0Ho8Zj8db3i5Vsw+ik6wC91TVuyY/fwb4n6r6TJIbgEuq6sbTbFvz\n7EPLJQnTJ3BhmR8jzoU2k4Sqysz7zXrAJLkTGAFvAdaBQ8DfAX8LvA14HvhwVb18mu2NuE5huKac\nC21mYRFfwECMuE5huKacC21m3oh7xaYkNWbEJakxIy5JjRlxSWrMiO+hlZVVkpCElZXVoYcj6Rzg\nu1P2kO9CmHIuppwLbcZ3p0jSEjDiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS\n1JgRl6TGjLgkNWbEJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEjLkmNGXFJ\nasyIS1JjRlySGjPiks4aKyurJCEJKyurQw+nhVTV7u4gqd3eRxdJgBNzEZZ5XpyLKediyrmYSkJV\nZdb9PBKXpMaMuCQ1ZsQlqbF9O9k4yXPA/wKvAT+oqvctYlCSpPnsKOJsxHtUVS8tYjCSpK3Z6emU\nLOBvSJK2aacBLuBrSR5J8slFDEiSNL+dnk45WFUvJvkxNmJ+tKoeeOOd1tbWXl8ejUaMRqMd7laS\nzi3j8ZjxeLzl7RZ2sU+SQ8D3qurmN6z3Yp8JL2SYci6mnIsp52Jq1y/2SXJBkgsny28GrgCe2u7f\nkyRt3U5Op1wGfDFJTf7OX1XVfYsZliRpHn52yh7yqeKUczHlXEw5F1N+dookLQEjLkmNGXFJasyI\nS1JjRlySGjPiktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIaM+KS1JgRl6TGjLgkNWbE\nJakxIy5JjRlxSWrMiEtSY0Zckhoz4pLUmBGXpMaMuCQ1ZsQlqTEjLkmNGXFJasyIS1JjRlySGjPi\nktTYvr3YyYc+dA0AV1wx4tprP7EXu5SkpZCq2t0dJAV/Afwb73znvTz77D/v6v7OZkmAE/Mddnvu\nz2bOxZRzMbXMc7Gyssr6+vP/b11VZdZ2e3Q65RrgN/ZmV5LU0EbA66TbfDwnLkmNGXFJamxHEU9y\nVZJvJfl2khsWNShJ0ny2HfEk5wF/BlwJ/Azw0SQ/taiBSZJm28mR+PuAZ6vq+ar6AfAF4OrFDEuS\nNI+dRPytwL+f9PMLk3WSpD2yJxf7XHTRb/Lqqy9z/vk/tBe7k6SlsZOI/wfw9pN+3j9Zd4pXXrkX\ngCNHTryZf5lN//2dC+diyrmYWua52Pq/77av2EzyJuAZ4FeAF4GHgY9W1dFt/UFJ0pZt+0i8ql5N\nch1wHxvn1m8z4JK0t3b9s1MkSbtn167Y9EKgqSS3JVlP8sTQYxlSkv1J7k9yJMmTST419JiGkuT8\nJA8leXwyH3889JiGluS8JI8luXvosQwpyXNJ/mXy2Hh45v1340h8ciHQt9k4X/6fwCPAR6rqWwvf\nWQNJfgE4Dny+qt499HiGkmQFWKmqw0kuBB4Frl7ix8UFVfX9yetLDwJ/UFUPDj2uoST5NPBe4KKq\n+sDQ4xlKkn8F3ltVL81z/906EvdCoJNU1QPAXP9BzmVVdayqDk+WjwNHWeJrC6rq+5PF89n4f3Fp\nHyNJ9gO/Btw69FjOAmELbd6tiHshkM4oySpwAHho2JEMZ3L64HHgGDCuqqeHHtOA/hT4Q7byGazn\nrgK+luSRJJ+cdWc/xVB7bnIq5S7g+skR+VKqqteq6mfZuMbiF5P80tBjGkKSXwfWJ8/SwnbeLH1u\nOVhV72HjmcnvTU7HntZuRXzuC4G0XJLsYyPgf1lVXxp6PGeDqnoF+Hvg54Yey0AOAh+YnAv+a+D9\nST4/8JgGU1UvTv75XeCLbJyePq3divgjwDuTvCPJDwMfAZb6FWc8wjjhduDpqrpl6IEMKcmPJrl4\nsvwjwK8Ch4cd1TCq6qaqentVXc5GK+6vqt8ZelxDSHLB5JkqSd4MXAE8daZtdiXiVfUqcOJCoCPA\nF5b5QqAkdwL/BPxEku8kWcpvi05yEPg48MuTt089luSqocc1kB8HvjE5J/5N4O6q+oeBx6ThXQY8\ncNLj4p6quu9MG3ixjyQ15gubktSYEZekxoy4JDVmxCWpMSMuSY0ZcUlqzIhLUmNGXJIa+z8zS46e\nC5aY7wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(list(map(get_number_of_activated_antenna, valid_networks)), bins=100)\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "6.208433278626436" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean(list(map(get_number_of_activated_antenna, valid_networks)))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "■-□-■-□-■\n", + "□-■-□-■-□\n", + "■-□-■-□-■\n", + "□-■-□-■-□\n", + "■-□-■-□-■\n", + "\n" + ] + } + ], + "source": [ + "print(pp_network(valid_networks[-1], box_size))" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "63" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(valid_networks)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "344239" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sum(list(map(get_number_of_activated_antenna, valid_networks)))" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "random.choice?" + ] + } + ], + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/fuck_me.ipynb b/fuck_me.ipynb new file mode 100644 index 0000000..474410e --- /dev/null +++ b/fuck_me.ipynb @@ -0,0 +1,112 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# algorithme\n", + "\n", + "\n", + "1. Posons $X^{(0)}$ un reseau quelconque.\n", + "\n", + "2. copie $X^{(0)}$ dans $X^{(1)}$\n", + "\n", + "3. pour change antenne dans $X^{(1)}$\n", + " \n", + " 1. si `antenne` peut etre allume et $U(0, 1) < 0.5$\n", + " 1. On allume `antenne` dans $X^{(1)}$ \n", + " 2. sinon\n", + " 1. On eteint `antenne` dans $X^{(1)}$ " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import (is_valid_network,\n", + " get_number_of_activated_antenna,\n", + " pp_network,\n", + " switch_antenna_in_network,\n", + " turn_off_antenna,\n", + " can_antenna_be_turned_on,\n", + " turn_on_antenna)\n", + "from numpy import mean\n", + "from random import Random\n", + "from itertools import cycle" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "random = Random(666)\n", + "\n", + "def get_all_antennas(box_size):\n", + " return range(box_size * box_size)\n", + " \n", + "def get_new_network_with_gibbs(network, box_size):\n", + " for antenna in get_all_antennas(box_size):\n", + " if can_antenna_be_turned_on(antenna, network, box_size) and random.choice([True, False]):\n", + " network = turn_on_antenna(antenna, network, box_size)\n", + " else:\n", + " network = turn_off_antenna(antenna, network)\n", + " return network" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "networks = [0]\n", + "box_size = 5\n", + "\n", + "for i in range(100):\n", + " last_network = networks[-1]\n", + " networks.append(get_new_network_with_gibbs(last_network, box_size))\n", + "\n" + ] + } + ], + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/gibbs_sampler.ipynb b/gibbs_sampler.ipynb new file mode 100644 index 0000000..9d2b5a1 --- /dev/null +++ b/gibbs_sampler.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import (is_valid_network,\n", + " get_number_of_activated_antenna,\n", + " can_antenna_be_turned_on,\n", + " pp_network,\n", + " turn_off_antenna,\n", + " turn_on_antenna)\n", + "from numpy import mean\n", + "from random import Random\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "random = Random(666)\n", + "\n", + "def create_new_network(network, box_size):\n", + " number_of_square = box_size * box_size\n", + " for antenna in range(number_of_square):\n", + " if can_antenna_be_turned_on(antenna, network, box_size) and random.choice([True, False]):\n", + " network = turn_on_antenna(antenna, network)\n", + " else:\n", + " network = turn_off_antenna(antenna, network)\n", + " return network\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "box_size = 8\n", + "n = 15000\n", + "\n", + "valid_networks = [0]\n", + "\n", + "for random_antenna in range(n):\n", + " network = valid_networks[-1]\n", + " new_network = create_new_network(network, box_size)\n", + " valid_networks.append(new_network)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEhRJREFUeJzt3VusXGd5xvH/k1ppCqRRRJtMZac4KAQcVBRCMa1SiaG0\nOfQijrgIhkpACRUiCaByg82N3aoqpBIoVFW44BQHgSKDBAltmpg0jBAXEBcSHLBJLFUOsYs3tBxK\nVKlKmrcXs2xPnG17nzxrtr//Txp57Xevb+ab5TX7mfWtU6oKSVKbzuq7A5Kk/hgCktQwQ0CSGmYI\nSFLDDAFJapghIEkNO2UIJPn1JN9O8nCSHyT5u65+fpJdSR5Lcn+S8ybabE2yP8m+JFdN1K9IsifJ\n40luOz1vSZK0UKcMgar6X+ANVfVq4FXAHye5EtgCPFBVLwceBLYCJLkMuAHYAFwL3J4k3dN9Arix\nqi4FLk1y9Uq/IUnSwi1oOKiq/qeb/PWuzc+BTcCOrr4DuL6bvg64q6qeqaoDwH5gY5IBcG5V7e7m\nu3OijSSpBwsKgSRnJXkYOAyMqmovcGFVzQFU1WHggm72tcCTE80PdbW1wMGJ+sGuJknqyZqFzFRV\nzwKvTvKbwP1JhsDx15vw+hOStMosKASOqKr/TnIv8PvAXJILq2quG+r5STfbIeCiiWbrutqJ6s+T\nxECRpCWoqpx6rmMWcnTQbx058ifJbwB/CjwM3AO8o5vt7cDd3fQ9wOYkZye5GLgEeKgbMvplko3d\njuK3TbSZ7434qGLbtm2992FWHi4Ll4XL4uSPpVjIlsDvADu6P9xnAZ+rqn/t9hHsTPJO4AnGRwRR\nVXuT7AT2Ak8DN9Wx3t0M3AGcA9xbVfctqdeSpBVxyhCoqkeBK+ap/wz4kxO0+TDw4Xnq3wF+b/Hd\nlCSdDp4xPOOGw2HfXZgZLotjXBbHuCyWJ0sdRzqdktQs9kuSZlkSaqV3DEuSzlyGgCQ1zBCQpIYZ\nApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEg\nSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIK0ig8F6kpCEwWB9393RGSBV1Xcf\nnidJzWK/pL4lAY58NoKfE01KQlVlMW3cEpCkhp0yBJKsS/Jgkh8keTTJe7v6tiQHk3y3e1wz0WZr\nkv1J9iW5aqJ+RZI9SR5PctvpeUuSpIU65XBQkgEwqKpHkrwI+A6wCXgz8Kuq+thx828AvgC8FlgH\nPAC8rKoqybeBW6pqd5J7gY9X1f3zvKbDQdI8HA7SyZyW4aCqOlxVj3TTTwH7gLVHXnOeJpuAu6rq\nmao6AOwHNnZhcm5V7e7muxO4fjGdlSStrEXtE0iyHrgc+HZXuiXJI0k+leS8rrYWeHKi2aGuthY4\nOFE/yLEwkST1YMEh0A0FfQl4f7dFcDvw0qq6HDgMfPT0dFGSdLqsWchMSdYwDoDPVdXdAFX104lZ\nPgl8tZs+BFw08bt1Xe1E9Xlt37796PRwOGQ4HC6kq5LUjNFoxGg0WtZzLOg8gSR3Av9ZVR+YqA2q\n6nA3/VfAa6vqrUkuAz4PvI7xcM/XOLZj+FvA+4DdwD8D/1BV983zeu4YlubhjmGdzFJ2DJ9ySyDJ\nlcCfA48meZjxGvgh4K1JLgeeBQ4A7waoqr1JdgJ7gaeBmyb+ot8M3AGcA9w7XwBIkqbHM4alVcQt\nAZ2MZwxLkhbFEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZAlJPBoP1\nJCEJg8H6vrujRnntIKknS7kOkNcO0sl47SBJz+MWh07GLQGpJ9PaEnDroR1uCUiSFsUQkKSGGQKS\n1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkN\nO2UIJFmX5MEkP0jyaJL3dfXzk+xK8liS+5OcN9Fma5L9SfYluWqifkWSPUkeT3Lb6XlLkqSFWsiW\nwDPAB6rqlcAfAjcneQWwBXigql4OPAhsBUhyGXADsAG4Frg947taAHwCuLGqLgUuTXL1ir4bSdKi\nnDIEqupwVT3STT8F7APWAZuAHd1sO4Dru+nrgLuq6pmqOgDsBzYmGQDnVtXubr47J9pIknqwqH0C\nSdYDlwPfAi6sqjkYBwVwQTfbWuDJiWaHutpa4OBE/WBXkyT1ZM1CZ0zyIuBLwPur6qkkx9+odEVv\nXLp9+/aj08PhkOFwuJJPL0mr3mg0YjQaLes5FnSj+SRrgH8C/qWqPt7V9gHDqprrhnq+XlUbkmwB\nqqpu7ea7D9gGPHFknq6+GXh9Vb1nntfzRvM643mjea2003mj+c8Ae48EQOce4B3d9NuBuyfqm5Oc\nneRi4BLgoW7I6JdJNnY7it820UaS1INTbgkkuRL4BvAo468TBXwIeAjYCVzE+Fv+DVX1i67NVuBG\n4GnGw0e7uvprgDuAc4B7q+r9J3hNtwR0xnNLQCttKVsCCxoOmjZDQC0wBLTSTudwkCTpDGQISFLD\nDAFJapghIEkNMwQkqWGGgKTnGQzWk4QkDAbr++6OTiMPEZV6MsuHiHpY6erkIaKSpEUxBCSpYYaA\nJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwSkFeBlFrRaedkIaQXM8uUcvGxEO7xshCRp\nUQwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlq2ClDIMmnk8wl2TNR25bk\nYJLvdo9rJn63Ncn+JPuSXDVRvyLJniSPJ7lt5d+KJGmxFrIl8Fng6nnqH6uqK7rHfQBJNgA3ABuA\na4HbM74ICcAngBur6lLg0iTzPackaYpOGQJV9U3g5/P8ar6LFG0C7qqqZ6rqALAf2JhkAJxbVbu7\n+e4Erl9alyVJK2U5+wRuSfJIkk8lOa+rrQWenJjnUFdbCxycqB/sapKkHq1ZYrvbgb+pqkryt8BH\ngXetXLdg+/btR6eHwyHD4XAln16SVr3RaMRoNFrWcyzofgJJXgJ8tapedbLfJdkCVFXd2v3uPmAb\n8ATw9ara0NU3A6+vqvec4PW8n4BWlVm+zr/3E2jH6byfQJjYB9CN8R/xJuD73fQ9wOYkZye5GLgE\neKiqDgO/TLKx21H8NuDuxXRUkrTyTjkclOQLwBB4cZIfMf5m/4YklwPPAgeAdwNU1d4kO4G9wNPA\nTRNf6W8G7gDOAe49ckSRJKk/3l5SWgGzPEzjcFA7vL2kJGlRDAFJapghIEkNMwQkqWGGgCQ1zBCQ\npIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlq\nmCEgacUNButJQhIGg/V9d0cn4T2GpRUwy/f+7eMew96juB/eY1iStCiGgCQ1zBCQpIYZApLUMENA\nOs7kkS0e3aIznUcHScd57pEtsJCjW2b5qB2PDmqHRwdJkhbFEJCkhp0yBJJ8Oslckj0TtfOT7Ery\nWJL7k5w38butSfYn2Zfkqon6FUn2JHk8yW0r/1YkSYu1kC2BzwJXH1fbAjxQVS8HHgS2AiS5DLgB\n2ABcC9ye8eAgwCeAG6vqUuDSJMc/pyRpyk4ZAlX1TeDnx5U3ATu66R3A9d30dcBdVfVMVR0A9gMb\nkwyAc6tqdzffnRNtJEk9Weo+gQuqag6gqg4DF3T1tcCTE/Md6mprgYMT9YNdTZLUozUr9DwrfvzX\n9u3bj04Ph0OGw+FKv4QkrWqj0YjRaLSs51jQeQJJXgJ8tape1f28DxhW1Vw31PP1qtqQZAtQVXVr\nN999wDbgiSPzdPXNwOur6j0neD3PE1BvPE/A8wRWq9N5nkC6xxH3AO/opt8O3D1R35zk7CQXA5cA\nD3VDRr9MsrHbUfy2iTaSpJ6ccjgoyReAIfDiJD9i/M3+I8AXk7yT8bf8GwCqam+SncBe4Gngpomv\n9DcDdwDnAPdW1X0r+1YkSYvlZSOk4zgc5HDQauVlIyRJi2IISFLDDAFJapghIEkNMwQkqWGGgCQ1\nzBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMM\nAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAGdcQaD9SQ5+hgM1vfdJS3A5P+b\n/2fTk6rquw/Pk6RmsV9aHZIAk+tPWMz6tJT2z22zsNc709r0/ZoaL8OqymLaLGtLIMmBJN9L8nCS\nh7ra+Ul2JXksyf1JzpuYf2uS/Un2JblqOa8tSVq+5Q4HPQsMq+rVVbWxq20BHqiqlwMPAlsBklwG\n3ABsAK4Fbs84+iVJPVluCGSe59gE7OimdwDXd9PXAXdV1TNVdQDYD2xEktSb5YZAAV9LsjvJu7ra\nhVU1B1BVh4ELuvpa4MmJtoe6miSpJ2uW2f7Kqvpxkt8GdiV5jOfuUWOenyVJM2JZIVBVP+7+/WmS\nrzAe3plLcmFVzSUZAD/pZj8EXDTRfF1Xm9f27duPTg+HQ4bD4XK6KklnnNFoxGg0WtZzLPkQ0SQv\nAM6qqqeSvBDYBfw18EbgZ1V1a5IPAudX1ZZux/DngdcxHgb6GvCy+Y4F9RBRLYeHiPbTpu/X1NIO\nEV3OlsCFwJeTVPc8n6+qXUn+DdiZ5J3AE4yPCKKq9ibZCewFngZu8i+9JPXLk8V0xnFLoJ82fb+m\nejhZTJK0uhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhS\nwwwBSavSYLCeJEcfg8H6vru0KnkVUZ1xvIpoP22m/ZrL/X8+E3kVUUnSohgCktQwQ0AzzXFf6fRy\nn4Bm2vLH5xfWZmVfc7bG6t0n0A73CUiSFsUQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSp\nYYaApGZ4BvrzGQKaqskPoR9ATdvc3BOMzzIeP8Y/t23qIZDkmiQ/TPJ4kg9O+/XVr8kPoR9AqX9T\nDYEkZwH/CFwNvBJ4S5JXTLMPq81oNOq7CzNk1HcHZsio7w7MkFHfHVjVpr0lsBHYX1VPVNXTwF3A\npin3YVUxBCaN+u7ADBn13YEZMuq7A6vatENgLfDkxM8Hu5okzaQzfT/WzO4YPrLQ9+zZ03dXdAJn\n+odDgqXtx1pNn42p3k8gyR8A26vqmu7nLUBV1a3Hzdf2RcElaYkWez+BaYfArwGPAW8Efgw8BLyl\nqvZNrROSpKPWTPPFqur/ktwC7GI8FPVpA0CS+jOTt5eUJE3HTO0Y9kSyY5IcSPK9JA8neajv/kxT\nkk8nmUuyZ6J2fpJdSR5Lcn+S8/rs47ScYFlsS3IwyXe7xzV99nFakqxL8mCSHyR5NMn7unpz68Y8\ny+K9XX3R68bMbAl0J5I9znh/wX8Au4HNVfXDXjvWkyT/Drymqn7ed1+mLckfAU8Bd1bVq7rarcB/\nVdXfd18Qzq+qLX32cxpOsCy2Ab+qqo/12rkpSzIABlX1SJIXAd9hfJ7RX9DYunGSZfFmFrluzNKW\ngCeSPVeYrf+fqamqbwLHh98mYEc3vQO4fqqd6skJlgWM14+mVNXhqnqkm34K2Aeso8F14wTL4sg5\nV4taN2bpj4wnkj1XAV9LsjvJX/bdmRlwQVXNwfgDAFzQc3/6dkuSR5J8qoXhj+MlWQ9cDnwLuLDl\ndWNiWXy7Ky1q3ZilENBzXVlVVwB/BtzcDQvomNkYx+zH7cBLq+py4DDQ2rDQi4AvAe/vvgUfvy40\ns27MsywWvW7MUggcAn534ud1Xa1JVfXj7t+fAl9mPFzWsrkkF8LR8dCf9Nyf3lTVT+vYzrxPAq/t\nsz/TlGQN4z96n6uqu7tyk+vGfMtiKevGLIXAbuCSJC9JcjawGbin5z71IskLuoQnyQuBq4Dv99ur\nqQvPHdu8B3hHN/124O7jG5zBnrMsuj90R7yJttaNzwB7q+rjE7VW143nLYulrBszc3QQjA8RBT7O\nsRPJPtJzl3qR5GLG3/6L8Ql9n29pWST5AjAEXgzMAduArwBfBC4CngBuqKpf9NXHaTnBsngD4zHg\nZ4EDwLuPjImfyZJcCXwDeJRjd4b5EOMrD+ykoXXjJMvirSxy3ZipEJAkTdcsDQdJkqbMEJCkhhkC\nktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWH/DzLuXOMKmYfLAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(list(map(get_number_of_activated_antenna, valid_networks)), bins=100)\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "15.072261849210053" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean(list(map(get_number_of_activated_antenna, valid_networks)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "□-■-□-□-■\n", + "□-□-■-□-□\n", + "□-□-□-□-□\n", + "□-■-□-■-□\n", + "■-□-□-□-■\n", + "\n", + "□-□-□-□-■\n", + "□-□-□-■-□\n", + "■-□-■-□-■\n", + "□-□-□-□-□\n", + "■-□-■-□-□\n", + "\n", + "□-□-■-□-□\n", + "□-■-□-□-□\n", + "□-□-■-□-■\n", + "□-□-□-■-□\n", + "■-□-□-□-□\n", + "\n", + "□-□-□-□-□\n", + "□-□-□-□-□\n", + "■-□-■-□-■\n", + "□-□-□-■-□\n", + "□-□-■-□-■\n", + "\n", + "■-□-□-□-□\n", + "□-□-□-■-□\n", + "■-□-□-□-□\n", + "□-■-□-■-□\n", + "■-□-□-□-■\n", + "\n", + "□-■-□-□-■\n", + "□-□-□-□-□\n", + "■-□-□-□-□\n", + "□-■-□-■-□\n", + "■-□-■-□-□\n", + "\n", + "□-□-□-□-■\n", + "□-□-□-□-□\n", + "■-□-□-□-□\n", + "□-■-□-□-□\n", + "□-□-□-□-□\n", + "\n", + "■-□-□-□-□\n", + "□-■-□-□-□\n", + "□-□-■-□-■\n", + "□-■-□-■-□\n", + "□-□-□-□-■\n", + "\n", + "□-□-□-■-□\n", + "□-■-□-□-□\n", + "■-□-□-□-■\n", + "□-□-□-□-□\n", + "□-■-□-□-■\n", + "\n", + "■-□-□-□-□\n", + "□-■-□-■-□\n", + "■-□-■-□-■\n", + "□-□-□-■-□\n", + "□-□-□-□-□\n", + "\n" + ] + } + ], + "source": [ + "for network in valid_networks[-10:]:\n", + " print(pp_network(network, box_size))" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "ValueError", + "evalue": "Cannot take a larger sample than population when 'replace=False'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mchoice\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m111\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreplace\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;32mmtrand.pyx\u001b[0m in \u001b[0;36mmtrand.RandomState.choice (numpy/random/mtrand/mtrand.c:11010)\u001b[1;34m()\u001b[0m\n", + "\u001b[1;31mValueError\u001b[0m: Cannot take a larger sample than population when 'replace=False'" + ] + } + ], + "source": [ + "random.choice([1,2,3,4], 111, replace=False)" + ] + }, + { + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/importance_sampling.ipynb b/importance_sampling.ipynb new file mode 100644 index 0000000..cf7f74f --- /dev/null +++ b/importance_sampling.ipynb @@ -0,0 +1,179 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline\n", + "from network.network import (is_valid_network,\n", + " create_network,\n", + " get_number_of_activated_antenna,\n", + " get_activated_antennas,\n", + " can_antenna_be_turned_on,\n", + " pp_network,\n", + " turn_off_antenna,\n", + " turn_on_antenna)\n", + "\n", + "from network.impSamp import AntennaGenerator\n", + "from numpy import mean\n", + "from random import Random" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "random = Random(666)\n", + "\n", + "def get_random_network(box_size, random):\n", + " antenna_generator = AntennaGenerator(box_size, random)\n", + " antennas = []\n", + " antenna, network_prob = antenna_generator.get_random_antenna()\n", + "\n", + " while antenna != -1:\n", + " antennas.append(antenna)\n", + " antenna, p = antenna_generator.get_random_antenna()\n", + " network_prob *= p\n", + " \n", + " return create_network(antennas), network_prob\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "box_size = 4\n", + "antannas = [0,1,2,3]\n", + "\n", + "for i in range(1000):\n", + " network = get_random_network(box_size)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "box_size = 5\n", + "for i in range(100000):\n", + " network, p = get_random_network(box_size)\n", + "# print(\"network\", network)\n", + "# print(pp_network(network, box_size))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds = list(reversed(range(10)))" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "all_networks = set()\n", + "\n", + "\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "5.582214353409138" + ] + }, + "execution_count": 45, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "get_estimation(1000000, random)" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "25318" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(all_networks)" + ] + } + ], + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/makefile b/makefile new file mode 100644 index 0000000..eeab73d --- /dev/null +++ b/makefile @@ -0,0 +1,2 @@ +test: + nosetests tests diff --git a/network/__init__.py b/network/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/network/acceptance_rejection.py b/network/acceptance_rejection.py new file mode 100644 index 0000000..400543f --- /dev/null +++ b/network/acceptance_rejection.py @@ -0,0 +1,22 @@ +from network.network import (get_number_of_activated_antenna, + is_valid_network) + + +def _get_random_networks(box_size, random): + '''yields random networks, broken or not''' + number_of_square = box_size * box_size + while True: + yield random.randint(0, 2 ** number_of_square - 1) + + +def get_random_networks(box_size, random): + for network in _get_random_networks(box_size, random): + if is_valid_network(network, box_size): + yield network + + +def estimate_number_of_activated_antenna(n, box_size, random): + estimate = 0 + for _, network in zip(range(n), get_random_networks(box_size, random)): + estimate += get_number_of_activated_antenna(network) + return estimate / n diff --git a/network/brute_force.py b/network/brute_force.py new file mode 100644 index 0000000..b9a493f --- /dev/null +++ b/network/brute_force.py @@ -0,0 +1,27 @@ +from network.network import is_valid_network, get_number_of_activated_antenna + +def _get_all_possible_network(box_size): + '''returns all networks, broken or not''' + return range( 2**(box_size*box_size) ) + + +def get_all_possible_network(box_size): + if box_size < 1: + return + for network in _get_all_possible_network(box_size): + if is_valid_network(network, box_size): + yield network + + +def _get_random_networks(n, box_size, random): + all_networks = list(get_all_possible_network(box_size)) + n = min(n, len(all_networks)) + return random.sample(all_networks, n), n + + +def estimate_number_of_activated_antenna(n, box_size, random): + networks, n = _get_random_networks(n, box_size, random) + estimate = 0 + for network in networks: + estimate += get_number_of_activated_antenna(network) + return estimate / n diff --git a/network/gibbs.py b/network/gibbs.py new file mode 100644 index 0000000..151d033 --- /dev/null +++ b/network/gibbs.py @@ -0,0 +1,36 @@ +from network.network import (get_number_of_activated_antenna, + can_antenna_be_turned_on, + turn_off_antenna, + turn_on_antenna) + +def _flip_coin(random): + return random.random() < 0.5 + +def create_new_network(network, box_size, random): + number_of_square = box_size * box_size + for antenna in range(number_of_square): + if can_antenna_be_turned_on(antenna, network, box_size) and _flip_coin(random): + network = turn_on_antenna(antenna, network) + else: + network = turn_off_antenna(antenna, network) + return network + + +def get_pseudo_random_networks(box_size, random): + network = 0 + while True: + network = create_new_network(network, box_size, random) + yield network + +def get_random_networks(box_size, random): + time_between_sample = 20 + for i, network in enumerate(get_pseudo_random_networks(box_size, random), 1): + if i % time_between_sample == 0: + yield network + + +def estimate_number_of_activated_antenna(n, box_size, random): + estimate = 0 + for _, network in zip(range(n), get_random_networks(box_size, random)): + estimate += get_number_of_activated_antenna(network) + return estimate / n diff --git a/network/importance_sampling.py b/network/importance_sampling.py new file mode 100644 index 0000000..65c3c58 --- /dev/null +++ b/network/importance_sampling.py @@ -0,0 +1,92 @@ +from random import Random +from network.network import create_network, get_number_of_activated_antenna + + +class AntennaGenerator(object): + + def __init__(self, box_size, random): + self.box_size = box_size + self._antennas = list(range(box_size * box_size)) + self.random = random + + def get_random_antenna(self): + idx, prob = self.get_random_index_and_its_probability() + if idx == -1: + self._antennas = [] + return -1, prob + choosen_antenna = self._antennas[idx] + self._remove_antenna(idx) + return choosen_antenna, prob + + +class SuccesiveAntennaGenerator(AntennaGenerator): + def _remove_antenna(self, idx): + ant = self._antennas[idx] + if ant == -1: + self._antennas = [] + return ant + if (ant + 1) % self.box_size: + self._antennas = self._antennas[idx + 2:] + else: + self._antennas = self._antennas[idx + 1:] + + if ant + self.box_size < self.box_size * self.box_size: + self._antennas.remove(ant + self.box_size) + + def get_random_index_and_its_probability(self): + idx = self.random.randint(-1, len(self._antennas)-1) + prob = (len(self._antennas) + 1) ** -1 + return idx, prob + + +class BadAntennaGenerator(SuccesiveAntennaGenerator): + '''this antenna generator does not return 'real' probabilities + and thus cannot be used in importance_sampling''' + def _remove_antenna(self, idx): + ant = self._antennas[idx] + to_remove = {ant, ant + self.box_size, ant - self.box_size} + if (ant + 1) % self.box_size: + to_remove.add(ant + 1) + + if ant % self.box_size: + to_remove.add(ant - 1) + self._antennas = [a for a in self._antennas if a not in to_remove] + + +class BetterSuccesiveAntennaGenerator(SuccesiveAntennaGenerator): + + def _get_random_index(self, max_number_inclusive): + if max_number_inclusive < 1: + return self.random.randint(-1, max_number_inclusive) + mean_ = max_number_inclusive/2 + x = int(self.random.expovariate(mean_**-1)) + if x > max_number_inclusive: + return -1 + return x + + def get_random_index_and_its_probability(self): + idx = self._get_random_index(len(self._antennas)-1) + prob = (len(self._antennas) + 1) ** -1 + return idx, prob + + +def get_random_network(box_size, random, Generator): + antenna_generator = Generator(box_size, random) + antennas = [] + antenna, network_prob = antenna_generator.get_random_antenna() + + while antenna != -1: + antennas.append(antenna) + antenna, p = antenna_generator.get_random_antenna() + network_prob *= p + + return create_network(antennas), network_prob + + +def estimate_number_of_activated_antenna(n, box_size, random, Generator): + ret, sum_of_probs = 0, 0 + for i in range(n): + network, p = get_random_network(box_size, random, Generator) + ret += get_number_of_activated_antenna(network) / p + sum_of_probs += p**-1 + return ret / sum_of_probs diff --git a/network/network.py b/network/network.py new file mode 100644 index 0000000..c0a3080 --- /dev/null +++ b/network/network.py @@ -0,0 +1,97 @@ +from collections import deque + + +def _get_print_symbol(network, position): + black_square = '\u25A0' + white_square = '\u25A1' + activated_antennas = set(get_activated_antennas(network)) + if position in activated_antennas: + return black_square + return white_square + + +def pp_network(network, box_size): + ret = "" + for line in range(box_size): + for column in range(box_size): + symbol = _get_print_symbol(network, line * box_size + column) + ret += symbol + '-' + ret = ret[:-1] + '\n' + return ret + + +def get_activated_antennas(network): + i = 0 + while network > 0: + if network % 2: + yield i + network = network >> 1 + i += 1 + + +def create_network(activated_antennas): + return sum(2**activated_antenna for activated_antenna in activated_antennas) + + +def _are_two_consicutive_horizontaly_antennas_activated(network, box_size): + predecessor = -2 + for antenna in get_activated_antennas(network): + if antenna == predecessor + 1: + return True + if (antenna + 1) % box_size: + predecessor = antenna + return False + + +def _are_two_consicutive_verticaly_antennas_activated(network, box_size): + invalid_entennas = set() + for antenna in get_activated_antennas(network): + if antenna in invalid_entennas: + return True + invalid_entennas.add(antenna + box_size) + + return False + + +def is_valid_network(network, box_size): + if _are_two_consicutive_horizontaly_antennas_activated(network, box_size): + return False + + if _are_two_consicutive_verticaly_antennas_activated(network, box_size): + return False + return True + + +def has_antenna(antenna, network): + return (network >> antenna) % 2 + + +def turn_on_antenna(antenna, network): # error with borders + if has_antenna(antenna, network): + return network + return network + 2**antenna + + +def turn_off_antenna(antenna, network): + if has_antenna(antenna, network): + return network - 2**antenna + return network + + +def can_antenna_be_turned_on(antenna, network, box_size): + for ant in get_activated_antennas(network): + if ant == antenna + 1 and ant % box_size: + return False + if ant == antenna - 1 and antenna % box_size: + return False + if ant == antenna + box_size: + return False + if ant == antenna - box_size: + return False + if ant > antenna + box_size: + return True + return True + + +def get_number_of_activated_antenna(network): + return len(list(get_activated_antennas(network))) diff --git a/play_around.ipynb b/play_around.ipynb new file mode 100644 index 0000000..6567959 --- /dev/null +++ b/play_around.ipynb @@ -0,0 +1,220 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from network.gibbs import *\n", + "from network.network import pp_network, can_antenna_be_turned_on\n", + "from collections import defaultdict\n", + "from random import Random\n", + "from numpy import matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "mc = defaultdict(dict)\n", + "random = Random(66)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "prec_network = 0\n", + "\n", + "for i, network in zip(range(10**5), get_random_networks(2, random)):\n", + " if network not in mc[prec_network]:\n", + " mc[prec_network][network] = 0\n", + " mc[prec_network][network] += 1\n", + " prec_network = network" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "defaultdict(dict,\n", + " {0: {0: 964, 1: 3671, 2: 1844, 4: 1759, 6: 1835, 8: 938, 9: 3387},\n", + " 1: {0: 891, 1: 3646, 2: 1832, 4: 1839, 6: 1729, 8: 938, 9: 3663},\n", + " 2: {0: 1810, 2: 3548, 4: 3517, 6: 3611, 8: 1811},\n", + " 4: {0: 1779, 2: 3519, 4: 3490, 6: 3486, 8: 1796},\n", + " 6: {0: 1885, 2: 3554, 4: 3465, 6: 3475, 8: 1757},\n", + " 8: {0: 3633, 1: 3553, 8: 3654, 9: 3572},\n", + " 9: {0: 3436, 1: 3668, 8: 3518, 9: 3527}})" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mc" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def get_prob(from_, to):\n", + " return mc[from_][to]/ sum(mc[from_].values())\n", + "\n", + "def get_matrix(mc):\n", + " n = len(mc)\n", + " mat = [[] for i in range(n)]\n", + " for i, from_ in enumerate(mc):\n", + " for j in enumerate(mc):\n", + " pass\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for key in mc:\n", + " print(key)\n", + " print( pp_network(key, 2))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "□-□\n", + "■-□\n", + "\n", + "1 True\n" + ] + } + ], + "source": [ + "ant = 1\n", + "print(get_prob(0, 2))\n", + "print(ant, can_antenna_be_turned_on(ant, 4, 2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "0.5**4" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "m = matrix([[1, 2], [3, 4]])\n", + "v = matrix([[1, 3]])" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from numpy import std" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "42.73458667121028" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "std([1432, 1376, 1461, 1436, 1436, 1363, 1496])" + ] + }, + { + "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.4.3+" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..14f517d --- /dev/null +++ b/requirements.txt @@ -0,0 +1,35 @@ +backports-abc==0.4 +cycler==0.9.0 +decorator==4.0.6 +ipykernel==4.2.1 +ipython==4.0.1 +ipython-genutils==0.1.0 +Jinja2==2.8 +jsonschema==2.5.1 +jupyter-client==4.1.1 +jupyter-core==4.0.6 +MarkupSafe==0.23 +matplotlib==1.5.0 +mistune==0.7.1 +nbconvert==4.1.0 +nbformat==4.0.1 +nose==1.3.7 +notebook==4.0.6 +numpy==1.10.2 +path.py==8.1.2 +pexpect==4.0.1 +pickleshare==0.5 +ptyprocess==0.5 +py==1.4.31 +Pygments==2.0.2 +pyparsing==2.0.6 +pytest==2.8.5 +python-dateutil==2.4.2 +pytz==2015.7 +pyzmq==15.1.0 +simplegeneric==0.8.1 +six==1.10.0 +terminado==0.5 +tornado==4.3 +traitlets==4.0.0 +wheel==0.24.0 diff --git a/tests/test_brute_force.py b/tests/test_brute_force.py new file mode 100644 index 0000000..97155e0 --- /dev/null +++ b/tests/test_brute_force.py @@ -0,0 +1,26 @@ +from unittest import TestCase +from network.brute_force import (get_all_possible_network, + estimate_number_of_activated_antenna) +from random import Random + + +class TestBruteForce(TestCase): + + def test_get_all_possible_network(self): + self.assertEqual(list(get_all_possible_network(0)), []) + + self.assertEqual(list(get_all_possible_network(1)), [0, 1]) + + expt = [0, 1, 2, 4, 6, 8, 9] + self.assertEqual(list(get_all_possible_network(2)), expt) + + def test_estimate_number_of_activated_antenna(self): + random = Random(666) + est = estimate_number_of_activated_antenna(7, 2, random) + self.assertEqual(est, (0 + 1 + 1 + 1 + 1 + 2 + 2) / 7) + + est = estimate_number_of_activated_antenna(100, 2, random) + self.assertEqual(est, (0 + 1 + 1 + 1 + 1 + 2 + 2) / 7) + + est = estimate_number_of_activated_antenna(3, 2, random) + self.assertEqual(est, 1) diff --git a/tests/test_gibbs.py b/tests/test_gibbs.py new file mode 100644 index 0000000..a75d955 --- /dev/null +++ b/tests/test_gibbs.py @@ -0,0 +1,35 @@ +from unittest import TestCase +from random import Random +from math import sqrt +from collections import defaultdict +from network.gibbs import get_random_networks, estimate_number_of_activated_antenna + + +def standard_deviation(lst): + sum_x = 0 + sum_x_2 = 0 + for elm in lst: + sum_x += elm + sum_x_2 += elm * elm + + n = len(lst) + return sqrt(sum_x_2 / n - (sum_x / n)**2) + + +class TestGibbs(TestCase): + + def test_get_random_networks(self): + counter = defaultdict(int) + random, n = Random(555), 1000 + for _, network in zip(range(n), get_random_networks(2, random)): + counter[network] += 1 + + self.assertSetEqual(set(counter.keys()), {0, 1, 2, 4, 6, 8, 9}) + self.assertLess(standard_deviation(counter.values()) / n, 0.1) + + def test_estimate_number_of_activated_antenna(self): + random = Random(666) + + est = estimate_number_of_activated_antenna(100, 2, random) + real = (0 + 1 + 1 + 1 + 1 + 2 + 2) / 7 + assert abs(est - real) < 0.1 diff --git a/tests/test_importance_sampling.py b/tests/test_importance_sampling.py new file mode 100644 index 0000000..dbf86e2 --- /dev/null +++ b/tests/test_importance_sampling.py @@ -0,0 +1,121 @@ +from unittest import TestCase +from random import Random +from network.importance_sampling import (AntennaGenerator, + SuccesiveAntennaGenerator, + BadAntennaGenerator, + BetterSuccesiveAntennaGenerator, + estimate_number_of_activated_antenna) + + +class TestSuccesiveAntennaGenerator(TestCase): + def test_get_random_antenna(self): + self.antenna_genetator = SuccesiveAntennaGenerator(5, Random(678)) + + ant, _ = self.antenna_genetator.get_random_antenna() + self.assertEqual(ant, 2) + expt = [4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24] + self.assertListEqual(self.antenna_genetator._antennas, expt) + + ant, _ = self.antenna_genetator.get_random_antenna() + self.assertEqual(ant, 14) + expt = [15,16,17,18,20,21,22,23,24] + self.assertListEqual(self.antenna_genetator._antennas, expt) + + ant, _ = self.antenna_genetator.get_random_antenna() + self.assertEqual(ant, -1) + expt = [] + self.assertListEqual(self.antenna_genetator._antennas, expt) + + def test_get_random_antenna_bottom_line(self): + self.antenna_genetator = SuccesiveAntennaGenerator(5, Random(607)) + ant, _ = self.antenna_genetator.get_random_antenna() + self.assertEqual(ant, 21) + expt = [23, 24, ] + self.assertListEqual(self.antenna_genetator._antennas, expt) + + def test_get_random_antenna_right_line(self): + self.antenna_genetator = SuccesiveAntennaGenerator(5, Random(637)) + ant, _ = self.antenna_genetator.get_random_antenna() + self.assertEqual(ant, 14) + expt = [15, 16, 17, 18, 20, 21, 22, 23, 24] + self.assertListEqual(self.antenna_genetator._antennas, expt) + + def test_get_random_antenna_last_box(self): + self.antenna_genetator = SuccesiveAntennaGenerator(5, Random(631)) + ant, _ = self.antenna_genetator.get_random_antenna() + self.assertEqual(24, ant) + expt = [] + self.assertListEqual(self.antenna_genetator._antennas, expt) + + def test_probability(self): + self.antenna_genetator = SuccesiveAntennaGenerator(5, Random(678)) + ant, prob = self.antenna_genetator.get_random_antenna() + self.assertAlmostEqual(prob, 26**-1) + self.assertEqual(ant, 2) + + _, prob = self.antenna_genetator.get_random_antenna() + self.assertAlmostEqual(prob, 21**-1) + +class TestImportanceSampling(TestCase): + + def test_SuccesiveAntennaGenerator(self): + random = Random(666) + + est = estimate_number_of_activated_antenna(150, 2, random, SuccesiveAntennaGenerator) + real = (0 + 1 + 1 + 1 + 1 + 2 + 2) / 7 + print(real) + print(est) + assert abs(est - real) > 0.1 + + def test_BadAntennaGenerator(self): + random = Random(666) + + est = estimate_number_of_activated_antenna(2500, 2, random, BadAntennaGenerator) + real = (0 + 1 + 1 + 1 + 1 + 2 + 2) / 7 + assert abs(est - real) > 0.1 + + def test_BetterSuccesiveAntennaGenerator(self): + random = Random(666) + + est = estimate_number_of_activated_antenna(150, 2, random, BetterSuccesiveAntennaGenerator) + real = (0 + 1 + 1 + 1 + 1 + 2 + 2) / 7 + print(real) + print(est) + assert abs(est - real) > 0.1 + + +class TestBadAntennaGenerator(TestCase): + def test_remove_antenna_easy(self): + gen = BadAntennaGenerator(3, Random(666)) + gen._remove_antenna(2) + expt = [0,3,4,6,7,8] + self.assertListEqual(gen._antennas, expt) + + gen = BadAntennaGenerator(3, Random(666)) + gen._remove_antenna(3) + expt = [1, 2, 5, 7, 8] + self.assertListEqual(gen._antennas, expt) + + gen = BadAntennaGenerator(3, Random(666)) + gen._remove_antenna(4) + expt = [0,2,6,8] + self.assertListEqual(gen._antennas, expt) + + gen = BadAntennaGenerator(3, Random(666)) + gen._remove_antenna(7) + expt = [0,1,2,3,5] + self.assertListEqual(gen._antennas, expt) + + def test_remove_antenna(self): + gen = BadAntennaGenerator(3, Random(666)) + gen._remove_antenna(gen._antennas.index(2)) + expt = [0,3,4,6,7,8] + self.assertListEqual(gen._antennas, expt) + + gen._remove_antenna(gen._antennas.index(0)) + expt = [4,6,7,8] + self.assertListEqual(gen._antennas, expt) + + gen._remove_antenna(gen._antennas.index(7)) + expt = [] + self.assertListEqual(gen._antennas, expt) diff --git a/tests/test_network.py b/tests/test_network.py new file mode 100644 index 0000000..cf27c15 --- /dev/null +++ b/tests/test_network.py @@ -0,0 +1,89 @@ +from unittest import TestCase +from network.network import (get_activated_antennas, + create_network, + pp_network, + _are_two_consicutive_horizontaly_antennas_activated, + _are_two_consicutive_verticaly_antennas_activated, + turn_on_antenna, + turn_off_antenna, + can_antenna_be_turned_on, + has_antenna) + + +class TestFunctions(TestCase): + def test_create_network(self): + assert create_network([]) == 0 + assert create_network([0,2,3,7,8,9]) == 909 + + def test_pp_network(self): + a = 111 + expt = """\ +■-■-■-■ +□-■-■-□ +□-□-□-□ +□-□-□-□ +""" + self.assertEqual( pp_network(a, 4), expt) + + + def test_get_activated_antennas(self): + self.assertListEqual(list(get_activated_antennas(0)), []) + self.assertListEqual(list(get_activated_antennas(1)), [0]) + self.assertListEqual(list(get_activated_antennas(2)), [1]) + self.assertListEqual(list(get_activated_antennas(2)), [1]) + + # 909 is written as 1110001101 in binary + expt = [0,2,3,7,8,9] + self.assertListEqual(list(get_activated_antennas(909)), expt) + + def test_are_two_consicutive_horizontaly_antennas_activated(self): + assert _are_two_consicutive_horizontaly_antennas_activated(160, 4) == False + assert _are_two_consicutive_horizontaly_antennas_activated(24, 4) == False + assert _are_two_consicutive_horizontaly_antennas_activated(96, 4) == True + + def test_are_two_consicutive_verticaly_antennas_activated(self): + assert _are_two_consicutive_verticaly_antennas_activated(24, 4) == False + assert _are_two_consicutive_verticaly_antennas_activated(16384, 4) == False + assert _are_two_consicutive_verticaly_antennas_activated(17408, 4) == True + assert _are_two_consicutive_verticaly_antennas_activated(136, 4) == True + assert _are_two_consicutive_verticaly_antennas_activated(1104, 4) == True + + def test_turn_on_antenna(self): + self.assertEqual(turn_on_antenna(5, 111), 111) + self.assertEqual(turn_on_antenna(4, 111), 111+2**4) + + def test_has_antenna(self): + network = create_network([0,2,3,7,8,9]) + assert has_antenna(0, network) + assert has_antenna(3, network) + assert has_antenna(9, network) + + assert not has_antenna(4, network) + assert not has_antenna(10, network) + + def test_turn_off_antenna(self): + network = create_network([0,2,3,7,8,9]) + assert turn_off_antenna(2, network) == 905 + assert turn_off_antenna(7, network) == 781 + assert turn_off_antenna(10, network) == 909 + assert turn_off_antenna(4, network) == 909 + + + def test_can_antenna_be_turned_on(self): + # □-□-□-□-■ + # □-□-□-□-□ + # □-□-■-□-□ + # □-□-□-□-■ + # □-□-□-□-□ + + network = create_network([4, 12, 20]) + assert can_antenna_be_turned_on(0, network, 5) + assert can_antenna_be_turned_on(1, network, 5) + assert can_antenna_be_turned_on(5, network, 5) + assert can_antenna_be_turned_on(19, network, 5) + + assert not can_antenna_be_turned_on(3, network, 5) + assert not can_antenna_be_turned_on(7, network, 5) + assert not can_antenna_be_turned_on(11, network, 5) + assert not can_antenna_be_turned_on(13, network, 5) + assert not can_antenna_be_turned_on(17, network, 5)