def _check(self, lblev=37.0,
            blev=9596.3, brlev=9500.0, brsvd1=9800.0,
            bhlev=0.35, bhrlev=0.31, brsvd2=0.39,
            dim=None):
     lbvc = 65
     lbcode = _lbcode(0)  # unused
     stash = STASH(1, 1, 1)  # unused
     coords_and_dims, factories = _convert_vertical_coords(
         lbcode=lbcode, lbvc=lbvc, blev=blev, lblev=lblev, stash=stash,
         bhlev=bhlev, bhrlev=bhrlev, brsvd1=brsvd1, brsvd2=brsvd2,
         brlev=brlev, dim=dim)
     expect_coords_and_dims = [
         (DimCoord(lblev,
                   standard_name='model_level_number',
                   attributes={'positive': 'up'}), dim)]
     brlev = np.atleast_1d(brlev)
     brsvd1 = np.atleast_1d(brsvd1)
     expect_coords_and_dims.append(
         (DimCoord(blev,
                   long_name='level_height', units='m',
                   bounds=np.vstack((brlev, brsvd1)).T,
                   attributes={'positive': 'up'}), dim))
     bhrlev = np.atleast_1d(bhrlev)
     brsvd2 = np.atleast_1d(brsvd2)
     expect_coords_and_dims.append(
         (AuxCoord(bhlev,
                   long_name='sigma',
                   bounds=np.vstack((bhrlev, brsvd2)).T), dim))
     expect_factories = [(HybridHeightFactory,
                          [{'long_name': 'level_height'},
                           {'long_name': 'sigma'},
                           Reference('orography')])]
     self.assertCoordsAndDimsListsMatch(coords_and_dims,
                                        expect_coords_and_dims)
     self.assertEqual(factories, expect_factories)

 def _make_segment(self,x,y,threshold=None):
     if threshold is None:
         threshold = self._segment_threshold
     x,y=np.atleast_1d(x),np.atleast_1d(y)
     d2 = np.sqrt((np.roll(x,1)-x)**2+(np.roll(y,1)-y)**2)
     w=np.where(d2 > threshold)[0]
     #w=w[w!=0]
     xx=[]
     yy=[]
     if len(w) == 1:
         x=np.roll(x,-w[0])
         y=np.roll(y,-w[0])
         xx.append(x)
         yy.append(y)
     elif len(w) >= 2:
         xx.append(x[0:w[0]])
         yy.append(y[0:w[0]])
         for i in xrange(len(w)-1):
             xx.append(x[w[i]:w[i+1]])
             yy.append(y[w[i]:w[i+1]])
         xx.append(x[w[-1]:])
         yy.append(y[w[-1]:])
     else:
         xx.append(x)
         yy.append(y)
     return xx,yy

    def mm_to_galvo(self, x, y):
        """ Given one or many points in mm space, map them to galvo space.
            e.g.,
            >>> Printer.mm_to_galvo(0, 0) # -> galvo ticks for middle of build area.
            >>> Printer.mm_to_galvo([[0, 1, 2], [0, 0, 0]]) # -> A three-segment line along the x axis.
            The returned array is 2xN, where N is the number of source points
        """
        xshape = np.shape(x)
        if self._grid_table is None:
            grid = np.array(self.read_grid_table())
            assert grid.shape == (5, 5, 2)

            pts_mm = np.linspace(-64, 64, 5) # Grid positions in mm

            # Interpolators for X and Y values (mm to galvo ticks)
            fit_x = scipy.interpolate.interp2d(pts_mm, pts_mm, grid[:,:,0])
            fit_y = scipy.interpolate.interp2d(pts_mm, pts_mm, grid[:,:,1])
            self._grid_table = (fit_x, fit_y)

        if np.shape(x) != np.shape(y):
            raise TypeError('x and y shapes must match. Got x.shape: {}, y.shape: {}'.format(np.shape(x), np.shape(y)))

        x = np.atleast_1d(x)
        y = np.atleast_1d(y)

        x_ = [self._grid_table[0](a, b) for a, b in zip(x, y)]
        y_ = [self._grid_table[1](a, b) for a, b in zip(x, y)]

        result = np.hstack([x_, y_]).T
        if xshape == (): # If it's called with scalars, return a flat result.
            return result.flatten()
        return result

 def record_values(self, updateTimeOnly=False):
    '''Records the current variable values of the space craft into a private list (to be plotted later)
    '''
    # Saves variable values!
    # Append the time variable:
    # O: FIX THESE
    self.__values["time"].append(self.__time)
    self.__values["cellid"].append(self.__cellid)
    self.__values["xcoordinates"].append(self.__coordinates[0])
    self.__values["ycoordinates"].append(self.__coordinates[1])
    self.__values["zcoordinates"].append(self.__coordinates[2])
    self.__values["timeindexes"].append(self.__time_index)
    #self.__timevariable.append(self.__time)
    #self.__cellidlist.append(self.__cellid)
    #self.__coordinateslist.append(self.__coordinates)
    # Set other variables if updateTimeOnly is false:
    if updateTimeOnly == False:
       for i in np.atleast_1d(self.__variables):
          # Append the value at spacecraft's position:
          self.__values[i].append(np.array(self.__vlsvfile.read_variable(i, self.__cellid)))
    else:
       for i in np.atleast_1d(self.__variables):
          # Append the value at spacecraft's position:
          # Note: Appends the latest value again so if list is [[5,2], [0, 3]] then after appending it looks like this: [[5,2,2], [0,3,3]]
          self.__values[i].append(np.array(self.__values[i][len(self.__values[i])-1]))

    def xy2lonlat(self, x, y):
        """Calculate x,y in own projection from given lon,lat (scalars/arrays).
        """
        if self.projected is True:
            if self.proj.is_latlong():
                return x, y
            else:
                if 'ob_tran' in self.proj4:
                    logging.info('NB: Converting degrees to radians ' +
                                 'due to ob_tran srs')
                    x = np.radians(np.array(x))
                    y = np.radians(np.array(y))
                return self.proj(x, y, inverse=True)
        else:
            np.seterr(invalid='ignore')  # Disable warnings for nan-values
            y = np.atleast_1d(np.array(y))
            x = np.atleast_1d(np.array(x))

            # NB: mask coordinates outside domain
            x[x < self.xmin] = np.nan
            x[x > self.xmax] = np.nan
            y[y < self.ymin] = np.nan
            y[y < self.ymin] = np.nan

            lon = map_coordinates(self.lon, [y, x], order=1,
                                  cval=np.nan, mode='nearest')
            lat = map_coordinates(self.lat, [y, x], order=1,
                                  cval=np.nan, mode='nearest')
            return (lon, lat)

def logit(prop, max_events=None):
    """Convert proportion (expressed in the range [0, 1]) to logit.

    Parameters
    ----------
    prop : float | array-like
        the occurrence proportion.
    max_events : int | array-like | None
        the number of events used to calculate ``prop``. Used in a correction
        factor for cases when ``prop`` is 0 or 1, to prevent returning ``inf``.
        If ``None``, no correction is done, and ``inf`` or ``-inf`` may result.

    Returns
    -------
    lgt : ``numpy.ndarray``, with shape matching ``numpy.array(prop).shape``.
    """
    prop = np.atleast_1d(prop).astype(float)
    if np.any([prop > 1, prop < 0]):
        raise ValueError('Proportions must be in the range [0, 1].')
    if max_events is not None:
        # add equivalent of half an event to 0s, and subtract same from 1s
        max_events = np.atleast_1d(max_events) * np.ones_like(prop)
        corr_factor = 0.5 / max_events
        for loc in zip(*np.where(prop == 0)):
            prop[loc] = corr_factor[loc]
        for loc in zip(*np.where(prop == 1)):
            prop[loc] = 1 - corr_factor[loc]
    return np.log(prop / (np.ones_like(prop) - prop))

    def __call__(self,x,y,dx=0,dy=0):
        """Interpolate the function.

        Parameters
        ----------
        x : 1D array
            x-coordinates of the mesh on which to interpolate.
        y : 1D array
            y-coordinates of the mesh on which to interpolate.
        dx : int >= 0, < kx
            Order of partial derivatives in x.
        dy : int >= 0, < ky
            Order of partial derivatives in y.

        Returns
        -------
        z : 2D array with shape (len(y), len(x))
            The interpolated values.

        """

        x = atleast_1d(x)
        y = atleast_1d(y)
        z = fitpack.bisplev(x, y, self.tck, dx, dy)
        z = atleast_2d(z)
        z = transpose(z)
        if len(z)==1:
            z = z[0]
        return array(z)

 def __init__(self, wsize, tstep, n_coefs):  # noqa: D102
     self.wsize = np.atleast_1d(wsize)
     self.tstep = np.atleast_1d(tstep)
     self.n_coefs = np.atleast_1d(n_coefs)
     self.n_dicts = len(tstep)
     self.n_freqs = wsize // 2 + 1
     self.n_steps = self.n_coefs // self.n_freqs

    def headalongline(self, x, y, t, layers=None):
        """Head along line or curve
        
        Parameters
        ----------
        x : array
            x values of line
        y : array
            y values of line
        t : list or array
            times for which grid is returned
        layers : integer, list or array, optional
            layers for which grid is returned
        
        Returns
        -------
        h : array size `nlayers, ntimes, nx`

        """
        
        xg = np.atleast_1d(x)
        yg = np.atleast_1d(y)
        if layers is None:
            Nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq
        else:
            Nlayers = len(np.atleast_1d(layers))
        nx = len(xg)
        if len(yg) == 1:
            yg = yg * np.ones(nx)
        t = np.atleast_1d(t)
        h = np.zeros( (Nlayers,len(t),nx) )
        for i in range(nx):
            h[:,:,i] = self.head(xg[i], yg[i], t, layers)
        return h

def coord_transform(x, y, z, affine):
    """ Convert the x, y, z coordinates from one image space to another
        space.

        Parameters
        ----------
        x : number or ndarray
            The x coordinates in the input space
        y : number or ndarray
            The y coordinates in the input space
        z : number or ndarray
            The z coordinates in the input space
        affine : 2D 4x4 ndarray
            affine that maps from input to output space.

        Returns
        -------
        x : number or ndarray
            The x coordinates in the output space
        y : number or ndarray
            The y coordinates in the output space
        z : number or ndarray
            The z coordinates in the output space

        Warning: The x, y and z have their Talairach ordering, not 3D
        numy image ordering.
    """
    coords = np.c_[np.atleast_1d(x).flat,
                   np.atleast_1d(y).flat,
                   np.atleast_1d(z).flat,
                   np.ones_like(np.atleast_1d(z).flat)].T
    x, y, z, _ = np.dot(affine, coords)
    return x.squeeze(), y.squeeze(), z.squeeze()

 def _compute_model(self, pset):
     """Computes a model and inserts results into the Mongo collection."""
     nBands = fsps.driver.get_n_bands()
     nLambda = fsps.driver.get_n_lambda()
     nAges = fsps.driver.get_n_ages()
     fsps.driver.comp_sp(pset['dust_type'], pset['zmet'], pset['sfh'],
         pset['tau'], pset['const'], pset['fburst'], pset['tburst'],
         pset['dust_tesc'], pset['dust1'], pset['dust2'],
         pset['dust_clumps'], pset['frac_nodust'], pset['dust_index'],
         pset['mwr'], pset['wgp1'], pset['wgp2'], pset['wgp3'],
         pset['duste_gamma'], pset['duste_umin'], pset['duste_qpah'],
         pset['tage'])
     if pset['tage'] == 0.:
         # SFH over all ages is returned
         mags = fsps.driver.get_csp_mags(nBands, nAges)
         specs = fsps.driver.get_csp_specs(nLambda, nAges)
         age, mass, lbol, sfr, dust_mass = fsps.driver.get_csp_stats(nAges)
     else:
         # get only a single age, stored in first age bin
         # arrays must be re-formated to appear like one-age versions of
         # the outputs from get_csp_mags, etc.
         mags = fsps.driver.get_csp_mags_at_age(1, nBands)
         specs = fsps.driver.get_csp_specs_at_age(1, nLambda)
         age, mass, lbol, sfr, dust_mass \
                 = fsps.driver.get_csp_stats_at_age(1)
         age = np.atleast_1d(age)
         mass = np.atleast_1d(mass)
         lbol = np.atleast_1d(lbol)
         sfr = np.atleast_1d(sfr)
         dust_mass = np.atleast_1d(dust_mass)
         mags = np.atleast_2d(mags)
         specs = np.atleast_2d(specs)
     dataArray = self._splice_mag_spec_arrays(age, mass, lbol, sfr,
             dust_mass, mags, specs, nLambda)
     self._insert_model(pset.name, dataArray)

def dataqc_gradienttest_wrapper(dat, x, ddatdx, mindx, startdat, toldat, strict_validation=False):
    if is_none(ddatdx) or is_fill(ddatdx) or is_none(mindx) or is_fill(mindx) or is_none(startdat) or is_fill(startdat) or is_none(toldat) or is_fill(toldat):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    outqc = dataqc_gradienttest(dat, x, [-np.atleast_1d(ddatdx)[-1], np.atleast_1d(ddatdx)[-1]], np.atleast_1d(mindx)[-1], np.atleast_1d(startdat)[-1], np.atleast_1d(toldat)[-1], strict_validation=strict_validation)
    return outqc

    def _make_tuples(self, key):
        print('Populating', key)
        spikes = np.vstack([s.squeeze() for s in (preprocess.Spikes.RateTrace() & key).fetch('rate_trace')])
        s = spikes.sum(axis=0)
        nans = np.isnan(s)

        key['leading_nans'] = int(nans[0])
        key['trailing_nans'] = int(nans[1])

        t = (preprocess.Sync() & key).fetch1('frame_times')  # does not need to be unique

        flip_first = (vis.Trial() * preprocess.Sync().proj('psy_id', trial_idx='first_trial') & key).fetch1('flip_times')
        flip_last = (vis.Trial() * preprocess.Sync().proj('psy_id', trial_idx='last_trial') & key).fetch1('flip_times')

        # (vis.Trial() * preprocess.Sync() & 'trial_idx between first_trial and last_trial')
        fro = np.atleast_1d(flip_first.squeeze())[0]
        to = np.atleast_1d(flip_last.squeeze())[
            -1]  # not necessarily where the stimulus stopped, just last presentation
        idx_fro = np.argmin(np.abs(t - fro))
        idx_to = np.argmin(np.abs(t - to)) + 1
        key['stimulus_nans'] = int(np.any(nans[idx_fro:idx_to]))
        if np.any(nans):
            key['nan_idx'] = nans
        key['stimulus_start'] = idx_fro + 1
        key['stimulus_end'] = idx_to

        self.insert1(key)

    def __new__(cls, val1, val2, scale, precision,
                in_subfmt, out_subfmt, from_jd=False):
        """
        Use __new__ instead of __init__ to output a class instance that
        is the same as the class of the first Time object in the list.
        """
        val1_0 = val1.flat[0]
        if not (isinstance(val1_0, Time) and all(type(val) is type(val1_0)
                                                 for val in val1.flat)):
            raise TypeError('Input values for {0} class must all be same '
                            'astropy Time type.'.format(cls.name))

        if scale is None:
            scale = val1_0.scale
        if val1.shape:
            vals = [getattr(val, scale)._time for val in val1]
            jd1 = np.concatenate([np.atleast_1d(val.jd1) for val in vals])
            jd2 = np.concatenate([np.atleast_1d(val.jd2) for val in vals])
        else:
            val = getattr(val1_0, scale)._time
            jd1, jd2 = val.jd1, val.jd2

        OutTimeFormat = val1_0._time.__class__
        self = OutTimeFormat(jd1, jd2, scale, precision, in_subfmt, out_subfmt,
                             from_jd=True)

        return self

def day_of_year_to_cal(year, N, gregorian=True):
    """Convert a day of year number to a month and day in the Julian or
    Gregorian calendars.

    Arguments:
      - `year`      : year
      - `N`         : day of year, 1..365 (or 366 for leap years)

    Keywords:
      - `gregorian` : If True, use Gregorian calendar, else use Julian calendar
        (default: True)

    Return:
      - (month, day) : (tuple)

    """
    year = np.atleast_1d(year)
    N = np.atleast_1d(N)
    year, N = np.broadcast_arrays(year, N)
    K = np.ones_like(N)
    K[:] = 2
    K[np.atleast_1d(is_leap_year(year, gregorian))] = 1
    mon = (9 * (K + N) / 275.0 + 0.98).astype(np.int64)
    mon[N < 32] = 1
    day = (N - (275 * mon / 9.0).astype(np.int64) +
           K * ((mon + 9) / 12.0).astype(np.int64) + 30).astype(np.int64)
    return _scalar_if_one(mon), _scalar_if_one(day)

    def test_prs_bottilt_ccmp(self):
        """
        Test prs_bottilt_ccmp function.

        Values based on those described in DPS as available on Alfresco:

        OOI (2012). Data Product Specification for Seafloor High-Resolution
            Tilt. Document Control Number 1341-00060.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00060_Data_Product_SPEC_BOTTILT_OOI.pdf)

        Note, DPS does not specify a test for this function. Values from CCMP
        lookup table in DPS used to create this test by the function and
        test_function author.

        Implemented by Christopher Wingard, July 2013
        """
        # set known inputs
        scmp = np.atleast_1d([row[2] for row in self.lily])
        snum = np.atleast_1d([row[3] for row in self.lily])

        # set known output for the corrected compass direction
        ccmp = np.atleast_1d([row[4] for row in self.lily])

        # calculate the corrected compass direction
        out = prsfunc.prs_bottilt_ccmp(scmp, snum)

        # How'd we do?
        np.testing.assert_array_equal(out, ccmp)

def frac_yr_to_jd(year, gregorian=True):
    """Convert a date in the Julian or Gregorian fractional year to the
    Julian Day Number (Meeus 7.1).

    Arguments:
      - `year` : (int, float)  year

    Keywords:
      - `gregorian` : (bool, default=True) If True, use Gregorian calendar,
        else use Julian calendar

    Returns:
      - (float)

    """
    year = np.atleast_1d(year)
    day = np.atleast_1d(0.0).astype(np.float64)
    year, day = list(map(np.array, np.broadcast_arrays(year, day)))
    # For float years abuse the day variable
    fyear = year - year.astype('i')
    mask = fyear > 0
    if np.any(mask):
        year = year.astype('i')
        days_in_year = cal_to_jd(year[mask] + 1) - cal_to_jd(year[mask])
        day[mask] = days_in_year*fyear[mask]
        return _scalar_if_one(cal_to_jd(year) + day)
    return _scalar_if_one(cal_to_jd(year))

    def test_prs_bottilt_tmag(self):
        """
        Test prs_bottilt_tmag function.

        Note, DPS specifies a test data set that does not vary the X-tilt and
        Y-tilt inputs to this function. They are 330 and -330 microradians,
        respectively, for the entirety of the dataset. Test values below were
        created by function and test_function author to ensure function
        operates as expected over a range of potential values.

        OOI (2012). Data Product Specification for Seafloor High-Resolution
            Tilt. Document Control Number 1341-00060.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00060_Data_Product_SPEC_BOTTILT_OOI.pdf)

        Implemented by Christopher Wingard, July 2013
        """
        # set inputs
        xtilt = np.atleast_1d([row[0] for row in self.lily])
        ytilt = np.atleast_1d([row[1] for row in self.lily])

        # set known output for the tilt magnitude
        tmag = np.atleast_1d([row[5] for row in self.lily])

        # calculate the tilt magnitude
        out = prsfunc.prs_bottilt_tmag(xtilt, ytilt)

        # How'd we do?
        np.testing.assert_allclose(out, tmag, rtol=0.1, atol=0.1)

def isnumeric(dat):
    """
    isnumeric - Determine whether input is numeric
    Syntax
    tf = isnumeric(A)
    Description
    tf = isnumeric(A) returns logical 1 (true) if A is a numeric array and logical 0 (false)
    otherwise. For example, sparse arrays and double-precision arrays are numeric, while strings,
    cell arrays, and structure arrays and logicals are not.
    Examples
    Given the following cell array,
    C{1,1} = pi;                 % double
    C{1,2} = 'John Doe';         % char array
    C{1,3} = 2 + 4i;             % complex double
    C{1,4} = ispc;               % logical
    C{1,5} = magic(3)            % double array
    C =
       [3.1416] 'John Doe' [2.0000+ 4.0000i] [1][3x3 double]
    isnumeric shows that all but C{1,2} and C{1,4} are numeric arrays.
    for k = 1:5
    x(k) = isnumeric(C{1,k});
    end
    x
    x =
         1     0     1     0     1

    """

    return np.array([np.atleast_1d(d).dtype.kind in NUMERIC_KINDS for d in np.nditer(np.atleast_1d(dat))]).astype('int8')

def hyp0f1(v, z):
    r"""Confluent hypergeometric limit function 0F1.

    Parameters
    ----------
    v, z : array_like
        Input values.

    Returns
    -------
    hyp0f1 : ndarray
        The confluent hypergeometric limit function.

    Notes
    -----
    This function is defined as:

    .. math:: _0F_1(v,z) = \sum_{k=0}^{\inf}\frac{z^k}{(v)_k k!}.

    It's also the limit as q -> infinity of ``1F1(q;v;z/q)``, and satisfies
    the differential equation :math:``f''(z) + vf'(z) = f(z)`.
    """
    v = atleast_1d(v)
    z = atleast_1d(z)
    v, z = np.broadcast_arrays(v, z)
    arg = 2 * sqrt(abs(z))
    old_err = np.seterr(all='ignore')  # for z=0, a<1 and num=inf, next lines
    num = where(z.real >= 0, iv(v - 1, arg), jv(v - 1, arg))
    den = abs(z)**((v - 1.0) / 2)
    num *= gamma(v)
    np.seterr(**old_err)
    num[z == 0] = 1
    den[z == 0] = 1
    return num / den