Staeckelecc - documentation and no_median

AUTHORS.txt

@@ -19,6 +19,7 @@ Wilma Trick (KuzminKutuzovPotential and various other contributions)
 Anna Juranova (FerrersPotential)
 Semyeong Oh (FerrersPotential)
 Jack Hong (SpiralArmsPotential)
+Ted Mackereth (addition to estimateDeltaStaeckel function, and fast orbit characterisation docs)
 
 General help:
 ==============

HISTORY.txt

@@ -61,6 +61,10 @@ v1.3 (2017-XX-XX)
 - Changed default method for computing actions, frequencies, and
   angles for Orbit instances to be the Staeckel approximation with an
   automatically-estimated delta parameter.
+
+- Added an option to the estimateDeltaStaeckel function to facilitate the
+  return of an estimated delta parameter at every phase space point passed,
+  rather than returning a median of the estimate at each point.
 
 - Generalized actionAngleStaeckel to allow for different focal lengths
   delta for different phase-space points. Also allowed the order of

doc/source/actionAngle.rst

@@ -1,3 +1,5 @@
+.. _actionangle:
+
 Action-angle coordinates
 =========================
 
@@ -379,6 +381,8 @@ vertical action to approximately five percent.
 .. WARNING::
    Frequencies and angles using the adiabatic approximation are not implemented at this time.
 
+.. _actionanglestaeckel:
+
 Action-angle coordinates using the Staeckel approximation
 -----------------------------------------------------------
 

doc/source/examples/dierickx_eccentricities.py

@@ -0,0 +1,191 @@
+import os, os.path
+import subprocess
+from astropy.io import fits, ascii
+from astropy import units
+import numpy as np
+from optparse import OptionParser
+import sys
+import tempfile
+from ftplib import FTP
+import shutil
+import pickle
+from tqdm import tqdm
+from galpy.potential import LogarithmicHaloPotential
+from galpy.potential import evaluatePotentials as evalPot
+from galpy.orbit import Orbit
+from galpy.actionAngle import estimateDeltaStaeckel, actionAngleStaeckel, UnboundError
+from galpy.util import bovy_coords
+import matplotlib.pyplot as plt
+import numpy
+
+_ERASESTR= "                                                                                "
+
+def calc_eccentricity(args, options):
+    table = os.path.join(args[0],'table2.dat')
+    readme = os.path.join(args[0],'ReadMe')
+    dierickx = ascii.read(table, readme=readme)
+    vxvv = np.dstack([dierickx['RAdeg'], dierickx['DEdeg'], dierickx['Dist']/1e3, dierickx['pmRA'], dierickx['pmDE'], dierickx['HRV']])[0]
+    ro, vo, zo = 8., 220., 0.025
+    ra, dec= vxvv[:,0], vxvv[:,1]
+    lb= bovy_coords.radec_to_lb(ra,dec,degree=True)
+    pmra, pmdec= vxvv[:,3], vxvv[:,4]
+    pmllpmbb= bovy_coords.pmrapmdec_to_pmllpmbb(pmra,pmdec,ra,dec,degree=True)
+    d, vlos= vxvv[:,2], vxvv[:,5]
+    rectgal= bovy_coords.sphergal_to_rectgal(lb[:,0],lb[:,1],d,vlos,pmllpmbb[:,0], pmllpmbb[:,1],degree=True)
+    vsolar= np.array([-10.1,4.0,6.7])
+    vsun= np.array([0.,1.,0.,])+vsolar/vo
+    X = rectgal[:,0]/ro
+    Y = rectgal[:,1]/ro
+    Z = rectgal[:,2]/ro
+    vx = rectgal[:,3]/vo
+    vy = rectgal[:,4]/vo
+    vz = rectgal[:,5]/vo
+    vsun= np.array([0.,1.,0.,])+vsolar/vo
+    Rphiz= bovy_coords.XYZ_to_galcencyl(X,Y,Z,Zsun=zo/ro)
+    vRvTvz= bovy_coords.vxvyvz_to_galcencyl(vx,vy,vz,Rphiz[:,0],Rphiz[:,1],Rphiz[:,2],vsun=vsun,Xsun=1.,Zsun=zo/ro,galcen=True)
+    #do the integration and individual analytic estimate for each object
+    ts= np.linspace(0.,20.,10000)
+    lp= LogarithmicHaloPotential(normalize=1.)
+    e_ana = numpy.zeros(len(vxvv))
+    e_int = numpy.zeros(len(vxvv))
+    print('Performing orbit integration and analytic parameter estimates for Dierickx et al. sample...')
+    for i in tqdm(range(len(vxvv))):
+        try:
+            orbit = Orbit(vxvv[i], radec=True, vo=220., ro=8.)
+            e_ana[i] = orbit.e(analytic=True, pot=lp, c=True)
+        except UnboundError:
+            e_ana[i] = np.nan
+        orbit.integrate(ts, lp)
+        e_int[i] = orbit.e(analytic=False)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.scatter(e_int, e_ana,  s=1, color='Black', lw=0.)
+    plt.xlabel(r'$\mathrm{galpy\ integrated}\ e$')
+    plt.ylabel(r'$\mathrm{galpy\ analytic}\ e$')
+    plt.xlim(0.,1.)
+    plt.ylim(0.,1.)
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0],'dierickx-integratedeanalytice.png'), format='png', dpi=200)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.hist(e_int, bins=30)
+    plt.xlim(0.,1.)
+    plt.xlabel(r'$\mathrm{galpy}\ e$')
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0], 'dierickx-integratedehist.png'), format='png', dpi=200)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.scatter(dierickx['e'], e_int,  s=1, color='Black', lw=0.)
+    plt.xlabel(r'$\mathrm{Dierickx\ et\ al.}\ e$')
+    plt.ylabel(r'$\mathrm{galpy\ integrated}\ e$')
+    plt.xlim(0.,1.)
+    plt.ylim(0.,1.)
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0],'dierickx-integratedee.png'), format='png', dpi=200)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.scatter(dierickx['e'], e_ana,  s=1, color='Black', lw=0.)
+    plt.xlabel(r'$\mathrm{Dierickx\ et\ al.}\ e$')
+    plt.ylabel(r'$\mathrm{galpy\ estimated}\ e$')
+    plt.xlim(0.,1.)
+    plt.ylim(0.,1.)
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0],'dierickx-analyticee.png'), format='png', dpi=200)
+    arr = numpy.recarray(len(e_ana), dtype=[('analytic_e', float), ('integrated_e', float)])
+    arr['analytic_e'] = e_ana
+    arr['integrated_e'] = e_int
+    with open(os.path.join(args[0],'eccentricities.dat'), 'w') as file:
+        pickle.dump(arr, file)
+        file.close()
+
+def get_table(args,options):
+    cat = 'J/ApJ/725/L186/'
+    tab2name = 'table2.dat.gz'
+    tab2readme = 'ReadMe'
+    out = args[0]
+    ensure_dir(os.path.join(out,tab2name))
+    vizier(cat, os.path.join(out,tab2name), os.path.join(out,tab2readme), catalogname=tab2name, readmename=tab2readme)
+    subprocess.call(['gunzip', os.path.join(out,tab2name)])
+
+def vizier(cat,filePath,ReadMePath,
+           catalogname='catalog.dat',readmename='ReadMe'):
+    """
+    NAME:
+       vizier
+    PURPOSE:
+       download a catalog and its associated ReadMe from Vizier
+    INPUT:
+       cat - name of the catalog (e.g., 'III/272' for RAVE, or J/A+A/... for journal-specific catalogs)
+       filePath - path of the file where you want to store the catalog (note: you need to keep the name of the file the same as the catalogname to be able to read the file with astropy.io.ascii)
+       ReadMePath - path of the file where you want to store the ReadMe file
+       catalogname= (catalog.dat) name of the catalog on the Vizier server
+       readmename= (ReadMe) name of the ReadMe file on the Vizier server
+    OUTPUT:
+       (nothing, just downloads)
+    HISTORY:
+       2016-09-12 - Written - Bovy (UofT)
+    """
+    _download_file_vizier(cat,filePath,catalogname=catalogname)
+    _download_file_vizier(cat,ReadMePath,catalogname=readmename)
+    return None
+
+
+def _download_file_vizier(cat,filePath,catalogname='catalog.dat'):
+    '''
+    Stolen from Jo Bovy's gaia_tools package!
+    '''
+    sys.stdout.write('\r'+"Downloading file %s ...\r" \
+                         % (os.path.basename(filePath)))
+    sys.stdout.flush()
+    try:
+        # make all intermediate directories
+        os.makedirs(os.path.dirname(filePath)) 
+    except OSError: pass
+    # Safe way of downloading
+    downloading= True
+    interrupted= False
+    file, tmp_savefilename= tempfile.mkstemp()
+    os.close(file) #Easier this way
+    ntries= 1
+    while downloading:
+        try:
+            ftp= FTP('cdsarc.u-strasbg.fr')
+            ftp.login('anonymous', 'test')
+            ftp.cwd(os.path.join('pub','cats',cat))
+            with open(tmp_savefilename,'wb') as savefile:
+                ftp.retrbinary('RETR %s' % catalogname,savefile.write)
+            shutil.move(tmp_savefilename,filePath)
+            downloading= False
+            if interrupted:
+                raise KeyboardInterrupt
+        except:
+            raise
+            if not downloading: #Assume KeyboardInterrupt
+                raise
+            elif ntries > _MAX_NTRIES:
+                raise IOError('File %s does not appear to exist on the server ...' % (os.path.basename(filePath)))
+        finally:
+            if os.path.exists(tmp_savefilename):
+                os.remove(tmp_savefilename)
+        ntries+= 1
+    sys.stdout.write('\r'+_ERASESTR+'\r')
+    sys.stdout.flush()        
+    return None
+
+def ensure_dir(f):
+    """ Ensure a a file exists and if not make the relevant path """
+    d = os.path.dirname(f)
+    if not os.path.exists(d):
+        os.makedirs(d)
+
+def get_options():
+    #no options yet - probably none needed?
+    usage = "usage: %prog [options] <outpath>"
+    parser = OptionParser(usage=usage)
+    return parser
+
+if __name__ == '__main__':
+    parser = get_options()
+    options, args= parser.parse_args()
+    get_table(args,options)
+    calc_eccentricity(args, options)

doc/source/orbit.rst

@@ -346,6 +346,79 @@ behavior
 
 .. image:: images/lp-orbit-integration-EzJz.png
 
+.. _fastchar:
+
+**NEW in v1.3** Fast orbit characterization
+--------------------------------------------
+
+It is also possible to use galpy for the fast estimation of orbit parameters as demonstrated
+in Mackereth & Bovy (2018, in prep.) via the Staeckel approximation (originally used by `Binney (2012) <http://adsabs.harvard.edu/abs/2012MNRAS.426.1324B>`_
+for the appoximation of actions in axisymmetric potentials), without performing any orbit integration. 
+The method uses the geometry of the orbit tori to estimate the orbit parameters. After initialising 
+an ``Orbit`` instance, the method is applied by specifying ``analytic=True`` and 
+selecting ``type='staeckel'``.
+
+>>> o.e(analytic=True, type='staeckel')
+
+if running the above without integrating the orbit, the potential should also be specified
+in the usual way
+
+>>> o.e(analytic=True, type='staeckel', pot=mp)
+
+This interface automatically estimates the necessary delta parameter based on the initial 
+condition of the ``Orbit`` object.
+
+While this is useful and fast for individual ``Orbit`` objects, it is likely that users will
+want to rapidly evaluate the orbit parameters of large numbers of objects. It is possible
+to perform the orbital parameter estimation above through the :ref:`actionAngle <actionangle>` 
+interface. To do this, we need arrays of the phase-space points ``R``, ``vR``, ``vT``, ``z``, ``vz``, and 
+``phi`` for the objects.  The orbit parameters are then calculated by first 
+specifying an ``actionAngleStaeckel`` instance (this requires a single ``delta`` focal-length parameter, see :ref:`the documentation of the actionAngleStaeckel class <actionanglestaeckel>`), then using the 
+``EccZmaxRperiRap`` method with the data points:
+
+>>> aAS = actionAngleStaeckel(pot=mp, delta=0.4)
+>>> e, Zmax, rperi, rap = aAS.EccZmaxRperiRap(R, vR, vT, z, vz, phi)
+
+Alternatively, you can specify an array for ``delta`` when calling ``aAS.EccZmaxRperiRap``, for example by first estimating good ``delta`` parameters as follows:
+
+>>> from galpy.actionAngle import estimateDeltaStaeckel
+>>> delta = estimateDeltaStaeckel(mp, R, z, no_median=True)
+
+where ``no_median=True`` specifies that the function return the delta parameter at each given point
+rather than the median of the calculated deltas (which is the default option). Then one can compute the eccetrncity etc. using individual delta values as:
+
+>>> e, Zmax, rperi, rap = aAS.EccZmaxRperiRap(R, vR, vT, z, vz, phi, delta=delta)
+
+Th ``EccZmaxRperiRap`` method also exists for the ``actionAngleIsochrone``, 
+``actionAngleSpherical``, and ``actionAngleAdiabatic`` modules. 
+
+We can test the speed of this method in iPython by finding the parameters at 100000 steps 
+along an orbit in MWPotential2014, like this
+
+>>> o= Orbit(vxvv=[1.,0.1,1.1,0.,0.1,0.])
+>>> ts = numpy.linspace(0,100,100000)
+>>> o.integrate(ts,MWPotential2014)
+>>> aAS = actionAngleStaeckel(pot=MWPotential2014,delta=0.3) 
+>>> R, vR, vT, z, vz, phi = o.getOrbit().T
+>>> delta = estimateDeltaStaeckel(MWPotential2014, R, z, no_median=True)
+>>> %timeit -n 10 es, zms, rps, ras = aAS.EccZmaxRperiRap(R,vR,vT,z,vz,phi,delta=delta)
+#10 loops, best of 3: 899 ms per loop
+
+you can see that in this potential, each phase space point is calculated in roughly 9µs.
+further speed-ups can be gained by using the ``actionAngleStaeckelGrid`` module, which first
+calculates the parameters using a grid-based interpolation
+
+>>> from galpy.actionAngle import actionAngleStaeckelGrid
+>>> aASG= actionAngleStaeckelGrid(pot=mp,delta=0.4,nE=51,npsi=51,nLz=61,c=True,interpecc=True)
+>>> %timeit -n 10 es, zms, rps, ras = aASG.EccZmaxRperiRap(R,vR,vT,z,vz,phi)
+#10 loops, best of 3: 587 ms per loop
+
+where ``interpecc=True`` is required to perform the interpolation of the orbit parameter grid.
+Looking at how the eccentricity estimation varies along the orbit, and comparing to the calculation
+using the orbit integration, we see that the estimation good job
+
+.. image:: images/lp-orbit-integration-et.png
+	:scale: 40 % 
 
 Accessing the raw orbit
 -----------------------
@@ -507,32 +580,74 @@ a set of thick disk stars. We start by downloading the sample of SDSS
 SEGUE (`2009AJ....137.4377Y
 <http://adsabs.harvard.edu/abs/2009AJ....137.4377Y>`_) thick disk
 stars compiled by Dierickx et al. (`2010arXiv1009.1616D
-<http://adsabs.harvard.edu/abs/2010arXiv1009.1616D>`_) at
-
-http://www.mpia-hd.mpg.de/homes/rix/Data/Dierickx-etal-tab2.txt
+<http://adsabs.harvard.edu/abs/2010arXiv1009.1616D>`_) from CDS at `this 
+link <http://vizier.cfa.harvard.edu/viz-bin/Cat?cat=J%2FApJ%2F725%2FL186&target=http&>`_.
+Downloading the table and the ReadMe will allow you to read in the data using ``astropy.io.ascii``
+like so
+
+>>> from astropy.io import ascii
+>>> dierickx = ascii.read('table2.dat', readme='ReadMe')
+>>> vxvv = numpy.dstack([dierickx['RAdeg'], dierickx['DEdeg'], dierickx['Dist']/1e3, dierickx['pmRA'], dierickx['pmDE'], dierickx['HRV']])[0]
 
 After reading in the data (RA,Dec,distance,pmRA,pmDec,vlos; see above)
 as a vector ``vxvv`` with dimensions [6,ndata] we (a) define the
 potential in which we want to integrate the orbits, and (b) integrate
-each orbit and save its eccentricity (running this for all 30,000-ish
+each orbit and save its eccentricity as calculated analytically following the :ref:`Staeckel 
+approximation method <fastchar>` and by orbit integration (running this for all 30,000-ish
 stars will take about half an hour)
 
+>>> from galpy.actionAngle import UnboundError
+>>> ts= np.linspace(0.,20.,10000)
 >>> lp= LogarithmicHaloPotential(normalize=1.)
->>> ts= nu.linspace(0.,20.,10000)
->>> mye= nu.zeros(ndata)
->>> for ii in range(len(e)):
-...         o= Orbit(vxvv[ii,:],radec=True,vo=220.,ro=8.) #Initialize
-...         o.integrate(ts,lp) #Integrate
-...         mye[ii]= o.e() #Calculate eccentricity
+>>> e_ana = numpy.zeros(len(vxvv))
+>>> e_int = numpy.zeros(len(vxvv))
+>>> for i in range(len(vxvv)):
+...	#calculate analytic e estimate, catch any 'unbound' orbits
+...     try:
+...         orbit = Orbit(vxvv[i], radec=True, vo=220., ro=8.)
+...         e_ana[i] = orbit.e(analytic=True, pot=lp, c=True)
+...     except UnboundError:
+...         #parameters cannot be estimated analytically
+...         e_ana[i] = np.nan
+...     #integrate the orbit and return the numerical e value
+...     orbit.integrate(ts, lp)
+...     e_int[i] = orbit.e(analytic=False)
+
+We then find the following eccentricity distribution (from the numerical eccentricities)
+
+.. image:: images/dierickx-integratedehist.png
+	:scale: 40 %
+
+The eccentricity calculated by integration in galpy compare well with those
+calculated by Dierickx et al., except for a few objects
 
-We then find the following eccentricity distribution
+.. image:: images/dierickx-integratedee.png
+	:scale: 40 %
 
-.. image:: images/dierickx-myehist.png
+and the analytical estimates are equally as good:
 
-The eccentricity calculated by galpy compare well with those
-calculated by Dierickx et al., except for a few objects
+.. image:: images/dierickx-analyticee.png
+	:scale: 40 %
+
+In comparing the analytic and integrated eccentricity estimates - one can see that in this case
+the estimation is almost exact, due to the spherical symmetry of the chosen potential:
+
+.. image:: images/dierickx-integratedeanalytice.png
+	:scale: 40 %
+
+A script that calculates and plots everything can be downloaded
+:download:`here <examples/dierickx_eccentricities.py>`. To generate the plots just run::
+
+    python dierickx_eccentricities.py ../path/to/folder
+
+specifiying the location you want to put the plots and data.
+
+Alternatively - one can transform the observed coordinates into spherical coordinates and perform 
+the estimations in one batch using the ``actionAngle`` interface, which takes considerably less time:
 
-.. image:: images/dierickx-myee.png
+>>> from galpy import actionAngle
+>>> deltas = actionAngle.estimateDeltaStaeckel(lp, Rphiz[:,0], Rphiz[:,2], no_median=True)
+>>> aAS = actionAngleStaeckel(pot=lp, delta=0.)
+>>> par = aAS.EccZmaxRperiRap(Rphiz[:,0], vRvTvz[:,0], vRvTvz[:,1], Rphiz[:,2], vRvTvz[:,2], Rphiz[:,1], delta=deltas)
 
-The script that calculates and plots everything can be downloaded
-:download:`here <examples/dierickx-edist.py>`.
+The above code calculates the parameters in roughly 100ms on a single core.