Skip to content

Instantly share code, notes, and snippets.

@adrn

adrn/Untitled1.ipynb

Created October 16, 2020 14:46
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save adrn/d64eda67fa073fadd4b40b05f081f609 to your computer and use it in GitHub Desktop.
Save adrn/d64eda67fa073fadd4b40b05f081f609 to your computer and use it in GitHub Desktop.
Download ZIP
projects/thejoker-notebooks/Untitled1.ipynb
Display the source blob
Display the rendered blob
Raw
Untitled1.ipynb
{
"cells": [
{
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-16T14:33:26.984232Z",
"start_time": "2020-10-16T14:33:25.655723Z"
},
"trusted": true
},
"cell_type": "code",
"source": "from astropy.io import ascii\nfrom astropy.time import Time\nimport astropy.units as u\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\n\nimport pymc3 as pm\nimport exoplanet.units as xu\nimport exoplanet as xo\nfrom exoplanet.distributions import Angle\nimport corner\n\nimport thejoker as tj",
"execution_count": 2,
"outputs": []
},
{
"metadata": {
"ExecuteTime": {
"start_time": "2020-10-16T14:40:07.383608Z",
"end_time": "2020-10-16T14:40:07.397488Z"
},
"trusted": true
},
"cell_type": "code",
"source": "tbl = ascii.read(\n \"\"\"2457945.5419639 -119.5 6.5\n 2457945.6301728 -16.9 7.2\n 2457956.5212405 -123.1 6.8\n 2457956.5352468 -107.9 6.4\n 2457956.5546942 -95.8 8.1\n 2457956.5687017 -65.6 8.2\n 2457956.5948590 -57.2 7.0\n 2457956.6088705 -38.0 6.0\n 2457966.4981876 -27.2 7.6\n 2457966.5122271 -18.4 20.1\n 2457966.5399208 38.7 7.2\n 2457966.5608956 49.6 8.7\n 2457966.5914866 115.1 10.2\n 2457967.5998999 -1.5 12.8\n 2457996.4786139 231.8 8.4\n 2457996.4931540 238.7 8.0\n 2457996.5105239 228.8 9.4\n 2458001.4795192 -144.0 11.4\n 2458001.4935388 -124.6 10.2\n 2458139.7564130 -136.8 7.7\n 2458139.7704090 -117.1 6.8\n 2458140.7400266 -194.5 9.8\n 2458140.7540222 -166.4 8.3\n 2458140.7770050 -150.8 7.4\n 2458140.7910007 -158.5 8.2\n 2458161.6721983 -121.2 7.9\n 2458161.6872977 -107.7 7.3\n 2458202.5822163 177.8 8.5\n 2458202.5962046 197.7 9.6\n 2458202.8350917 190.9 7.4\n 2458202.8490793 176.1 7.6\n 2458223.5928781 231.9 7.7\n 2458223.8317433 41.2 8.4\n 2458223.8456941 32.1 14.9\n 2458243.6779211 -87.6 8.6\n 2458243.6919415 -112.3 10.5\n 2458243.7115353 -125.5 9.5\n 2458243.7263475 -141.9 10.2\n 2458243.7459698 -130.3 10.9\n 2458247.5062024 131.9 11.5\n 2458247.5435496 160.5 14.2\n 2458278.5472619 197.1 15.9\n 2458278.5613912 183.7 15.7\"\"\",\n names=['BJD', 'rv', 'rv_err'])\ntbl['rv'].unit = u.km/u.s\ntbl['rv_err'].unit = u.km/u.s",
"execution_count": 7,
"outputs": []
},
{
"metadata": {
"ExecuteTime": {
"start_time": "2020-10-16T14:40:15.641982Z",
"end_time": "2020-10-16T14:40:15.647804Z"
},
"trusted": true
},
"cell_type": "code",
"source": "data = tj.RVData(\n t=Time(tbl['BJD'], format='jd', scale='tcb'),\n rv=u.Quantity(tbl['rv']),\n rv_err=u.Quantity(tbl['rv_err']))",
"execution_count": 8,
"outputs": []
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-16T14:33:30.305200Z",
"start_time": "2020-10-16T14:33:30.302436Z"
},
"trusted": true
},
"cell_type": "code",
"source": "# set up a random number generator to ensure reproducibility\nseed = 42\nrnd = np.random.default_rng(seed=seed)",
"execution_count": 3,
"outputs": []
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-16T14:33:47.163338Z",
"start_time": "2020-10-16T14:33:40.695860Z"
},
"trusted": true
},
"cell_type": "code",
"source": "with pm.Model() as model:\n # Allow extra error to account for under-estimated error bars\n e = xu.with_unit(pm.Constant('e', 0),\n u.one)\n\n prior = tj.JokerPrior.default(\n P_min=0.1*u.day, P_max=100*u.day, # Range of periods to consider\n sigma_K0=50*u.km/u.s, P0=1*u.year, # scale of the prior on semiamplitude, K\n sigma_v=50*u.km/u.s, # std dev of the prior on the systemic velocity, v0\n pars={'e': e}\n )",
"execution_count": 4,
"outputs": []
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-16T14:33:51.790850Z",
"start_time": "2020-10-16T14:33:51.109700Z"
},
"trusted": true
},
"cell_type": "code",
"source": "prior.sample(size=100)",
"execution_count": 5,
"outputs": [
{
"data": {
"text/plain": "<JokerSamples [P, e, omega, M0, s] (100 samples)>"
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"ExecuteTime": {
"start_time": "2020-10-16T14:40:43.718537Z",
"end_time": "2020-10-16T14:40:43.898792Z"
},
"trusted": true
},
"cell_type": "code",
"source": "joker = tj.TheJoker(prior, random_state=rnd)\nsamples = joker.rejection_sample(data,\n prior_samples=1_000,\n max_posterior_samples=256)\nsamples",
"execution_count": 10,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 10,
"data": {
"text/plain": "<JokerSamples [P, e, omega, M0, s, K, v0] (1 samples)>"
},
"metadata": {}
}
]
},
{
"metadata": {
"ExecuteTime": {
"start_time": "2020-10-16T14:44:46.053985Z",
"end_time": "2020-10-16T14:44:46.057363Z"
},
"trusted": true
},
"cell_type": "code",
"source": "import os\nimport pathlib\n\ndef tj_rvs(data, prior_samples_filename, joker_samples_filename, \n overwrite=False):\n # TODO: Construct the prior \n \n if not os.path.exists(prior_samples_filename) or overwrite:\n prior_samples = prior.sample(...)\n prior_samples.write(prior_samples_filename)\n \n prior_samples = JokerSamples.read(prior_samples_filename)\n \n # TODO: construct TheJoker instance\n \n if not os.path.exists(joker_samples_filename) or overwrite:\n joker_samples, logprobs = thejoker.rejection_sample(...)\n joker_samples.write(joker_samples_filename)\n \n # TODO: also save logprobs\n \n return joker_samples, logprobs",
"execution_count": 11,
"outputs": []
},
{
"metadata": {
"ExecuteTime": {
"start_time": "2020-10-16T14:45:44.941482Z",
"end_time": "2020-10-16T14:45:44.954959Z"
},
"trusted": true
},
"cell_type": "code",
"source": "project_root = pathlib.WindowsPath(r'C:\\Users\\apricewhelan\\projects\\thejoker')\nproject_root / 'my_joker_samples.hdf5'",
"execution_count": 15,
"outputs": [
{
"output_type": "error",
"ename": "NotImplementedError",
"evalue": "cannot instantiate 'WindowsPath' on your system",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-15-13885e13d9f5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mproject_root\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpathlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mWindowsPath\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr'C:\\Users\\apricewhelan\\projects\\thejoker'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mproject_root\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'my_joker_samples.hdf5'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda/lib/python3.8/pathlib.py\u001b[0m in \u001b[0;36m__new__\u001b[0;34m(cls, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1038\u001b[0m \u001b[0mself\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_from_parts\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1039\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flavour\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_supported\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1040\u001b[0;31m raise NotImplementedError(\"cannot instantiate %r on your system\"\n\u001b[0m\u001b[1;32m 1041\u001b[0m % (cls.__name__,))\n\u001b[1;32m 1042\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_init\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNotImplementedError\u001b[0m: cannot instantiate 'WindowsPath' on your system"
]
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"kernelspec": {
"name": "conda-root-py",
"display_name": "Python [conda env:root] *",
"language": "python"
},
"language_info": {
"name": "python",
"version": "3.8.3",
"mimetype": "text/x-python",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"pygments_lexer": "ipython3",
"nbconvert_exporter": "python",
"file_extension": ".py"
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"base_numbering": 1,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
},
"gist": {
"id": "",
"data": {
"description": "projects/thejoker-notebooks/Untitled1.ipynb",
"public": true
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment