Merge pull request #73 from California-Planet-Search/next-release

Release 0.9.5

README.md

@@ -1,6 +1,6 @@
 # RadVel
 
-General Toolkit for Modeling Radial Velocities. Version 0.9.4
+General Toolkit for Modeling Radial Velocities.
 
 ## Attribution
 

radvel/__init__.py

@@ -11,7 +11,7 @@
 import os 
 import sys
 
-__version__ = '0.9.4'
+__version__ = '0.9.5'
 
 MODULEDIR, filename = os.path.split(__file__)
 DATADIR = os.path.join(sys.prefix,'radvel_example_data')

radvel/basis.py

@@ -1,4 +1,4 @@
-import numpy as np        
+import numpy as np
 import copy
 import pandas as pd
 from collections import OrderedDict
@@ -11,7 +11,8 @@
 'per tc secosw sesinw logk',
 'per tc secosw sesinw k',
 'per tc ecosw esinw k',
-'per tc e w k'
+'per tc e w k',
+'logper tc secosw sesinw k'
 
 ]
 
@@ -48,7 +49,8 @@ class Basis(object):
         'per tc secosw sesinw logk'  \n 
         'per tc secosw sesinw k'  \n
         'per tc ecosw esinw k'  \n
-        'per tc e w k'
+        'per tc e w k' \n
+        'logper tc secosw sesinw logk'
     """
 
     cps_params = 'per tp e w k'.split()
@@ -149,6 +151,20 @@ def _delpar(key):
                 w = np.arctan2(sesinw , secosw)
                 tp = timetrans_to_timeperi(tc, per, e, w)
 
+            if basis_name=='logper tc secosw sesinw k':
+                # pull out parameters
+                logper = _getpar('logper')
+                tc = _getpar('tc')
+                secosw = _getpar('secosw')
+                sesinw = _getpar('sesinw')
+                k = _getpar('k')
+
+                # transform into CPS basis
+                per = np.exp(logper)
+                e = secosw**2 + sesinw**2
+                w = np.arctan2(sesinw , secosw)
+                tp = timetrans_to_timeperi(tc, per, e, w)
+
             if basis_name=='per tc ecosw esinw k':
                 # pull out parameters
                 per = _getpar('per')
@@ -275,6 +291,31 @@ def _delpar(key):
                 self.name = newbasis
                 self.params = newbasis.split()
 
+            if newbasis == 'logper tc secosw sesinw k':
+                per = _getpar('per')
+                e = _getpar('e')
+                w = _getpar('w')
+                k = _getpar('k')
+                try:
+                    tp = _getpar('tp')
+                except KeyError:
+                    tp = timetrans_to_timeperi(_getpar('tc'), per, e, w)
+                    _setpar('tp', tp)
+                _setpar('logper', np.log(per))
+                _setpar('secosw', np.sqrt(e)*np.cos(w) )
+                _setpar('sesinw', np.sqrt(e)*np.sin(w) )
+                _setpar('k', k )
+                _setpar('tc', timeperi_to_timetrans(tp, per, e, w) )
+
+                if not kwargs.get('keep', True):
+                    _delpar('per')
+                    _delpar('tp')
+                    _delpar('e')
+                    _delpar('w')
+
+                self.name = newbasis
+                self.params = newbasis.split()
+
             if newbasis == 'per tc ecosw esinw k':
                 per = _getpar('per')
                 e = _getpar('e')

radvel/cli.py

@@ -14,7 +14,12 @@
 
 def main():
     psr = ArgumentParser(
-        description="RadVel: The Radial Velocity Toolkit", prog='radvel'
+        description="RadVel: The Radial Velocity Toolkit", prog='RadVel'
+    )
+    psr.add_argument('--version',
+        action='version',
+        version="%(prog)s {}".format(radvel.__version__),
+        help="Print version number and exit."
     )
 
     subpsr = psr.add_subparsers(title="subcommands", dest='subcommand')
@@ -38,6 +43,8 @@ def main():
         help="Include decorrelation in likelihood."
     )
 
+
+
     # Fitting    
     psr_fit = subpsr.add_parser(
         'fit', parents=[psr_parent],

radvel/mcmc.py

@@ -4,7 +4,6 @@
 import copy
 import pp
 import threading
-from multiprocessing.pool import Pool
 from scipy import optimize
 import sys
 import time
@@ -147,6 +146,9 @@ def _crunch(sampler, ipos, checkinterval):
         if par.startswith('per'):
             pscale = np.abs(val * 1e-5*np.log10(val))
             pscale_per = pscale
+        elif par.startswith('logper'):
+            pscale = np.abs(1e-5 * val)
+            pscale_per = pscale
         elif par.startswith('tc'):
             pscale = 0.1
         else:
@@ -204,14 +206,6 @@ def _crunch(sampler, ipos, checkinterval):
         ch = CheckThread(convergence_check, statevars.server, statevars.samplers)
         ch.start()
 
-        # Use multiprocessing
-        # result = pool.apply_async(convergence_check,
-        #                 (statevars.server, statevars.samplers))
-
-        # ch = CheckThread(status_message, statevars)
-        # ch.start()
-
-        #convergence_check(statevars.server, statevars.samplers)
         # Burn-in complete after maximum G-R statistic first reaches burnGR
         # reset samplers
         if not statevars.burn_complete and statevars.maxgr <= burnGR:

radvel/model.py

@@ -9,6 +9,7 @@
 
 texdict = {
     'per': 'P',
+    'logper': '\\ln{P}',
     'logk': '\\ln{K}',
     'tc': 'T\\rm{conj}',
     'secosw': '\\sqrt{e}\\cos{\\omega}',

radvel/posterior.py

@@ -59,8 +59,10 @@ def bic(self):
             float: BIC
         """
 
-        return -2.0 * self.logprob() + len(self.likelihood.get_vary_params()) + np.log(len(self.likelihood.y))
-
+        n = len(self.likelihood.y)
+        k = len(self.likelihood.get_vary_params())
+        _bic = np.log(n) * k - 2.0 * self.logprob() 
+        return _bic
 
     def logprob_array(self, params_array):
         """Log probability for parameter vector

tests/K2-24_Fitting+MCMC.ipynb

@@ -361,64 +361,6 @@
     "print post"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Run the MCMC again"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "df = radvel.mcmc(post,threads=1,nwalkers=20,nrun=1000)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Convert into more physical parameters and make another corner plot"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "df['k1'] = np.exp(df['logk1'])\n",
-    "df['k2'] = np.exp(df['logk2'])\n",
-    "df['jit'] = np.exp(df['jit'])\n",
-    "df['e1'] = df['secosw1']**2 + df['sesinw1']**2\n",
-    "df['e2'] = df['secosw2']**2 + df['sesinw2']**2\n",
-    "\n",
-    "df['ecosw1'] = df['secosw1']*np.sqrt(df['e1'])\n",
-    "df['esinw1'] = df['sesinw1']*np.sqrt(df['e1'])\n",
-    "\n",
-    "df['ecosw2'] = df['secosw2']*np.sqrt(df['e2'])\n",
-    "df['esinw2'] = df['sesinw2']*np.sqrt(df['e2'])\n",
-    "\n",
-    "\n",
-    "labels = 'k1 k2 jit e1 e2'.split()\n",
-    "fig = corner.corner(\n",
-    "    df[labels],\n",
-    "    labels=labels,\n",
-    "    levels=[0.68,0.95],\n",
-    "    plot_datapoints=False,\n",
-    "    smooth=True,\n",
-    "    bins=20,\n",
-    "    quantiles=[.14,.5,.84]\n",
-    "    )"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -463,7 +405,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython2",
-   "version": "2.7.11"
+   "version": "2.7.13"
   }
  },
  "nbformat": 4,