#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Toolbox of various functions and generic utilities.

Authors: Steve Morley, Jon Niehof, Brian Larsen, Josef Koller, Dan Welling
Institution: Los Alamos National Laboratory
Contact: smorley@lanl.gov, jniehof@lanl.gov, balarsen@lanl.gov, jkoller@lanl.gov, dwelling@lanl.gov
Los Alamos National Laboratory

Copyright 2010 Los Alamos National Security, LLC.
"""

#If you add functions here, be sure to:
#1) add to the __all__ list
#2) add to functions in Doc/source/toolbox.rst so it goes in the docs
from __future__ import absolute_import
from __future__ import division

import calendar
import datetime
import glob
try:
    import html.parser
except ImportError: #python2
    import HTMLParser as html
    html.parser = html
try:
    import http.client
except ImportError: # Python2
    import httplib as http
    http.client = http
import numbers
import os
import os.path
import re
try:
    import urllib.request
except ImportError: #python2
    import urllib2 as urllib
    urllib.request = urllib
import select
import shutil
import socket
import subprocess
import sys
import tempfile
import time
import warnings
import zipfile

try:
    import cPickle as pickle
except:
    import pickle

import numpy as np

try:
    from spacepy import help
except ImportError:
    pass
except:
    pass

import spacepy
from spacepy import time as spt

#Try to pull in the C version. Assumption is that if you import this module,
#you want to do some association analysis, so the overhead in the import
#is OK.
from spacepy import lib
if lib.have_libspacepy:
    import ctypes

#Py3k compatibility renamings
try:
    xrange
except NameError:
    xrange = range

__all__ = ['tOverlap', 'tOverlapHalf', 'tCommon', 'loadpickle', 'savepickle', 'assemble',
           'human_sort', 'dictree', 'update', 'progressbar',
           'windowMean', 'medAbsDev', 'binHisto', 'bootHisto',
           'logspace', 'geomspace', 'linspace', 'arraybin', 'mlt2rad',
           'rad2mlt', 'pmm', 'getNamedPath', 'query_yes_no',
           'interpol', 'normalize', 'intsolve', 'dist_to_list',
           'bin_center_to_edges', 'bin_edges_to_center', 'thread_job', 'thread_map',
           'eventTimer', 'isview', 'interweave', 'indsFromXrange', 'hypot',
           'do_with_timeout', 'TimeoutError', 'timeout_check_call',
           'poisson_fit', 'unique_columns']

__contact__ = 'Brian Larsen: balarsen@lanl.gov'

def unique_columns(inval, axis=0):
    """
    Given a multidimensional input return the unique rows or columns along the given axis.
    Based largely on http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array
    axis=0 is unique rows, axis=1 is unique columns

    Parameters
    ==========
    inval :  array-like
        array to find unique columns or rows of

    Optional Parameters
    ===================
    axis : int
        The axis to find unique over, default: 0

    Returns
    =======
    out : array
        N-dimensional array of the unique values along the axis

    Examples
    ========
    """
    # this is a nice trick taking advantage of structed arrays where each row or column
    #   is the value, so returl unique works
    # np.ascontiguousarray() is to be really sure it will work
    if axis == 0:
        val = np.ascontiguousarray(np.transpose(inval))
    else:
        val = np.ascontiguousarray(inval)
    b = val.view(np.dtype((np.void, val.dtype.itemsize * val.shape[1])))
    unique_a = np.unique(b).view(val.dtype).reshape(-1, val.shape[1])
    return unique_a

    
def hypot(*args):
    """
    compute the N-dimensional hypot of an iterable or many arguments

    Parameters
    ==========
    args : many numbers or array-like
        array like or many inputs to compute from

    Returns
    =======
    out : float
        N-dimensional hypot of a number

    Notes
    =====
    This function has a complicated speed function.
     - if a numpy array of floats is input this is passed off to C
     - if iterables are passed in they are made into numpy arrays and comptaton is done local
     - if many scalar agruments are passed in calculation is done in a loop
    For max speed:
     - <20 elements expand them into scalars
         >>> tb.hypot(*vals)
         >>> tb.hypot(vals[0], vals[1]...) #alternate
     - >20 elements premake them into a numpy array of doubles

    Examples
    ========
    >>> from spacepy import toolbox as tb
    >>> print tb.hypot([3,4])
    5.0
    >>> print tb.hypot(3,4)
    5.0
    >>> # Benchmark ####
    >>> from spacepy import toolbox as tb
    >>> import numpy as np
    >>> import timeit
    >>> num_list = []
    >>> num_np = []
    >>> num_np_double = []
    >>> num_scalar = []
    >>> tot = 500
    >>> for num in tb.logspace(1, tot, 10):
    >>>     print num
    >>>     num_list.append(timeit.timeit(stmt='tb.hypot(a)',
                            setup='from spacepy import toolbox as tb;
                            import numpy as np; a = [3]*{0}'.format(int(num)), number=10000))
    >>>     num_np.append(timeit.timeit(stmt='tb.hypot(a)',
                          setup='from spacepy import toolbox as tb;
                          import numpy as np; a = np.asarray([3]*{0})'.format(int(num)), number=10000))
    >>>     num_scalar.append(timeit.timeit(stmt='tb.hypot(*a)',
                              setup='from spacepy import toolbox as tb;
                              import numpy as np; a = [3]*{0}'.format(int(num)), number=10000))
    >>> from pylab import *
    >>> loglog(tb.logspace(1, tot, 10),  num_list, lw=2, label='list')
    >>> loglog(tb.logspace(1, tot, 10),  num_np, lw=2, label='numpy->ctypes')
    >>> loglog(tb.logspace(1, tot, 10),  num_scalar, lw=2, label='scalar')
    >>> legend(shadow=True, fancybox=1, loc='upper left')
    >>> title('Different hypot times for 10000 runs')
    >>> ylabel('Time [s]')
    >>> xlabel('Size')

    .. image:: ../../source/images/hypot_no_extension_speeds_3cases.png
    """
    if lib.have_libspacepy:
        if len(args) == 1 and isinstance(args[0], np.ndarray):  # it is an array
            # make sure everything is C-ready
            ans = lib.hypot_tb(
                np.require(args[0], dtype=np.double, requirements='C'),
                np.product(args[0].shape))
            return ans
    ans = 0.0
    for arg in args:
        if hasattr(arg, '__iter__'):
            tmp = np.asanyarray(arg)
            ans += np.sum(tmp**2)
        else:
            ans += arg**2
    return np.sqrt(ans)

def tOverlap(ts1, ts2, *args, **kwargs):
    """
    Finds the overlapping elements in two lists of datetime objects

    Parameters
    ==========
    ts1 : datetime
        first set of datetime object
    ts2 : datetime
        datatime object
    args :
        additional arguments passed to tOverlapHalf

    Returns
    =======
    out : list
        indices of ts1 within interval of ts2, & vice versa

    Examples
    ========
    Given two series of datetime objects, event_dates and omni['Time']:

    >>> import spacepy.toolbox as tb
    >>> from spacepy import omni
    >>> import datetime
    >>> event_dates = st.tickrange(datetime.datetime(2000, 1, 1), datetime.datetime(2000, 10, 1), deltadays=3)
    >>> onni_dates = st.tickrange(datetime.datetime(2000, 1, 1), datetime.datetime(2000, 10, 1), deltadays=0.5)
    >>> omni = omni.get_omni(onni_dates)
    >>> [einds,oinds] = tb.tOverlap(event_dates, omni['ticks'])
    >>> omni_time = omni['ticks'][oinds[0]:oinds[-1]+1]
    >>> print omni_time
    [datetime.datetime(2000, 1, 1, 0, 0), datetime.datetime(2000, 1, 1, 12, 0),
    ... , datetime.datetime(2000, 9, 30, 0, 0)]

    See Also
    ========
    tOverlapHalf
    tCommon
    """
    idx_1in2 = tOverlapHalf(ts2, ts1, *args, **kwargs)
    idx_2in1 = tOverlapHalf(ts1, ts2, *args, **kwargs)
    if len(idx_2in1) == 0:
        idx_2in1 = None
    if len(idx_1in2) == 0:
        idx_1in2 = None
    return idx_1in2, idx_2in1

def tOverlapHalf(ts1, ts2, presort=False):
    """
    Find overlapping elements in two lists of datetime objects

    This is one-half of tOverlap, i.e. it finds only occurrences where
    ts2 exists within the bounds of ts1, or the second element
    returned by tOverlap.

    Parameters
    ==========
    ts1 : list
        first set of datetime object
    ts2 : list
        datatime object
    presort : bool
        Set to use a faster algorithm which assumes ts1 and
                   ts2 are both sorted in ascending order. This speeds up
                   the overlap comparison by about 50x, so it is worth sorting
                   the list if one sort can be done for many calls to tOverlap

    Returns
    =======
    out : list
        indices of ts2 within interval of ts1

        **note:** Returns empty list if no overlap found

    See Also
    ========
    tOverlap
    tCommon
    """
    if presort:
        import bisect
        t_lower, t_upper = ts1[0], ts1[-1]
        return xrange(bisect.bisect_left(ts2, t_lower),
                     bisect.bisect_right(ts2, t_upper))
    else:
        t_lower, t_upper = min(ts1), max(ts1)
        return [i for i in range(len(ts2))
                if ts2[i] >= t_lower and ts2[i] <= t_upper]

def tCommon(ts1, ts2, mask_only=True):
    """
    Finds the elements in a list of datetime objects present in another

    Parameters
    ==========
    ts1 : list or array-like
        first set of datetime objects
    ts2 : list or array-like
        second set of datetime objects

    Returns
    =======
    out : tuple
        Two element tuple of truth tables (of 1 present in 2, & vice versa)

    See Also
    ========
    tOverlapHalf
    tOverlap

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> import numpy as np
    >>> import datetime as dt
    >>> ts1 = np.array([dt.datetime(2001,3,10)+dt.timedelta(hours=a) for a in range(20)])
    >>> ts2 = np.array([dt.datetime(2001,3,10,2)+dt.timedelta(hours=a*0.5) for a in range(20)])
    >>> common_inds = tb.tCommon(ts1, ts2)
    >>> common_inds[0] #mask of values in ts1 common with ts2
    array([False, False,  True,  True,  True,  True,  True,  True,  True,
            True,  True,  True, False, False, False, False, False, False,
           False, False], dtype=bool)
    >>> ts2[common_inds[1]] #values of ts2 also in ts1

    The latter can be found more simply by setting the mask_only keyword to False

    >>> common_vals = tb.tCommon(ts1, ts2, mask_only=False)
    >>> common_vals[1]
    array([2001-03-10 02:00:00, 2001-03-10 03:00:00, 2001-03-10 04:00:00,
           2001-03-10 05:00:00, 2001-03-10 06:00:00, 2001-03-10 07:00:00,
           2001-03-10 08:00:00, 2001-03-10 09:00:00, 2001-03-10 10:00:00,
           2001-03-10 11:00:00], dtype=object)
    """
    from matplotlib.dates import date2num, num2date

    tn1, tn2 = date2num(ts1), date2num(ts2)

    el1in2 = np.in1d(tn1, tn2, assume_unique=True)  #makes mask of present/absent
    el1in2 = np.in1d(tn1, tn2, assume_unique=True)  #makes mask of present/absent
    el2in1 = np.in1d(tn2, tn1, assume_unique=True)

    if mask_only:
        return el1in2, el2in1
    else:
        truemask1 = np.abs(np.array(el1in2)-1)
        truemask2 = np.abs(np.array(el2in1)-1)
        time1 = np.ma.masked_array(tn1, mask=truemask1)
        time2 = np.ma.masked_array(tn2, mask=truemask2)
        dum1 = num2date(time1.compressed())
        dum2 = num2date(time2.compressed())
        #        dum1 = [val.replace(tzinfo=None) for val in dum1] # this hits ValueError: microsecond must be in 0..999999
        #        dum2 = [val.replace(tzinfo=None) for val in dum2]
        dum11 = []
        dum22 = []
        for v1, v2 in zip(dum1, dum2):
            try:
                dum11.append(v1.replace(tzinfo=None))
            except ValueError: # ValueError: microsecond must be in 0..999999
                dum11.append((datetime.datetime(v1.year, v1.month, v1.day, v1.hour, v1.minute, v1.second, 0) +
                             datetime.timedelta(seconds=1)).replace(tzinfo=None))
            try:
                dum22.append(v2.replace(tzinfo=None))
            except ValueError: # ValueError: microsecond must be in 0..999999
                dum22.append((datetime.datetime(v2.year, v2.month, v2.day, v2.hour, v2.minute, v2.second, 0) +
                             datetime.timedelta(seconds=1)).replace(tzinfo=None))
        dum1 = dum11
        dum2 = dum22
        if type(ts1)==np.ndarray or type(ts2)==np.ndarray:
            dum1 = np.array(dum1)
            dum2 = np.array(dum2)
        return dum1, dum2

def loadpickle(fln):
    """
    load a pickle and return content as dictionary

    Parameters
    ==========
    fln : string
        filename

    Returns
    =======
    out : dict
        dictionary with content from file

    See Also
    ========
    savepickle

    Examples
    ========
    **note**: If fln is not found, but the same filename with '.gz'
           is found, will attempt to open the .gz as a gzipped file.

    >>> d = loadpickle('test.pbin')
    """
    if not os.path.exists(fln) and os.path.exists(fln + '.gz'):
        gzip = True
        fln += '.gz'
    else:
        try:
            with open(fln, 'rb') as fh:
                try: #Py3k
                    return pickle.load(fh, encoding='latin1')
                except TypeError:
                    return pickle.load(fh)
        except pickle.UnpicklingError: #maybe it's a gzip?
            gzip = True
        else:
            gzip = False
    if gzip:
        try:
            import zlib
            with open(fln, 'rb') as fh:
                stream = zlib.decompress(fh.read(), 16 + zlib.MAX_WBITS) 
                try: #Py3k
                    return pickle.loads(stream, encoding='latin1')
                except TypeError:
                    return pickle.loads(stream)
        except MemoryError:
            import gzip
            with open(fln) as fh:
                gzh = gzip.GzipFile(fileobj=fh)
                try: #Py3k
                    contents = pickle.load(gzh, encoding='latin1')
                except TypeError:
                    contents = pickle.load(gzh)
                gzh.close()
            return contents


# -----------------------------------------------
def savepickle(fln, dict, compress=None):
    """
    save dictionary variable dict to a pickle with filename fln

    Parameters
    ----------
    fln : string
        filename
    dict : dict
        container with stuff
    compress : bool
        write as a gzip-compressed file
                     (.gz will be added to ``fln``).
                     If not specified, defaults to uncompressed, unless the
                     compressed file exists and the uncompressed does not.

    See Also
    ========
    loadpickle

    Examples
    ========
    >>> d = {'grade':[1,2,3], 'name':['Mary', 'John', 'Chris']}
    >>> savepickle('test.pbin', d)
    """
    if compress == None: # Guess at compression
        # Assume compressed if compressed already exists (and no uncompressed)
        compress = not os.path.exists(fln) and os.path.exists(fln + '.gz')
    if compress:
        import gzip
        with open(fln + '.gz', 'wb') as fh:
            gzh = gzip.GzipFile(fln, 'wb', compresslevel=3, fileobj=fh)
            pickle.dump(dict, gzh, 2)
            gzh.close()
    else:
        with open(fln, 'wb') as fh:
            pickle.dump(dict, fh, 2) # 2 ... fast binary


# -----------------------------------------------
def assemble(fln_pattern, outfln, sortkey='ticks', verbose=True):
    """
    assembles all pickled files matching fln_pattern into single file and
    save as outfln. Pattern may contain simple shell-style wildcards \*? a la fnmatch
    file will be assembled along time axis given by Ticktock (key: 'ticks') in dictionary
    If sortkey = None, then nothing will be sorted

    Parameters
    ----------
    fln_pattern : string
        pattern to match filenames
    outfln : string
        filename to save combined files to

    Returns
    =======
    out : dict
        dictionary with combined values

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> a, b, c = {'ticks':[1,2,3]}, {'ticks':[4,5,6]}, {'ticks':[7,8,9]}
    >>> tb.savepickle('input_files_2001.pkl', a)
    >>> tb.savepickle('input_files_2002.pkl', b)
    >>> tb.savepickle('input_files_2004.pkl', c)
    >>> a = tb.assemble('input_files_*.pkl', 'combined_input.pkl')
    ('adding ', 'input_files_2001.pkl')
    ('adding ', 'input_files_2002.pkl')
    ('adding ', 'input_files_2004.pkl')
    ('\\n writing: ', 'combined_input.pkl')
    >>> print(a)
    {'ticks': array([1, 2, 3, 4, 5, 6, 7, 8, 9])}
    """
    # done this way so it works before install
    from spacepy import coordinates as c

    filelist = glob.glob(fln_pattern)
    filelist = human_sort(filelist)
    # read all files
    d = {}
    for fln in filelist:
        if verbose: print("adding ", fln)
        d[fln] = loadpickle(fln)

    # combine them
    dcomb = d[filelist[0]]  # copy the first file over
    for fln in filelist[1:]:
       # check if sortkey is actually available
        assert (sortkey in d[fln] or sortkey==None), 'provided sortkey ='+sortkey+' is not available'

        if sortkey:
            TAIcount = len(d[fln][sortkey])
        else:
            TAIcount = len(d[fln][ list(d[fln].keys())[0] ])
        for key in d[fln]:
            #print fln, key
            dim = np.array(np.shape(d[fln][key]))
            ax = np.where(dim==TAIcount)[0]
            if len(ax) == 1: # then match with TAI length is given (jump over otherwise like for 'parameters')
                if isinstance(dcomb[key], spt.Ticktock):
                    dcomb[key] = dcomb[key].append(d[fln][key])
                elif isinstance(dcomb[key], c.Coords):
                    dcomb[key] = dcomb[key].append(d[fln][key])
                else:
                    dcomb[key] = np.append(dcomb[key], d[fln][key], axis=ax[0])

    if sortkey:    #  then sort
        if isinstance(dcomb[sortkey], spt.Ticktock):
            idx = np.argsort(dcomb[sortkey].RDT)
        else:
            idx = np.argsort(dcomb[sortkey])
        TAIcount = len(dcomb[sortkey])
        for key in dcomb: # iterates over keys by default
            dim = np.array(np.shape(dcomb[key]))
            ax = np.where(dim==TAIcount)[0]
            if len(ax) == 1: # then match with length of TAI
                dcomb[key] = dcomb[key][idx] # resort
    else:
        # do nothing
        pass

    if verbose: print('\n writing: ', outfln)
    savepickle(outfln, dcomb)
    return dcomb

def human_sort( l ):
    """ Sort the given list in the way that humans expect.
    http://www.codinghorror.com/blog/2007/12/sorting-for-humans-natural-sort-order.html

    Parameters
    ----------
    l : list
        list of objects to human sort

    Returns
    -------
    out : list
        sorted list

    Examples
    --------
    >>> import spacepy.toolbox as tb
    >>> dat = ['r1.txt', 'r10.txt', 'r2.txt']
    >>> dat.sort()
    >>> print dat
    ['r1.txt', 'r10.txt', 'r2.txt']
    >>> tb.human_sort(dat)
    ['r1.txt', 'r2.txt', 'r10.txt']
    """
    convert = lambda text: int(text) if text.isdigit() else text
    alphanum_key = lambda key: [ convert(c) for c in re.split(r'([0-9]+)', key) ]
    try:
        l.sort( key=alphanum_key )
    except TypeError:
        l.sort()
    return l


def dictree(in_dict, verbose=False, spaces=None, levels=True, attrs=False, **kwargs):
    """
    pretty print a dictionary tree

    Parameters
    ----------
    in_dict : dict
        a complex dictionary (with substructures)
    verbose : boolean (optional)
        print more info
    spaces : string (optional)
        string will added for every line
    levels : integer (optional)
        number of levels to recurse through (True means all)
    attrs : boolean (optional)
        display information for attributes

    Examples
    --------
    >>> import spacepy.toolbox as tb
    >>> d = {'grade':{'level1':[4,5,6], 'level2':[2,3,4]}, 'name':['Mary', 'John', 'Chris']}
    >>> tb.dictree(d)
    +
    |____grade
         |____level1
         |____level2
    |____name

    More complicated example using a datamodel:

    >>> from spacepy import datamodel
    >>> counts = datamodel.dmarray([2,4,6], attrs={'units': 'cts/s'})
    >>> data = {'counts': counts, 'PI': 'Dr Zog'}
    >>> tb.dictree(data)
    +
    |____PI
    |____counts
    >>> tb.dictree(data, attrs=True, verbose=True)
    +
    |____PI (str [6])
    |____counts (spacepy.datamodel.dmarray (3,))
        :|____units (str [5])

    Attributes of, e.g., a CDF or a datamodel type object (obj.attrs)
    are denoted by a colon.
    """
    try:
        assert hasattr(in_dict, 'keys')
    except AssertionError:
        try:
            assert hasattr(in_dict, 'attrs')
        except:
            raise TypeError('dictree: Input must be dictionary-like')

    if not spaces:
        spaces = ''
        print('+')

    if 'toplev' in kwargs:
        toplev = kwargs['toplev']
    else:
        toplev = True
    try:
        if toplev and attrs:
            dictree(in_dict.attrs, spaces = ':', verbose = verbose, levels = levels, attrs=attrs, toplev=True)
            toplev = False
    except:
        pass

    # TODO, if levels is True why check again?
    if levels:
        try:
            assert levels is True
        except AssertionError:
            levels -= 1
            if levels == 0:
                levels = None

    try:
        for key in sorted(in_dict.keys()):
            bar = '|____' + str(key)
            if verbose:
                typestr = str(type(in_dict[key])).split("'")[1]
                #check entry for dict-like OR .attrs dict
                try:
                    dimstr = in_dict[key].shape
                    dimstr = ' ' + str(dimstr)
                except AttributeError:
                    try:
                        dimstr = len(in_dict[key])
                        dimstr = ' [' + str(dimstr) + ']'
                    except:
                        dimstr = ''
                print(spaces + bar + ' ('+ typestr + dimstr + ')')
            else:
                print(spaces + bar)
            if hasattr(in_dict[key], 'attrs') and attrs:
                dictree(in_dict[key].attrs, spaces = spaces + '    :', verbose = verbose, levels = levels, attrs=attrs, toplev=False)
            if hasattr(in_dict[key], 'keys') and levels:
                dictree(in_dict[key], spaces = spaces + '     ', verbose = verbose, levels = levels, attrs=attrs, toplev=False)
    except:
        pass
    return None


def _crawl_yearly(base_url, pattern, datadir, name=None, cached=True,
                  startyear=None):
    """Crawl files in a directory-by-year structure

    Parameters
    ==========
    base_url : str
        Base of the data. This URL should point to a directory containing
        yearly directories (YYYY).
    pattern : str
        Regular expression to match filenames. Will download files in each
        yearly directory that match the pattern.
    datadir : str
        Directory to store downloaded files. A mirror will be maintained in
        this directory. Note this is a "flat" mirror without the year
        directories.
    name : str (optional)
        Name of the data set, used only in status messages.
    cached : boolean (optional)
        Only update files if timestamp on server is newer than
        timestamp on local file (default). Set False to always
        download files.
    startyear : int (optional)
        First year to crawl, as four-digit number or four-character string.
        If not specified, will download all years. If specified, will delete
        years prior to this which have already been downloaded.

    Returns
    =======
    list
        All the filenames that were mirrored, in order; or None if there
        were no updates. If there are any updates, all filenames are
        included.
    """
    name = '' if name is None else '{} '.format(name)
    #Find all the files to download
    print("Finding {}files to download ...".format(name))
    progressbar(0, 1, 1, text='Listing files')
    conn = None
    if spacepy.config['keepalive']:
        try:
            data, conn = get_url(base_url, keepalive=True)
        except socket.error: #Give up on keepalives
            pass
    if conn is None:
        data = get_url(base_url)
    if str is not bytes:
        data = data.decode('utf-8')
    p = LinkExtracter()
    p.feed(data)
    p.close()
    yearlist = [y[0:4] for y in p.links if re.match(r'\d{4}/', y)
                and (startyear is None or y[0:4] >= str(startyear))]
    downloadme = {}
    for i, y in enumerate(yearlist):
        yearurl = '{}{}/'.format(base_url, y)
        if conn is None:
            data = get_url(yearurl)
        else:
            data, conn = get_url(yearurl, keepalive=True, conn=conn)
        if str is not bytes:
            data = data.decode('utf-8')
        p = LinkExtracter()
        p.feed(data)
        p.close()
        for f in p.links:
            if not re.match(pattern, f):
                continue
            downloadme[f] = yearurl + f
        progressbar(i + 1, 1, len(yearlist), text='Listing files')
    print("Retrieving {}files ...".format(name))
    filenames = sorted(list(downloadme.keys()))
    if not os.path.exists(datadir):
        os.makedirs(datadir)
    newdata = False
    #Check for existing files (delete them if no longer on server)
    have_files = os.listdir(datadir)
    for f in have_files:
        if not f in filenames:
            os.remove(os.path.join(datadir, f))
            #File was removed, so need to reparse even if no new downloads
            newdata = True
    #Download
    for i, fname in enumerate(filenames):
        if conn is None:
            res = get_url(downloadme[fname], os.path.join(datadir, fname),
                          cached=cached)
        else:
            res, conn = get_url(downloadme[fname], os.path.join(datadir, fname),
                                cached=cached, keepalive=True, conn=conn)
        if res is not None:
            newdata = True
        progressbar(i + 1, 1, len(filenames))
    if conn is not None:
        conn.close()
    return filenames if newdata else None


def _get_qindenton_daily(qd_daily_url=None, cached=True, startyear=None):
    """Download the Qin-Denton OMNI-like daily files
    
    Parameters
    ==========
    qd_daily_url : str (optional)
        Base of the Qin-Denton data, in hourly JSON-headed ASCII. This URL
        should point to the directory containing the yearly directories.
        Default from ``qd_daily_url`` in config file.
    cached : boolean (optional)
        Only update files if timestamp on server is newer than
        timestamp on local file (default). Set False to always
        download files.
    startyear : int (optional)
        If specified, start downloading files from the year given,
        rather than all years. This will delete older files!

    Returns
    =======
    SpaceData
        The data extracted from the Q-D dataset, fully processed for saving
        as SpacePy HDF5 OMNI data.
    """
    if qd_daily_url is None:
        qd_daily_url = spacepy.config['qd_daily_url']
    datadir = os.path.join(spacepy.DOT_FLN, 'data', 'qindenton_daily_files')
    _crawl_yearly(qd_daily_url, r'QinDenton_\d{8}_hour.txt',
                  datadir, name='Q-D daily', cached=cached, startyear=startyear)
    #Read and process
    print("Processing Q-D daily files ...")
    return _assemble_qindenton_daily(datadir)


def _assemble_qindenton_daily(qd_daily_dir):
    """Assemble Qin-Denton daily files into OMNI structure

    Parameters
    ==========
    qd_daily_dir : str
        Directory with Qin-Denton daily files.

    Returns
    =======
    SpaceData
        The data extracted from the Q-D dataset, fully processed for saving
        as SpacePy HDF5 OMNI data.
    """
    import spacepy.datamodel
    filelist = sorted(glob.glob(os.path.join(qd_daily_dir, '*_hour.txt')))
    data = [spacepy.datamodel.readJSONheadedASCII(f)
            for f in filelist]
    omnidata = spacepy.datamodel.SpaceData()
    for k in data[0].keys():
        if k in ('DateTime', 'Minute', 'OriginFile', 'Second'):
            continue
        omnidata[k] = spacepy.datamodel.dmarray(
            np.concatenate([d[k] for d in data]), dtype=np.float32)
    del data
    ntimes = set([len(v) for v in omnidata.values()])
    if len(ntimes) != 1:
        raise ValueError(
            'Input Q-D daily file has different size for different variables')
    ntimes = ntimes.pop()
    # Renaming from the names in new file to old file names
    for oldname, newname in (('dens', 'Den_P'),
                             ('velo', 'Vsw')):
        omnidata[oldname] = omnidata[newname]
        del omnidata[newname]
        if '{}_status'.format(newname) in omnidata:
            omnidata['{}_status'.format(oldname)] \
                = omnidata['{}_status'.format(newname)]
            del omnidata['{}_status'.format(newname)]
    # Quality flags
    qbits = spacepy.datamodel.SpaceData()
    for k in list(omnidata.keys()): # Edit while iterate
        if not k.endswith('_status'):
            continue
        v = spacepy.datamodel.dmarray(omnidata[k], dtype=np.int8)
        basename = k.split('_')[0]
        if basename in ('G', 'W'): # Arrays
            for i in range(omnidata[k].shape[1]):
                qbits['{}{:d}'.format(basename, i + 1)] = v[:, i]
        else:
            qbits[basename] = v
        del omnidata[k]
    omnidata['Qbits'] = qbits
    # Reformat some arrays to multi-variable
    for name in ('Bz', 'G', 'W'):
        for i in range(omnidata[name].shape[1]):
            omnidata['{}{:d}'.format(name, i + 1)] = omnidata[name][:, i]
        del omnidata[name]
    # Make integers of integers
    for k in ('Dst', 'Year', 'Month', 'Day', 'Hour'):
        omnidata[k] = spacepy.datamodel.dmarray(omnidata[k], dtype=np.int16)
    # Process time formats
    omnidata['UTC'] = spacepy.datamodel.dmarray([
        datetime.datetime(omnidata['Year'][i],
                          omnidata['Month'][i],
                          omnidata['Day'][i],
                          omnidata['Hour'][i])
        for i in range(len(omnidata['Year']))])
    omnidata['DOY'] = spacepy.datamodel.dmarray([
        dt.timetuple().tm_yday for dt in omnidata['UTC']], dtype=np.int16)
    omnidata['RDT'] = spt.Ticktock(omnidata['UTC'], 'UTC').RDT
    for k in ('Year', 'Hour', 'Month', 'Day'):
        del omnidata[k]
    return omnidata


def _get_cdaweb_omni2(omni2url=None):
    """Download the OMNI2 data from SPDF

    Parameters
    ==========
    omni2url : str (optional)
        Base of the OMNI2 data at SPDF, in hourly CDF form. This URL
        should point to the directory containing the yearly directories.
        Default from ``omni2_url`` in config file.

    Returns
    =======
    SpaceData
        The data extracted from the OMNI2 dataset, with variables renamed
        to match the old ViRBO combined OMNI2 CDF. Returns ``None``
        if there are no new data.
    """
    import spacepy.pycdf
    import spacepy.pycdf.istp
    if omni2url is None:
        omni2url = spacepy.config['omni2_url']
    datadir = os.path.join(spacepy.DOT_FLN, 'data', 'omni2cdfs')
    filenames = _crawl_yearly(omni2url, r'omni2_h0_mrg1hr_\d{8}_v\d+\.cdf',
                              datadir, name='OMNI2')
    if filenames is None:
        return None
    #Read and process
    print("Reading OMNI2 files ...")
    data = spacepy.pycdf.concatCDF([spacepy.pycdf.CDF(os.path.join(datadir, f))
                                    for f in filenames])
    #Map keyed by the variable names as in the OMNI2 data from CDAWeb,
    #valued by the variable names as in the OMNI2 data from ViRBO
    #i.e. this is from: to
    #None means not in VirBO, so delete.
    keymap = {
        'ABS_B': 'Field_mag_ave', #This is average of the magnitude
        'AE': 'AE_index',
        'AL_INDEX': 'AL_index',
        'AP_INDEX': 'Ap_index',
        'AU_INDEX': 'AU_index',
        'BX_GSE': 'Bx_GSE',
        'BY_GSE': 'By_GSE',
        'BZ_GSE': 'Bz_GSE',
        'BY_GSM': 'By_GSM',
        'BZ_GSM': 'Bz_GSM',
        'Beta': 'Plasma_beta',
        'Day': 'Decimal_Day',
        'DST': 'Dst_index',
        'E': 'Electric_field_GSM',
        #'Epoch': 'Epoch', #same, do not change
        'F': 'Mag_of_ave_field_vect', #Magnitude of average
        'F10_INDEX': 'f10_7_index',
        'HR': 'Hour',
        'IMF': 'IMF_spacecraft_ID',
        'IMF_PTS': 'Num_pts_in_IMF_aves',
        'KP': 'Kp_index',
        'Mach_num': 'Alfven_mach_number',
        'MFLX': 'Proton_flux_flag',
        'Mgs_mach_num': None, #magnetosonic mach number
        'N': 'Ion_density',
        'PC_N_INDEX': 'PC_N_index',
        'PHI-V': 'Plasma_bulk_flow_long_angle',
        'PHI_AV': 'Long_angle_ave_field_vector',
        'PLS': 'SW_plasma_spacecraft_ID',
        'PLS_PTS': 'Num_pts_in_plasma_aves',
        'PR-FLX_1': 'Proton_flux_gt_1_MeV',
        'PR-FLX_10': 'Proton_flux_gt_10_MeV',
        'PR-FLX_2': 'Proton_flux_gt_2_MeV',
        'PR-FLX_30': 'Proton_flux_gt_30_MeV',
        'PR-FLX_4': 'Proton_flux_gt_4_MeV',
        'PR-FLX_60': 'Proton_flux_gt_60_MeV',
        'Pressure': 'Flow_pressure',
        'R': 'Sunspot_number',
        'Ratio': 'Na_over_Np',
        'Rot#': 'Bartels_rotation_num',
        'SIGMA-ABS_B': 'sigma_field_mag_ave',
        'SIGMA-B': 'sigma_mag_of_ave_field_vect',
        'SIGMA-Bx': 'sigma_Bx_GSE',
        'SIGMA-By': 'sigma_By_GSE',
        'SIGMA-Bz': 'sigma_Bz_GSE',
        'SIGMA-N': 'sigma_ion_density',
        'SIGMA-PHI-V': 'sigma_plasma_flow_long_angle',
        'SIGMA-T': 'sigma_temp',
        'SIGMA-THETA-V': 'sigma_plasma_flow_lat_angle',
        'SIGMA-V': 'sigma_plasma_bulk_speed',
        'SIGMA-ratio': 'sigma_ratio',
        'Solar_Lyman_alpha': None,
        'T': 'Plasma_temp',
        'THETA-V': 'Plasma_bulk_flow_lat_angle',
        'THETA_AV': 'Lat_angle_ave_field_vector',
        'V': 'Plasma_bulk_speed',
        'YR': 'Year',
        }
    for k, v in keymap.items():
        if v is not None:
            data[v] = data[k]
        del data[k]
    #The CDAWeb Epochs are erroneously tagged as START of collection,
    #not midpoint. So fix that
    data['Epoch'] = data['Epoch'] + datetime.timedelta(minutes=30)
    #Castings to match ViRBO: everything is an int, and fill turns to NaN
    for k, v in data.items():
        if v.dtype == object: #Skip epoch
            continue
        if v.dtype == np.int32:
            data[k] = v.astype(np.float32)
            v = data[k]
        spacepy.pycdf.istp.nanfill(v)
    return data


class LinkExtracter(html.parser.HTMLParser):
    """Finds all links in a HTML page, useful for crawling.

    After HTML has been parsed, the ``links`` attribute contains
    a list of link targets.
    """

    def reset(self, *args, **kwargs):
        if isinstance(LinkExtracter, type):
            super(LinkExtracter, self).reset(*args, **kwargs)
        else: #py2k, old-style class
            html.parser.HTMLParser.reset(self, *args, **kwargs)
        self.links = []
        """List of link targets found in the page"""

    def handle_starttag(self, tag, attrs):
        if tag != 'a':
            return
        for name, value in attrs:
            if name != 'href':
                continue
            self.links.append(value)


def get_url(url, outfile=None, reporthook=None, cached=False,
            keepalive=False, conn=None):
    """Read data from a URL

    Open an HTTP URL, honoring the user agent as specified in the
    SpacePy config file. Returns the data, optionally also writing
    out to a file.

    This is similar to the deprecated ``urlretrieve``.

    Parameters
    ==========
    url : str
        The URL to open
    outfile : str (optional)
        Full path to file to write data to
    reporthook : callable (optional)
        Function for reporting progress; takes arguments of block
        count, block size, and total size.
    cached : bool (optional)
        Compare modification time of the URL to the modification time
        of ``outfile``; do not retrieve (and return None) unless
        the URL is newer than the file.
    keepalive : bool (optional)
        Attempt to keep the connection open to retrieve more URLs.
        The return becomes a tuple of (data, conn) to return the
        connection used so it can be used again. This mode does not
        support proxies. (Default False)
    conn : http.client.HTTPConnection (optional)
        An established http connection (HTTPS is also okay) to use with
        ``keepalive``. If not provided, will attempt to make a connection.

    Returns
    =======
    bytes
        The HTTP data from the server.

    See Also
    ========
    progressbar

    Notes
    =====
    This function honors proxy settings as described in
    :func:`urllib.request.getproxies`. Cryptic error messages (such as
    ``Network is unreachable``) may indicate that proxy settings
    should be defined as appropriate for your environment (e.g. with
    ``HTTP_PROXY`` or ``HTTPS_PROXY`` environment variables).
    """
    if not keepalive and conn is not None:
        raise ValueError('Cannot specify connection without keepalive')
    if keepalive:
        scheme, _, host, path = url.split('/', 3)
        path = '/' + path # Explicitly root on the server
        if conn is not None and conn.sock is not None:
            readable, writeable, _ = select.select(
                [conn.sock], [conn.sock], [], 0)
            # Make sure no stale data to read on socket, and can write to it
            if readable or len(writeable) != 1:
                conn.close()
                conn = None
        if conn is None:
            ctype = http.client.HTTPConnection if scheme == 'http:'\
                    else http.client.HTTPSConnection
            conn = ctype(host)
        clheaders = {
            "Connection": "keep-alive",
        }
        if spacepy.config.get('user_agent', ''):
            clheaders['User-Agent'] =  spacepy.config['user_agent']
        conn.request('HEAD' if cached else 'GET', path, headers=clheaders)
        def checkresponse(conn):
            """Get the response on a connection, return response and headers"""
            r = conn.getresponse()
            if r.status >= 400:
                raise RuntimeError(
                    'HTTP status {} {}'.format(r.status, r.reason))
            headers = dict(((k.title(), v) for k, v in r.getheaders()))
            return r, headers
        r, headers = checkresponse(conn)
    else:
        r = urllib.request.Request(url)
        if spacepy.config.get('user_agent', ''):
            r.add_header('User-Agent', spacepy.config['user_agent'])
        r = urllib.request.urlopen(r)
        if r.getcode() >= 400:
            r.close()
            raise RuntimeError('HTTP status {} {}'.format(r.code, r.msg))
        headers = r.info()
    modified = headers.get('Last-Modified', None)
    if modified is not None:
        # strptime is affected by locale (including the month name) but the
        # header is a constant format, so massage
        modified = modified.split()[1:5] # Get rid of day of week and 'GMT'
        modified[1] = str(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
                           'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']\
                          .index(modified[1]) + 1) # Make month numerical
        modified = datetime.datetime.strptime(
            ' '.join(modified), "%d %m %Y %H:%M:%S")
        modified = calendar.timegm(modified.timetuple())
    if cached:
        if outfile is None:
            if not keepalive:
                r.close()
            raise RuntimeError('Must specify outfile if cached is True')
        if os.path.exists(outfile) and modified is not None:
            #Timestamp is truncated to second, so do same for local
            local_mod = int(os.path.getmtime(outfile))
            if modified <= local_mod:
                if keepalive:
                    r.read()
                else:
                    r.close()
                return (None, conn) if keepalive else None
        if keepalive: # Replace previous header request with full get
            r.read()
            conn.request('GET', path, headers=clheaders)
            r, headers = checkresponse(conn)
    size = int(headers.get('Content-Length', 0))
    blocksize = 1024
    count = 0
    data = []
    while True:
        newdata = r.read(blocksize)
        if not newdata:
            break
        data.append(newdata)
        count += 1
        if reporthook:
            reporthook(count, blocksize, size)
    if not keepalive:
        r.close()
    if outfile:
        with open(outfile, 'wb') as f:
            for d in data:
                f.write(d)
        if modified is not None: #Copy web mtime to file
            os.utime(outfile, (int(time.time()), modified))
    data = b''.join(data)
    return (data, conn) if keepalive else data


def update(all=True, QDomni=False, omni=False, omni2=False, leapsecs=False,
           PSDdata=False, cached=True):
    """
    Download and update local database for omni, leapsecs etc

    Web access is via :func:`get_url`; notes there may be helpful in
    debugging errors. See also the ``keepalive`` configuration option.

    Parameters
    ==========
    all : boolean (optional)
        if True, update OMNI2, Qin-Denton and leapsecs
    omni : boolean (optional)
        if True. update only omni (Qin-Denton)
    omni2 : boolean (optional)
        if True, update only original OMNI2
    QDomni : boolean (optional)
        if True, update OMNI2 and Qin-Denton
    leapsecs : boolean (optional)
        if True, update only leapseconds
    cached : boolean (optional)
        Only update files if timestamp on server is newer than
        timestamp on local file (default). Set False to always
        download files.

    Returns
    =======
    out : string
        data directory where things are saved

    See Also
    ========
    get_url

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.update(omni=True)
    """
    from spacepy.datamodel import SpaceData, dmarray, fromCDF, toHDF5
    from spacepy import DOT_FLN, config

    datadir = os.path.join(DOT_FLN, 'data')
    if not os.path.exists(datadir):
        os.mkdir(datadir)
        os.chmod(datadir, 0o777)

    #leapsec_url ='ftp://maia.usno.navy.mil/ser7/tai-utc.dat'
    leapsec_fname = os.path.join(datadir, 'tai-utc.dat')

    # define location for getting omni
    #omni_url = 'ftp://virbo.org/QinDenton/hour/merged/latest/WGhour-latest.d.zip'
    omni_fname_zip = os.path.join(datadir, 'WGhour-latest.d.zip')
    omni2_fname_zip = os.path.join(datadir, 'omni2-latest.cdf.zip')
    omni_fname_pkl = os.path.join(datadir, 'omnidata.pkl')
    omni_fname_json = os.path.join(datadir, 'omnidata.txt')
    omni_fname_h5 = os.path.join(datadir, 'omnidata.h5')
    omni2_fname_h5 = os.path.join(datadir, 'omni2data.h5')

    PSDdata_fname = os.path.join('psd_dat.sqlite')

    if (omni or omni2 or QDomni or leapsecs or PSDdata):
        all = False #if an option is explicitly selected, turn 'all' off

    if all == True:
        omni = True
        omni2 = True
        leapsecs = True

    if QDomni == True:
        omni = True
        omni2 = True

    if omni == True:
        # retrieve omni, unzip and save as table
        print("Retrieving initial Qin-Denton file ...")
        get_url(config['qindenton_url'], omni_fname_zip, progressbar,
                cached=cached)
        fh_zip = zipfile.ZipFile(omni_fname_zip)
        data = fh_zip.read(fh_zip.namelist()[0])
        fh_zip.close()
        if not str is bytes:
            data = data.decode('ascii')
        A = np.array(data.split('\n'))
        print("Processing initial Qin-Denton file ...")

        # create a keylist
        keys = A[0].split()
        keys.remove('8')
        keys.remove('6')
        keys[keys.index('status')] = '8_status'
        keys[keys.index('stat')] = '6_status'
        keys[keys.index('dst')] = 'Dst'
        keys[keys.index('kp')] = 'Kp'
        #keys[keys.index('Hr')] = 'Hr'
        keys[keys.index('V_SW')] = 'velo'
        keys[keys.index('Den_P')] = 'dens'
        keys[keys.index('Day')] = 'DOY'
        keys[keys.index('Year')] = 'Year'

        # remove keyword lines and empty lines as well
        idx = np.where(A != '')[0]
        # put it into a 2D table
        tab = [val.split() for val in A[idx[1:]]]
        stat8 = [val[11] for val in tab]
        stat6 = [val[27] for val in tab]

        tab = np.array(tab, dtype='float32')
        # take out where Dst not available ( = 99999) or year == 0
        idx = np.where((tab[:,12] !=99.0) & (tab[:,0] != 0))[0]
        tab = tab[idx,:]
        stat8 = np.array(stat8)[idx]
        stat6 = np.array(stat6)[idx]

        omnidata = SpaceData()
        # sort through and make an omni dictionary
        # extract keys from line above
        for ikey, i  in zip(keys,range(len(keys))):
            if ikey in ('Year', 'DOY', 'Hr', 'Dst'):
                omnidata[ikey] = dmarray(tab[:, i], dtype='int16')
            else:
                omnidata[ikey] = dmarray(tab[:,i])

        # add TAI to omnidata
        nTAI = len(omnidata['DOY'])

        # add interpolation quality flags
        omnidata['Qbits'] = SpaceData()
        arr = dmarray(stat8.view(stat8.dtype.kind + '1'),
                      dtype=np.byte).reshape((8, nTAI))
        for ik, key in enumerate(['ByIMF', 'BzIMF', 'velo', 'dens', 'Pdyn', 'G1', 'G2', 'G3']):
            omnidata['Qbits'][key] = arr[ik,:]
        if stat6.dtype.str[1:] == 'U6':
            stat6 = np.require(stat6, dtype='|S6')
        arr = dmarray(stat6.view(stat6.dtype.kind + '1'),
                      dtype=np.byte).reshape((6, nTAI))
        for ik, key in enumerate(['W1', 'W2', 'W3', 'W4', 'W5', 'W6']):
            omnidata['Qbits'][key] = arr[ik,:]

        #remove string status keys
        foo = omnidata.pop('6_status')
        foo = omnidata.pop('8_status')

        # add time information to omni pickle (long loop)
        omnidata['UTC'] = dmarray([datetime.datetime(int(omnidata['Year'][i]), 1, 1) +
                 datetime.timedelta(days=int(omnidata['DOY'][i]) - 1,
                                    hours=int(omnidata['Hr'][i]))
                 for i in range(nTAI)])

        omnidata['ticks'] = spt.Ticktock(omnidata['UTC'], 'UTC')
        omnidata['RDT'] = omnidata['ticks'].RDT
        del omnidata['ticks'] #Can be quickly regenerated on import
        startyear = omnidata['Year'][-1]
        del omnidata['Year']
        del omnidata['Hr']

        # Supplement with daily files
        print('Supplementing with latest Q-D daily files,'
              ' this will take a while...')
        dailyomnidata = _get_qindenton_daily(cached=cached, startyear=startyear)
        # Find where new files start
        idx = np.searchsorted(omnidata['UTC'], dailyomnidata['UTC'][0])
        for k in sorted(omnidata.keys()):
            if k == 'Qbits':
                for qk in sorted(omnidata[k].keys()):
                    omnidata[k][qk] = spacepy.dmarray(np.concatenate((
                        omnidata[k][qk][:idx, ...], dailyomnidata[k][qk])))
            else:
                omnidata[k] = spacepy.dmarray(np.concatenate((
                    omnidata[k][:idx, ...], dailyomnidata[k])))

        print("Now saving... ")
        ##for now, make one file -- think about whether monthly/annual files makes sense
        toHDF5(omni_fname_h5, omnidata)
        print('Complete.')

    if omni2 == True:
        omni2_url = config['omni2_url']
        if omni2_url.endswith('.zip'):
            # adding missing values from original omni2
            print("Retrieving OMNI2 file ...")
            get_url(omni2_url, omni2_fname_zip, progressbar, cached=cached)
            fh_zip = zipfile.ZipFile(omni2_fname_zip)
            flist = fh_zip.namelist()
            if len(flist) != 1:
                raise RuntimeError('Unable to find OMNI2 file in zip file.')
            file_to_read = flist[0]
            if os.path.dirname(file_to_read):
                raise RuntimeError('Unexpected contents of OMNI2 zip file.')
            td = tempfile.mkdtemp()
            try:
                fh_zip.extract(file_to_read, td)
                omnicdf = fromCDF(os.path.join(td, file_to_read))
            finally:
                fh_zip.close()
                os.remove(omni2_fname_zip)
                shutil.rmtree(td)
        else:
            omnicdf = _get_cdaweb_omni2(omni2_url)

        if omnicdf is not None: #If no new data, skip all this
            #add RDT
            omnicdf['RDT'] = spt.Ticktock(omnicdf['Epoch'],'UTC').RDT
            #remove keys that get in the way
            del omnicdf['Hour']
            del omnicdf['Year']
            del omnicdf['Decimal_Day']

            # save as HDF5
            toHDF5(omni2_fname_h5, omnicdf)

    if leapsecs == True:
        print("Retrieving leapseconds file ... ")
        get_url(config['leapsec_url'], leapsec_fname, progressbar,
                cached=cached)
        # Reload leap seconds if they've already been used.
        if 'spacepy.time' in sys.modules\
           and hasattr(sys.modules['spacepy.time'], 'TAIleaps'):
            sys.modules['spacepy.time']._read_leaps()

    if PSDdata == True:
        print("Retrieving PSD sql database")
        get_url(config['psddata_url'], PSDdata_fname, progressbar,
                cached=cached)
    return datadir

def indsFromXrange(inxrange):
    '''return the start and end indices implied by an xrange, useful when xrange is zero-length

    Parameters
    ==========
    inxrange : xrange
        input xrange object to parse

    Returns
    =======
    list of int
       List of start, stop indices in the xrange. The return value is not
       defined if a stride is specified or if stop is before start (but
       will work when stop equals start).

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> foo = xrange(23, 39)
    >>> foo[0]
    23
    >>> tb.indsFromXrange(foo)
    [23, 39]
    >>> foo1 = xrange(23, 23)
    >>> tb.indsFromXrange(foo) #indexing won't work in this case
    [23, 23]
    '''
    if not isinstance(inxrange, xrange): return None
    valstr = inxrange.__str__()
    if ',' not in valstr:
        res = re.search(r'(\d+)', valstr)
        retval = [int(0), int(res.group(1))]
    else:
        res = re.search(r'(\d+), (\d+)', valstr)
        retval = [int(res.group(1)), int(res.group(2))]
    return retval

def progressbar(count, blocksize, totalsize, text='Download Progress'):
    """
    print a progress bar with urllib.urlretrieve reporthook functionality

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> import urllib
    >>> urllib.urlretrieve(config['psddata_url'], PSDdata_fname, reporthook=tb.progressbar)
    """
    percent = count * blocksize * 100. / totalsize
    if percent > 100:
        percent = 100.
    sys.stdout.write("\r{} ...{:.0f}%".format(text, percent))
    if percent >= 100:
        sys.stdout.write("\n")
    sys.stdout.flush()

def windowMean(data, time=[], winsize=0, overlap=0, st_time=None, op=np.mean):
    """
    Windowing mean function, window overlap is user defined

    Parameters
    ==========
    data : array_like
        1D series of points
    time : list (optional)
        series of timestamps, optional (format as numeric or datetime)
        For non-overlapping windows set overlap to zero.
    winsize : integer or datetime.timedelta (optional)
        window size
    overlap : integer or datetime.timedelta (optional)
        amount of window overlap
    st_time : datetime.datetime (optional)
        for time-based averaging, a start-time other than the first
        point can be specified
    op : callable (optional)
        the operator to be called, default numpy.mean

    Returns
    =======
    out : tuple
        the windowed mean of the data, and an associated reference time vector

    Examples
    ========
    For non-overlapping windows set overlap to zero.
    e.g. (time-based averaging)
    Given a data set of 100 points at hourly resolution (with the time tick in
    the middle of the sample), the daily average of this, with half-overlapping
    windows is calculated:

    >>> import spacepy.toolbox as tb
    >>> from datetime import datetime, timedelta
    >>> wsize = datetime.timedelta(days=1)
    >>> olap = datetime.timedelta(hours=12)
    >>> data = [10, 20]*50
    >>> time = [datetime.datetime(2001,1,1) + datetime.timedelta(hours=n, minutes = 30) for n in range(100)]
    >>> outdata, outtime = tb.windowMean(data, time, winsize=wsize, overlap=olap, st_time=datetime.datetime(2001,1,1))
    >>> outdata, outtime
    ([15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0],
     [datetime.datetime(2001, 1, 1, 12, 0),
      datetime.datetime(2001, 1, 2, 0, 0),
      datetime.datetime(2001, 1, 2, 12, 0),
      datetime.datetime(2001, 1, 3, 0, 0),
      datetime.datetime(2001, 1, 3, 12, 0),
      datetime.datetime(2001, 1, 4, 0, 0),
      datetime.datetime(2001, 1, 4, 12, 0)])

    When using time-based averaging, ensure that the time tick corresponds to
    the middle of the time-bin to which the data apply. That is, if the data
    are hourly, say for 00:00-01:00, then the time applied should be 00:30.
    If this is not done, unexpected behaviour can result.

    e.g. (pointwise averaging),

    >>> outdata, outtime = tb.windowMean(data, winsize=24, overlap=12)
    >>> outdata, outtime
    ([15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0], [12.0, 24.0, 36.0, 48.0, 60.0, 72.0, 84.0])

    where winsize and overlap are numeric,
    in this example the window size is 24 points (as the data are hourly) and
    the overlap is 12 points (a half day). The output vectors start at
    winsize/2 and end at N-(winsize/2), the output time vector
    is basically a reference to the nth point in the original series.

    **note** This is a quick and dirty function - it is NOT optimized, at all.
    """
    #check inputs and initialize
    #Set resolution to 1 if no times supplied
    if len(time) == 0:
        startpt, res = 0, 1
        time = list(range(len(data)))
        pts = True
    else:
        if len(data) != len(time):
            raise ValueError('windowmean error: data and time must have same length')
        #First check if datetime objects
        try:
            assert type(winsize) == datetime.timedelta
            assert type(overlap) == datetime.timedelta
        except AssertionError:
            raise TypeError('windowmean error: winsize/overlap must be timedeltas')
        pts = False #force time-based averaging
        if (type(time[0]) != datetime.datetime):
            startpt = time[0]

    #now actually do windowing mean
    outdata, outtime = [], []
    data = np.array(data)
    if pts:
        #loop for fixed number of points in window
        try:
            inttypes = (int, long)
        except NameError:
            inttypes = (int,)
        if not isinstance(winsize, inttypes):
            winsize = int(round(winsize))
            warnings.warn('windowmean: non-integer windowsize, rounding to %d' \
            % winsize)
        if winsize < 1:
            winsize = 1
            warnings.warn('windowmean: window length < 1, defaulting to 1')
        if overlap >= winsize:
            overlap = winsize - 1
            warnings.warn('''windowmean: overlap longer than window, truncated to
            %d''' % overlap)
        lastpt = winsize-1 #set last point to end of window size
        while lastpt < len(data):
            datwin = np.ma.masked_where(np.isnan(data[startpt:startpt+winsize]), \
                data[startpt:startpt+winsize])
            getmean = op(datwin.compressed()) #mean of window, excl. NaNs
            gettime = (time[startpt+winsize] - time[startpt])/2. \
                + time[startpt]#new timestamp
            startpt = startpt+winsize-overlap
            lastpt = startpt+winsize
            outdata.append(getmean) #construct output arrays
            outtime.append(gettime)
    else:
        #loop with time-based window
        lastpt = time[0] + winsize
        delta = datetime.timedelta(microseconds=1) #TODO: replace this with an explicit check for times on the boundary?
        if st_time:
            startpt = st_time
        else:
            startpt = time[0]
        if overlap >= winsize:
            raise ValueError('Overlap requested greater than size of window')
        while startpt < time[-1]:
            getinds = tOverlapHalf([startpt,startpt+winsize-delta], time, presort=True)
            if getinds: #if not None
                getdata = np.ma.masked_where(np.isnan(data[getinds[0]:getinds[-1]+1]), data[getinds[0]:getinds[-1]+1])
                getmean = op(getdata.compressed()) #find mean excluding NaNs
            else:
                getmean = np.nan
            gettime = startpt + winsize//2 #new timestamp -floordiv req'd with future division
            startpt = startpt + winsize - overlap #advance window start
            lastpt = startpt + winsize
            outdata.append(getmean) #construct output arrays
            outtime.append(gettime)

    return outdata, outtime

def medAbsDev(series, scale=False):
    """
    Calculate median absolute deviation of a given input series

    Median absolute deviation (MAD) is a robust and resistant measure of
    the spread of a sample (same purpose as standard deviation). The
    MAD is preferred to the inter-quartile range as the inter-quartile
    range only shows 50% of the data whereas the MAD uses all data but
    remains robust and resistant. See e.g. Wilks, Statistical methods
    for the Atmospheric Sciences, 1995, Ch. 3. For additional details on 
    the scaling, see Rousseeuw and Croux, J. Amer. Stat. Assoc., 88 (424),
    pp. 1273-1283, 1993.

    Parameters
    ==========
    series : array_like
        the input data series

    Other Parameters
    ================
    scale : bool
        if True (default: False), scale to standard deviation of a normal distribution

    Returns
    =======
    out : float
        the median absolute deviation

    Examples
    ========
    Find the median absolute deviation of a data set. Here we use the log-
    normal distribution fitted to the population of sawtooth intervals, see
    Morley and Henderson, Comment, Geophysical Research Letters, 2009.

    >>> import numpy
    >>> import spacepy.toolbox as tb
    >>> numpy.random.seed(8675301)
    >>> data = numpy.random.lognormal(mean=5.1458, sigma=0.302313, size=30)
    >>> print data
    array([ 181.28078923,  131.18152745, ... , 141.15455416, 160.88972791])
    >>> tb.medAbsDev(data)
    28.346646721370192

    **note** This implementation is robust to presence of NaNs
    """
    #ensure input is numpy array (and make 1-D)
    series = (np.array(series, dtype=float)).ravel()
    #mask for NaNs
    series = np.ma.masked_where(np.isnan(series),series)
    #get median absolute deviation of unmasked elements
    perc50 = np.median(series.compressed())
    mad = np.median(abs(series.compressed()-perc50))
    if scale:
        mad *= 1.4826 #scale so that MAD is same as SD for normal distr.
    return mad

def binHisto(data, verbose=False):
    """
    Calculates bin width and number of bins for histogram using Freedman-Diaconis rule, if rule fails, defaults to square-root method

    The Freedman-Diaconis method is detailed in:
        Freedman, D., and P. Diaconis (1981), On the histogram as a density estimator: L2 theory, Z. Wahrscheinlichkeitstheor. Verw. Geb., 57, 453–476

    and is also described by:
        Wilks, D. S. (2006), Statistical Methods in the Atmospheric Sciences, 2nd ed.

    Parameters
    ==========
    data : array_like
        list/array of data values
    verbose : boolean (optional)
        print out some more information

    Returns
    =======
    out : tuple
        calculated width of bins using F-D rule, number of bins (nearest integer) to use for histogram

    Examples
    ========
    >>> import numpy, spacepy
    >>> import matplotlib.pyplot as plt
    >>> numpy.random.seed(8675301)
    >>> data = numpy.random.randn(1000)
    >>> binw, nbins = spacepy.toolbox.binHisto(data)
    >>> print(nbins)
    19
    >>> p = plt.hist(data, bins=nbins, histtype='step', density=True)

    See Also
    ========
    matplotlib.pyplot.hist
    """
    pul = np.percentile(data, (25, 75)) #get confidence interval
    ql, qu = pul[0], pul[1]
    iqr = qu-ql
    binw = 2.*iqr/(len(data)**(1./3.))
    if binw != 0:
        nbins = int(round((max(data)-min(data))/binw))
    # if nbins is 0, NaN or inf don't use the F-D rule just use sqrt(num) rule
    if binw == 0 or nbins == 0 or not np.isfinite(nbins) or nbins >= len(data)/2.:
        nbins = int(round(np.sqrt(len(data))))
        binw = len(data)/nbins
        if verbose:
            print("Used sqrt rule")
    else:
        if verbose:
            print("Used F-D rule")
    return (binw, nbins)

def bootHisto(data, inter=90., n=1000, seed=None,
              plot=False, target=None, figsize=None, loc=None, **kwargs):
    """Bootstrap confidence intervals for a histogram.

    All other keyword arguments are passed to :func:`numpy.histogram`
    or :func:`matplotlib.pyplot.bar`.

    .. versionchanged:: 0.2.3
       This argument pass-through did not work in earlier versions of
       SpacePy.

    Parameters
    ==========

    data : array_like
        list/array of data values
    inter : float (optional; default 90)
        percentage confidence interval to return. Default 90% (i.e.
        lower CI will be 5% and upper will be 95%)
    n : int (optional; default 1000)
        number of bootstrap iterations
    seed : int (optional)
        Optional seed for the random number generator. If not
        specified; numpy generator will not be reseeded.
    plot : bool (optional)
        Plot the result. Plots if True or ``target``, ``figsize``,
        or ``loc`` specified.
    target : (optional)
        Target on which to plot the figure (figure or axes). See
        :func:`spacepy.plot.utils.set_target` for details.
    figsize : tuple (optional)
        Passed to :func:`spacepy.plot.utils.set_target`.
    loc : int (optional)
        Passed to :func:`spacepy.plot.utils.set_target`.

    Returns
    =======
    out : tuple
      tuple of bin_edges, low, high, sample[, bars]. Where
      ``bin_edges`` is the edges of the bins used; ``low`` is the
      histogram with the value for each bin from the bottom of that
      bin's confidence interval; ``high`` similarly for the top;
      ``sample`` is the histogram of the input sample without
      resampling. If plotting, also returned is ``bars``, the
      container object returned from matplotlib.

    Notes
    =====
    .. versionadded:: 0.2.1

    The confidence intervals are calculated for each bin individually and thus
    the resulting low/high histograms may not have actually occurred in the
    calculation from the surrogates. If using a probability density histogram,
    this can have "interesting" implications for interpretation.

    Examples
    ========
    .. plot::
        :include-source:

        >>> import numpy.random
        >>> import spacepy.toolbox
        >>> numpy.random.seed(0)
        >>> data = numpy.random.randn(1000)
        >>> bin_edges, low, high, sample, bars = spacepy.toolbox.bootHisto(
        ...     data, plot=True)

    See Also
    ========
    binHisto
    plot.utils.set_target
    numpy.histogram
    matplotlib.pyplot.hist
    """
    import spacepy.poppy
    histogram_allowed_kwargs = (
        'bins', 'range', 'normed', 'weights', 'density')
    histogram_kwargs = {k: v for k, v in kwargs.items()
                        if k in histogram_allowed_kwargs}
    bar_kwargs = {k: v for k, v in kwargs.items()
                  if k not in histogram_allowed_kwargs}
    sample, bin_edges = np.histogram(data, **histogram_kwargs)
    histogram_kwargs['bins'] = bin_edges
    ci_low, ci_high = spacepy.poppy.boots_ci(
        data, n, inter,
        lambda x: np.histogram(x, **histogram_kwargs)[0], seed=seed,
        nretvals=len(bin_edges) - 1)
    if not plot and all([x is None for x in (target, figsize, loc)]):
        return bin_edges, ci_low, ci_high, sample
    import spacepy.plot.utils
    _, ax = spacepy.plot.utils.set_target(
        target, figsize=figsize, loc=(111 if loc is None else loc))
    if 'ecolor' not in bar_kwargs:
        bar_kwargs['ecolor'] = 'k'
    bars = ax.bar(
        bin_edges[:-1], height=sample, width=np.diff(bin_edges),
        align='edge', yerr=np.stack((sample - ci_low, ci_high - sample)),
        **bar_kwargs)
    return bin_edges, ci_low, ci_high, sample, bars


def logspace(min, max, num, **kwargs):
    """
    Returns log-spaced bins. Same as numpy.logspace except the min and max are the min and max
    not log10(min) and log10(max)

    Parameters
    ==========
    min : float
        minimum value
    max : float
        maximum value
    num : integer
        number of log spaced bins

    Other Parameters
    ================
    kwargs : dict
        additional keywords passed into matplotlib.dates.num2date

    Returns
    =======
    out : array
        log-spaced bins from min to max in a numpy array

    Notes
    =====
    This function works on both numbers and datetime objects

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.logspace(1, 100, 5)
    array([   1.        ,    3.16227766,   10.        ,   31.6227766 ,  100.        ])

    See Also
    ========
    geomspace
    linspace
    """
    if isinstance(min, datetime.datetime):
        from matplotlib.dates import date2num, num2date
        ans = num2date(np.logspace(np.log10(date2num(min)), np.log10(date2num(max)), num, **kwargs))
        ans = spt.no_tzinfo(ans)
        return np.array(ans)
    else:
        return np.logspace(np.log10(min), np.log10(max), num, **kwargs)

def linspace(min, max, num, **kwargs):
    """
    Returns linear-spaced bins. Same as numpy.linspace except works with datetime
    and is faster

    Parameters
    ==========
    min : float, datetime
        minimum value
    max : float, datetime
        maximum value
    num : integer
        number of linear spaced bins

    Other Parameters
    ================
    kwargs : dict
        additional keywords passed into matplotlib.dates.num2date

    Returns
    =======
    out : array
        linear-spaced bins from min to max in a numpy array

    Notes
    =====
    This function works on both numbers and datetime objects

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.linspace(1, 10, 4)
    array([  1.,   4.,   7.,  10.])

    See Also
    ========
    geomspace
    logspace
    """
    if hasattr(min, 'shape') and min.shape == ():
        min = min.item()
    if hasattr(max, 'shape') and max.shape == ():
        max = max.item()
    if isinstance(min, datetime.datetime):
        from matplotlib.dates import date2num, num2date
        ans = num2date(np.linspace(date2num(min), date2num(max), num, **kwargs))
        ans = spt.no_tzinfo(ans)
        return np.array(ans)
    else:
        return np.linspace(min, max, num, **kwargs)

def geomspace(start, ratio=None, stop=False, num=50):
    """
    Returns geometrically spaced numbers.

    Parameters
    ==========
    start : float
        The starting value of the sequence.
    ratio : float (optional)
        The ratio between subsequent points
    stop : float (optional)
        End value, if this is selected `num` is overridden
    num : int (optional)
        Number of samples to generate. Default is 50.

    Returns
    =======
    seq : array
        geometrically spaced sequence

    See Also
    ========
    linspace
    logspace

    Examples
    ========
    To get a geometric progression between 0.01 and 3 in 10 steps

    >>> import spacepy.toolbox as tb
    >>> tb.geomspace(0.01, stop=3, num=10)
    [0.01,
     0.018846716378431192,
     0.035519871824902655,
     0.066943295008216955,
     0.12616612944575134,
     0.23778172582285118,
     0.44814047465571644,
     0.84459764235318191,
     1.5917892219322083,
     2.9999999999999996]

    To get a geometric progression with a specified ratio, say 10

    >>> import spacepy.toolbox as tb
    >>> tb.geomspace(0.01, ratio=10, num=5)
    [0.01, 0.10000000000000001, 1.0, 10.0, 100.0]
    """
    if not ratio and stop != False:
        ratio = (stop/start)**(1/(num-1))
    seq = []
    seq.append(start)
    if stop == False:
        for j in range(1, num):
            seq.append(seq[j-1]*ratio)
        return seq
    else:
        val, j = start, 1
        while val <= stop or np.allclose(val, stop, ):
            val = seq[j-1]*ratio
            seq.append(val)
            j+=1
        return seq[:-1]

def arraybin(array, bins):
    """
    Split a sequence into subsequences based on value.

    Given a sequence of values and a sequence of values representing the
    division between bins, return the indices grouped by bin.

    Parameters
    ==========
    array : array_like
        the input sequence to slice, must be sorted in ascending order
    bins : array_like
        dividing lines between bins. Number of bins is len(bins)+1,
            value that exactly equal a dividing value are assigned
            to the higher bin

    Returns
    =======
    out : list
        indices for each bin (list of lists)

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.arraybin(range(10), [4.2])
    [[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]]
    """
    bin_it = lambda value: (i for i in range(len(array)) if array[i] >= value)
    splits = [next(bin_it(value), len(array)) for value in bins]
    return [list(range(start_idx, stop_idx)) for (start_idx, stop_idx)
            in zip([0] + splits, splits + [len(array)])]

def mlt2rad(mlt, midnight = False):
    """
    Convert mlt values to radians for polar plotting
    transform mlt angles to radians from -pi to pi
    referenced from noon by default

    Parameters
    ==========
    mlt : numpy array
        array of mlt values
    midnight : boolean (optional)
        reference to midnight instead of noon

    Returns
    =======
    out : numpy array
        array of radians

    Examples
    ========
    >>> from numpy import array
    >>> mlt2rad(array([3,6,9,14,22]))
    array([-2.35619449, -1.57079633, -0.78539816,  0.52359878,  2.61799388])

    See Also
    ========
    rad2mlt
    """
    if midnight:
        try:
            mlt_arr =  [val + 12 for val in mlt]
        except TypeError:
            mlt_arr = mlt + 12
    else:
        mlt_arr = mlt
    try:
        rad_arr = [(val-12)*np.pi/12. for val in mlt_arr]
    except TypeError:
        rad_arr = (mlt_arr-12)*np.pi/12.
    return rad_arr

def rad2mlt(rad, midnight=False):
    """
    Convert radians values to mlt
    transform radians from -pi to pi to mlt
    referenced from noon by default

    Parameters
    ==========
    rad : numpy array
        array of radian values
    midnight : boolean (optional)
        reference to midnight instead of noon

    Returns
    =======
    out : numpy array
        array of mlt values

    Examples
    ========
    >>> rad2mlt(array([0,pi, pi/2.]))
    array([ 12.,  24.,  18.])

    See Also
    ========
    mlt2rad
    """
    if midnight:
        rad_arr = rad + np.pi
    else:
        rad_arr = rad
    mlt_arr=rad_arr*(12/np.pi) + 12
    return mlt_arr

def pmm(*args):
    """
    print min and max of input arrays

    Parameters
    ==========
    a : array-like
        arbitrary number of input arrays (or lists)

    Returns
    =======
    out : list
        list of min, max for each array

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> from numpy import arange
    >>> tb.pmm(arange(10), arange(10)+3)
    [[0, 9], [3, 12]]
    """
    ans = []
    for a in args:
        try:
            ind = np.isfinite(a)
        except TypeError:
            ind = np.arange(len(a)).astype(int)
            import warnings
            warnings.warn('pmm: Unable to exclude non-finite values, results may be incorrect', RuntimeWarning)
        try:
            ans.append([np.min(a[ind]), np.max(a[ind])])
        except TypeError:
            a_tmp = np.asarray(a)
            if a_tmp.dtype.type in [np.dtype('S').type, np.dtype('U').type]:
                a_tmp = np.require(a_tmp, dtype=object)
            ans.append([np.min(a_tmp[ind]), np.max(a_tmp[ind])])
    return ans

def getNamedPath(name):
    """
    Return the full path of a parent directory with name as the leaf

    Parameters
    ==========
    name : string
        the name of the parent directory to locate

    Examples
    ========
    Run from a directory
    /mnt/projects/dream/bin/Ephem
    with 'dream' as the name, this function
    would return '/mnt/projects/dream'
    """
    def findNamed(path):
        pp = os.path.split(path)
        if pp[-1] == '':
            return None
        if pp[-1] != name:
            path = findNamed(pp[0])
        return path
    return findNamed(os.getcwd())

def query_yes_no(question, default="yes"):
    """
    Ask a yes/no question via raw_input() and return their answer.

    "question" is a string that is presented to the user.
    "default" is the presumed answer if the user just hits <Enter>.
    It must be "yes" (the default), "no" or None (meaning
    an answer is required of the user).

    The "answer" return value is one of "yes" or "no".

    Parameters
    ==========
    question : string
        the question to ask
    default : string (optional)

    Returns
    =======
    out : string
        answer ('yes' or 'no')

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.query_yes_no('Ready to go?')
    Ready to go? [Y/n] y
    'yes'
    """
    valid = {"yes":"yes",   "y":"yes",  "ye":"yes",
             "no":"no",     "n":"no"}
    if default == None:
        prompt = " [y/n] "
    elif default == "yes":
        prompt = " [Y/n] "
    elif default == "no":
        prompt = " [y/N] "
    else:
        raise ValueError("invalid default answer: '%s'" % default)
    while 1:
        sys.stdout.write(question + prompt)
        if sys.version_info[0]==2:
            choice = raw_input().lower()
        elif sys.version_info[0]>2:
            choice = input().lower()
        if default is not None and choice == '':
            return default
        elif choice in valid:
            return valid[choice]
        else:
            sys.stdout.write("Please respond with 'yes' or 'no' "\
                             "(or 'y' or 'n').\n")

def interpol(newx, x, y, wrap=None, **kwargs):
    """
    1-D linear interpolation with interpolation of hours/longitude

    Parameters
    ==========
    newx : array_like
        x values where we want the interpolated values
    x : array_like
        x values of the original data (must be monotonically increasing or wrapping)
    y : array_like
        y values of the original data
    wrap : string, optional
        for continuous x data that wraps in y at 'hours' (24), 'longitude' (360),
        or arbitrary value (int, float)
    kwargs : dict
        additional keywords, currently accepts baddata that sets baddata for
        masked arrays

    Returns
    =======
    out : numpy.masked_array
        interpolated data values for new abscissa values

    Examples
    ========
    For a simple interpolation

    >>> import spacepy.toolbox as tb
    >>> import numpy
    >>> x = numpy.arange(10)
    >>> y = numpy.arange(10)
    >>> tb.interpol(numpy.arange(5)+0.5, x, y)
    array([ 0.5,  1.5,  2.5,  3.5,  4.5])

    To use the wrap functionality, without the wrap keyword you get the wrong answer

    >>> y = range(24)*2
    >>> x = range(len(y))
    >>> tb.interpol([1.5, 10.5, 23.5], x, y, wrap='hour').compressed() # compress removed the masked array
    array([  1.5,  10.5,  23.5])
    >>> tb.interpol([1.5, 10.5, 23.5], x, y)
    array([  1.5,  10.5,  11.5])
    """
    if 'baddata' in kwargs:
        y = np.ma.masked_equal(y, kwargs['baddata'])
        x = np.ma.masked_array(x)
        x.mask = y.mask
        kwargs.__delitem__('baddata')
    else:
        tst = np.ma.core.MaskedArray
        if type(x)!=tst or type(y)!=tst or type(newx)!=tst:
            x = np.ma.masked_array(x)
            y = np.ma.masked_array(y)
            newx = np.ma.masked_array(newx)

    def wrap_interp(xout, xin, yin, sect):
        dpsect=360/sect
        yc = np.cos(np.deg2rad(y*dpsect))
        ys = np.sin(np.deg2rad(y*dpsect))
        new_yc = np.interp(newx, x.compressed(), yc.compressed(), **kwargs)
        new_ys = np.interp(newx, x.compressed(), ys.compressed(), **kwargs)
        try:
            new_bad = np.interp(newx, x, y.mask)
        except ValueError:
            new_bad = np.zeros((len(newx)))
        newy = np.rad2deg(np.arctan(new_ys/new_yc))/dpsect
        #1st quadrant is O.K
        #2nd quadrant
        idx = [n for n in range(len(new_yc)) if new_yc[n]<0 and new_ys[n]>0]
        newy[idx] = sect/2 + newy[idx]
        #print('Sector 2 inds: %s' % idx)
        #3rd quadrant
        idx = [n for n in range(len(new_yc)) if new_yc[n]<0 and new_ys[n]<0]
        newy[idx] = sect/2 + newy[idx]
        #print('Sector 3 inds: %s' % idx)
        #4th quadrant
        idx = [n for n in range(len(new_yc)) if new_yc[n]>0 and new_ys[n]<0]
        newy[idx] = sect + newy[idx]
        #print('Sector 4 inds: %s' % idx)
        new_bad = np.ma.make_mask(new_bad)
        newy = np.ma.masked_array(newy, mask=new_bad)
        return newy

    if wrap=='hour':
        newy = wrap_interp(newx, x.compressed(), y.compressed(), 24)
    elif wrap=='lon':
        newy = wrap_interp(newx, x.compressed(), y.compressed(), 360)
    elif isinstance(wrap, numbers.Real):
        newy = wrap_interp(newx, x.compressed(), y.compressed(), wrap)
    else:
        newy = np.interp(newx, x.compressed(), y.compressed(), **kwargs)
    return newy

# -----------------------------------------------

def normalize(vec, low=0.0, high=1.0):
    """
    Given an input vector normalize the vector to a given range

    Parameters
    ==========
    vec : array_like
        input vector to normalize
    low : float
        minimum value to scale to, default 0.0
    high : float
        maximum value to scale to, default 1.0

    Returns
    =======
    out : array_like
        normalized vector

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.normalize([1,2,3])
    [0.0, 0.5, 1.0]
    """
    return np.interp(vec, (np.nanmin(vec), np.nanmax(vec)), (low, high))


def intsolve(func, value, start=None, stop=None, maxit=1000):
    """
    Find the function input such that definite integral is desired value.

    Given a function, integrate from an (optional) start point until the
    integral reached a desired value, and return the end point of the
    integration.

    Parameters
    ==========
    func : callable
        function to integrate, must take single parameter
    value : float
        desired final value of the integral
    start : float (optional)
        value at which to start integration, default -Infinity
    stop : float (optional)
        value at which to stop integration, default +Infinity
    maxit : integer
        maximum number of iterations

    Returns
    =======
    out : float
        x such that the integral of L{func} from L{start} to x is L{value}

    **Note:** Assumes func is everywhere positive, otherwise solution may
           be multi-valued.
    """
    from scipy import inf
    from scipy.integrate import quad
    from warnings import warn
    if start is None:
        start = -inf
    if stop is None:
        stop = inf
    lower_bound = start
    upper_bound = stop
    it = 0
    while it < maxit:
        it += 1
        if upper_bound == inf:
            if lower_bound == -inf:
                test_bound = 0
            elif lower_bound < 1.0:
                test_bound = 1.0
            else:
                test_bound = lower_bound * 2
        else:
            if lower_bound == -inf:
                if upper_bound > -1.0:
                    test_bound = -1.0
                else:
                    test_bound = upper_bound * 2
            else:
                test_bound = (lower_bound + upper_bound) / 2.0
        (test_value, err) = quad(func, start, test_bound)
        if abs(value - test_value) <= err:
            break
        elif value < test_value:
            upper_bound = test_bound
        else:
            lower_bound = test_bound

    if abs(value - test_value) > err:
        warn('Difference between desired value and actual is ' +
             str(abs(value - test_value)) +
             ', greater than integral error ' +
             str(err), UserWarning, stacklevel=2)
    return test_bound

def dist_to_list(func, length, min=None, max=None):
    """
    Convert a probability distribution function to a list of values

    This is a deterministic way to produce a known-length list of values
    matching a certain probability distribution. It is likely to be a closer
    match to the distribution function than a random sampling from the
    distribution.

    Parameters
    ==========
    func : callable
        function to call for each possible value, returning
            probability density at that value (does not need to be
            normalized.)
    length : int
        number of elements to return
    min : float
        minimum value to possibly include
    max : float
        maximum value to possibly include

    Examples
    ========
    >>> import matplotlib
    >>> import numpy
    >>> import spacepy.toolbox as tb
    >>> gauss = lambda x: math.exp(-(x ** 2) / (2 * 5 ** 2)) / (5 * math.sqrt(2 * math.pi))
    >>> vals = tb.dist_to_list(gauss, 1000, -numpy.inf, numpy.inf)
    >>> print vals[0]
    -16.45263...
    >>> p1 = matplotlib.pyplot.hist(vals, bins=[i - 10 for i in range(21)], facecolor='green')
    >>> matplotlib.pyplot.hold(True)
    >>> x = [i / 100.0 - 10.0 for i in range(2001)]
    >>> p2 = matplotlib.pyplot.plot(x, [gauss(i) * 1000 for i in x], 'red')
    >>> matplotlib.pyplot.draw()
    """
    from scipy import inf
    from scipy.integrate import quad
    if min is None:
        min = -inf
    if max is None:
        max = inf
    total = quad(func, min, max)[0]
    step = float(total) / length
    return [intsolve(func, (0.5 + i) * step, min, max) for i in range(length)]

def bin_center_to_edges(centers):
    """
    Convert a list of bin centers to their edges

    Given a list of center values for a set of bins, finds the start and
    end value for each bin. (start of bin n+1 is assumed to be end of bin n).
    Useful for e.g. matplotlib.pyplot.pcolor.

    Edge between bins n and n+1 is arithmetic mean of the center of n
    and n+1; edge below bin 0 and above last bin are established to make
    these bins symmetric about their center value.

    Parameters
    ==========
    centers : list
        list of center values for bins

    Returns
    =======
    out : list
        list of edges for bins

    **note:** returned list will be one element longer than centers

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.bin_center_to_edges([1,2,3])
    [0.5, 1.5, 2.5, 3.5]
    """
    edges = bin_edges_to_center(centers)
    edges = np.append(centers[0]-(edges[0]-centers[0]), edges)
    edges = np.append(edges, centers[-1]+(centers[-1]-edges[-1]))
    return edges

def bin_edges_to_center(edges):
    """
    Convert a list of bin edges to their centers

    Given a list of edge values for a set of bins, finds the center of each bin.
    (start of bin n+1 is assumed to be end of bin n).

    Center of bin n is arithmetic mean of the edges of the adjacent bins.

    Parameters
    ==========
    edges : list
        list of edge values for bins

    Returns
    =======
    out : numpy.ndarray
        array of centers for bins

    **note:** returned array will be one element shorter than edges

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.bin_center_to_edges([1,2,3])
    [0.5, 1.5, 2.5, 3.5]
    """
    df = np.diff(edges)
    if isinstance(df[0], datetime.timedelta) and sys.version_info[0:2]<=(3,2):
        return edges[:-1] + df//2
    else:
        return edges[:-1] + df/2

def thread_job(job_size, thread_count, target, *args, **kwargs):
    """
    Split a job into subjobs and run a thread for each

    Each thread spawned will call L{target} to handle a slice of the job.

    This is only useful if a job:
      1. Can be split into completely independent subjobs
      2. Relies heavily on code that does not use the Python GIL, e.g.
         numpy or ctypes code
      3. Does not return a value. Either pass in a list/array to hold the
         result, or see L{thread_map}

    Examples
    ========
    squaring 100 million numbers:

    >>> import numpy
    >>> import spacepy.toolbox as tb
    >>> numpy.random.seed(8675301)
    >>> a = numpy.random.randint(0, 100, [100000000])
    >>> b = numpy.empty([100000000], dtype='int64')
    >>> def targ(in_array, out_array, start, count): \
             out_array[start:start + count] = in_array[start:start + count] ** 2
    >>> tb.thread_job(len(a), 0, targ, a, b)
    >>> print(b[0:5])
    [2704 7225  196 1521   36]

    This example:
      - Defines a target function, which will be called for each thread.
        It is usually necessary to define a simple "wrapper" function
        like this to provide the correct call signature.
      - The target function receives inputs C{in_array} and C{out_array},
        which are not touched directly by C{thread_job} but are passed
        through in the call. In this case, C{a} gets passed as
        C{in_array} and C{b} as C{out_array}
      - The target function also receives the start and number of elements
        it needs to process. For each thread where the target is called,
        these numbers are different.

    Parameters
    ==========
    job_size : int
        Total size of the job. Often this is an array size.
    thread_count : int
        Number of threads to spawn. If =0 or None, will
            spawn as many threads as there are cores available on
            the system. (Each hyperthreading core counts as 2.)
            Generally this is the Right Thing to do.
            If NEGATIVE, will spawn abs(thread_count) threads,
            but will run them sequentially rather than in
            parallel; useful for debugging.
    target : callable
        Python callable (generally a function, may also be an
            imported ctypes function) to run in each thread. The
            *last* two positional arguments passed in will be a
            "start" and a "subjob size," respectively;
            frequently this will be the start index and the number
            of elements to process in an array.
    args : sequence
        Arguments to pass to L{target}. If L{target} is an instance
            method, self must be explicitly passed in. start and
            subjob_size will be appended.
    kwargs : dict
        keyword arguments to pass to L{target}.
    """
    try:
        import threading
    except:
        return (target(*(args +  (0, job_size)), **kwargs), )
    if thread_count == None or thread_count == 0:
        try:
            import multiprocessing
            thread_count = multiprocessing.cpu_count()
        except: #Do something not too stupid
            thread_count = 8
    if thread_count < 0:
        thread_count *= -1
        seq = True
    else:
        seq = False
    count = float(job_size) / thread_count
    starts = [int(count * i + 0.5) for i in range(thread_count)]
    subsize = [(starts[i + 1] if i < thread_count - 1 else job_size) -
               starts[i]
               for i in range(thread_count)]
    threads = []
    for i in range(thread_count):
        t = threading.Thread(
            None, target, None,
            args + (starts[i], subsize[i]), kwargs)
        t.start()
        if seq:
            t.join()
        else:
            threads.append(t)
    if not seq:
        for t in threads:
            t.join()

def thread_map(target, iterable, thread_count=None, *args, **kwargs):
    """
    Apply a function to every element of a list, in separate threads

    Interface is similar to multiprocessing.map, except it runs in threads

    This is made largely obsolete in python3 by from concurrent import futures 
    
    Examples
    ========
    find totals of several arrays

    >>> import numpy
    >>> from spacepy import toolbox
    >>> inputs = range(100)
    >>> totals = toolbox.thread_map(numpy.sum, inputs)
    >>> print(totals[0], totals[50], totals[99])
    (0, 50, 99)

    >>> # in python3
    >>> from concurrent import futures
    >>> with futures.ThreadPoolExecutor(max_workers=4) as executor:
    ...:     for ans in executor.map(numpy.sum, [0,50,99]):
    ...:         print ans
    #0
    #50
    #99
   
    Parameters
    ==========
    target : callable
        Python callable to run on each element of iterable.
            For each call, an element of iterable is appended to
            args and both args and kwargs are passed through.
            Note that this means the iterable element is always the
            *last* positional argument; this allows the specification
            of self as the first argument for method calls.
    iterable : iterable
        elements to pass to each call of L{target}
    args : sequence
        arguments to pass to target before each element of
            iterable
    thread_count : integer
        Number of threads to spawn; see L{thread_job}.
    kwargs : dict
        keyword arguments to pass to L{target}.

    Returns
    =======
    out : list
        return values of L{target} for each item from L{iterable}
    """
    try:
        jobsize = len(iterable)
    except TypeError:
        iterable = list(iterable)
        jobsize = len(iterable)
    def array_targ(function, it, retvals, arglist, kwarglist, start, size):
        for i in range(start, start + size):
            retvals[i] = function(*(arglist + (it[i],)), **kwarglist)
    retvals = [None] * jobsize
    thread_job(jobsize, thread_count, array_targ,
               target, iterable, retvals, args, kwargs)
    return retvals

def eventTimer(Event, Time1):
    """
    Times an event then prints out the time and the name of the event,
    nice for debugging and seeing that the code is progressing

    Parameters
    ==========
    Event : str
        Name of the event, string is printed out by function
    Time1 : time.time
        the time to difference in the function

    Returns
    =======
    Time2 : time.time
        the new time for the next call to EventTimer

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> import time
    >>> t1 = time.time()
    >>> t1 = tb.eventTimer('Test event finished', t1)
    ('4.40', 'Test event finished')
    """
    Time2 = time.time()
    print(("%4.2f" % (Time2 - Time1), Event))
    return Time2


def isview(array1, array2=None):
    """
    Returns if an object is a view of another object.  More precisely if one array argument is specified
    True is returned is the arrays owns its data.  If two arrays arguments are specified a tuple is returned
    of if the first array owns its data and the the second if they point at the same memory location

    Parameters
    ==========
    array1 : numpy.ndarray
        array to query if it owns its data

    Other Parameters
    ================
    array2 : object (optional)
        array to query if array1 is a view of this object at the specified memory location

    Returns
    =======
    out : bool or tuple
        If one array is specified bool is returned, True is the array owns its data.  If two arrays
        are specified a tuple where the second element is a bool of if the array point at the same
        memory location

    Examples
    ========
    import numpy
    import spacepy.toolbox as tb
    a = numpy.arange(100)
    b = a[0:10]
    tb.isview(a)
    # False
    tb.isview(b)
    # True
    tb.isview(b, a)
    # (True, True)
    tb.isview(b, b)
    # (True, True)  # the conditions are met and numpy cannot tell this
    """
    # deal with the one input case first
    if array2 is None:
        try:
            if array1.base is None:
                return False
            return True
        except AttributeError:
            return False # if it is not an array then it is not a view
    # there are two arrays input
    try:
        if array1.base is None:
            return (False, False)
        return (True, array1.base is array2)
    except AttributeError:
            return (False, False) # if it is not an array then it is not a view

def interweave(a, b):
    """
    given two array-like variables interweave them together.
    Discussed here: http://stackoverflow.com/questions/5347065/interweaving-two-numpy-arrays

    Parameters
    ==========
    a : array-like
        first array

    b : array-like
        second array

    Returns
    =======
    out : numpy.ndarray
        interweaved array
    """
    a = np.asanyarray(a)
    b = np.asanyarray(b)
    ans = np.empty((a.size + b.size), dtype=a.dtype)
    ans[0::2] = a
    ans[1::2] = b
    return ans


class TimeoutError(Exception):
    """Raised when a time-limited process times out"""
    pass


def do_with_timeout(timeout, target, *args, **kwargs):
    """
    Execute a function (or method) with a timeout.

    Call the function (or method) ``target``, with arguments ``args`` and
    keyword arguments ``kwargs``. Normally return the return value
    from ``target``, but if ``target`` takes more than ``timeout`` seconds
    to execute, raises ``TimeoutError``.

    .. note::
        This is, at best, a blunt instrument. Exceptions from ``target`` may
        not propagate properly (tracebacks will be hard to follow.) The
        function which failed to time out may continue to execute until the
        interpreter exits; trapping the TimeoutError and continuing normally
        is not recommended.

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> import time
    >>> def time_me_out():
    ...     time.sleep(5)
    >>> tb.do_with_timeout(0.5, time_me_out) #raises TimeoutError

    Parameters
    ==========
    timeout : float
        Timeout, in seconds.
    target : callable
        Python callable (generally a function, may also be an
            imported ctypes function) to run.
    args : sequence
        Arguments to pass to ``target``.
    kwargs : dict
        keyword arguments to pass to ``target``.

    Raises
    ======
    TimeoutError : If ``target`` does not return in ``timeout`` seconds.

    Returns
    =======
    out :
        return value of ``target``
    """
    import threading
    class ReturningThread(threading.Thread):
        def __init__(self, group=None, target=None, name=None,
                     args=(), kwargs={}):
            #we're handling target, args, kwargs
            super(ReturningThread, self).__init__(group, name=name)
            self._target = target
            self._args = args
            self._kwargs = kwargs
            self._retval = None
            self._exception = None
            
        def run(self):
            try:
                self._retval = self._target(*self._args, **self._kwargs)
            except:
                self._exception = sys.exc_info()

        def join(self, *args, **kwargs):
            super(ReturningThread, self).join(*args, **kwargs)
            if not self._exception is None:
                try:
                    raise self._exception[1].with_traceback(self._exception[2])
                except AttributeError:
                    raise (self._exception[1], None, self._exception[2])
            return self._retval
            
    t = ReturningThread(None, target, None, args, kwargs)
    t.start()
    retval = t.join(timeout)
    if t.is_alive():
        raise TimeoutError()
    return retval


def timeout_check_call(timeout, *args, **kwargs):
    """
    Call a subprocess with a timeout.

    Like :func:`subprocess.check_call`, but will terminate the process and
    raise :exc:`TimeoutError` if it runs for too long.

    This will only terminate the single process started; any child processes
    will remain running (this has implications for, say, spawing shells.)

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> tb.timeout_check_call(1, 'sleep 30', shell=True) #raises TimeoutError

    Parameters
    ==========
    timeout : float
        Timeout, in seconds. Fractions are acceptable but the resolution is of
        order 100ms.
    args : sequence
        Arguments passed through to :class:`subprocess.Popen`
    kwargs : dict
        keyword arguments to pass to :class:`subprocess.Popen`

    Raises
    ======
    TimeoutError : If subprocess does not return in ``timeout`` seconds.
    CalledProcessError : if command has non-zero exit status

    Returns
    =======
    out : int
        0 on successful completion
    """
    resolution = 0.1
    pro = subprocess.Popen(*args, **kwargs)
    starttime = time.time()
    while pro.poll() is None:
        time.sleep(resolution)
        if time.time() - starttime > timeout:
            to = time.time() - starttime
            pro.terminate()
            time.sleep(resolution)
            if pro.poll() is None:
                pro.kill()
            raise TimeoutError('Timed out after {0:.1f} seconds'.format(to))
    if pro.returncode:
        raise subprocess.CalledProcessError(
            pro.returncode, kwargs['args'] if 'args' in kwargs else args[0])
    return 0


def poisson_fit(data, initial=None, method='Powell'):
    """
    Fit a Poisson distribution to data using the method and initial guess provided.

    Parameters
    ==========
    data : array-like
        Data to fit a Poisson distribution to.
    initial : int or None
        initial guess for the fit, if None np.median(data) is used
    method : str
        method passed to scipy.optimize.minimize, default='Powell'

    Examples
    ========
    >>> import spacepy.toolbox as tb
    >>> from scipy.stats import poisson
    >>> import matplotlib.pyplot as plt
    >>> import numpy as np
    >>> data = poisson.rvs(20, size=1000)
    >>> res = tb.poisson_fit(data)
    >>> print(res.x)
    19.718000038769095
    >>> xvals = np.arange(0, np.max(data)+5)
    >>> plt.hist(data, bins=xvals, normed=True)
    >>> plt.plot(xvals, poisson.pmf(xvals, np.round(res.x)))

    Returns
    =======
    result : scipy.optimize.optimize.OptimizeResult
        Resulting fit results from scipy.optimize, answer is result.x,
        user should likely round.

    """
    from scipy.optimize import minimize
    from scipy.stats import poisson

    def negLogLikelihood(params, data):
        """ the negative log-Likelihood-Function"""
        return (- np.sum(np.log(poisson.pmf(data, params[0]))))

    if initial is None:
        initial = np.median(data)

    ans = minimize(negLogLikelihood,  # function to minimize
                   x0=initial,  # start value
                   args=(data,),  # additional arguments for function
                   method=method,  # minimization method, see docs
                   )
    return ans
