From ae7277647cc47a107581badd0d4891ccab5ace66 Mon Sep 17 00:00:00 2001 From: Brett Morris Date: Sat, 1 Jul 2017 12:07:19 -0700 Subject: [PATCH 1/2] Py3 compatibility, PEP8 compliance --- demo.ipynb | 61 +++--- func.py | 112 +++++------ mr_forecast.py | 495 +++++++++++++++++++++++++------------------------ 3 files changed, 344 insertions(+), 324 deletions(-) diff --git a/demo.ipynb b/demo.ipynb index d0e1d8b..1b0215f 100644 --- a/demo.ipynb +++ b/demo.ipynb @@ -21,16 +21,18 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ + "%matplotlib inline\n", + "from __future__ import (print_function, absolute_import,\n", + " division, unicode_literals)\n", "import numpy as np\n", "import mr_forecast as mr\n", - "import matplotlib.pyplot as plt\n", - "%matplotlib inline" + "import matplotlib.pyplot as plt" ] }, { @@ -44,7 +46,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -63,7 +65,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -72,12 +74,12 @@ "name": "stdout", "output_type": "stream", "text": [ - "R = 1.00 (+ 0.12 - 0.10) REarth\n" + "R = 1.04 (+ 0.08 - 0.10) REarth\n" ] } ], "source": [ - "print 'R = %.2f (+ %.2f - %.2f) REarth' % (Rmedian, Rplus, Rminus)" + "print('R = %.2f (+ %.2f - %.2f) REarth' % (Rmedian, Rplus, Rminus))" ] }, { @@ -91,7 +93,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -111,16 +113,16 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEVCAYAAADQC4MUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2UHHWd7/H3F1AejDpGJQ8KTPYqohskPqF3h5XGSxSY\niOi6qEs0Z292lyOKDysJWaMmLAYwQa/HuwcOxw2Kj5H1gYVkFQbMKIYrLi6ZBGOMejMoKkElrXLd\nZVf93j+qOlPT0z1d3V39q6ruz+ucOl3VXVX97Z6kvv17qN/P3B0REZFQDss7ABERGSxKPCIiEpQS\nj4iIBKXEIyIiQSnxiIhIUEo8IiISVKESj5nNNbMxM9tnZreZ2VCT/c4ys71m9n0zuzTx/Hoze8DM\n7o2Xs8JFLyIiaRQq8QBrgDF3PxG4I96exswOB/4BOAt4DvAGM3t2/LIDH3L358XLVwLFLSIiKRUt\n8ZwL3BCv3wCc12CfU4EfuPuku/8XsAV4VeJ1622IIiLSjaIlnnnufiBePwDMa7DP04AfJ7YfiJ+r\nudjMJsxsc7OqOhERyU/wxBO34exusJyb3M+jsXwajecz2xg/1wKLgCXAz4APZha4iIhk4ojQb+ju\nS5u9ZmYHzGy+uz9oZguAhxrs9hPguMT2cUSlHtz90P5m9o/ALU3eRwPUiYh0wN27bs4oWlXbzcCK\neH0FcFODfe4Bnmlmw2b2WOB18XHEyarm1cDuZm/k7oVf1q1bl3sMilNxljVGxZn9kpWiJZ6rgKVm\ntg94WbyNmS00s20A7v474K3ArcAe4HPu/t34+A+Y2S4zmwBOB94Z+gOIiMjsgle1zcbdHwbObPD8\nT4HRxPaXgS832O9NPQ1QRES6VrQSjyRUKpW8Q0hFcWarDHGWIUZQnEVlWdbblYWZ+SB+bhGRbpgZ\n3oedC0REpM8p8YiISFBKPCIiEpQSj4iIBKXEIyIiQSnxiIhIUEo8IiISlBKPiIgEpcQjIiJBKfGI\niEhQSjwiIhKUEo+IiASlxCMiIkEp8YiISFBKPCIiEpQSj4iIBKXEIyKSk23boFqd/ly1Gj3fzwqV\neMxsrpmNmdk+M7vNzIaa7He9mR0ws92dHC8iUgQjI7B27VTyqVaj7ZGRfOPqtUIlHmANMObuJwJ3\nxNuNfAw4q4vjRURyNzQEGzZEyWZyMnrcsCF6vp+Zu+cdwyFmthc43d0PmNl8YNzdT2qy7zBwi7uf\n3O7xZuZF+twiMtgmJ2HRIti/H4aH846mOTPD3a3b8xStxDPP3Q/E6weAeYGPFxEJqlqFTZuipLNp\n08w2n350ROg3NLMxYH6Dl9YmN9zdzazjYkmr49evX39ovVKpUKlUOn0rEZGO1Np0atVrtWq3olS3\njY+PMz4+nvl5i1jVVnH3B81sAbC9g6q2lserqk1EimDbtqgjQTLJVKuwYweMjuYXVzP9WtV2M7Ai\nXl8B3BT4eBGRYEZHZ5ZshoaKmXSyVLQSz1zgRuB4YBI4392rZrYQ+Ki7j8b7fRY4HXgy8BDwPnf/\nWLPjG7yPSjwiIm3KqsRTqMQTihKPiEj7+rWqTURE+pwSj4iIBKXEIyIiQSnxiIhIUEo8IiISlBKP\niIgEpcQjIiJBKfGIiEhQSjwiIhKUEo+IiASlxCMiIkEp8YiISFBKPCIiEpQSj4iIBKXEIyIiQSnx\niIhIUEo8IiISlBKPiIgEpcQjIiJBFSrxmNlcMxszs31mdpuZDTXZ73ozO2Bmu+ueX29mD5jZvfFy\nVpjIRWQQbNsG1er056rV6HlJr1CJB1gDjLn7icAd8XYjHwMaJRUHPuTuz4uXr/QoThEZQCMjsHbt\nVPKpVqPtkZF84yqboiWec4Eb4vUbgPMa7eTudwIHm5zDehCXiAhDQ7BhQ5RsJiejxw0bouclvSPy\nDqDOPHc/EK8fAOZ1cI6LzexNwD3Au9y92uoAEZG0hoZg1SpYtAj271fS6UTwxGNmY8D8Bi+tTW64\nu5uZt3n6a4G/j9cvBz4IrGy04/r16w+tVyoVKpVKm28lIoOoWoVNm6Kks2lTf5d4xsfHGR8fz/y8\n5t7utb13zGwvUHH3B81sAbDd3U9qsu8wcIu7n9zu62bmRfrcIlIOtTadWrKp3+53Zoa7d92cUbQ2\nnpuBFfH6CuCmdg6Ok1XNq4HdzfYVEWnXjh3Tk0ytzWfHjnzjKpuilXjmAjcCxwOTwPnuXjWzhcBH\n3X003u+zwOnAk4GHgPe5+8fM7BPAEqLebfuBCxNtRsn3UYlHRKRNWZV4CpV4QlHiERFpX79WtYmI\nSJ9T4hERkaCUeEREJCglHhERCUqJR0REglLiERGRoJR4REQkKCUeEREJSolHRESCUuIREamjmUZ7\nS4lHRKSOZhrtLSUeEUltUEoCzWYa3bFjMD5/rynxiEhqg1QSSM40umpVtD1In7+XlHhEJLVmJYF+\nnAStfqbRanWwPn8vaVoEEWnb9dfDypXRRXl4OHquWo2qokZHcw0tE61mGp2cjEpCyc8/CDQtgojk\nolqFu+6C5cvh/e+PttNUOZWpfWi2mUYblYQ6VabvJFPuPnBL9LFFpF0HD7pfdFH0ePCg+8qV7suX\nR48HD6Y/ttF2GWT9Gcr2ncTXzu6vwVmcpGyLEo9IZ7ZunX5R3L8/uops3pzu+NqFdf/+Yl9gm6n/\n/O7R9tatnZ+zTN9JVolHbTwi0pFa9dqqVVGVU9pG9kFtH5lNWb6TvmzjMbO5ZjZmZvvM7DYzm/HP\n2MyOM7PtZvYdM7vPzN7WzvEi0r1kY/vw8FRPr1btHVm2j/SLgfxOsig2ZbUAG4HV8fqlwFUN9pkP\nLInX5wDfA05Ke7yrqk2ka51UOWXdntGLaq/QBrWNp1BVbWa2Fzjd3Q+Y2Xxg3N1PanHMTcD/dvc7\n0h6vqjaR8LZti3q9JavjuumC3arLcxlk/Z30WlZVbUVLPAfd/UnxugEP17ab7D8MfA34Y3d/JO3x\nSjwi/aHTdqa0ypYYei2rxHNEFsG0w8zGiKrL6q1Nbri7m1nT7GBmc4DPA29390fqX291/Pr16w+t\nVyoVKpVKy9hFpFiSw9rs3599Sac2RE6jUtUgGB8fZ3x8PPPzFq3EsxeouPuDZrYA2N6kquwxwFbg\ny+7+4Q6OV4lHpA/0usQT6j3Koi97tQE3Ayvi9RXATfU7xFVom4E9yaST9ngR6Q+d9qxrV6PBQqU7\nRUs8VwFLzWwf8LJ4GzNbaGa1QSRGgOXAGWZ2b7ycNdvxItJ/ZhvWZjbtDlMzkN2de6xQVW2hqKpN\nZHC10xuuH3rOZakwvdrM7HB3/323gYSkxDN41DtJktK22+jfzXS5Jx4zOwpYAjwN+Alwr7s/2m1A\nISjxDB79cpV6ZRmmpkiK0LngJHf/JvCD+PE53QYj0iuawEuS1G6Tr25KPGcA/wk8A/gBcLS7355h\nbD2jEs/g0q9cUem3c0Uo8TxAlHB2AvuAyW6DEeml0L9yizzJV1FjCxFXp73hJEPdDPQGnARUgGdm\nMXBcqAUNEjpw8hiMscgDQBY1tqLGJRE0EZwSj6SX10jGRZ7kK3Rsaf8GRf7OBl1WiadlG4+ZHQ0s\nBb7h7g8nnj8WmOfuu3tQEOsptfFISMl2pe98p1jdc0O2ebXTtlLGtrhB6Hodso3nQ8BfAjeZ2dFm\ndriZHe3uDwHP7jYAkX5W3660ePH0YV1qF9+Rkfxj63WbV9qehWXtcVYbULQIf9vCa1UkAt4RP84D\nriSahmA/8Cjw+SyKXaEXVNUmATRrr5iczL8qKc+2lP373SF6zDquvCeH6/dqQkK18QB/lVh/c2L9\nyCwCyGNR4pEQZrsIznbxDaGobV7dxlWEzgl5/217KWTi+T7wAeAc4DV1rz01iyBCL0o8kqd+/lWc\nTBy19VriOHjQfeVK9y1botfTJoV2k1Ge328//23dwyaedwNnEk3U9iXg68CN8fYNWQQRelHikbzk\n8Ys8ZOkm+XlqiWblymh9y5ap9XbiaPWdNfp8ExMevNRRhNJWrwVLPA0PghOA84E7sggi9KLEI3nJ\no4or9AUx+au/lni6LQHMVpKo/zyTk+6LF0fJJ+SFP+/2pRByTTyHDoaXZhFE6EWJRwZN6CqgZDtH\nVm0ejc6TrM676KIo2Tz72e6bN0ev96o6b1AFSTzAsyjZqASpPrQSjwygUI3eIUs8ycRS+3znn599\ndZ5Esko8re7j+SFwgpm9xczebGYvTNlLW0QKJNS9Mc1uCk3ew9PqvZPjtW3bBvffD5dcAi99aXQz\n6erVsHx5tE/tvJdcAu99b/T84x8//XxDQ61v4Gxn9PKijnNXKu1kKeBU4M3AW4g6HByRRfYLvaAS\njwyQkL/mZ+vVVnvvdkofBw+6L13qvnz51PY550RValu3RsvkZFTKqZXmJifdR0c7+3xpSoWDXDoi\n7zYeomq4C4GLgdcAj8sioBCLEo8MkjK2X8xWXVe7AbeWiCoV9xNPjJJRbd/JyfY/XzvtYP3ebbqZ\nYIkHOBo4F5hb9/yxwMnx+kLgdV0HA3OBMaJpFm4DhhrscxywHfgOcB/wtsRr64mma7g3Xs5q8j6Z\n/SFEpDdm66BQ35lg+fLpSardRNBJKaafbxRtJmTiuZap+3eOBg4nmvQN4Pwsgki810Zgdbx+KXBV\ng33mA0vi9TnA94hmQwVYB/xtivfJ5I8gIr3RrMRzzjlRacZ96sJ/551RL7ZaEuikNFemm1TzlFXi\nSTNI6Pfc/dXAnwPvA74K7DGzR4nu5cnSucAN8foNwHn1O7j7g+6+M15/BPgu8LTELl2PnCois8uy\ngb3+XNXqVGeC+g4K11wDy5bBrl1RJ4mJCfibv4Ht26c6TkD7o0GPjs7sSNCsU0KyA8XwcPpOEzIl\nTeJ5BMDdDwA/cvfT3X0R8AR3f23G8cyL3wfgANHApE2Z2TDwPODuxNMXm9mEmW02M01kK6VTxF5T\n9TGNjETJ4XOfi7a7GYm5flTnW2+NHl/ximhKgauvjpYdO+CEE+DTn4Zzz4ULL4SPfARe8AI48sj2\nes51QzOYdi/NfDzfB75INCr1Ue7+xcRrT3X3n7f1hmZjRNVl9WpD8Dwpse/D7j63yXnmAOPA+939\npvi5Y4FaPJcDC9x9ZYNjfd26dYe2K5UKlUqlnY8h0jPtzFuTZ0yXXBK99p73RCWNbuKrnX/Vqubn\nqs13s2MHHHccnHJKlHiOPXYqSY2O9t8cOHkaHx9nfHz80PZll12GZzAfT5p2l2BjtQF7gfnx+gJg\nb5P9HgPcSjxlQ5N9hoHdTV5rv3JTJKAitiE0iinLBvZW56q9f61X28RENDROrc1Heo9QbTzufoW7\n3+7uG9z91e7+UmAV0ajVT+868013M7AiXl8B3FS/g5kZsBnY4+4frnttQWLz1UDpZkeVwdOoag2i\nKqRFi6JSQF4lnaShoSiWWkyQ3U2p9Te4fu5zzb+TZcuiarbrroOtW2HjRrWvlE43WYuMx2oj6k59\nO3XdqYm6a2+L108D/gDspK7bNPAJYBcwQZS05jV5n24Tv0hmGnXlzWqomSw16mmWxU2Us33++vNv\n2TJz5Omi35PUT8j7BtIyL0o8UjS9uqhnHV8thk6nOGikWVfmLVtmVu3VV/dt2VK+m2PLLKvE07Jz\nQT8yMx/Ezy3FNjkZVWNt3gyvec306rW8G8xrDfuhY6p9J/v3R+/drIPD1VcXpyNGPzMzPIPOBWm6\nUzcL4Hgze5GZHd9tECKDLtnG8e1vR88l235q95Tk1a26nftcslLf7nPrrTO7MV99NSxdmm5wTymO\njko8ZnYhcCTRPT5DwB+8rqG/yFTikSJp1n169eqo4bxI3apDabdLebJkNDwcOtrBkVWJp9M2kjPr\nts/Iot4v1ILaeKRAZhuupZtu1WUcHLSmnVGus+p6XubvKxRynvr6VOBq4BrgSuBPswgm1KLEI2XS\n6b0y/TJ8/2yfo9FryfHckufQZHDdyzXxlH1R4pGy6PbXfB43ovai5NDsczR6r/r5eNpJIEW8cbdI\n8i7xLEwsTwOWZxFMqEWJR8ogq1/goYfv71XJoZ3P0U0CGcTpDtLKO/G8imj8tnXx8pksggm1KPFI\nWnnW+3fz3sl2kdHRqUnStmxp7zydyrrk0Mn5OkkgKvHMLveqNuIx1eL1Y7MIJtSixCNplbXeP3n3\n/+Rk9Lh8+dR2iM+QVcmhk79BJwkkr791mTo15JJ4gOuADwGvBRZmEUAeixKPtKOsv4JrowskR0OY\nmIga33v9GbL8zjqdpK3dBJJXAijTj5usEk9b9/GY2QqiqalfDJweP+4G1rv7T1OfKGe6j0faVdb7\nRJJxQ5jPkOe0Dtu2wSOPRNMk1N6rWo1uPp0zp7hTJaSZFqIIcrmPB3gPMCex/WfAE4FVWWTBUAsq\n8UgbylriaTT+W4h7XfKsOipT6aFeGTo1kFNV20Lgn4mmL9gEXBU/f14WwYRalHgkrSJeyNJc2Ovv\ndUkOPNrNZyji91GvjD8UyhJzLonHpy7cJwBLgMOJpqf+WBbBhFqUeCStIjb8prn4t3Pnf6fvX+SL\nZBlKDzVlSOY1WSWedtt4ng1cBBwEPunu3++8ki8/auORssu7TaDIbV55fzftymvk707kNTr1KHAt\n8H+ANWZ2drcBiEj76mcDDXlhrR81utvZPxvNwNrpKNzJjgzDw9Hj2rXFnqE0j5G/89Zu4vm5u+9x\n9y+7+0rg2F4EJSKzy/ri3877Zn1hHxmZfo7ae4yMtH+uHTtmTp2wYUP0vBRHu1Vty4DlwKeBHwEv\nd/dNPYqtZ1TVJmVWrcIFF8A118AJJ0yfRuG++3r7S7lX1UJlqx4bVFlVtbVMPGZ2OfBN4G53/4WZ\nnQisAM4G3urud3UbRGhKPMVRpvrtkGb7XgAWL54+V8/998Nb3gKf+lR5L9hFbjeSSMg2nqOB44FN\nZrYNeB9wAHgrcFq3ASSZ2VwzGzOzfWZ2m5nN+C9kZkeZ2d1mttPM9pjZle0cL8WSZTVLaM3aJtav\n777NYrbvZXQ0KunUqrkmJ6MklFXSaafNJav2mbyqDiUn7XaDI7ph9EzgUuC1WXStS5x7I7A6Xr+U\n+D6hBvsdEz8eQVQaG2nz+HZ7EUoPlaF7biPNusHWj4XW6ZAttftvNm9ufnwvug230703i67AZepO\nPOjIe5DQXizAXmBevD4f2Nti/2OAfwWe087xSjzFU6b7LpKaJc2sBqlcvrz599LLhN3OubuNo4j3\nSklj/Zp4DibWLbldt99hwE7gN8DGDo7v6suXbJW1xFPTLGl2Oyz/bMPchCgltBN/WX84SHtKm3iI\nBhnd3WA5tz5RAA+3ONcT46q2itclntmOB3zdunWHlu3bt3f795AOlb2aJcsST03tIl6pTE3hnKzG\n27q196WEkCWetNatmzml9eRk9Lz0xvbt26ddK0ubeGYNJqoqmx+vL2hV1Rbv917gXe0crxJPcZS5\nmiXrNp7kvps3u7/xjVPjq7nPnNK5V0K38aQ1Oem+ePFU8qnflt7r18SzEbg0Xl/TqHMA8BRgKF4/\nGvg68D/SHu9KPJKRZklz3brOkmmji/jLXx618yRLE71OzO38GOjVvs3Uks2ddyrp5KFfE89c4HZg\nH3BbIsEsBLbF688F/i1u49lFYkqGZsc3eJ8s/gYimdq6NZq8LXlxnpx0P+00P9R+UraqyKSsSkd3\n3hl9H3femX2MMrusEk9bIxf0C91AKr2Qxc2w9ZOm3X8/nH02nHwyPPaxcOSRcPXV5b1JtNsRCu6/\nH5Ytg2uvhTe/GbZuje5pAt2MHEIuE8H1y4JKPJKh5LQDyeqwLVs6+0VfO8/ExFR1Uq3DwfLlM6c7\nqD+26O1jnfaAa9XGU/aOKmVAP1a1hVqUeCRL9Qln5cooQSQ7BrSrdnGemJjZxXrLlpnv22i7iLrp\nAZemV1vZu+YXnRKPEo80kFcpIHnBm+2mz3bPlZw5NPlaFt22QwuVKHVPUe8o8SjxSAN5lgKS1WGd\nJoL6eLdsmVlyqk+kZbnQhvhRUKZEXEZKPEo80kQeF59GVWydJL12b5LUhXZKGasey0aJR4lHZhGy\nFFC7wNW6Qte3+bTzi76oN2/2SpaloLJ2tigTJR4lHmki61JAqwta1he8tPH3w4W2H5LnIFHiUeKR\nBlpdyDq5WOdxcSxLu00WVF1YHko8SjzSQKvE0mkSCXlxHMQL8SAl2jJT4lHikQ51emEPcXFMkxj7\noYotaRATbVkp8SjxSBfaTSKhLo5pkko/tYv002cZBEo8SjzSoXaTSBEvjrN9hjKViLIe4Vt6S4lH\niUc60EkSKeqFvFmprdFNqEuXTr8/6ODB6Pm8P0MzRUz2kl3i0ejUMlD6ZQTjVqM8V6uwfDlccQV8\n5CPw6KPRyNbvfS9885swNhbtV+SRrrsdyVqyl9Xo1Eo8IiVTP3VC/XbNrl1wyikwMQHHHw8XXwzf\n+hYsWQKPf3yxk07N5CQsWgT798PwcN7RSFaJ57AsghGRcHbsmJ5khoai7R07pvapVuG666Kkc8EF\n8KMfRc/v2wc33gjveU/xk061GpV09u+PHqvVvCOSrKjEI9Jn6ktAtZLPeefBnDnRPkWfUC5tqU7C\nUlVbF5R4pJ8l27GqVbjkEvjlL+Gee+Ab34AnPjF6DoqbfPqlLa7fKPF0QYlnsKS9iPXbxa5WSnjp\nS6Ptl7wENm6MSg0At94alYDK+NkkH3059TUwFxgD9gG3AUMN9jkKuBvYCewBrky8th54ALg3Xs5q\n8j7t9SGUUkvbNbffuvAWtRu4lBf92J3azDYCv3D3jWZ2KfAkd1/TYL9j3P23ZnYE8A3gXe6+w8zW\nAb9x9w+1eB8v0ueW3kvbNVddeEWa68uqNjPbC5zu7gfMbD4w7u4nzbL/McDXgBXuvidOPI+4+wdb\nvI8SzwBK2zVXXXhFGuvX7tTz3P1AvH4AmNdoJzM7zMx2xvtsd/c9iZcvNrMJM9tsZvqtKkD6rrnq\nwivSe8FLPGY2Bsxv8NJa4AZ3f1Ji34fdfe4s53oicCuwxt3HzexY4Ofxy5cDC9x9ZYPjfN26dYe2\nK5UKlUqlk48jJZC2a26/deHtt84SEt74+Djj4+OHti+77LK+rWqruPuDZraAqDTTtKotPua9wL+7\n+9V1zw8Dt7j7yQ2OUVXbABn0Xm39kkglf/3axrMR+KW7f8DM1hD1altTt89TgN+5e9XMjiYq8Vzm\n7neY2QJ3/1m83zuBF7n7XzR4HyUeGQjqLCFZ6tfEMxe4ETgemATOjxPMQuCj7j5qZs8FPk7UPnUY\n8El33xQf/wlgCeDAfuDCRJtR8n2UeGRg9LqzRL+VFKW5vuxc4O4Pu/uZ7n6iu7/c3avx8z9199F4\nfZe7P9/dl7j7c2tJJ37tTfFzp7j7eY2Sjsig2LYN7r9/emeJ+++Pns/SyEhUqqp1xKiVskZGsn2f\nbdtmdvaoVrP/PNJ7hUo8IpKdxYth2TJYvToq6axeHW0vXpzt+9QGKV27Nipd9aodKVSCk95T4pFS\n0q/f1u67D7ZujYbJmZyMHrdujZ5PK+33PDQUtSMtWhQ99qIdKVSCkwCyGP6gbAsaMqf0+m14m15q\nNlNpGu0ON5R2OvFudPN5pDto6mslnkEX8mJXVll8R63OEfJHgP7m+VLiUeIR16/f2WSZEGb7nkMN\nRqpSbv6ySjxq45HSKvrwNnm3Q6WZqTSNVt/z6OjMdpahoey7Un/4w1EHieTnWb06el5KJovsVbYF\nlXhKrwy/fssQYytF+gxFimVQoao2JZ5BVpa5ZsreJlG077nR91m0GPtZVomnUCMXhKKRCyQkTbOQ\nrfrvU2PShdOXIxeI9Juit0OVTaPvU/f3lI9KPNI3ijZmmH6JZ6vV96mSZe+pxCNSp2hDqmTVq0wi\ns32fKlmWi0o80lc0DcDgUckynL6cFiEUJZ7+piqXwVK0KtZ+pqo2kQZU5TJ4Qt3AKtlR4uljre6c\nz/vO+qwlq1iGh6d6Oin5iBSLEk8fa9XYXrTG+G6pMV+kHNTG0+daNbarMV5E0lLngi4MUuKB1o3t\naowXkTT6snOBmc01szEz22dmt5lZ09/eZna4md1rZrd0cvygaNXYrsZ4EQmtUIkHWAOMufuJwB3x\ndjNvB/YAyaJLO8f3vVaN7WqMF5E8FKqqzcz2Aqe7+wEzmw+Mu/tJDfZ7OvBxYAPwt+7+yjaPH4iq\ntlb3N+j+BxFpR1+28ZjZQXd/UrxuwMO17br9/gm4AngCcEki8aQ9fiASj4hIlrJKPEdkEUw7zGwM\nmN/gpbXJDXd3M5uRHcxsGfCQu99rZpVm79Ps+Jr169cfWq9UKlQqTU8lIjKQxsfHGR8fz/y8RSvx\n7AUq7v6gmS0AttdXlZnZFcAbgd8BRxGVer7g7m9Kc3x8DpV4RETa1Je92oCbgRXx+grgpvod3P3d\n7n6cuy8CXg981d3flPZ4ERHJV9ESz1XAUjPbB7ws3sbMFppZs4FckkWXhseLiEhxFKqqLRRVtYmI\ntK9fq9pERKTPKfGIiEhQSjwiIhKUEo+IiASlxCMiIkEp8YiISFBKPCIiEpQSj4iIBKXEIyIiQSnx\niIhIUEo8IiISlBKPiIgEpcQjIiJBKfHkYNs2qFanP1etRs+LiPQ7JZ4cjIzA2rVTyadajbZHRvKN\nS0QkBM3Hk5NqFS64AK68Eq67DjZsgKGh6PkdO2B0NNfwRERm0Hw8Xcq7qmtoKEo6p5wCF144lXRU\n8hGRfjewiSfvqq5qNSrpTExEJZ9du6IYaiUfEZF+VaiqNjObC3wOOAGYBM5392qTfQ8H7gEecPdX\nxs+tB/4K+Hm829+5+1caHOsHDzpr18KqVbBpU9gLfi3R1d5z166o5DMxAc99bpgYRETa1a9VbWuA\nMXc/Ebgj3m7m7cAeIJk5HfiQuz8vXmYknZqhoSjpLFoUPYYsZezYMb1Np1byefe7p1cBjo+Phwuq\nC4ozW2WIswwxguIsqqIlnnOBG+L1G4DzGu1kZk8HzgH+EajPvqmycbUalXT2748e69t8eml0dHqb\nzoYNUUkav0UsAAAIcElEQVTnU5+aXgVYln+MijNbZYizDDGC4iyqoiWeee5+IF4/AMxrst//AlYB\nf2jw2sVmNmFmm82saTmmdsEfHo4ekxf8UJIlH4geN2yInhcR6VfBE4+ZjZnZ7gbLucn94v7OMxqg\nzGwZ8JC738vM0s21wCJgCfAz4IPN4ijCBb9W8kkaGlJXahHpb0XrXLAXqLj7g2a2ANju7ifV7XMF\n8Ebgd8BRwBOAL7j7m+r2GwZucfeTG7xPcT60iEiJZNG5oGiJZyPwS3f/gJmtAYbcvWkHAzM7Hbgk\n0attgbv/LF5/J/Aid/+LELGLiEg6RWvjuQpYamb7gJfF25jZQjNrdntnMnN+wMx2mdkEcDrwzp5G\nKyIibStUiUdERPpf0Uo8PWFmf25m3zGz35vZ82fZbzIuMd1rZt8KGWP8/qnijPc9PI7zllDxJd67\nZZxmdpSZ3W1mO81sj5ldWdA4jzOz7fF+95nZ24oWY7zf9WZ2wMx2h4wv8f5p4zzLzPaa2ffN7NKQ\nMcbvPzfuwLTPzG5r1rPVzN4ed2q6z8zeXuA4/y7+3neb2WfM7MiixWlmz4qvRbXlV63+Hw1E4gF2\nA68Gvt5iPyfq3PA8dz+192HNkDZOaHwDbSgt43T3/wDOcPclwHOBM8zstEDx1aT5Pv8LeKe7/zHw\nEuAtZvbsEMHF0v7NPwac1ftwmmoZZzyayD8Qxfkc4A2Bv0tIcRO6mS0mGuHkRcApwDIz+29Bo0wX\n5zDw18Dz405ShwOvDxgjpIjT3b9Xu2kfeAHwW+BLs510IBKPu+91930pd++6x0an0sbZ4gbanksb\np7v/Nl59LNF/mod7GtjM928Zp7s/6O474/VHgO8CC0PEF79n2u/yTuBggJCavX+aOE8FfuDuk+7+\nX8AW4FW9j26aNDehnwTc7e7/4e6/B74GvCZQfDVp4vw10Q+jY8zsCOAY4Cdhwjsk1U39CWcCP3T3\nH8+200AknjY4cLuZ3WNmf513MLOY7QbawjCzw8xsJ9HNwNvdfU/eMc0m/oX5PODufCMpracByQvO\nA/FzIaW5Cf0+4E/jaqRjgFHg6aECjLWM090fJroX8UfAT4Gqu98eLkQg/U39Na8HPtPqpEd0G1VR\nmNkYML/BS+9297TtICPu/jMzeyowZmZ741+amek2zuQNtGZWyTK2uvfp+vt09z8AS8zsicCtZlZx\n9/EMw8zq746ZzQE+D7w9LvlkJqsYey2DOINU+84S59ppwbh7o3v23H2vmX0AuA34f8C99OBHXLdx\nxtV/7wCGgV8B/2RmF7j7p4sUZ+I8jwVeCbRs2+ubxOPuSzM4x8/ix5+b2ZeIqg4yTTwZxPknwLlm\ndg7xDbRm9on6G2i7lcX3mTjXr+Lu8C8ExrM6b3zuruM0s8cAXwA+5e43dR/VdFl+l72UQZw/AY5L\nbB9HVOrJ1Gxxxh0w5iduQn+oyTmuB66Pj7mCqFRRtDhfCNzl7r+Mj/ki0f//TBNPFt9n7Gzg2+7+\n81n2AQazqq1hm4iZHWNmj4/XHwe8nKhBNS8N43T3d7v7ce6+iKhY+9Wsk06bmn2fT6n1gDGzo4Gl\nRL8s89IsTgM2A3vc/cNhQ5oZTs7vn1azOO8Bnmlmw/Gv39cBN4cLC+L3WxGvrwAa/pAws2Pjx+OJ\nOk20rB7KWJo49wIvMbOj43+nZxJ1KAop1fcZewPw2VRndfe+X4j+Yf0Y+HfgQeDL8fMLgW3x+h8B\nO+PlPqK5fAoXZ93+pwM3FzFOop5s/xZ/n7uAVQWN8zSiapadRInxXuCsIsUYb3+WqJ7/0Xj/vyza\ndxlvnw18D/hBTv+H5gK3A/uIqtKGmsT5deA78d/9jALHuTqOczdR4/5jChrn44BfAI9Pc17dQCoi\nIkENYlWbiIjkSIlHRESCUuIREZGglHhERCQoJR4REQlKiUdERIJS4hERkaCUeERKxMyOMLNn5R2H\nSDeUeERSMrOrzezyHp37IjP7tZk9ue75G83s44l5bSokBrRs4ziRwlDiEUnvh8A3e3TubwFfJhqJ\nGIB4ps85wOXu/t346We5+/c7OE6kMJR4RNI7ld7N1XMC0UjoxyeemwMc6+4/TDxXP3x/2uNECqNv\npkUQCeBYd/9FPCfSk4GnEg2U+F0zWwj8T6KBNP/E3S9s89wWHzsMYGb/Hfi/JIahN7NTgX9t9ziR\nolGJRySFeDK7g2Z2IrDc3W8A/gW4KN7lWuDDwBjR5GKd+DFwXDw/kANLiKrSal7g7vd0cJxIoSjx\niKTzIqKL+QqmJuI6gSgZnQCYRzOXvhi4C8DMXhWXhDCzU83sTDObURIysycAB4kSyPHAS9z9m8ys\n2jusw+NECkWJRySdFxBNdHY2U7NVvhb4JPAkojloIJoj6S4zm8/UBFoA57v77cCR8eRjSS9iaubG\nPwJ+Ez9/qGot7kL9vXaPqzGzZ5rZGWa2xMyOauuTi2RMiUcknR8STRr3VuDlZrYC+Hzcw2w38Hsz\n+zPgxe7+U3d/EJhIHH9M/PgIMK/2pJmdBlwJLIuf+oa77zSzi4iqzE6Ln6+QmDa8jeMws5OAXwP/\n6e4749dFcqPOBSIpuPvnE5t31b38FHdfE7cDvaLutdo00b+KH4eAA4nzfoOohFLbfkf8eA1wTeI8\nj3X333VwHMB84EjgGWb2CPA0Mzvc3X/f/BOL9I4Sj0j3rjKzfwaeCawHMLNjgWcBZwCfAv7FzM4A\n/uDuP2p2okbidqKfdBHfT4AqMMfdJ8zsaCUdyZOmvhYpODN7HbDV3TvtLYeZPYOoy/UvgT3u/mhG\n4Ym0TYlHRESCUucCEREJSolHRESCUuIREZGglHhERCQoJR4REQlKiUdERIJS4hERkaCUeEREJCgl\nHhERCer/A6IO23M4mdSJAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEOCAYAAAC0BAELAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2UHHWd7/H3V0IiiQYI4SEIY/SIV1AgykRx3YgjZFeQ\nuzwcDShm41HvYNC7rA/EcNiZu5fohUTC7qq7KOJiwsWFIIuw+ESSGXXdRU2QQECMAdcHIBd84ly4\n7rKi3/vHr9rp6XTPVHVXdf2q+/M6p053V1dXfbszqW/9Hur3M3dHREQkb88qOwAREelNSjAiIlII\nJRgRESmEEoyIiBRCCUZERAqhBCMiIoVQghERkUIowYiISCGUYEREpBAzyg6gTPPnz/eFCxeWHYaI\nSKXcddddP3f3g6fbrq8TzMKFC9m+fXvZYYiIVIqZ/TjNdqoiExGRQijBiIhIIZRgRESkEEowIiJS\nCCUYEREphBKMiPS1detgfHzyuvHxsF46owQjIn1t8WJYtmwiyYyPh9eLF5cbVy/o6/tgRESGhmDT\nppBUVq6Eq64Kr4eGyo6s+lSCEZG+NzQUksuaNeFRySUfUSQYM5tnZpvNbHfyeGCL7b5iZk+Y2e0N\n6z9rZv9mZjuSZVF3IheRXjA+HkouIyPhsbFNRtoTRYIBVgNb3f0oYGvyupmPAstbvHeRuy9Klh1F\nBCkivafW5rJpE1x66UR1mZJM52JJMGcAG5LnG4Azm23k7luBJ7sVlIj0vm3bJre51Npktm0rN65e\nEEsj/6HuvgfA3feY2SFt7OMjZjZKUgJy96dzjVBEetKqVXuvGxpSO0weupZgzGwLcFiTty7JYfcX\nA/8HmAlcDXwIuLRFHMPAMMDAwEAOhxYRkWa6lmDc/ZRW75nZY2a2ICm9LAAez7jvPcnTp83sWuCD\nU2x7NSEJMTg46FmOIyIi6cXSBnMbsCJ5vgK4NcuHk6SEmRmh/ea+XKMTEZHMYkkwlwNLzWw3sDR5\njZkNmtk1tY3M7J+Bm4CTzexhM/vj5K3rzWwnsBOYD3y4q9GLSOVpyJj8RdHI7+6/AE5usn478K66\n10tafP71xUUnIv2gNmRMrUdZffdlaU8UCUZEpGwaMiZ/sVSRiYiUTkPG5EsJRkQkoSFj8qUEIyKC\nhowpghKMiKTWyz2tNGRM/pRgRCS1Xp6ca9WqvdtchoaaDyUj6agXmYikpp5WkoVKMCKSiXpaSVpK\nMCKSSb/0tOrl9qZuUYIRkdT6qadVL7c3dYsSjIik1k89rerbm0ZHJw8jI+mYe/+OWD84OOjbt28v\nOwwRidjoaGhvGhkJpTYBM7vL3Qen204lGBGRFvqlvakoSjAiIk30U3tTUZRgRESaiLW9qUq925Rg\nRESaiPXO/ir1btOd/CIiFVKl0RRUghERqZiqjKagBCMikrNm7STnnx+Weu22nVSld5sSjIhIzpq1\nk9x4I9xwQ+dtJ1Xq3aY2GBGRnDVrJ7nllvBep20nU/Vui62qLIoSjJnNM7PNZrY7eTywyTaLzOxO\nM7vfzO41s3Pq3nuBmX07+fyNZjazu99ARKqoyC6/zdpJ8mg7ibV3WzNRJBhgNbDV3Y8CtiavG/0a\n+FN3fynwBuCvzeyA5L21wF8ln/8V8M4uxCwiFVdkl99m7SRVaTvJjbuXvgC7gAXJ8wXArhSfuQc4\nCjDg58CMZP2rga+mOe4JJ5zgItLfxsbc5893HxkJj2Nj+e2ztq+xMff993efO3fyuryO123Adk9x\njo2lDeZQd98D4O57zOyQqTY2s1cCM4GHgIOAJ9z9meTth4HnFRmsiPSO+mqrkZF82jGatZOcc87E\n89pjrG0neenaaMpmtgU4rMlblwAb3P2Aum1/5e57tcMk7y0AvgascPdvmdnBwJ3u/qLk/SOBL7n7\nsS0+PwwMAwwMDJzw4x//uINvJSJVV6sWi/2mxZikHU25ayUYdz+l1Xtm9piZLUhKLwuAx1tsNxf4\nIvAX7v6tZPXPgQPMbEZSijkCeHSKOK4GroYwXH9730ZEekF9l99aI7zmfclPLI38twErkucrgFsb\nN0h6ht0CbHT3m2rrk/rAceBNU31eRKRRrANa9oooJhwzs4OATcAA8BPgze7+SzMbBN7t7u8ys7cB\n1wL313307e6+w8xeCNwAzAPuBt7m7k9Pd1xNOCYikl3aKrIoEkxZlGBERLLTjJYiIhVRpTleslCC\nEZG29eqJsduqNMdLFkowItK2xhPj+efDWWdNPjG2m3BiSl5Fx1I/dtnoaP492cr6LZVgRKRtjSfG\nG2+E+mbdTq7EY7qq70YsRc7xsngxnH46XHlleF2Lf8aMgpNMmtv9e3XRUDEi+RgZcYfwmOfQK0UM\n4xJrLEXvf/16dzP35cvD/tevb/84pBwqpvSTfJmLEoxI55qdGOsTTqc62dfatXufQMfGwvpuxzKV\nZmOXFZFkli8P8S9Z0tn+lWCUYEQKN9WgjjGUYPI8cRdZwsg7ETZTi3/JknDmX768/X0pwSjBiBSu\n8cQ4NhaSy/DwxOtOT+idJoc8EkO3ShhFqcVbqxZbvjxUl61f397+0iYYNfKL9LgiexA1Tn61bRt8\n4QvwqU+F150MvZLXMC55NJ5XfUiZbdvg4ovhsstC3Bs3whVXhI4Zhc5JkyYL9eqiEozEpoiqkqpf\nfXcqpo4CZcrzbwtVkSnBSPUUlQzyPMl2o70gr2P3e3ItStoEoyoykYgUdcNdnvdYlHl/StZjV71q\nq/LSZKFeXVSCkVjl3R0272qiMqudVOVVPlRFpgQj1VRUMsi7mqioe0JiP7akTzCqIhOJSP0Mi5de\nOlFd1klPnyKqicbHw/TCIyPhsdCeSBEdO08xjbVWmDRZqFcXlWAkNmU2oKdVZsN5J8eO7betcgcE\nUpZgNOGYJhwTyWTdutCoXt9RYHw8lIhWrYr32LXS4dlnw7nnhnW10iJ0J/5WMa1cGUpjeY6gXCTN\naJmCEkzvKPOkJ3FI8zcwPg5nngnPPBNGEv7CF8L6vIfHz2J0NPTuGxkJ1aJVoBktpa/k3XW2L+rH\nu6zo3/Shh8JcNPV/A2edFdbXDA3BhRfCr38dkkx9m1cZyaVX2pNaSlOP1quL2mB6SxHDxFexfrwI\nebRfFP2b1sZB23//8DdQG3Szcay0+fPdTz7ZfdYsn9QTrdvtMVP9HrG1FzVC3ZSVYPpRnt1X++l+\ni+lOaEUMPDlnjvvpp++dAIaH2z+Rjo2577df+BvYb7/myWVsbGJulFmzQiLqZG6Udk31m8d+gaME\nowTTd4pICP1yv0WaE1pev2/tN12+fKLEMTY2uQTS7r7Hxtxnzw77nz178n5qJ+76kYXnznV/+cun\nHlm4rNJEzBc4lUowwDxgM7A7eTywyTaLgDuB+4F7gXPq3vss8G/AjmRZlOa4SjC9o4juq8PD8f4H\nL0KaE1qnCbfxGLWT/OzZocTRaXKpn4umVbI69dSJZFL7PkuXhvVTxVxGaSLWC5yqJZh1wOrk+Wpg\nbZNtXgwclTw/HNgDHOATCeZNWY+rBNM7OrnKbHYCaay/j62KohNT/VZTndA6vaJudaKuzbLY6Yl0\neHjvf7P6uWka48gyN0oZpQmVYPJLMLuABcnzBcCuFJ+5py7hKMFIRxr/Mw8Px93I2olWJ/raCbfZ\nCS2Pq/hmiW39+tAOkkcJJstFRjvz03ezNKE2mHwTzBMNr381zfavBB4AnuUTCWZXUnX2V8CsKT47\nDGwHtg8MDOT0c0sviLU6ogjNqqqmOqEVNU9N3m0waa1dO1FyStOLrNulCfUiy55EtgD3NVnOyJJg\naiUc4MSGdQbMAjYAo2liUglGamKujihKfUIt44S2du3eJcVOe5GlleXfO/bSRBmiSzBTBpGyigyY\nC3wXePMU+3odcHua4yrBiHt/nkD6MaHWZP33ru8UULN+fetOAf0gbYKJ5U7+24AVyfMVwK2NG5jZ\nTOAWYKO739Tw3oLk0YAzCSUjkVT6bVKqIkZsrpKs/94XXRTmsq8fIeCyy8J6mVoUY5GZ2UHAJmAA\n+AmhhPJLMxsE3u3u7zKztwHXErop17zd3XeY2RhwMKGabEfymaemO67GIpN+pHHbptbs97nyyjCc\nywc+UK1BKYuiwS5TUIKRXqUk0r7G8clqr089Fa67rlqDUhZFg12KdEGsg2JmHfwz1u9RhlqV2bJl\nYaTjZcvg4ovhy1/u4UEpi5KmoaZXFzXyS6di7iCgnlKdqR/SJq/fpszux3kemyr1IitrUYKRPMTc\nIyvLvT0xf49uqD8B136LV7zCfcaMyb3IOulKXWYiz/PYSjBKMNJFMd6kWZ8wZs/eu6tts6vXLN8j\n9psBs2oc0aA2TtqcORM3fzYOR1P7XDvTFpSRyPM6thKMEox0SYxX/o1Xp7WhUWpJptnVa9bvUZVq\ntSyJcGwsJOP6arHaCAMnn7z3SAPtfufpEnmRyTuPi6FCEwwwAzgemNfO52NZlGCkU7GeZFuN+zVn\nTr5jjcWYXBtl/W7NRjioXzfV75glnqk+X9TfVfQlGMKd9tcDnwSuAM7Kuo9YFiUY6VTVqolaXb12\n8j1irB5slPbE2rhdrZqscZbMxnHMssaRJnHknbwr0QYDfJQwHMs48FLg1qz7iGVRgpF+UkRpowol\nmJrpEmGzE3CtDaZ+jplam0w73zlrIs8zeVeiFxnwt8A7CHfMDwNfyrqPWBYlGOkXRVS5xFo92Eya\nRNjsBDw8HJb6Lsv77z8xx0wR37lx5s1ayalxXpsyFZlgXkwYsmUD8H5gZdZ9xLIowUi/KKIqryrV\ng50mwnZ643Uab1nTGKSVNsG0NVRMMnbYa4EfuPv9020fKw0VI9L7Ohk2p9WwMUWPRXb++XDDDXDh\nhRNjn0E8Q/3kNhaZmS0HrgSeBi5x9w1mdiJwOnCqu5+QR8BlUIIRmaDxy/ZW5m8yOgpr1sQ59lme\nY5GNAqcBi4AXmNlm4CZgX+DPO4pSoqcxqvpH1vHL+sGqVXuXVIaGik8u4+Oh5FL5sc+mq0MD7q57\nbsDjwAFp6t9iX9QGM70qNeRK56rUKywvsbUlVeH/HDlOOHaYmQ2b2UnAocDD7v5EUQlP4tJsZNl+\nnwujlw0NwcqVoWpm5cqp/517pXQbW8mtpybAmy4DEboifwL4OvBL4DfAFsL9MG9Nk8ViXVSCSa8K\nN9P1mjKurLN0562/sh4bC91oY7vSTqsfS26doMBuykcQ2mQ+BFyX9fMxLUow6eg/Xzm6XVWS9njD\nw3t3oZ01KyxV/tvQRVR6uSUY4NUkvc16bVGCmV4V6oN7WTeTe9oSU/19GSMjIbHUbkLsZhx5yvo7\nx9Zu0215JphPAncDNwBvBw5Ls+MqLEow0+v3/0gxiPHKemzMfb/9QlzNJuXqdN8xltzKjDE2uVeR\nAS8B3gd8BbgT+F+Emy33SbuP2BYlGIldrNWTY2MTJZda1VieJ9kYS25lxhibwtpgwr7ZL2mH+Xja\nA8W4KMFIzGK9Sh4bC9Vjz352KLk0DmuSV+k2xpJboyrEWIRCE0wRCzAP2AzsTh4PbLLN84G7CANt\n3g+8u+69E4CdwIPAx9K0GynBSMxirJ5cu9b9jW+cPKvj+vUh2eQ5GGOr0kHZv0mzaZWXLw8jLJed\n+LspzzaY5cDPgIeBFcm6E4EPA3elOUiqQGAdsDp5vhpY22SbmcCs5PlzgB8Bhyevv1PrkAB8mTCM\njRKMSI7GxiYP+Fg7ya5fn99JfqqSW5mlurVrJ0+nPH+++8qVE79HDKXLbskzwewGFgPzgb9MShc/\nTRLCkjQHSRUI7AIWJM8XALum2f4g4CfA4cn236977y3Ap6Y7phKMtKvsK+kyYym67WG671NW20d9\nMp09233p0r2noe6Xzi95JpiuDBUDPNHw+lcttjsSuBf4NfCeZN0gsKVumyXA7S0+PwxsB7YPDAzk\n+JNLP4mpfaSMWMpueyjr+LXfdskSz7VrdtXkmWD2JCflk4DDgO+m2XGLfW0B7muynJE2wdS9f3hS\nLXZoUsJqTDD/NF08KsFIJ2LqRZQ2ljxKO2V/77KPX5suecmS8v/dy5JngunKUDFZq8iS7a4F3qQq\nMilL2VfyWWPptLTT7dJSY0Ks9WArckbJqaxfH6rFavf99FvbS01uCWavDxQ0VEySsOob+de1OPZ+\nyfMDgR8AxyavtxE6H9Qa+U+b7phKMNKJIq+ks5Y0ssTSSdxltffUjjk8PLkH23THzxJvmrafOXOK\n7eBQFYUlmKKWpNF+K6FTwVZgXrJ+ELgmeb40aX+5J3kcrvv8YFLd9lBS4lI3ZSlM0VfyWfbfTiwx\nlbym00lCzPN3jKljR9kql2DKWJRgpF3dONkU1a5SdhtGK1N9j04SYrdKd/1ECUYJRnpA3iWNmHq/\nNWoVW62do5OTfpbfsUqlu7IowSjBSMUVcTUdazVP4xwzIyOhreX00ztPiCrB5K/IRv63EkZWvh74\nHPCWrPuIZVGCkVjVn0jr7yCvP9GWnRTyVP99ayWI2bNDo34nCbHotqx+lTbBpJkyudFJ7n6uu5/n\n7m8F/rCNfYhEr8wpgeunzV28GC67DC6+OKzPOqVv/feoPa//HjFMc1ybFviss+CKK2D2bJgxA849\nd+9pm4eGYNWqdPvNMv1wT01VHIs0Wah+AT4LvBE4jtBd+e+z7iOWRSUYmUpMV7R59aSq3UdS6+ob\n01X62NjEHDMjI3HFJpNRYAnmAsI9KKcRRkB+Ty6ZTiQytSvYZctgdDQ81l/hdjuWlSthzZrwmCWG\n+u8xPh6mCDObKAmV9Z0a3XADzJwJIyNw1VVhnUoQFZcmC021AB/qdB9lLSrBSBox9CrKo/G59j2W\nL5/8nWJoz8laWiyqs0KsnSBiQ4GN/JvqlpuA3Vn3EcuiBCPTiaFXUR5VdbXP1MbRmjMnfKf66rKs\n8jwZt3svT97VlzFVi8asyARzTcPrq7LuI5ZFCUamEsvJpv7kW9+dt34Ik6lO6o1tMHPm+O9LMvWz\nUWZV9u9TVPKP4aIidkUmmBc0vJ6XdR+xLEowMpVmV9XDw3vP3NjNKpR2TurNElStJNNYRVa1UQGK\nqL5cu3by7+OuarJGhSWYXlqUYCSrsq/a64/Z7kl9qs+38/3Knpsl7+SmEZOnl3uCIQxC+Y/A/yDM\n37Iw7WdjXZRgpB1lX7W7t39ST5NAqnDne9FtMLWksnz55FkrJSgiwXwQ+DvgvcAngaeAncClwL5p\n9xPTogQj7erkqr3TxvFOTuppj92NuWU60Y1eZPW97lQ9NlkRCWZHw+tFwN8kiefjafcT06IEI+3I\nq4qqnRNzN07qab9fL3fpjaGUGrMiEszXgeMa1n07eWx7GuUyFyUYySqvE3y7J7CiT+oxtDGVTb/B\n9NImmCx38r8buNbMPmNm/93MPgH8LnlvZob9iFRWXuNVtXtn/qpVnY3NNR2Nx6XfIE8WklHKjc32\nAc4mjEP2C+A64D+A97n7hwuJsECDg4O+ffv2ssOQPlQbpmXlyjAsSizDtYikYWZ3ufvgdNvNyLDD\necD7gEOA7wEb3f1XyduVSy4iZWkcA2xoKK4xwUTykqWK7AbgSeCfgNnAN83slYVEJdLDVAUj/SJ1\nFZmZ7XT3Y+tevxD4nLufWFRwRVMVmYhIdmmryLKUYH5pZsfVXrj7DwklGRERkb2kboMBhoGbzeyf\nCTdYvhR4KI8gkvadG4GFwI+AZXXtO7Vtnk8YSWAfYF/CvTefTN77GrAA+Pdk8z9y98fziE1ERNoz\nbQnGzDaa2fuB5wGvB8aBg4G7gbfkFMdqYKu7HwVsTV432gP8gbsvAl4FrDazw+veP8/dFyWLkouI\nZFLmFNm9Kk0V2YbkcQVwB3A5sJhQ2vivOcVxRt1xNgBnNm7g7v/p7k8nL2eRrXpPRGRKixdPzPoJ\nE739Fi8uN64qm7aKzN23EkoVAJjZDOAY4HhCSeKmHOI41N33JMfbY2aHNNvIzI4Evgi8CLjI3R+t\ne/taM/stcDPwYc9yg4+I9L36qaV1f1I+srTBAODuzwD3Jst1aT9nZluAw5q8dUmGY/8UOC6pGvuC\nmX3e3R8jVI89YmbPJSSY5cDGFnEME9qTGBgYSHtoEekD9SMsjIwouXSqa9VM7n6Ku7+syXIr8JiZ\nLQBIHqdsQ0lKLvcDS5LXjySPTwKfA1ren+PuV7v7oLsPHnzwwfl8ORHpCePjoeQyMhIeG9tkJJtY\n2jFuI7TxkDze2riBmR1hZvslzw8EXgPsMrMZZjY/Wb8vcDpwX1eiFpGeUT/CwqWXTlSXKcm0L5YE\nczmw1Mx2A0uT15jZoJldk2xzNPBtM7uHMLLzFe6+k9Dg/1UzuxfYATwCfLrbX0CkG9TTqTgaYSF/\nmQa77DW6k1+qpnEcs8bXvWDdutBzq/77jI+HE31eo0ZLZ4q4k19ESlbf02l0tPeSC6i7cC9RghGp\nmHbnkqmKfkii/UIJRqRi+qGnU68n0X6hBCNSIf3S06kfkmg/UIIRqZB+6OnUL0m0HyjBiHRRp92M\nV63au7poaKi3ele1m0TVhTs+SjAiXaQeUtNrN4nqt41P5rHIRKR9GlCxOPpt46MSjEiXqYdUcfTb\nxkUJRqTL1EOqOPpt46IEI9JF6iFVHP228VGCEemifuhmXBb9tvHRYJca7FJ6mAaOlCJosEuRksVw\nX4a67kqZlGBEChLDyV0DR0qZlGBEChLLyV1dd6UsSjAiBYrh5K6uu1IWJRiRApV9clfXXSmTEoxI\nQWI4uavrrpRJ3ZTVTVkKoi7C0qvSdlNWgsmQYHTCEBGp4H0wZjbPzDab2e7k8cAptp1rZo+Y2Sfq\n1p1gZjvN7EEz+5iZWd4xxtDtVESkKqJJMMBqYKu7HwVsTV63sgb4esO6q4Bh4KhkeUPeAcbS7VRE\npApiSjBnABuS5xuAM5ttZGYnAIcCd9StWwDMdfc7PdT5bWz1+U7F0O1URKQKYkowh7r7HoDk8ZDG\nDczsWcB64KKGt54HPFz3+uFkXe7K7nYqIlIVXZ3R0sy2AIc1eeuSlLu4APiSu/+0oYmlWXtL094L\nZjZMqEpjYGAg5WGD+m6nQ0NhUTWZiEhzXU0w7n5Kq/fM7DEzW+Due5Iqr8ebbPZqYImZXQA8B5hp\nZk8BfwMcUbfdEcCjLWK4GrgaQi+yLPFPdU+BEoyIyGRdTTDTuA1YAVyePN7auIG7n1d7bmZvBwbd\nfXXy+kkzOxH4NvCnwMfzDrBZV+RaSUZERCaLqQ3mcmCpme0GliavMbNBM7smxedXAtcADwIPAV8u\nKlAREZmebrTUnfwiIplU7kZLqZ4YJtQSkXgpwUjbNLKBiEwlpkZ+qZj6kQ1Wrgz3BanLtojUqAQj\nHdHIBiLSihKMdEQjG4hIK0ow0rYYJtQSkXgpwUjbNFuiiExF98HoPhgRkUx0H4yIiJRKCUZERAqh\nBCMiIoVQghERkUIowYiISCGUYEREpBBKMCIiUgglGBERKYQSjIiIFEIJRkRECqEEIyIihVCCERGR\nQijBiIhIIaJIMGY2z8w2m9nu5PHAKbada2aPmNkn6tZ9zcx2mdmOZDmkO5GLiEgrUSQYYDWw1d2P\nArYmr1tZA3y9yfrz3H1RsjxeRJAiIpJeLAnmDGBD8nwDcGazjczsBOBQ4I4uxSUiIm2KJcEc6u57\nAJLHvaq4zOxZwHrgohb7uDapHhsxMysuVBERSWNGtw5kZluAw5q8dUnKXVwAfMndf9okf5zn7o+Y\n2XOBm4HlwMYWcQwDwwADAwMpDy0iIll1LcG4+ymt3jOzx8xsgbvvMbMFQLM2lFcDS8zsAuA5wEwz\ne8rdV7v7I8kxnjSzzwGvpEWCcfergashTJnc2bfqvnXrYPFiGBqaWDc+Dtu2wapV5cUlItIoliqy\n24AVyfMVwK2NG7j7ee4+4O4LgQ8CG919tZnNMLP5AGa2L3A6cF93wu6+xYth2bKQVCA8LlsW1ouI\nxCSWBHM5sNTMdgNLk9eY2aCZXTPNZ2cBXzWze4EdwCPAp4sMtkxDQ7BpU0gqo6PhcdOmySUaEZEY\nmHvlaolyMzg46Nu3by87jLaMjsKaNTAyApdeWnY0ItJPzOwudx+cbrtYSjCSwfg4XHVVSC5XXTVR\nXSYiEhMlmIqptbls2hRKLrXqMiUZEYmNEkzFbNs2uc2l1iazbVu5cYmINFIbTEXbYEREyqI2mB63\nbt3e1WLj42G9iEgMlGAqSvfDiEjsunYnv+Sr/n6YlStDbzLdDyMiMVEJpsKGhkJyWbMmPCq5iEhM\nlGAqTPfDiEjMlGAqSvfDiEjslGAqSvfDiEjsdB+M7oMREclE98GIiEiplGBERKQQSjAiIlIIJRgR\nESmEEoyIiBSir3uRmdnPgB9n+Mh84OcFhVOEKsVbpVihWvFWKVZQvEXKK9bnu/vB023U1wkmKzPb\nnqZrXiyqFG+VYoVqxVulWEHxFqnbsaqKTERECqEEIyIihVCCyebqsgPIqErxVilWqFa8VYoVFG+R\nuhqr2mBERKQQKsGIiEghlGCmYGZvNrP7zex3Ztay54WZ/cjMdprZDjMrbfTMtPEm2+5jZneb2e3d\niq/h+NPGambPNrPvmNk9ybb/s9tx1sWSJt4jzWzczB5Itr2w23EmcaT9u/17M3vczO7rZnxN4kgb\n7xvMbJeZPWhmq7sZY0Mc88xss5ntTh4PbLHdWjO7L1nO6XacSQxpY12X/Bs8YGYfMzPL4/hKMFO7\nDzgb+EaKbYfcfVHJ3RWzxHsh8ECx4UwpTaxPA6939+OBRcAbzOzEbgTXRJp4nwE+4O5HAycC7zGz\nY7oRXIO0fwefBd5QeDTTmzZeM9sH+FvgVOAY4C0l/bYAq4Gt7n4UsDV5PYmZvRF4BeHv9lXARWY2\nt6tRBmli/QPgNcBxwMuAxcBJeRxcCWYK7v6Au+8qO4600sZrZkcAbwSuKT6q5tLE6sFTyct9k6WU\nRsOU8e5x9+8mz58kJPDndSO+hjhS/R24+zeAX3YhpOniSBPvK4EH3f2H7v6fwA3AGcVH19QZwIbk\n+QbgzCZBcr1MAAAFbklEQVTbHAN83d2fcff/B9xDOck8TawOPBuYCcwi/D97LI+DK8Hkw4E7zOwu\nMxsuO5gU/hpYBfyu7ECmk1Tl7QAeBza7+7fLjikNM1sIvByoRLwV8Dzgp3WvH6aE5J041N33QLio\nAA5pss09wKlmNtvM5gNDwJFdjLFm2ljd/U5gHNiTLF9191xqN2bksZMqM7MtwGFN3rrE3W9NuZvX\nuPujZnYIsNnMvp9cHeau03jN7HTgcXe/y8xel3d8Dcfq+Ld1998Ci8zsAOAWM3uZuxfSZpDT3wJm\n9hzgZuDP3f3/5hVfwzFyibVbcoi3WZtAYaXZqeJN83l3v8PMFgP/CvwMuJNQhZq7TmM1sxcBRwNH\nJKs2m9lr8ziH9X2CcfdTctjHo8nj42Z2C6E4X0iCySHe1wB/YmanEYrFc83sf7v72zqPbrI8ftu6\nfT1hZl8jVDMUkmDyiNfM9iUkl+vd/R87j6q5PH/bbsgh3oeZXAI4Ani0w322NFW8ZvaYmS1w9z1m\ntoBQum62j48AH0k+8zlgd6SxngV8q1YdbWZfJrQhdnwOUxVZh8xsjpk9t/Yc+CMKOgHmwd0vdvcj\n3H0hcC4wVkRyyYOZHZyUXDCz/YBTgO+XG1VrSc+bzwAPuPuVZcfTY7YBR5nZC8xsJuFv97aSYrkN\nWJE8XwHsVQJLqnYPSp4fR2hAv6NrEU6YNlbgJ8BJZjYjuUA6ibw6ALm7lhYLIbM/TOjN9BihbhLg\ncOBLyfMXEupb7wHuJxT5o423YfvXAbfHGivhP+XdwL2EpD0a828L/CGh2uZeYEeynBZjrMnrfyDU\nuf8m2f6dsf62yevTgB8AD5X8/+wgQo+s3cnjvGT9IHBN8vzZwPeS5VvAoohj3Qf4FCGpfA+4Mq/j\n605+EREphKrIRESkEEowIiJSCCUYEREphBKMiIgUQglGREQKoQQjIiKFUIIREZFCKMGIZGBmp5jZ\ndV061vlm5mZ2Ut269ybrTqlb9ykze03a7UW6RQlGJJvjCaM2dMNxhFEBjgYws9nAOwmDJ+6s2+5V\nhLvF024v0hVKMCLZHA/sMLOXmNk3klkAtyRDsmNmRyfr7zWzi8zswQ6OdSxhOJeXJK//DLgJ+J27\nP1Y7HvADD6NOT7u9SDcpwYhkczyhNHAzcKG7vxTYDLzPzGYA1yfrjyOMU9fJwKdHA5uAl5jZ/sA5\nhOHf6/d5KvCVDNuLdI0SjEhKyUizcwmDhH7T3e9O3voeYSKns4F7Gtbfk3z2hWb2GTP7fPJ6jplt\nMLNPm9l5TY51JPALd/9hsu9VwMeBFxOqwWr+GPhK2u3N7L+Y2d+Z2XozO7zzX0WkNSUYkfSOIYw4\newyT2zSOJSST4wgjKNe8rPbaw1S/76x772zg8+7+34A/aXKs4+qO8SRhHpwNybF2wu/bWA7wMB9R\nmu3nAxcD3wVuB65Ihr4XKYQSjEh6xxMSxiOEJIOZvRBYDmwEfkEoMWBmi4C30bpDwBFMTAH82ybv\n/z4xAB8F3lvXzlIrkQwRprpNu/1rCf/n/wz4S+Ap4KXTfmuRNvX9jJYiGRwPfIcwidNpZrYT+Hfg\nHe7+i6T78hfNbBthitwfJVVWzTxMSDI7aH6hdyyhnQd3v71u/TGE0hKE9pfPZ9h+NyHh/SvwTeAd\nTJ7nXiRXmg9GJCdm9hyfmHb2ImB/d/+L5PVBhOlzlwLXAB8DPgH8B6E95/o2jvdd4FXu/psMn7mA\nUCX3G+Cz7n5z1uOKpKUEI5ITMxshTOX7G+BfgPe7+9PlRiVSHiUYEREphBr5RUSkEEowIiJSCCUY\nEREphBKMiIgUQglGREQKoQQjIiKFUIIREZFCKMGIiEgh/j9nBQFBAjwPGwAAAABJRU5ErkJggg==\n", "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -145,7 +147,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -154,7 +156,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Terran 69.0 %, Neptunian 31.0 %, Jovian 0.0 %, Star 0.0 %\n" + "Terran 62.0 %, Neptunian 38.0 %, Jovian 0.0 %, Star 0.0 %\n" ] } ], @@ -164,7 +166,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -178,7 +180,7 @@ } ], "source": [ - "print 'M = %.3f (+ %.3f - %.3f) MEarth' % (Mmedian, Mplus, Mminus)" + "print('M = %.3f (+ %.3f - %.3f) MEarth' % (Mmedian, Mplus, Mminus))" ] }, { @@ -192,7 +194,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -201,7 +203,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Terran 72.0 %, Neptunian 28.0 %, Jovian 0.0 %, Star 0.0 %\n" + "Terran 73.0 %, Neptunian 27.0 %, Jovian 0.0 %, Star 0.0 %\n" ] } ], @@ -212,16 +214,16 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEVCAYAAAASFwXVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEdBJREFUeJzt3X+sZGV9x/H3BxatYOuK4O6qKG1SfyUmoBGt0jr+LDYE\nMbUYm6bbai2JRa1JW1Abva1/+KNqTGP0H9FurLUhWikoVBbDRA0BhLIICKIWKlZ2wYoKpUWUb/+Y\ns8vlevfe+XFn5s7D+5VM7pkz55znOfvc/cy5z/nxpKqQJC22Q+ZdAUnS5AxzSWqAYS5JDTDMJakB\nhrkkNcAwl6QGrBnmSX4pyeVJ9iT5RpJ3d/OPTLI7yU1JLkqydTbVlSStJutdZ57k8Kq6J8kW4KvA\nXwCnAD+oqvclORN4dFWdNf3qSpJWs243S1Xd000+DDgUuJNBmO/q5u8CTp1K7SRJQ1k3zJMckmQP\nsA+4pKquB7ZV1b5ukX3AtinWUZK0ji3rLVBV9wPHJXkU8MUkL1zxeSXxmQCSNEfrhvl+VfXjJF8A\nngXsS7K9qvYm2QHcvnJ5A16SxlNVGXWd9a5mOWr/lSpJHgG8FLgaOA/Y2S22Ezj3IBVq9vXOd75z\n7nVw/9y/h+L+tbxvVeMfA693ZL4D2JXkEAbB/8mq+lKSq4FzkrwOuAU4bewaSJImtmaYV9W1wDNX\nmf9D4CXTqpQkaTTeATqmXq837ypMlfu32Frev5b3bRLr3jQ09oaTmta2JalVSaiNPgEqSVoMhrkk\nNcAwl6QGGOaS1ADDXJIaYJhLUgMMc0lqgGEuSQ0wzCWpAYa5JDXAMJekBhjmktQAw1ySGmCYS1ID\nDHNJaoBhLkkNWG8MUEkPMcnI4yJMxEFsNoZhLmkVswrY2X5xtMxuFklqgGEuSQ0wzCWpAYa5JDXA\nMJekBhjmktQAw1ySGmCYS1ID1gzzJMckuSTJ9UmuS/Kmbv5Sku8lubp7nTSb6kqSVpO1bqVNsh3Y\nXlV7kjwSuAo4FTgNuKuqPrjGuuVtutLiGdzOP7s7QM2JB0tCVY18a+yat/NX1V5gbzd9d5IbgMfv\nL3PkWkqSpmLoPvMkxwLHA5d1s96Y5JokZyfZOoW6SZKGNFSYd10snwHeXFV3Ax8FfhU4DrgN+MDU\naihJWte6T01MchjwWeAfq+pcgKq6fdnnHwPOX23dpaWlA9O9Xo9erzdZbSWpMf1+n36/P/F21jsB\nGmAX8N9V9ZZl83dU1W3d9FuAZ1fV769Y1xOg0gLyBOh8jXsCdL0wPxH4MvB1HmjdtwGvYdDFUsDN\nwOlVtW/Fuoa5tIAM8/maSphPwjCXFpNhPl/jhrl3gEpSAwxzSWqAYS5JDTDMJakBhrkkNcAwl6QG\nGOaS1ADDXJIaYJhLUgMMc0lqgGEuSQ0wzCWpAYa5JDXAMJekBhjmktQAw1ySGmCYS1IDDHNJaoBh\nLkkNMMwlqQGGuSQ1wDCXpAYY5pLUAMNckhpgmEtSAwxzSWqAYS5JDTDMJakBhrkkNWDNME9yTJJL\nklyf5Lokb+rmH5lkd5KbklyUZOtsqitJWk2q6uAfJtuB7VW1J8kjgauAU4E/Bn5QVe9Lcibw6Ko6\na8W6tda2JW1OSYBZ/d8N5sSDJaGqMup6ax6ZV9XeqtrTTd8N3AA8HjgF2NUttotBwEuS5mToPvMk\nxwLHA5cD26pqX/fRPmDbhtdMkjS0LcMs1HWxfBZ4c1XdNfgzbKCqKsmqfyctLS0dmO71evR6vUnq\nKknN6ff79Pv9ibezZp85QJLDgM8DF1bVh7p5NwK9qtqbZAdwSVU9dcV69plLC8g+8/maSp95Bq16\nNvCN/UHeOQ/Y2U3vBM4dtWBJ0sZZ72qWE4EvA1/nga/qtwJXAOcATwRuAU6rqh+tWNcjc2kBeWQ+\nX+Mema/bzTIuw1xaTIb5fE2lm0WStBgMc0lqgGEuSQ0wzCWpAYa5JDXAMJekBhjmktQAw1ySGmCY\nS1IDDHNJaoBhLkkNMMwlqQGGuSQ1wDCXpAYY5pLUAMNckhpgmEtSA7bMuwKShjMYAUhanWEuLZRZ\nDLHml8YisptFkhpgmEtSAwxzSWqAYS5JDTDMJakBhrkkNcBLE6UJef23NgPDXNoQXv+t+bKbRZIa\nsG6YJ/l4kn1Jrl02bynJ95Jc3b1Omm41JUlrGebI/BPAyrAu4INVdXz3+reNr5okaVjrhnlVfQW4\nc5WP7MCTpE1ikj7zNya5JsnZSbZuWI0kSSMb92qWjwJ/202/C/gA8LqVCy0tLR2Y7vV69Hq9MYuT\npDb1+336/f7E20nV+pdUJTkWOL+qnjHsZ0lqmG1Li25wnfmsLk1sqZxBWebEgyWhqkbuxh6rmyXJ\njmVvXwlce7BlJUnTt243S5JPAy8AjkpyK/BOoJfkOAZf3zcDp0+1lpKkNQ3VzTLWhu1m0UOE3SyT\nlWVOPNhMu1kkSZuLYS5JDTDMJakBhrkkNcAwl6QGGOaS1ADDXJIaYJhLUgMMc0lqgGEuSQ1wQGdJ\nczV4HML0tf7YAMNc0pzN6nkzbbObRZIaYJhLUgMMc0lqgGEuSQ0wzCWpAYa5JDXAMJekBhjmktQA\nw1ySGmCYS1IDDHNJaoBhLkkNMMwlqQE+NVHNmtWjVaXNwDBX43y8qh4a7GaRpAasG+ZJPp5kX5Jr\nl807MsnuJDcluSjJ1ulWU5K0lmGOzD8BnLRi3lnA7qp6MvCl7r0kaU7WDfOq+gpw54rZpwC7uuld\nwKkbXC9J0gjG7TPfVlX7uul9wLYNqo8kaQwTX81SVZVk1UsGlpaWDkz3ej16vd6kxUlSU/r9Pv1+\nf+LtpGr9S7eSHAucX1XP6N7fCPSqam+SHcAlVfXUFevUMNuWpmVwnfmsLk20nM1dVliUPEpCVY18\nveu43SznATu76Z3AuWNuR5K0AdY9Mk/yaeAFwFEM+sffAfwrcA7wROAW4LSq+tGK9Twy11x5ZL7Z\ny5llWe0fmQ/VzTIOw1zzZphv9nJmWVb7Ye4doJLUAMNckhpgmEtSAwxzSWqAYS5JDTDMJakBhrkk\nNcAwl6QGGOaS1ADDXJIaYJhLUgMMc0lqgGEuSQ0wzCWpAYa5JDXAMJekBhjmktQAw1ySGrBl3hXQ\nQ89gODdJG8kw15zMaixL6aHBbhZJaoBhLkkNMMwlqQGGuSQ1wDCXpAYY5pLUAMNckhpgmEtSAya6\naSjJLcBPgJ8D91XVCRtRKUnSaCa9A7SAXlX9cCMqI0kaz0Z0s3jPtCTN2aRhXsDFSa5M8vqNqJAk\naXSTdrM8v6puS3I0sDvJjVX1lY2omCRpeBOFeVXd1v28I8nngBOAA2G+tLR0YNler0ev15ukOElq\nTr/fp9/vT7ydVI33KNIkhwOHVtVdSY4ALgL+pqou6j6vcbettg2eZz6rR+BazuYtZ5ZlhUXJoyRU\n1cjnIic5Mt8GfK4baGAL8Kn9QS5Jmq2xj8zX3bBH5joIj8wtZ/ZltX9k7h2gktQAw1ySGmCYS1ID\nHNBZwP5+bKlds/odn1ffvGGuZWZ50kuatVmdPJ4Pu1kkqQGGuSQ1wDCXpAYY5pLUAMNckhpgmEtS\nAwxzSWqAYS5JDTDMJakBhrkkNcAwl6QGGOaS1ADDXJIaYJhLUgMMc0lqgGEuSQ0wzCWpAYa5JDXA\nMJekBjgG6Bhuvvlm7rrrrpmUdfTRR7Njx46ZlCVpcWVaI0knqXmNUj1tL3rRyVx++XVs2fIrUy3n\n3ntv54wzdvL+9793quXA/pHLZzmg86wG17WczVvOLMuaXTmT5l4SqmrkkaE9Mh/DfffBPfd8GDh5\nyiW9l6ofTrkMSS2wz1ySGjB2mCc5KcmNSb6V5MyNrJQkaTRjhXmSQ4EPAycBTwdek+RpG1mxze/a\neVdgyvrzrsCU9eddgSnrz7sCU9SfdwU2pXGPzE8Avl1Vt1TVfcA/A6/YuGotAsN8sfXnXYEp68+7\nAlPUn3cFNqVxw/zxwK3L3n+vmydJmoNxr2Zp85rDIR16KGzZ8iWOOOLFUy3n3ntvIXnVVMuQ1Iax\nrjNP8lxgqapO6t6/Fbi/qt67bJmHdOBL0rjGuc583DDfAnwTeDHwfeAK4DVVdcPIG5MkTWysbpaq\n+lmSM4AvAocCZxvkkjQ/U7udX5I0Oxt2B2iSv0tyQ5JrkvxLkkcdZLmFvNkoye8luT7Jz5M8c43l\nbkny9SRXJ7lilnWcxAj7t6jtd2SS3UluSnJRkq0HWW5h2m+Ytkjy993n1yQ5ftZ1nMR6+5ekl+TH\nXVtdneSv51HPcST5eJJ9SQ56jfPIbVdVG/ICXgoc0k2/B3jPKsscCnwbOBY4DNgDPG2j6jDNF/BU\n4MnAJcAz11juZuDIedd3Gvu34O33PuCvuukzV/v9XKT2G6YtgN8BLuimnwNcNu96b/D+9YDz5l3X\nMffvN4HjgWsP8vnIbbdhR+ZVtbuq7u/eXg48YZXFFvZmo6q6sapuGnLxkc9Ez9uQ+7ew7QecAuzq\npncBp66x7CK03zBtcWCfq+pyYGuSbbOt5tiG/V1bhLb6BVX1FeDONRYZue2m9aCt1wIXrDL/oXCz\nUQEXJ7kyyevnXZkNtsjtt62q9nXT+4CD/cdYlPYbpi1WW2a1g6zNaJj9K+B5XTfEBUmePrPaTd/I\nbTfS1SxJdgPbV/nobVV1frfM24GfVtU/rbLcpj7bOsz+DeH5VXVbkqOB3Ulu7L6F524D9m9R2+/t\ny99UVa1xH8Smbb8Vhm2LlUeum7oNlxmmnv8OHFNV9yR5OXAug67CVozUdiOFeVW9dM2Skz9i0Ndz\nsFsj/ws4Ztn7Yxh842wK6+3fkNu4rft5R5LPMfhzcVOEwQbs38K2X3eyaXtV7U2yA7j9INvYtO23\nwjBtsXKZJ3TzFsG6+1dVdy2bvjDJR5IcWW0MAjBy223k1SwnAX8JvKKq/u8gi10J/HqSY5M8DHg1\ncN5G1WGGVu2nS3J4kl/upo8AXsZiPpHrYP2Qi9x+5wE7u+mdDI7iHmTB2m+YtjgP+EM4cNf2j5Z1\nNW126+5fkm0ZDJFFkhMYXGrdQpDDOG23gWdnvwX8J3B19/pIN/9xwBeWLfdyBnePfht467zPKo+w\nf69k0If1v8Be4MKV+wf8GoOz7nuA61rbvwVvvyOBi4GbgIuArYvefqu1BXA6cPqyZT7cfX4Na1yF\ntRlf6+0f8GddO+0BLgWeO+86j7Bvn2Zw9/xPu/93r5207bxpSJIa4LBxktQAw1ySGmCYS1IDDHNJ\naoBhLkkNMMwlqQGGuSQ1wDCXVkiyJclT5l0PaRSGuTaVJO9P8q4pbfsNSX6S5DEr5p+T5B+SPK2b\n1QPuH2M9aW4Mc2023wEum9K2rwAuZDDgAQDdqEqPBN5VD4xj+5Sq+tYY60lzY5hrszmBweAm0/Ak\nBk9AfOKyeY8EHltV31k2734ebNj1pLkZ6RG40gw8tqp+kORk4DHA0QwehHVDkscxeCDRrcDzqur0\nEbedbt1jAZL8BvAfLHscbvf0va+Nup40bx6Za9PoBgG/M8mTgT+oql0MRqx6Q7fIR4EPAbuB/xmz\nmFuBY5IcxuBh/8cx6EbZ71lVdeUY60lzZZhrM3k2g4DcCXyqm/ckBgH/JAbPq76bwQC3lwIkeUV3\nxE6SE5K8JMkvHLEn+RUGYy7eyqC75LlVdRm/2K1zyJjrSXNlmGszeRaDQQleDny3m/cq4JPAoxk8\n2xrgBcClSbbzwIATAKdV1cXAw5Ms79+GwRfFVVV1B4Pnlu8fpeZAt0p3OeI3R11vvyRPS3Jckmd3\nR/DSzBjm2ky+A5wInAG8LMlO4DPdlSXXAj9P8rvAc6rq+1W1l8GD+/c7vPt5N8sGbE5yIvBu4ORu\n1lerak+SNzDoLjmxm98D+mOsR5LjgDuqak9VfQ34rYn+JaQReQJUm0ZVfWbZ20tXfHxUVZ3V9av/\n9orP9g9z9+Pu51bgwBBbVfVVBkfS+9//effzI8BHlm3nYVX1szHWAziiO3H7MeAdDEZskmbGI3Mt\nivckORX4U2AJIMljgacAL+yWuSDJC4H7q+q7q27lILp+90kGO747ydFV9SdV9X3g4RNsSxqZw8ZJ\nQJJXA5+vqnGvkiHJ04FHMPiL96rlR/nStBnmktQAu1kkqQGGuSQ1wDCXpAYY5pLUAMNckhpgmEtS\nAwxzSWqAYS5JDTDMJakB/w/zNH7nMOt6igAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEOCAYAAACZ2uz0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADg5JREFUeJzt3X2MHHd9x/HPp5iCGiB14rNrIO41akpimthIp1DVEk3I\nQ/MgJSEtolGLXOHWVG1U2iLQ9UlFIFQjSFGfVedBuG2aQpNGcesoYAytRQuUM3ViBxOcRoYaW7aT\nVBAkmhLn2z/md3hx9ryzu7Ozt997v6TT7s7O7fxudvX2eHZn1hEhAMDk+75xDwAA0AyCDgBJEHQA\nSIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgiWVtLmzFihUxPT3d5iIBYOLt2bPnyYiY6jVfq0Gf\nnp7W3Nxcm4sEgIln+6t15mOXCwAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACTR\n6pGiABaP6dkdY1v2oS3Xj23ZmbGFDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIO\nAEkQdABIgqADQBIEHQCSIOgAkARBB4Akegbd9nm2P237gO1Hbb+jTD/H9k7bB8vl8tEPFwCwkDpb\n6M9JemdEXCTpJyT9mu21kmYl7YqICyTtKrcBAGPSM+gRcTQivliuPyPpgKRXSbpR0rYy2zZJN41q\nkACA3vrah257WtLrJH1e0qqIOCpV0Ze0sunBAQDqqx102y+TdJ+k34iIb/bxe5ttz9meO3HixCBj\nBADUUCvotl+sKuZ3R8Q/lsnHbK8u96+WdLzb70bE1oiYiYiZqampJsYMAOiizqdcLOlOSQci4o86\n7touaWO5vlHSA80PDwBQ17Ia82yQ9FZJ+2zvLdN+R9IWSR+zvUnS1yS9eTRDBADU0TPoEfEZSV7g\n7iuaHQ4AYFAcKQoASRB0AEiCoANAEgQdAJKo8ykXACM0Pbtj3ENAEmyhA0ASBB0AkiDoAJAEQQeA\nJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANA\nEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJDEsnEPAMDSMz27YyzL\nPbTl+rEsty1soQNAEgQdAJIg6ACQBEEHgCQIOgAk0TPotu+yfdz2/o5p77H9ddt7y891ox0mAKCX\nOlvoH5F0TZfpH46I9eXnwWaHBQDoV8+gR8RuSU+3MBYAwBCG2Yd+q+1Hyi6Z5Y2NCAAwkEGPFP1L\nSe+TFOXyNklv6zaj7c2SNkvSmjVrBlwcMFrjOnIRaNJAW+gRcSwiTkbE85Jul3TpGebdGhEzETEz\nNTU16DgBAD0MFHTbqztuvknS/oXmBQC0o+cuF9v3SLpM0grbhyX9gaTLbK9XtcvlkKS3j3CMAIAa\negY9Im7pMvnOEYwFADAEjhQFgCQIOgAkQdABIAmCDgBJ8BV0WFQ4wAcYHFvoAJAEQQeAJAg6ACRB\n0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg\n6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQ\ndABIgqADQBIEHQCSIOgAkETPoNu+y/Zx2/s7pp1je6ftg+Vy+WiHCQDopc4W+kckXXPatFlJuyLi\nAkm7ym0AwBj1DHpE7Jb09GmTb5S0rVzfJummhscFAOjToPvQV0XEUUkqlyubGxIAYBAjf1PU9mbb\nc7bnTpw4MerFAcCSNWjQj9leLUnl8vhCM0bE1oiYiYiZqampARcHAOhl0KBvl7SxXN8o6YFmhgMA\nGFSdjy3eI+mzkl5j+7DtTZK2SLrK9kFJV5XbAIAxWtZrhoi4ZYG7rmh4LACAIXCkKAAkQdABIAmC\nDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARB\nB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQxLJxDwCLz/TsjnEPAcAA2EIHgCQIOgAk\nQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBJDnT7X\n9iFJz0g6Kem5iJhpYlAAgP41cT70yyPiyQYeBwAwBHa5AEASwwY9JH3C9h7bm5sYEABgMMPuctkQ\nEUdsr5S00/aXI2J35wwl9Jslac2aNUMuDgCwkKG20CPiSLk8Lul+SZd2mWdrRMxExMzU1NQwiwMA\nnMHAQbd9lu2Xz1+XdLWk/U0NDADQn2F2uaySdL/t+cf5u4h4qJFRAQD6NnDQI+IJSesaHAsAYAh8\nbBEAkiDoAJAEQQeAJAg6ACTRxLlcMCLTszvGPQQAE4QtdABIgqADQBIEHQCSIOgAkARBB4AkCDoA\nJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASfAUd\ngCVjnF/reGjL9SNfBlvoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSmJgDi7IfEAAA\nw2ILHQCSIOgAkARBB4AkCDoAJEHQASCJoYJu+xrbj9l+3PZsU4MCAPRv4KDbfpGkP5d0raS1km6x\nvbapgQEA+jPMFvqlkh6PiCci4v8k/b2kG5sZFgCgX8ME/VWS/rvj9uEyDQAwBsMcKeou0+IFM9mb\nJW0uN79l+7GOu1dIenKIMbTCHxj5IiZiPbSEdVFhPZySYl0M2ZEfrjPTMEE/LOm8jtuvlnTk9Jki\nYqukrd0ewPZcRMwMMYYUWA+nsC4qrIdTWBf1DbPL5QuSLrD9I7a/X9LPSdrezLAAAP0aeAs9Ip6z\nfaukj0t6kaS7IuLRxkYGAOjLUGdbjIgHJT04xEN03RWzBLEeTmFdVFgPp7AuanLEC97HBABMIA79\nB4AkWg267TfbftT287YXfNfa9iHb+2zvtT3X5hjb0Md6SH9qBdvn2N5p+2C5XL7AfCfL62Gv7TRv\nvvd6jm2/xPZHy/2ftz3d/ihHr8Z6+EXbJzpeA780jnEudm1voe+XdLOk3TXmvTwi1if9uFLP9bCE\nTq0wK2lXRFwgaVe53c23y+thfUTc0N7wRqfmc7xJ0v9ExI9K+rCk0R8V0bI+Xusf7XgN3NHqICdE\nq0GPiAMR8VjvOXOruR6WyqkVbpS0rVzfJummMY6lbXWe4871c6+kK2x3O6hvki2V1/rILdZ96CHp\nE7b3lCNNl6KlcmqFVRFxVJLK5coF5nup7Tnbn7OdJfp1nuPvzhMRz0n6hqRzWxlde+q+1n/G9iO2\n77V9Xpf7l7zGvyTa9icl/VCXu343Ih6o+TAbIuKI7ZWSdtr+ckTU2U2zaDSwHmqdWmESnGld9PEw\na8pr4nxJn7K9LyL+q5kRjk2d5zjN6+AM6vyN/yTpnoh41vavqPpfyxtHPrIJ03jQI+LKBh7jSLk8\nbvt+Vf8lm6igN7Aeap1aYRKcaV3YPmZ7dUQctb1a0vEFHmP+NfGE7X+R9DpJkx70Os/x/DyHbS+T\ndLakp9sZXmt6roeIeKrj5u1K+F5CExbdLhfbZ9l++fx1SVerehNxqVkqp1bYLmljub5R0gv+92J7\nue2XlOsrJG2Q9KXWRjg6dZ7jzvXzs5I+FfkOHum5Hso/9vNukHSgxfFNjoho7UfSm1T9a/yspGOS\nPl6mv1LSg+X6+ZIeLj+PqtpF0eo4F8N6KLevk/QVVVui6dZD+RvPVfXploPl8pwyfUbSHeX6T0ra\nV14T+yRtGve4G/z7X/AcS3qvpBvK9ZdK+gdJj0v6D0nnj3vMY1oPf1h68LCkT0u6cNxjXow/HCkK\nAEksul0uAIDBEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0LHo2b7S9t+0tKy32w7bP9Ux7dYy\n7cqOaX9le0Pd+YE2EHRMgnWqjhBswyWSHpF0kSTZ/gFV5yQ/oeoo1Xmvl/S5PuYHRo6gYxKsk7TX\n9oW2d5dve/pkOa+LbF9Upj9i+122Hx9iWRdLukfSheX2r6s69P75iDg2vzxJX4mIk3XmB9pC0DEJ\n1qna2r1P0jsi4rWSdkr6zXIGwrvL9EtUnQtomJO5XSTpY5IutH22pLdI+vfTHvNaSQ/1MT/QCoKO\nRc32iyW9QtJlkj4TEf9Z7vqSqi/DuFnSw6dNf7j87vm277R9b7l9lu1ttm+3/fNdlnWepKci4ony\n2O+W9KeSfkzVbpV5Py3pobrz236N7b+wfZvtVw6/VoDuCDoWu7WqTpW6Vt+7T/piVfG+RNLejuk/\nPn87qq8029Rx382S7o2IX1Z1CtbTXdKxjGckXaPqixQunp9e9pH/YFTnZ68z/wpJvy3pPapi/6Fy\niligcQQdi906VYH+uqqoq3xr0Vsl/bWkp1RtEcv2ekm/oIXfQH21Tn3V2cku9383xJI+KOnWjv3k\n81vcl6s6fWvd+d8g6Y9VnTJ5lardRq/t/WcD/Wv8G4uAhq1TdR7w7ZKus71P0rclvS0iniofZ9xh\n+wuSPivpUNkF0s1hVVHfq+4bMxerCq4i4p87pq/VqS/UuFbVlzXXnf+gpKsj4jZJsv0WSf9a4+8G\n+sb50DHRbL8sIr5Vrr9L0tkR8Xvl9rmS3i/pKkl3SPoTSX8m6X9V7Y+/e4DlfVHS6yPiO338zq9K\nulLVd2f+bUTc1+9ygToIOiaa7d9X9ZVl35H0b5J+KyKeHe+ogPEg6ACQBG+KAkASBB0AkiDoAJAE\nQQeAJAg6ACRB0AEgCYIOAEkQdABI4v8BSVciV4cqFd8AAAAASUVORK5CYII=\n", "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -234,6 +236,15 @@ "plt.show()" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, @@ -246,21 +257,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.10" + "pygments_lexer": "ipython3", + "version": "3.5.2" } }, "nbformat": 4, diff --git a/func.py b/func.py index e3732e5..793e20d 100644 --- a/func.py +++ b/func.py @@ -1,79 +1,81 @@ +from __future__ import (print_function, absolute_import, + division, unicode_literals) import numpy as np -from scipy.stats import norm, truncnorm +from scipy.stats import norm ### fix the number of different populations n_pop = 4 def indicate(M, trans, i): - ''' - indicate which M belongs to population i given transition parameter - ''' - ts = np.insert(np.insert(trans, n_pop-1, np.inf), 0, -np.inf) - ind = (M>=ts[i]) & (M=ts[i]) & (M 3e5: - print 'Mass range out of model expectation. Returning None.' - return None - - ## convert to radius - sample_size = len(mass) - logm = np.log10(mass) - prob = np.random.random(sample_size) - logr = np.ones_like(logm) - - hyper_ind = np.random.randint(low = 0, high = np.shape(all_hyper)[0], size = sample_size) - hyper = all_hyper[hyper_ind,:] - - if classify == 'Yes': - classification(logm, hyper[:,-3:]) - - - for i in range(sample_size): - logr[i] = piece_linear(hyper[i], logm[i], prob[i]) - - radius_sample = 10.** logr - - ## convert to right unit - if unit == 'Jupiter': - radius = radius_sample / rearth2rjup - else: - radius = radius_sample - - return radius - - - -def Mstat2R(mean, std, unit='Earth', sample_size=1000, classify = 'No'): - """ - Forecast the mean and standard deviation of radius given the mena and standard deviation of the mass. - Assuming normal distribution with the mean and standard deviation truncated at the mass range limit of the model. - - Parameters - --------------- - mean: float - Mean (average) of mass. - std: float - Standard deviation of mass. - unit: string (optional) - Unit of the mass. Options are 'Earth' and 'Jupiter'. - sample_size: int (optional) - Number of mass samples to draw with the mean and std provided. - Returns - --------------- - mean: float - Predicted mean of radius in the input unit. - std: float - Predicted standard deviation of radius. - """ - - # unit - if unit == 'Earth': - pass - elif unit == 'Jupiter': - mean = mean * mearth2mjup - std = std * mearth2mjup - else: - print "Input unit must be 'Earth' or 'Jupiter'. Using 'Earth' as default." - - # draw samples - mass = truncnorm.rvs( (mlower-mean)/std, (mupper-mean)/std, loc=mean, scale=std, size=sample_size) - if classify == 'Yes': - radius = Mpost2R(mass, unit='Earth', classify='Yes') - else: - radius = Mpost2R(mass, unit='Earth') - - if unit == 'Jupiter': - radius = radius / rearth2rjup - - r_med = np.median(radius) - onesigma = 34.1 - r_up = np.percentile(radius, 50.+onesigma, interpolation='nearest') - r_down = np.percentile(radius, 50.-onesigma, interpolation='nearest') - - return r_med, r_up - r_med, r_med - r_down + """ + Forecast the Radius distribution given the mass distribution. + + Parameters + --------------- + mass: one dimensional array + The mass distribution. + unit: string (optional) + Unit of the mass. + Options are 'Earth' and 'Jupiter'. Default is 'Earth'. + classify: string (optional) + If you want the object to be classifed. + Options are 'Yes' and 'No'. Default is 'No'. + Result will be printed, not returned. + + Returns + --------------- + radius: one dimensional array + Predicted radius distribution in the input unit. + """ + + # mass input + mass = np.array(mass) + assert len(mass.shape) == 1, "Input mass must be 1-D." + + # unit input + if unit == 'Earth': + pass + elif unit == 'Jupiter': + mass = mass * mearth2mjup + else: + print("Input unit must be 'Earth' or 'Jupiter'. " + "Using 'Earth' as default.") + + # mass range + if np.min(mass) < 3e-4 or np.max(mass) > 3e5: + print('Mass range out of model expectation. Returning None.') + return None + + ## convert to radius + sample_size = len(mass) + logm = np.log10(mass) + prob = np.random.random(sample_size) + logr = np.ones_like(logm) + + hyper_ind = np.random.randint(low=0, high=np.shape(all_hyper)[0], + size=sample_size) + hyper = all_hyper[hyper_ind,:] + + if classify == 'Yes': + classification(logm, hyper[:,-3:]) + + for i in range(sample_size): + logr[i] = piece_linear(hyper[i], logm[i], prob[i]) + + radius_sample = 10.** logr + + ## convert to right unit + if unit == 'Jupiter': + radius = radius_sample / rearth2rjup + else: + radius = radius_sample + + return radius + + + +def Mstat2R(mean, std, unit='Earth', sample_size=1000, classify='No'): + """ + Forecast the mean and standard deviation of radius given the mean and + standard deviation of the mass. Assuming normal distribution with the mean + and standard deviation truncated at the mass range limit of the model. + + Parameters + ---------- + mean: float + Mean (average) of mass. + std: float + Standard deviation of mass. + unit: string (optional) + Unit of the mass. Options are 'Earth' and 'Jupiter'. + sample_size: int (optional) + Number of mass samples to draw with the mean and std provided. + Returns + ------- + mean: float + Predicted mean of radius in the input unit. + std: float + Predicted standard deviation of radius. + """ + + # unit + if unit == 'Earth': + pass + elif unit == 'Jupiter': + mean = mean * mearth2mjup + std = std * mearth2mjup + else: + print("Input unit must be 'Earth' or 'Jupiter'. " + "Using 'Earth' as default.") + + # draw samples + mass = truncnorm.rvs((mlower-mean)/std, (mupper-mean)/std, loc=mean, + scale=std, size=sample_size) + if classify == 'Yes': + radius = Mpost2R(mass, unit='Earth', classify='Yes') + else: + radius = Mpost2R(mass, unit='Earth') + + if unit == 'Jupiter': + radius = radius / rearth2rjup + + r_med = np.median(radius) + onesigma = 34.1 + r_up = np.percentile(radius, 50.+onesigma, interpolation='nearest') + r_down = np.percentile(radius, 50.-onesigma, interpolation='nearest') + + return r_med, r_up - r_med, r_med - r_down def Rpost2M(radius, unit='Earth', grid_size = 1e3, classify = 'No'): - """ - Forecast the mass distribution given the radius distribution. - - Parameters - --------------- - radius: one dimensional array - The radius distribution. - unit: string (optional) - Unit of the mass. Options are 'Earth' and 'Jupiter'. - grid_size: int (optional) - Number of grid in the mass axis when sampling mass from radius. - The more the better results, but slower process. - classify: string (optional) - If you want the object to be classifed. - Options are 'Yes' and 'No'. Default is 'No'. - Result will be printed, not returned. - - Returns - --------------- - mass: one dimensional array - Predicted mass distribution in the input unit. - """ - - # unit - if unit == 'Earth': - pass - elif unit == 'Jupiter': - radius = radius * rearth2rjup - else: - print "Input unit must be 'Earth' or 'Jupiter'. Using 'Earth' as default." - - - # radius range - if np.min(radius) < 1e-1 or np.max(radius) > 1e2: - print 'Radius range out of model expectation. Returning None.' - return None - - - - # sample_grid - if grid_size < 10: - print 'The sample grid is too sparse. Using 10 sample grid instead.' - grid_size = 10 - - ## convert to mass - sample_size = len(radius) - logr = np.log10(radius) - logm = np.ones_like(logr) - - hyper_ind = np.random.randint(low = 0, high = np.shape(all_hyper)[0], size = sample_size) - hyper = all_hyper[hyper_ind,:] - - logm_grid = np.linspace(-3.522, 5.477, grid_size) - - for i in range(sample_size): - prob = ProbRGivenM(logr[i], logm_grid, hyper[i,:]) - logm[i] = np.random.choice(logm_grid, size=1, p = prob) - - mass_sample = 10.** logm - - if classify == 'Yes': - classification(logm, hyper[:,-3:]) - - ## convert to right unit - if unit == 'Jupiter': - mass = mass_sample / mearth2mjup - else: - mass = mass_sample - - return mass - - - -def Rstat2M(mean, std, unit='Earth', sample_size=1e3, grid_size=1e3, classify = 'No'): - """ - Forecast the mean and standard deviation of mass given the mean and standard deviation of the radius. - - Parameters - --------------- - mean: float - Mean (average) of radius. - std: float - Standard deviation of radius. - unit: string (optional) - Unit of the radius. Options are 'Earth' and 'Jupiter'. - sample_size: int (optional) - Number of radius samples to draw with the mean and std provided. - grid_size: int (optional) - Number of grid in the mass axis when sampling mass from radius. - The more the better results, but slower process. - Returns - --------------- - mean: float - Predicted mean of mass in the input unit. - std: float - Predicted standard deviation of mass. - """ - # unit - if unit == 'Earth': - pass - elif unit == 'Jupiter': - mean = mean * rearth2rjup - std = std * rearth2rjup - else: - print "Input unit must be 'Earth' or 'Jupiter'. Using 'Earth' as default." - - # draw samples - radius = truncnorm.rvs( (0.-mean)/std, np.inf, loc=mean, scale=std, size=sample_size) - if classify == 'Yes': - mass = Rpost2M(radius, 'Earth', grid_size, classify='Yes') - else: - mass = Rpost2M(radius, 'Earth', grid_size) - - if mass is None: - return None - - if unit=='Jupiter': - mass = mass / mearth2mjup - - m_med = np.median(mass) - onesigma = 34.1 - m_up = np.percentile(mass, 50.+onesigma, interpolation='nearest') - m_down = np.percentile(mass, 50.-onesigma, interpolation='nearest') - - return m_med, m_up - m_med, m_med - m_down - - - \ No newline at end of file + """ + Forecast the mass distribution given the radius distribution. + + Parameters + --------------- + radius: one dimensional array + The radius distribution. + unit: string (optional) + Unit of the mass. Options are 'Earth' and 'Jupiter'. + grid_size: int (optional) + Number of grid in the mass axis when sampling mass from radius. + The more the better results, but slower process. + classify: string (optional) + If you want the object to be classifed. + Options are 'Yes' and 'No'. Default is 'No'. + Result will be printed, not returned. + + Returns + --------------- + mass: one dimensional array + Predicted mass distribution in the input unit. + """ + + # unit + if unit == 'Earth': + pass + elif unit == 'Jupiter': + radius = radius * rearth2rjup + else: + print("Input unit must be 'Earth' or 'Jupiter'. " + "Using 'Earth' as default.") + + + # radius range + if np.min(radius) < 1e-1 or np.max(radius) > 1e2: + print('Radius range out of model expectation. Returning None.') + return None + + + + # sample_grid + if grid_size < 10: + print('The sample grid is too sparse. Using 10 sample grid instead.') + grid_size = 10 + + ## convert to mass + sample_size = len(radius) + logr = np.log10(radius) + logm = np.ones_like(logr) + + hyper_ind = np.random.randint(low=0, high=np.shape(all_hyper)[0], + size=sample_size) + hyper = all_hyper[hyper_ind,:] + + logm_grid = np.linspace(-3.522, 5.477, grid_size) + + for i in range(sample_size): + prob = ProbRGivenM(logr[i], logm_grid, hyper[i,:]) + logm[i] = np.random.choice(logm_grid, size=1, p = prob) + + mass_sample = 10.** logm + + if classify == 'Yes': + classification(logm, hyper[:,-3:]) + + ## convert to right unit + if unit == 'Jupiter': + mass = mass_sample / mearth2mjup + else: + mass = mass_sample + + return mass + + + +def Rstat2M(mean, std, unit='Earth', sample_size=1e3, grid_size=1e3, classify = 'No'): + """ + Forecast the mean and standard deviation of mass given the mean and standard deviation of the radius. + + Parameters + --------------- + mean: float + Mean (average) of radius. + std: float + Standard deviation of radius. + unit: string (optional) + Unit of the radius. Options are 'Earth' and 'Jupiter'. + sample_size: int (optional) + Number of radius samples to draw with the mean and std provided. + grid_size: int (optional) + Number of grid in the mass axis when sampling mass from radius. + The more the better results, but slower process. + Returns + --------------- + mean: float + Predicted mean of mass in the input unit. + std: float + Predicted standard deviation of mass. + """ + # unit + if unit == 'Earth': + pass + elif unit == 'Jupiter': + mean = mean * rearth2rjup + std = std * rearth2rjup + else: + print("Input unit must be 'Earth' or 'Jupiter'. " + "Using 'Earth' as default.") + + # draw samples + radius = truncnorm.rvs( (0.-mean)/std, np.inf, loc=mean, scale=std, + size=sample_size) + if classify == 'Yes': + mass = Rpost2M(radius, 'Earth', grid_size, classify='Yes') + else: + mass = Rpost2M(radius, 'Earth', grid_size) + + if mass is None: + return None + + if unit=='Jupiter': + mass = mass / mearth2mjup + + m_med = np.median(mass) + onesigma = 34.1 + m_up = np.percentile(mass, 50.+onesigma, interpolation='nearest') + m_down = np.percentile(mass, 50.-onesigma, interpolation='nearest') + + return m_med, m_up - m_med, m_med - m_down From 57e187eaac1066159d05155018c8247ebba54625 Mon Sep 17 00:00:00 2001 From: Brett Morris Date: Sat, 1 Jul 2017 12:27:04 -0700 Subject: [PATCH 2/2] Adding explicit import path for hyperparams --- mr_forecast.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mr_forecast.py b/mr_forecast.py index 12c0205..04a9546 100644 --- a/mr_forecast.py +++ b/mr_forecast.py @@ -4,6 +4,7 @@ import numpy as np from scipy.stats import truncnorm import h5py +import os ## constant mearth2mjup = 317.828 @@ -19,7 +20,7 @@ n_pop = 4 ## read parameter file -hyper_file = 'fitting_parameters.h5' +hyper_file = os.path.join(os.path.dirname(__file__), 'fitting_parameters.h5') h5 = h5py.File(hyper_file, 'r') all_hyper = h5['hyper_posterior'][:] h5.close()