Merge pull request #15 from telegraphic/dev

PyGDSM 1.4.0

README.md

@@ -71,7 +71,7 @@ Q & A
      the two outputs are indeed comparable. 
 
 **Q. What's the difference between this and the [Zheng et. al. github repo](https://github.com/jeffzhen/gsm2016)?**
-     `pygdsm` provides two classes: `GlobalSkyModel2016()` and `GSMObserver2016()`, which once instantiated
+     `pygdsm` provides two classes: `GlobalSkyModel16()` and `GSMObserver16()`, which once instantiated
      provide methods for programatically generating sky models. The Zheng et. al. github repo is a 
      simple, low-dependency, command line tool. Have a look at the 
      [GSM2016 quickstart guide](http://nbviewer.ipython.org/github/telegraphic/PyGDSM/blob/master/docs/pygdsm2016_quickstart.ipynb)

docs/pygdsm2016_quickstart.ipynb → docs/pygdsm16_quickstart.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {
     "collapsed": false,
     "jupyter": {
@@ -37,13 +37,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {
     "tags": []
    },
    "outputs": [],
    "source": [
-    "from pygdsm import GlobalSkyModel2016\n",
+    "from pygdsm import GlobalSkyModel16\n",
     "from pygdsm import GlobalSkyModel"
    ]
   },
@@ -90,7 +90,7 @@
     }
    ],
    "source": [
-    "gsm_2016 = GlobalSkyModel2016(freq_unit='MHz')\n",
+    "gsm_2016 = GlobalSkyModel16(freq_unit='MHz')\n",
     "gsm_2016.generate(500)\n",
     "gsm_2016.view(logged=True)\n",
     "\n",
@@ -164,7 +164,7 @@
     }
    ],
    "source": [
-    "gsm = GlobalSkyModel2016(freq_unit='GHz', data_unit='MJysr', resolution='low')\n",
+    "gsm = GlobalSkyModel16(freq_unit='GHz', data_unit='MJysr', resolution='low')\n",
     "gsm.generate(30) # Generate at 30 GHz\n",
     "gsm.view(logged=False)"
    ]
@@ -224,7 +224,7 @@
     }
    ],
    "source": [
-    "gsm = GlobalSkyModel2016(freq_unit='GHz')\n",
+    "gsm = GlobalSkyModel16(freq_unit='GHz')\n",
     "freqs = np.linspace(0.1, 4, 10)\n",
     "map_cube = gsm.generate(freqs)\n",
     "\n",
@@ -262,12 +262,12 @@
    },
    "outputs": [],
    "source": [
-    "from pygdsm import GSMObserver2016, GSMObserver\n",
+    "from pygdsm import GSMObserver16, GSMObserver\n",
     "from datetime import datetime\n",
     "\n",
     "# Setup observatory location - in this case, Parkes Australia\n",
     "(latitude, longitude, elevation) = ('-32.998370', '148.263659', 100)\n",
-    "ov = GSMObserver2016()\n",
+    "ov = GSMObserver16()\n",
     "ov.lon = longitude\n",
     "ov.lat = latitude\n",
     "ov.elev = elevation\n",
@@ -350,6 +350,31 @@
    "source": [
     "d = ov.view_observed_gsm()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gsm = GlobalSkyModel16(freq_unit='GHz')\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "GlobalSKyModel16?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -368,7 +393,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,

pygdsm/__init__.py

@@ -1,9 +1,11 @@
 
-from .pygdsm import GlobalSkyModel
-from .pygdsm import GSMObserver
+from .gsm08 import GlobalSkyModel
+from .gsm08 import GSMObserver
 
-from .pygdsm2016 import GlobalSkyModel2016
-from .pygdsm2016 import GSMObserver2016
+GlobalSkyModel08 = GlobalSkyModel
+
+from .gsm16 import GlobalSkyModel16
+from .gsm16 import GSMObserver16
 
 from .lfsm import LowFrequencySkyModel
 from .lfsm import LFSMObserver

pygdsm/pygdsm2016.py → pygdsm/gsm16.py

@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.interpolate import interp1d, pchip
 import h5py
 from astropy import units
 import healpy as hp
@@ -44,11 +45,11 @@ def rotate_map(hmap, rot_theta, rot_phi, nest=True):
     return rot_map
 
 
-class GlobalSkyModel2016(BaseSkyModel):
+class GlobalSkyModel16(BaseSkyModel):
     """ Global sky model (GSM) class for generating sky models.
     """
 
-    def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot=0, phi_rot=0):
+    def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot=0, phi_rot=0, interpolation='pchip'):
         """ Global sky model (GSM) class for generating sky models.
 
         Upon initialization, the map PCA data are loaded into memory and interpolation
@@ -63,9 +64,19 @@ def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot
         resolution: 'hi' or 'low'
             Resolution of output map. Either 300 arcmin (low) or 24 arcmin (hi).
             For frequencies under 10 GHz, output is 48 arcmin.
+        interpolation: 'cubic' or 'pchip'
+            Choose whether to use cubic spline interpolation or
+            piecewise cubic hermitian interpolating polynomial (PCHIP).
+            PCHIP is designed to never locally overshoot data, whereas
+            splines are designed to have smooth first and second derivatives.
 
         Notes
         -----
+        The scipy `interp1d` function does not allow one to explicitly
+        set second derivatives to zero at the endpoints, as is done in
+        the original GSM. As such, results will differ. Further, we default
+        to use PCHIP interpolation.
+
 
         """
 
@@ -79,8 +90,9 @@ def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot
         else:
             raise RuntimeError("RESOLUTION ERROR: Must be either hi or low, not %s" % resolution)
 
-        super(GlobalSkyModel2016, self).__init__('GSM2016', GSM2016_FILEPATH, freq_unit, data_unit, basemap='')
+        super(GlobalSkyModel16, self).__init__('GSM2016', GSM2016_FILEPATH, freq_unit, data_unit, basemap='')
 
+        self.interpolation_method = interpolation
         self.resolution = resolution
 
         # Map data to load
@@ -137,27 +149,46 @@ def generate(self, freqs):
 
         spec_nf = self.spec_nf
         nfreq = spec_nf.shape[1]
+
+        # Now borrow code from the orignal GSM2008 model to do a sensible interpolation
+
+        pca_freqs_ghz = spec_nf[0]
+        pca_scaling   = spec_nf[1]
+        pca_comps     = spec_nf[2:]
+         # Interpolate to the desired frequency values
+        ln_pca_freqs = np.log(pca_freqs_ghz)
+        if self.interpolation_method == 'cubic':
+            spl_scaling = interp1d(ln_pca_freqs, np.log(pca_scaling), kind='cubic')
+            spl1 = interp1d(ln_pca_freqs,   pca_comps[0],   kind='cubic')
+            spl2 = interp1d(ln_pca_freqs,   pca_comps[1],   kind='cubic')
+            spl3 = interp1d(ln_pca_freqs,   pca_comps[2],   kind='cubic')
+            spl4 = interp1d(ln_pca_freqs,   pca_comps[3],   kind='cubic')
+            spl5 = interp1d(ln_pca_freqs,   pca_comps[4],   kind='cubic')
+            spl6 = interp1d(ln_pca_freqs,   pca_comps[5],   kind='cubic')
 
-        output = np.zeros((len(freqs_ghz), map_ni.shape[1]), dtype='float32')
+        else:
+            spl_scaling = pchip(ln_pca_freqs, np.log(pca_scaling))
+            spl1 = pchip(ln_pca_freqs,   pca_comps[0])
+            spl2 = pchip(ln_pca_freqs,   pca_comps[1])
+            spl3 = pchip(ln_pca_freqs,   pca_comps[2])
+            spl4 = pchip(ln_pca_freqs,   pca_comps[3])
+            spl5 = pchip(ln_pca_freqs,   pca_comps[4])
+            spl6 = pchip(ln_pca_freqs,   pca_comps[5])
+
+        self.interp_comps = (spl_scaling, spl1, spl2, spl3, spl4, spl5, spl6)
+
+        ln_freqs = np.log(freqs_ghz)
+        comps = np.row_stack((spl1(ln_freqs), spl2(ln_freqs), spl3(ln_freqs), spl4(ln_freqs), spl5(ln_freqs), spl6(ln_freqs)))
+        scaling = np.exp(spl_scaling(ln_freqs))
+
+        # Finally, compute the dot product via einsum (awesome function)
+        # c=comp, f=freq, p=pixel. We want to dot product over c for each freq
+        #print comps.shape, self.pca_map_data.shape, scaling.shape
+
+        output = np.single(np.einsum('cf,pc,f->fp', comps, map_ni.T, scaling))
+
         for ifreq, freq in enumerate(freqs_ghz):
 
-            left_index = -1
-            for i in range(nfreq - 1):
-                if spec_nf[0, i] <= freq <= spec_nf[0, i + 1]:
-                    left_index = i
-                    break
-
-            # Do the interpolation
-            interp_spec_nf = np.copy(spec_nf)
-            interp_spec_nf[0:2] = np.log10(interp_spec_nf[0:2])
-            x1 = interp_spec_nf[0, left_index]
-            x2 = interp_spec_nf[0, left_index + 1]
-            y1 = interp_spec_nf[1:, left_index]
-            y2 = interp_spec_nf[1:, left_index + 1]
-            x = np.log10(freq)
-            interpolated_vals = (x * (y2 - y1) + x2 * y1 - x1 * y2) / (x2 - x1)
-            output[ifreq] = np.sum(10.**interpolated_vals[0] * (interpolated_vals[1:, None] * map_ni), axis=0)
-
             output[ifreq] = hp.pixelfunc.reorder(output[ifreq], n2r=True)
 
             if self.data_unit == 'TCMB':
@@ -181,11 +212,10 @@ def generate(self, freqs):
         return output
 
 
-class GSMObserver2016(BaseObserver):
+class GSMObserver16(BaseObserver):
     def __init__(self):
         """ Initialize the Observer object.
 
         Calls ephem.Observer.__init__ function and adds on gsm
         """
-        super(GSMObserver2016, self).__init__(gsm=GlobalSkyModel2016)
-
+        super(GSMObserver16, self).__init__(gsm=GlobalSkyModel16)

setup.py

@@ -25,7 +25,7 @@
     # Versions should comply with PEP440.  For a discussion on single-sourcing
     # the version across setup.py and the project code, see
     # https://packaging.python.org/en/latest/single_source_version.html
-    version='1.3.1',
+    version='1.4.0',
 
     description='Python Global Sky Model of diffuse Galactic radio emission',
     long_description=long_description,
@@ -47,7 +47,7 @@
         #   3 - Alpha
         #   4 - Beta
         #   5 - Production/Stable
-        'Development Status :: 4 - Beta',
+        'Development Status :: 5 - Production/Stable',
 
         # Indicate who your project is intended for
         'Intended Audience :: Science/Research',
@@ -61,6 +61,10 @@
         'Programming Language :: Python :: 3',
         'Programming Language :: Python :: 3.6',
         'Programming Language :: Python :: 3.7',
+        'Programming Language :: Python :: 3.8',
+        'Programming Language :: Python :: 3.9',
+        'Programming Language :: Python :: 3.10',
+        'Programming Language :: Python :: 3.11',
     ],
 
     # What does your project relate to?
@@ -77,20 +81,4 @@
     install_requires=requirements,
     tests_require=test_requirements,
 
-    # List additional groups of dependencies here (e.g. development
-    # dependencies). You can install these using the following syntax,
-    # for example:
-    # $ pip install -e .[dev,test]
-    extras_require={
-    },
-
-
-    # To provide executable scripts, use entry points in preference to the
-    # "scripts" keyword. Entry points provide cross-platform support and allow
-    # pip to create the appropriate form of executable for the target platform.
-    #entry_points={
-    #    'console_scripts': [
-    #        'sample=sample:main',
-    #    ],
-    #},
 )

tests/test_gsm2016.py

@@ -3,17 +3,17 @@
 
 import pylab as plt
 import healpy as hp
+import numpy as np
 from datetime import datetime
 
 import pytest
 
-from pygdsm import GlobalSkyModel2016, GSMObserver
-from pygdsm import GlobalSkyModel
-from pygdsm import GSMObserver2016
+from pygdsm import GlobalSkyModel16, GSMObserver16
+from pygdsm import GlobalSkyModel, GSMObserver
 
 
 def test_compare_gsm_to_old():
-    g = GlobalSkyModel2016(freq_unit='GHz', resolution='hi', data_unit='TCMB')
+    g = GlobalSkyModel16(freq_unit='GHz', resolution='hi', data_unit='TCMB')
     d = g.generate(0.408)
     g.view()
 
@@ -26,7 +26,7 @@ def test_observer_test():
 
     # Setup observatory location - in this case, Parkes Australia
     (latitude, longitude, elevation) = ('-32.998370', '148.263659', 100)
-    ov = GSMObserver2016()
+    ov = GSMObserver16()
     ov.lon = longitude
     ov.lat = latitude
     ov.elev = elevation
@@ -45,26 +45,46 @@ def test_observer_test():
     d = ov.view(logged=True)
     plt.show()
 
+def test_interp():
+    f = np.arange(40, 80, 5)
+    for interp in ('pchip', 'cubic'):
+        for SkyModel in (GlobalSkyModel, GlobalSkyModel16):
+            name = str(SkyModel).strip("<>").split('.')[-1].strip("' ")
+            gsm = SkyModel(freq_unit='MHz', interpolation=interp)
+            d = gsm.generate(f)
+
+            sky_spec = d.mean(axis=1)
+            fit = np.poly1d(np.polyfit(f, sky_spec, 5))(f)
+
+            plt.plot(f, sky_spec - fit, label=f'{name}: {interp}')
+
+    plt.xlabel("Frequency [MHz]")
+    plt.ylabel("Residual [K]")
+    plt.legend()
+    plt.show()
+
+
 
 def test_gsm_opts():
-    g = GlobalSkyModel2016(freq_unit='GHz', resolution='hi', data_unit='TCMB')
+    g = GlobalSkyModel16(freq_unit='GHz', resolution='hi', data_unit='TCMB')
     d = g.generate(0.408)
-    g = GlobalSkyModel2016(freq_unit='GHz', resolution='lo', data_unit='MJysr')
+    g = GlobalSkyModel16(freq_unit='GHz', resolution='lo', data_unit='MJysr')
     d = g.generate(0.408)
-    g = GlobalSkyModel2016(freq_unit='MHz', resolution='lo', data_unit='TRJ')
+    g = GlobalSkyModel16(freq_unit='MHz', resolution='lo', data_unit='TRJ')
     d = g.generate(408)
 
     with pytest.raises(RuntimeError):
         g.generate(5e12)
 
     with pytest.raises(RuntimeError):
-        g = GlobalSkyModel2016(resolution='oh_hai')
+        g = GlobalSkyModel16(resolution='oh_hai')
 
     with pytest.raises(RuntimeError):
-        g = GlobalSkyModel2016(data_unit='furlongs/fortnight')
+        g = GlobalSkyModel16(data_unit='furlongs/fortnight')
 
 
 if __name__ == "__main__":
     test_compare_gsm_to_old()
     test_observer_test()
-    test_gsm_opts()
+    test_gsm_opts()
+    test_interp()