def cloud_statistics(file_name):
    """
    Return core duration, minimum core base, maximum core height, mean core 
    mass, formation time, dissipation time, maximum depth, depth evolution and 
    corresponding times for tracked clouds.
        
    Parameters
    ----------
    file_name : netCDF file name
        id_profile file for a tracked core with dimensions double t(t), 
        double z(z).
      
    Return
    ------
    tuple : cloud_id, lifetime, base, top, mass, l_min, l_max, depths,
        max_depth, times
    """
    
    # Read netCDF dataset
    data = Dataset(file_name)
    
    # Core ID
    cloud_id = int(file_name[-11:-3])
    
    # Core duration (seconds)
    times = data.variables['t'][...]
    lifetime = len(times)*mc.dt
    
    # Formation time, dissipation time (seconds)
    l_min = times.min()*mc.dt
    l_max = times.max()*mc.dt

    # Minimum core base, maximum core height, maximum depth, depth evolution 
    # (metres)
    area = ma.masked_invalid(data.variables['AREA'][...])
    z = data.variables['z'][...]
    z = z*np.ones(np.shape(area))
    z = ma.masked_array(z, ma.getmask(area)) 
    bases = z.min(axis=1)
    tops = z.max(axis=1)
    depths = tops - bases + mc.dz
    max_depth = depths.max()
    base = bases.min()
    top = tops.max()

    # Mean core mass mass (kilograms)
    qn = ma.masked_invalid(data.variables['QN'][...])
    rho = ma.masked_invalid(data.variables['RHO'][...])
    mass = np.mean(np.sum(area*rho*mc.dz, axis=1))

    # Remove missing values
    times = ma.masked_array(times, ma.getmask(depths))
    depths = depths[~depths.mask]
    times = times[~times.mask]
    
    data.close()
    
    return cloud_id, lifetime, base, top, mass, l_min, l_max, depths, \
        max_depth, times

 def outer (self, a, b):
     "Return the function applied to the outer product of a and b."
     a = _makeMaskedArg(a)
     b = _makeMaskedArg(b)
     ma = getmask(a)
     mb = getmask(b)
     if ma is nomask and mb is nomask:
         m = None
     else:
         ma = getmaskarray(a)
         mb = getmaskarray(b)
         m = logical_or.outer(ma, mb)
     d = numpy.maximum.outer(filled(a), filled(b))
     return TransientVariable(d, mask=m)

 def test_testCI(self):
     # Test of conversions and indexing
     x1 = np.array([1, 2, 4, 3])
     x2 = array(x1, mask=[1, 0, 0, 0])
     x3 = array(x1, mask=[0, 1, 0, 1])
     x4 = array(x1)
     # test conversion to strings
     str(x2)  # raises?
     repr(x2)  # raises?
     assert_(eq(np.sort(x1), sort(x2, fill_value=0)))
     # tests of indexing
     assert_(type(x2[1]) is type(x1[1]))
     assert_(x1[1] == x2[1])
     assert_(x2[0] is masked)
     assert_(eq(x1[2], x2[2]))
     assert_(eq(x1[2:5], x2[2:5]))
     assert_(eq(x1[:], x2[:]))
     assert_(eq(x1[1:], x3[1:]))
     x1[2] = 9
     x2[2] = 9
     assert_(eq(x1, x2))
     x1[1:3] = 99
     x2[1:3] = 99
     assert_(eq(x1, x2))
     x2[1] = masked
     assert_(eq(x1, x2))
     x2[1:3] = masked
     assert_(eq(x1, x2))
     x2[:] = x1
     x2[1] = masked
     assert_(allequal(getmask(x2), array([0, 1, 0, 0])))
     x3[:] = masked_array([1, 2, 3, 4], [0, 1, 1, 0])
     assert_(allequal(getmask(x3), array([0, 1, 1, 0])))
     x4[:] = masked_array([1, 2, 3, 4], [0, 1, 1, 0])
     assert_(allequal(getmask(x4), array([0, 1, 1, 0])))
     assert_(allequal(x4, array([1, 2, 3, 4])))
     x1 = np.arange(5) * 1.0
     x2 = masked_values(x1, 3.0)
     assert_(eq(x1, x2))
     assert_(allequal(array([0, 0, 0, 1, 0], MaskType), x2.mask))
     assert_(eq(3.0, x2.fill_value))
     x1 = array([1, 'hello', 2, 3], object)
     x2 = np.array([1, 'hello', 2, 3], object)
     s1 = x1[1]
     s2 = x2[1]
     assert_equal(type(s2), str)
     assert_equal(type(s1), str)
     assert_equal(s1, s2)
     assert_(x1[1:1].shape == (0,))

 def __call__(self, value, clip=None):
     if clip is None:
         clip = self.clip
     if cbook.iterable(value):
         vtype = 'array'
         val = ma.asarray(value).astype(np.float)
     else:
         vtype = 'scalar'
         val = ma.array([value]).astype(np.float)
     self.autoscale_None(val)
     vmin, vmax = self.vmin, self.vmax
     cmin, cmax = self.cmin * vmin, self.cmax * vmax
     if vmin > vmax:
         raise ValueError("minvalue must be less than or equal to maxvalue")
     elif vmin == vmax:
         result = 0.0 * val
     else:
         if clip:
             mask = ma.getmask(val)
             val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                             mask=mask)
         result = 0. * val + 0.5
         result[val > cmax] = (ma.log10(val[val > cmax]) - ma.log10(cmax)) / (np.log10(vmax) - np.log10(cmax)) / 2. + 0.5
         result[val < cmin] = -(ma.log10(-val[val < cmin]) - ma.log10(-cmin)) / (np.log10(-vmin) - np.log10(-cmin)) / 2. + 0.5
     if vtype == 'scalar':
         result = result[0]
     return result

def plot_spin_temp_lv(values, ax, key, label, max_clip=None):
    scale_vals = values[key]
    dataset = values[~ma.getmask(scale_vals)]
    if max_clip:
        scale_vals = np.clip(dataset[key], 0, max_clip)
    base = np.array([dataset['longitude'], dataset['Velocity'], scale_vals],
                    dtype=[('l', float), ('b', float), ('scale', float)])

    sample = base
    #vmax = np.max(sample[2,:])
    #cm = plt.cm.get_cmap('RdYlBu_r')
    cm = plt.cm.get_cmap('viridis')

    x = sample[0,:]
    y = sample[1,:]
    z = sample[2,:]
    idx = z.argsort()
    x, y, z = x[idx], y[idx], z[idx]

    # spin temp plot needs log, min 3
    sc = ax.scatter(x, y, c=z, s=25, cmap=cm, norm=matplotlib.colors.PowerNorm(0.5, vmin=3))
    formatter = LogFormatter(10, labelOnlyBase=False)
    ticks = [0, 5, 10, 20, 40, 70, 100, 150, 200, 300]
    cb = plt.colorbar(sc, ticks=ticks, format=formatter)

    ax.set_xlim(values['longitude'].max() + 5, values['longitude'].min() - 5)

    ax.set_xlabel('Galactic longitude (deg)')
    ax.set_ylabel('LSR Velocity (km/s)')
    cb.set_label(label)

    return None

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        result, is_scalar = self.process_value(value)
        self.autoscale_None(result)
        vmin, vmax = self.vmin, self.vmax

        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin == vmax:
            result.fill(0)
        else:
            if clip:
                mask = ma.getmask(result)
                result = ma.array(np.clip(result.filled(vmax), vmin, vmax),
                                  mask=mask)
            # in-place equivalent of above can be much faster
            resdat = self._transform(result.data)
            resdat -= self._lower
            resdat /= (self._upper - self._lower)

        if is_scalar:
            result = result[0]
        return result

 def __setslice__(self, i, j, value):
     "Sets the slice described by [i,j] to `value`."
     _localdict = self.__dict__
     d = self._data
     m = _localdict['_fieldmask']
     names = self.dtype.names
     if value is masked:
         for n in names:
             m[i:j][n] = True
     elif not self._hardmask:
         fval = filled(value)
         mval = getmaskarray(value)
         for n in names:
             d[n][i:j] = fval
             m[n][i:j] = mval
     else:
         mindx = getmaskarray(self)[i:j]
         dval = np.asarray(value)
         valmask = getmask(value)
         if valmask is nomask:
             for n in names:
                 mval = mask_or(m[n][i:j], valmask)
                 d[n][i:j][~mval] = value
         elif valmask.size > 1:
             for n in names:
                 mval = mask_or(m[n][i:j], valmask)
                 d[n][i:j][~mval] = dval[~mval]
                 m[n][i:j] = mask_or(m[n][i:j], mval)
     self._fieldmask = m

    def __call__(self, value, clip=None):
        """Map value to the interval [0, 1]. The clip argument is unused."""

        result, is_scalar = self.process_value(value)

        self.autoscale_None(result)
        vmin, vcenter, vmax = self.vmin, self.vcenter, self.vmax
        if vmin == vmax == vcenter:
            result.fill(0)
        elif not vmin <= vcenter <= vmax:
            raise ValueError("minvalue must be less than or equal to "
                             "centervalue which must be less than or "
                             "equal to maxvalue")
        else:
            vmin = float(vmin)
            vcenter = float(vcenter)
            vmax = float(vmax)
            # in degenerate cases, prefer the center value to the extremes
            degen = (result == vcenter) if vcenter == vmax else None

            x, y = [vmin, vcenter, vmax], [0, 0.5, 1]
            result = ma.masked_array(np.interp(result, x, y),
                                     mask=ma.getmask(result))
            if degen is not None:
                result[degen] = 0.5

        if is_scalar:
            result = np.atleast_1d(result)[0]
        return(result)

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        if cbook.iterable(value):
            vtype = 'array'
            val = ma.asarray(value).astype(numpy.float)
        else:
            vtype = 'scalar'
            val = ma.array([value]).astype(numpy.float)
        
        if self.staticrange is None:
            self.autoscale_None(val)
            vmin, vmax = self.vmin, self.vmax
        else:
            self.vmin, self.vmax = None, None
            self.autoscale_None(val)
            vmin, vmax = self.vmax - self.staticrange, self.vmax
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin==vmax:
            result = 0.0 * val
        else:
            vmin = float(vmin)
            vmax = float(vmax)
            rmin = float(self.rmin)
            rmax = float(self.rmax)
            if clip:
                mask = ma.getmask(val)
                val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                                mask=mask)
            result = (val-vmin) * ((rmax-rmin) / (vmax-vmin)) + rmin
        if vtype == 'scalar':
            result = result[0]
        return result

    def onMouseOver(self, event):
        """

        """

        if event.inaxes == self.ax1:
            if (event.xdata != None and event.ydata != None):

                xidx = (int(event.xdata/
                        (float(self.omf.parameters['xstepsize'])*1.0e10)))

                yidx = (int(event.ydata/
                        (float(self.omf.parameters['ystepsize'])*1.0e10)))

                value = self.angle[xidx,yidx]


                if ma.getmask(self.angle)[xidx,yidx]:
                    self.datavalue.SetStatusText('Angle Value: MASKED')

                else:
                    value = -degrees(value)
                    if (value < 0.0): value +=360
                    self.datavalue.SetStatusText('Angle Value: '
                                                 +str('%.2f'%value))
        else:

            self.datavalue.SetLabel('')

        return

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        result, is_scalar = self.process_value(value)

        result = ma.masked_less_equal(result, 0, copy=False)

        self.autoscale_None(result)
        vmin, vmax = self.vmin, self.vmax
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin <= 0:
            raise ValueError("values must all be positive")
        elif vmin == vmax:
            result.fill(0)
        else:
            if clip:
                mask = ma.getmask(result)
                result = ma.array(np.clip(result.filled(vmax), vmin, vmax), mask=mask)
            # in-place equivalent of above can be much faster
            resdat = result.data
            mask = result.mask
            if mask is np.ma.nomask:
                mask = resdat <= 0
            else:
                mask |= resdat <= 0
            cbook._putmask(resdat, mask, 1)
            np.log(resdat, resdat)
            resdat -= np.log(vmin)
            resdat /= np.log(vmax) - np.log(vmin)
            result = np.ma.array(resdat, mask=mask, copy=False)
        if is_scalar:
            result = result[0]
        return result

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        result, is_scalar = self.process_value(value)

        self.autoscale_None(result)
        gamma = self.gamma
        vmin, vmax = self.vmin, self.vmax
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin == vmax:
            result.fill(0)
        else:
            if clip:
                mask = ma.getmask(result)
                val = ma.array(np.clip(result.filled(vmax), vmin, vmax),
                                mask=mask)
            resdat = result.data
            resdat -= vmin
            np.power(resdat, gamma, resdat)
            resdat /= (vmax - vmin) ** gamma
            result = np.ma.array(resdat, mask=result.mask, copy=False)
            result[(value < 0)&~result.mask] = 0
        if is_scalar:
            result = result[0]
        return result

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        if cbook.iterable(value):
            vtype = 'array'
            val = ma.asarray(value).astype(np.float)
        else:
            vtype = 'scalar'
            val = ma.array([value]).astype(np.float)

        self.autoscale_None(val)
        vmin, vmax = self.vmin, self.vmax
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin<=0:
            raise ValueError("values must all be positive")
        elif vmin==vmax:
            return 0.0 * val
        else:
            if clip:
                mask = ma.getmask(val)
                val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                                mask=mask)
            result = (ma.log(val)-np.log(vmin))/(np.log(vmax)-np.log(vmin))
        if vtype == 'scalar':
            result = result[0]
        return result

    def __call__(self, value, clip=None, midpoint=None):


        if clip is None:
            clip = self.clip

        if cbook.iterable(value):
            vtype = 'array'
            val = ma.asarray(value).astype(np.float)
        else:
            vtype = 'scalar'
            val = ma.array([value]).astype(np.float)

        self.autoscale_None(val)
        vmin, vmax = self.vmin, self.vmax

        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin==vmax:
            return 0.0 * val
        else:
            if clip:
                mask = ma.getmask(val)
                val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                                mask=mask)
            result = (val-vmin) * (1.0/(vmax-vmin))
            #result = (ma.arcsinh(val)-np.arcsinh(vmin))/(np.arcsinh(vmax)-np.arcsinh(vmin))
            result = result**(1./self.nthroot)
        if vtype == 'scalar':
            result = result[0]
        return result

def _unscented_correct(cross_sigma, mu_pred, sigma2_pred, obs_mu_pred, obs_sigma2_pred, z):
    """Correct predicted state estimates with an observation

    Parameters
    ----------
    cross_sigma : [n_dim_state, n_dim_obs] array
        cross-covariance between the state at time t given all observations
        from timesteps [0, t-1] and the observation at time t
    mu_pred : [n_dim_state] array
        mean of state at time t given observations from timesteps [0, t-1]
    sigma2_pred : [n_dim_state, n_dim_state] array
        square root of covariance of state at time t given observations from
        timesteps [0, t-1]
    obs_mu_pred : [n_dim_obs] array
        mean of observation at time t given observations from times [0, t-1]
    obs_sigma2_pred : [n_dim_obs] array
        square root of covariance of observation at time t given observations
        from times [0, t-1]
    z : [n_dim_obs] array
        observation at time t

    Returns
    -------
    mu_filt : [n_dim_state] array
        mean of state at time t given observations from time steps [0, t]
    sigma2_filt : [n_dim_state, n_dim_state] array
        square root of covariance of state at time t given observations from
        time steps [0, t]
    """
    n_dim_state = len(mu_pred)
    n_dim_obs = len(obs_mu_pred)

    if not np.any(ma.getmask(z)):
        ##############################################
        # Same as this, but more stable (supposedly) #
        ##############################################
        # K = cross_sigma.dot(
        #     linalg.pinv(
        #         obs_sigma2_pred.T.dot(obs_sigma2_pred)
        #     )
        # )
        ##############################################

        # equivalent to this MATLAB code
        # K = (cross_sigma / obs_sigma2_pred.T) / obs_sigma2_pred
        K = linalg.lstsq(obs_sigma2_pred, cross_sigma.T)[0]
        K = linalg.lstsq(obs_sigma2_pred.T, K)[0]
        K = K.T

        # correct mu, sigma
        mu_filt = mu_pred + K.dot(z - obs_mu_pred)
        U = K.dot(obs_sigma2_pred)
        sigma2_filt = cholupdate(sigma2_pred, U.T, -1.0)
    else:
        # no corrections to be made
        mu_filt = mu_pred
        sigma2_filt = sigma2_pred
    return (mu_filt, sigma2_filt)

def approx (a, b, fill_value=True, rtol=1.e-5, atol=1.e-8):
    """Returns true if all components of a and b are equal subject to given tolerances.

If fill_value is True, masked values considered equal. Otherwise, masked values
are considered unequal.
The relative error rtol should be positive and << 1.0
The absolute error atol comes into play for those elements of b that are very
small or zero; it says how small a must be also.
    """
    m = mask_or(getmask(a), getmask(b))
    d1 = filled(a)
    d2 = filled(b)
    if d1.dtype.char == "O" or d2.dtype.char == "O":
        return np.equal(d1,d2).ravel()
    x = filled(masked_array(d1, copy=False, mask=m), fill_value).astype(float_)
    y = filled(masked_array(d2, copy=False, mask=m), 1).astype(float_)
    d = np.less_equal(umath.absolute(x-y), atol + rtol * umath.absolute(y))
    return d.ravel()

 def ask_chunks(self, count, price, inplace=False):
     assert 0 < count <= self.max_ask
     leftover = self.leftover - count * self.min_bet
     chunks = self.chunks if inplace else self.chunks.copy()
     for _ in range(count):
         idx = np.where(ma.getmask(chunks))[0][0]
         chunks[idx] = self.chunk_value(price)
         chunks.mask[idx] = False
     sort_chunks(chunks)
     return self._get_ret_val(leftover, chunks, inplace)

def forward_fill(marr, maxgap=None):
    """
    Forward fills masked values in a 1-d array when there are less ``maxgap``
    consecutive masked values.

    Parameters
    ----------
    marr : MaskedArray
        Series to fill
    maxgap : {int}, optional
        Maximum gap between consecutive masked values.
        If ``maxgap`` is not specified, all masked values are forward-filled.


    Examples
    --------
    >>> x = ma.arange(20)
    >>> x[(x%5)!=0] = ma.masked
    >>> print x
    [0 -- -- -- -- 5 -- -- -- -- 10 -- -- -- -- 15 -- -- -- --]
    >>> print forward_fill(x)
    [0 0 0 0 0 5 5 5 5 5 10 10 10 10 10 15 15 15 15 15]

    """
    # !!!: We should probably port that to C.
    # Initialization ..................
    if np.ndim(marr) > 1:
        raise ValueError,"The input array should be 1D only!"
    a = ma.array(marr, copy=True)
    amask = getmask(a)
    if amask is nomask or a.size == 0:
        return a
    #
    adata = getdata(a)
    # Get the indices of the masked values (except a[0])
    idxtofill = amask[1:].nonzero()[0] + 1
    currGap = 0
    if maxgap is not None:
        previdx = -1
        for i in idxtofill:
            if i != previdx + 1:
                currGap = 0
            currGap += 1
            if currGap <= maxgap and not amask[i-1]:
                adata[i] = adata[i-1]
                amask[i] = False
                previdx = i
            else:
                amask[i-maxgap:i] = True
    else:
        for i in idxtofill:
            if not amask[i-1]:
                adata[i] = adata[i-1]
                amask[i] = False
    return a

    def _numpy_interpolation(self, point_num, eval_points):
        """

        Parameters
        ----------
        point_num: int
            Index of class position in values list
        eval_points: ndarray
            Inputs used to evaluate class member function

        Returns
        -------
        ndarray: output from member function
        """
        is_masked = ma.is_masked(eval_points)

        shape = point_num.shape
        ev_shape = eval_points.shape

        vals = self.values[point_num.ravel()]
        eval_points = np.repeat(eval_points, shape[1], axis=0)
        it = np.arange(eval_points.shape[0])

        it = np.repeat(it, eval_points.shape[1], axis=0)

        eval_points = eval_points.reshape(
            eval_points.shape[0] * eval_points.shape[1],
            eval_points.shape[-1]
        )

        scaled_points = eval_points.T
        if is_masked:
            mask = np.invert(ma.getmask(scaled_points[0]))
        else:
            mask = np.ones_like(scaled_points[0], dtype=bool)

        it = ma.masked_array(it, mask)
        scaled_points[0] = (
            (scaled_points[0] - (self._bounds[0][0])) /
            (self._bounds[0][1] - self._bounds[0][0])
        ) * (vals.shape[-2] - 1)
        scaled_points[1] += (
            (scaled_points[1] - (self._bounds[1][0])) /
            (self._bounds[1][1] - self._bounds[1][0])
        ) * (vals.shape[-1] - 1)
        scaled_points = np.vstack((it, scaled_points))

        output = np.zeros(scaled_points.T.shape[:-1])
        output[mask] = map_coordinates(vals, scaled_points.T[mask].T, order=1)

        new_shape = (*shape, ev_shape[-2])
        output = output.reshape(new_shape)

        return ma.masked_array(output, mask=mask)

  def __call__(self, value, clip=None):
    if clip is None:
      clip = self.clip

    if cbook.iterable(value):
      vtype = 'array'
      val = ma.asarray(value).astype(np.float)
    else:
      vtype = 'scalar'
      val = ma.array([value]).astype(np.float)

    self.autoscale_None(val)
    vmin, vmax = self.vmin, self.vmax
    vin, cin = self.vin, self.cin
    if vmin > vmax:
      raise ValueError("minvalue must be less than or equal to maxvalue")
    elif vmin > 0:
      raise ValueError("minvalue must be less than 0")
    elif vmax < 0:
      raise ValueError("maxvalue must be greater than 0")
    elif vmin==vmax:
      result = 0.0 * val
    else:
      if clip:
        mask = ma.getmask(val)
        val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                        mask=mask)
      ipos = (val > vin)
      ineg = (val < -vin)
      izero = ~(ipos | ineg)

      result = ma.empty_like(val)
      result[izero] = 0.5 + cin * val[izero] / vin
      result[ipos] = 0.5 + cin + (0.5 - cin) * \
                    (ma.log(val[ipos]) - np.log(vin)) / (np.log(vmax) - np.log(vin))
      result[ineg] = 0.5 - cin - (0.5 - cin) * \
                    (ma.log(-val[ineg]) - np.log(vin)) / (np.log(-vmin) - np.log(vin))
      result.mask = ma.getmask(val)
    if vtype == 'scalar':
      result = result[0]
    return result