Merge pull request #22 from efiring/spectral

Rework periodogram; add Lomb-Scargle; add documentation; streamline API

doc/conf.py

@@ -24,7 +24,7 @@
 # -- General configuration ------------------------------------------------
 
 # If your documentation needs a minimal Sphinx version, state it here.
-#needs_sphinx = '1.0'
+needs_sphinx = '1.3'
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
@@ -35,6 +35,8 @@
     'sphinx.ext.todo',
     'sphinx.ext.mathjax',
     'sphinx.ext.viewcode',
+    'sphinx.ext.napoleon', # requires sphinx >= 1.3
+    #'numpydoc',           # napoleon handles numpydoc format
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -78,7 +80,7 @@
 
 # The reST default role (used for this markup: `text`) to use for all
 # documents.
-#default_role = None
+default_role = 'py:obj'
 
 # If true, '()' will be appended to :func: etc. cross-reference text.
 #add_function_parentheses = True
@@ -105,7 +107,15 @@
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'pyramid'
+html_theme = 'nature' #
+
+# 'nature' looks like one of the most readable options with numpydoc
+#       It's clean, and the width scales well.
+# 'sphinxdoc' works reasonably well.
+# alabaster is a little too pale, but it has comfortable numpydoc formats
+#'pyramid' is ugly with numpydoc
+# sphinx_rtd_theme also works badly with numpydoc
+# 'agogo' by default is too wide, otherwise so-so
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the

doc/index.rst

@@ -28,6 +28,7 @@ Contents:
    :maxdepth: 2
 
    strategy
+   internals/index.rst
 
 Indices and tables
 ==================

doc/internals/index.rst

@@ -0,0 +1,8 @@
+################
+Internal details
+################
+
+.. toctree::
+   :maxdepth: 1
+
+   periodogram.rst

doc/internals/periodogram.rst

@@ -0,0 +1,13 @@
+***********
+periodogram
+***********
+
+
+:mod:`utide.periodogram`
+========================
+
+.. automodule:: utide.periodogram
+   :members:
+   :private-members:
+   :undoc-members:
+   :show-inheritance:

doc/strategy.rst

@@ -5,7 +5,9 @@ In translating the algorithms from Matlab to Python, the
 first step is generating a set of Python functions that
 mimic their Matlab counterparts.  Gradually, however, the
 code evolves to become closer to what might have been done
-if the original implementation had been in Python.
+if the original implementation had been in Python.  This
+evolution will include extensive renaming of variables and
+functions to improve readability.
 
 Array dimension ordering
 ^^^^^^^^^^^^^^^^^^^^^^^^
@@ -43,9 +45,9 @@ has a leading underscore.
 There is an overwhelming number of options; we might be able
 to find ways of making this interface friendlier.
 
-Options are being held internally in a `dict`.  We might
-switch to using a `Bunch` or similar more flexible
-structure.
+Options are being held internally in a `dict`.  We will
+switch to using a `Bunch` so as to provide both dictionary
+and attribute access syntax.
 
 Time
 ^^^^
@@ -61,10 +63,10 @@ investigate handling whatever Pandas produces.
 
 Missing values
 ^^^^^^^^^^^^^^^
-We will add support for masked array inputs to the public
-functions; the degree to
-which masked arrays will be used internally is open to
-discussion, but most likely their internal use will be at
-most highly localized.
+The `t`, `u`, `v` inputs to `solve` now support any combination
+of nans and masked array inputs to indicate missing values.
+
+The degree to which masked arrays will be used internally is
+unclear, but most likely their use will be highly localized.
 
 

utide/_solve.py

@@ -51,8 +51,11 @@ def solve(tin, uin, vin=None, lat=None, **opts):
 def _solv1(tin, uin, vin, lat, **opts):
 
     print('solve: ')
+
+    # The following returns a possibly modified copy of tin (ndarray).
+    # t, u, v are fully edited ndarrays (unless v is None)
     packed = _slvinit(tin, uin, vin, lat, **opts)
-    t, u, v, tref, lor, elor, opt, tgd, uvgd = packed
+    tin, t, u, v, tref, lor, elor, opt = packed
     nt = len(t)
 
     # opt['cnstit'] = cnstit
@@ -151,7 +154,7 @@ def _solv1(tin, uin, vin, lat, **opts):
             coef['slope'] = np.real(m[-1])/lor
 
     if opt['conf_int'] is True:
-        coef = _confidence(coef, opt, t, e, tin, tgd, uvgd, elor, xraw, xmod,
+        coef = _confidence(coef, opt, t, e, tin, elor, xraw, xmod,
                            W, m, B, nc, Xu, Yu, Xv, Yv)
 
     # diagnostics
@@ -232,21 +235,36 @@ def _slvinit(tin, uin, vin, lat, **opts):
     opt = dict(twodim=(vin is not None))
 
     # Step 1: remove invalid times from tin, uin, vin
-    tgd = ~np.isnan(tin)
-    uin = uin[tgd]
-    tin = tin[tgd]
-    uvgd = ~np.isnan(uin)
+    tin = np.ma.masked_invalid(tin)
+    uin = np.ma.masked_invalid(uin)
     if vin is not None:
-        vin = vin[tgd]
-        uvgd &= ~np.isnan(vin)
+        vin = np.ma.masked_invalid(vin)
+    if np.ma.is_masked(tin):
+        mask = np.ma.getmaskarray(tin)
+        uin = uin.compress(mask)
+        if vin is not None:
+            vin = vin.compress(mask)
+
+    tin = tin.compressed() # no longer masked
 
     # Step 2: generate t, u, v from edited tin, uin, vin
-    t = tin[uvgd]
-    u = uin[uvgd]
     v = None
-    if vin is not None:
-        v = vin[uvgd]
+    if np.ma.is_masked(uin) or np.ma.is_masked(vin):
+        mask = np.ma.getmaskarray(uin)
+        if vin is not None:
+            mask = np.ma.mask_or(np.ma.getmaskarray(vin), mask)
+        t = tin.compress(mask)
+        u = uin.compress(mask).filled()
+        if vin is not None:
+            v = vin.compress(mask).filled()
+    else:
+        t = tin
+        u = uin.filled()
+        if vin is not None:
+            v = vin.filled()
 
+    # Now t, u, v, tin are clean ndarrays; uin and vin are masked,
+    # but don't necessarily have masked values.
 
     # Are the times equally spaced?
     eps = np.finfo(np.float64).eps
@@ -307,6 +325,6 @@ def _slvinit(tin, uin, vin, lat, **opts):
 
     opt['tunconst'] = opt['tunconst'] / opt['tunrdn']
 
-    return t, u, v, tref, lor, elor, opt, tgd, uvgd
+    return tin, t, u, v, tref, lor, elor, opt
 
 

utide/confidence.py

@@ -11,10 +11,10 @@
 import numpy as np
 import scipy.interpolate as sip
 
-from .periodogram import ut_pdgm
+from .periodogram import band_psd
 
 
-def _confidence(coef, opt, t, e, tin, tgd, uvgd, elor, xraw, xmod, W, m, B,
+def _confidence(coef, opt, t, e, tin, elor, xraw, xmod, W, m, B,
                 nc, Xu, Yu, Xv, Yv):
     """
     This confidence interval calculation does not correspond
@@ -33,26 +33,17 @@ def _confidence(coef, opt, t, e, tin, tgd, uvgd, elor, xraw, xmod, W, m, B,
     if not opt['white']:
         # Band-averaged (ba) spectral densities.
         if opt['equi']:
-            if np.sum(tgd) > np.sum(uvgd):
-                # efill = np.interp1(t, e, tin(tgd))
-                # efill = np.interp(t, e, tin[tgd])
-                eTemp = sip.interp1d(t, e)
-                efill = eTemp(tin[tgd])
-                # Fill start&/end nans w/ nearest good.
-                if np.any(np.isnan(efill)):
-                    ind = np.where(np.isnan(efill))[0]
-                    # ind2 = ind(ind<find(~isnan(efill),1,'first'))
-                    ind2 = ind[ind < np.where(~np.isnan(efill), 1, 'first')]
-                    efill[ind2] = efill[np.max(ind2) + 1]
-                    ind2 = ind[ind > np.where(~np.isnan(efill), 1, 'last')]
-                    efill[ind2] = efill[np.min(ind2) - 1]
-
-                ba = ut_pdgm(tin[tgd], efill, coef['aux']['frq'], 1, 0)
+            if len(tin) > len(t):
+                efill = np.interp(tin, t, e)
+                ba = band_psd(tin, efill, coef['aux']['frq'],
+                              equi=True)
             else:
-                ba = ut_pdgm(tin[tgd], e, coef['aux']['frq'], 1, 0)
+                ba = band_psd(tin, e, coef['aux']['frq'],
+                              equi=True)
 
         else:
-            ba = ut_pdgm(t, e, coef['aux']['frq'], 0, opt['lsfrqosmp'])
+            ba = band_psd(t, e, coef['aux']['frq'],
+                          equi=False, frqosamp=opt['lsfrqosmp'])
 
         # import pdb; pdb.set_trace()
         # power [ (e units)^2 ] from spectral density [ (e units)^2 / cph ]