def test_kendalltau(self):
     # Tests some computations of Kendall's tau
     x = ma.fix_invalid([5.05, 6.75, 3.21, 2.66, np.nan])
     y = ma.fix_invalid([1.65, 26.5, -5.93, 7.96, np.nan])
     z = ma.fix_invalid([1.65, 2.64, 2.64, 6.95, np.nan])
     assert_almost_equal(np.asarray(mstats.kendalltau(x, y)), [+0.3333333, 0.4969059])
     assert_almost_equal(np.asarray(mstats.kendalltau(x, z)), [-0.5477226, 0.2785987])
     #
     x = ma.fix_invalid([0, 0, 0, 0, 20, 20, 0, 60, 0, 20, 10, 10, 0, 40, 0, 20, 0, 0, 0, 0, 0, np.nan])
     y = ma.fix_invalid([0, 80, 80, 80, 10, 33, 60, 0, 67, 27, 25, 80, 80, 80, 80, 80, 80, 0, 10, 45, np.nan, 0])
     result = mstats.kendalltau(x, y)
     assert_almost_equal(np.asarray(result), [-0.1585188, 0.4128009])

def r2eff_LM63(r20=None, phi_ex=None, kex=None, cpmg_frqs=None, back_calc=None):
    """Calculate the R2eff values for the LM63 model.

    See the module docstring for details.


    @keyword r20:           The R20 parameter value (R2 with no exchange).
    @type r20:              numpy float array of rank [NE][NS][NM][NO][ND]
    @keyword phi_ex:        The phi_ex parameter value (pA * pB * delta_omega^2).
    @type phi_ex:           numpy float array of rank [NE][NS][NM][NO][ND]
    @keyword kex:           The kex parameter value (the exchange rate in rad/s).
    @type kex:              float
    @keyword cpmg_frqs:     The CPMG nu1 frequencies.
    @type cpmg_frqs:        numpy float array of rank [NE][NS][NM][NO][ND]
    @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies.
    @type back_calc:        numpy float array of rank [NE][NS][NM][NO][ND]
    """

    # Flag to tell if values should be replaced if phi_ex is zero.
    t_phi_ex_zero = False

    # Catch divide with zeros (to avoid pointless mathematical operations).
    if kex == 0.0:
        back_calc[:] = r20
        return

    # Catch zeros (to avoid pointless mathematical operations).
    # This will result in no exchange, returning flat lines.
    if min(phi_ex) == 0.0:
        t_phi_ex_zero = True
        mask_phi_ex_zero = masked_where(phi_ex == 0.0, phi_ex)

    # Repetitive calculations (to speed up calculations).
    rex = phi_ex / kex
    kex_4 = 4.0 / kex

    # Calculate R2eff.
    back_calc[:] = r20 + rex * (1.0 - kex_4 * cpmg_frqs * tanh(kex / (4.0 * cpmg_frqs)))

    # Replace data in array.
    # If phi_ex is zero.
    if t_phi_ex_zero:
        back_calc[mask_phi_ex_zero.mask] = r20[mask_phi_ex_zero.mask]

    # Catch errors, taking a sum over array is the fastest way to check for
    # +/- inf (infinity) and nan (not a number).
    if not isfinite(sum(back_calc)):
        # Replaces nan, inf, etc. with fill value.
        fix_invalid(back_calc, copy=False, fill_value=1e100)

    def test_spearmanr(self):
        # Tests some computations of Spearman's rho
        (x, y) = ([5.05, 6.75, 3.21, 2.66], [1.65, 2.64, 2.64, 6.95])
        assert_almost_equal(mstats.spearmanr(x, y)[0], -0.6324555)
        (x, y) = ([5.05, 6.75, 3.21, 2.66, np.nan], [1.65, 2.64, 2.64, 6.95, np.nan])
        (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y))
        assert_almost_equal(mstats.spearmanr(x, y)[0], -0.6324555)

        x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1, 1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7]
        y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6, 0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4]
        assert_almost_equal(mstats.spearmanr(x, y)[0], 0.6887299)
        x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1, 1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7, np.nan]
        y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6, 0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4, np.nan]
        (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y))
        assert_almost_equal(mstats.spearmanr(x, y)[0], 0.6887299)

 def test_zscore(self):
     # This is not in R, so tested by using:
     #     (testcase[i]-mean(testcase,axis=0)) / sqrt(var(testcase)*3/4)
     y = mstats.zscore(self.testcase)
     desired = ma.fix_invalid([-1.3416407864999, -0.44721359549996,
                               0.44721359549996, 1.3416407864999, np.nan])
     assert_almost_equal(desired, y, decimal=12)

def impute_missing_total_reads(total_reads, missing_variant_confidence):
  # Change NaNs to masked values via SciPy.
  masked_total_reads = ma.fix_invalid(total_reads)

  # Going forward, suppose you have v variants and s samples in a v*s matrix of
  # read counts. Missing values are masked.

  # Calculate geometric mean of variant read depth in each sample. Result: s*1
  sample_means = gmean(masked_total_reads, axis=0)
  assert np.sum(sample_means <= 0) == np.sum(np.isnan(sample_means)) == 0
  # Divide every variant's read count by its mean sample read depth to get read
  # depth enrichment relative to other variants in sample. Result: v*s
  normalized_to_sample = np.dot(masked_total_reads, np.diag(1./sample_means))
  # For each variant, calculate geometric mean of its read depth enrichment
  # across samples. Result: v*1
  variant_mean_reads = gmean(normalized_to_sample, axis=1)
  assert np.sum(variant_mean_reads <= 0) == np.sum(np.isnan(variant_mean_reads)) == 0

  # Convert 1D arrays to vectors to permit matrix multiplication.
  imputed_counts = np.dot(variant_mean_reads.reshape((-1, 1)), sample_means.reshape((1, -1)))
  nan_coords = np.where(np.isnan(total_reads))
  total_reads[nan_coords] = imputed_counts[nan_coords]
  assert np.sum(total_reads <= 0) == np.sum(np.isnan(total_reads)) == 0

  total_reads[nan_coords] *= missing_variant_confidence
  return np.floor(total_reads).astype(np.int)

def ccf(x, y, periodogram=True):
    """Computes the auto-correlation of the series x and y at different lags.
The computations are performed on anomalies (deviations from average).
Gaps in the series are filled first, anomalies are then computed and missing
values filled with 0.
If x and y are valid TimeSeries object, they are aligned so that their starting
and ending point match.

Parameters
----------
    x : sequence
        Input data.
    y : sequence
        Input data.
        If y is longer than x, it is truncated to match the length of x.
        If y is shorter than x, x is truncated.
    periodogram : {True, False} optional
        Whether to return a periodogram or a standard estimate of the autocovariance.

Returns
-------
    cvf : ma.array
        Cross-correlation at lags [0,1,...,n,n-1,...,-1]
    """
    ccf_ = cvf(x,y,periodogram)
    return ma.fix_invalid(ccf_/ccf_[0])

def _moving_func(data, cfunc, kwargs):

    data = ma.fix_invalid(data)
    data = ma.array(data.filled(0), mask=data._mask)

    if data.ndim == 1:
        kwargs['array'] = data
        result_dict = cfunc(**kwargs)
        return _process_result_dict(data, result_dict)

    elif data.ndim == 2:
        for i in range(data.shape[-1]):
            kwargs['array'] = data[:,i]
            result_dict = cfunc(**kwargs)

            if i == 0:
                rtype = result_dict['array'].dtype
                result = data.astype(rtype)
                print data.dtype, result.dtype

            rmask = result_dict.get('mask', ma.nomask)

            curr_col = marray(result_dict['array'], mask=rmask, copy=False)
            result[:,i] = curr_col

        return result

    else:
        raise ValueError, "Data should be at most 2D"

 def test_cov(self):
     "Tests the cov function."
     x = ma.array([[1,2,3],[4,5,6]], mask=[[1,0,0],[0,0,0]])
     c = mstats.cov(x[0])
     assert_equal(c, x[0].var(ddof=1))
     c = mstats.cov(x[1])
     assert_equal(c, x[1].var(ddof=1))
     c = mstats.cov(x)
     assert_equal(c[1,0], (x[0].anom()*x[1].anom()).sum())
     #
     x = [[nan,nan,  4,  2, 16, 26,  5,  1,  5,  1,  2,  3,  1],
          [  4,  3,  5,  3,  2,  7,  3,  1,  1,  2,  3,  5,  3],
          [  3,  2,  5,  6, 18,  4,  9,  1,  1,nan,  1,  1,nan],
          [nan,  6, 11,  4, 17,nan,  6,  1,  1,  2,  5,  1,  1]]
     x = ma.fix_invalid(x).T
     (winter,spring,summer,fall) = x.T
     #
     assert_almost_equal(mstats.cov(winter,winter,bias=True),
                         winter.var(ddof=0))
     assert_almost_equal(mstats.cov(winter,winter,bias=False),
                         winter.var(ddof=1))
     assert_almost_equal(mstats.cov(winter,spring)[0,1], 7.7)
     assert_almost_equal(mstats.cov(winter,spring)[1,0], 7.7)
     assert_almost_equal(mstats.cov(winter,summer)[0,1], 19.1111111, 7)
     assert_almost_equal(mstats.cov(winter,summer)[1,0], 19.1111111, 7)
     assert_almost_equal(mstats.cov(winter,fall)[0,1], 20)
     assert_almost_equal(mstats.cov(winter,fall)[1,0], 20)

def hdquantiles_sd(data, prob=list([.25,.5,.75]), axis=None):
    """
    The standard error of the Harrell-Davis quantile estimates by jackknife.

    Parameters
    ----------
    data : array_like
        Data array.
    prob : sequence, optional
        Sequence of quantiles to compute.
    axis : int, optional
        Axis along which to compute the quantiles. If None, use a flattened
        array.

    Returns
    -------
    hdquantiles_sd : MaskedArray
        Standard error of the Harrell-Davis quantile estimates.

    See Also
    --------
    hdquantiles

    """
    def _hdsd_1D(data, prob):
        "Computes the std error for 1D arrays."
        xsorted = np.sort(data.compressed())
        n = len(xsorted)

        hdsd = np.empty(len(prob), float_)
        if n < 2:
            hdsd.flat = np.nan

        vv = np.arange(n) / float(n-1)
        betacdf = beta.cdf

        for (i,p) in enumerate(prob):
            _w = betacdf(vv, (n+1)*p, (n+1)*(1-p))
            w = _w[1:] - _w[:-1]
            mx_ = np.fromiter([np.dot(w,xsorted[np.r_[list(range(0,k)),
                                                      list(range(k+1,n))].astype(int_)])
                                  for k in range(n)], dtype=float_)
            mx_var = np.array(mx_.var(), copy=False, ndmin=1) * n / float(n-1)
            hdsd[i] = float(n-1) * np.sqrt(np.diag(mx_var).diagonal() / float(n))
        return hdsd

    # Initialization & checks
    data = ma.array(data, copy=False, dtype=float_)
    p = np.array(prob, copy=False, ndmin=1)
    # Computes quantiles along axis (or globally)
    if (axis is None):
        result = _hdsd_1D(data, p)
    else:
        if data.ndim > 2:
            raise ValueError("Array 'data' must be at most two dimensional, "
                             "but got data.ndim = %d" % data.ndim)
        result = ma.apply_along_axis(_hdsd_1D, axis, data, p)

    return ma.fix_invalid(result, copy=False).ravel()

 def test_kendalltau_seasonal(self):
     # Tests the seasonal Kendall tau.
     x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1],
          [4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3],
          [3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan],
          [nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]]
     x = ma.fix_invalid(x).T
     output = mstats.kendalltau_seasonal(x)
     assert_almost_equal(output['global p-value (indep)'], 0.008, 3)
     assert_almost_equal(output['seasonal p-value'].round(2),
                         [0.18,0.53,0.20,0.04])

    def test_kendalltau(self):
        # Tests some computations of Kendall's tau
        x = ma.fix_invalid([5.05, 6.75, 3.21, 2.66,np.nan])
        y = ma.fix_invalid([1.65, 26.5, -5.93, 7.96, np.nan])
        z = ma.fix_invalid([1.65, 2.64, 2.64, 6.95, np.nan])
        assert_almost_equal(np.asarray(mstats.kendalltau(x,y)),
                            [+0.3333333,0.4969059])
        assert_almost_equal(np.asarray(mstats.kendalltau(x,z)),
                            [-0.5477226,0.2785987])
        #
        x = ma.fix_invalid([0, 0, 0, 0,20,20, 0,60, 0,20,
                            10,10, 0,40, 0,20, 0, 0, 0, 0, 0, np.nan])
        y = ma.fix_invalid([0,80,80,80,10,33,60, 0,67,27,
                            25,80,80,80,80,80,80, 0,10,45, np.nan, 0])
        result = mstats.kendalltau(x,y)
        assert_almost_equal(np.asarray(result), [-0.1585188, 0.4128009])

        # test for namedtuple attributes
        res = mstats.kendalltau(x, y)
        attributes = ('correlation', 'pvalue')
        check_named_results(res, attributes, ma=True)

def hdquantiles_sd(data, prob=list([0.25, 0.5, 0.75]), axis=None):
    """Computes the standard error of the Harrell-Davis quantile estimates by jackknife.


Parameters
----------
    data: ndarray
        Data array.
    prob: sequence
        Sequence of quantiles to compute.
    axis : int
        Axis along which to compute the quantiles. If None, use a flattened array.

Notes
-----
    The function is restricted to 2D arrays.

    """

    def _hdsd_1D(data, prob):
        "Computes the std error for 1D arrays."
        xsorted = np.sort(data.compressed())
        n = len(xsorted)
        # .........
        hdsd = np.empty(len(prob), float_)
        if n < 2:
            hdsd.flat = np.nan
        # .........
        vv = np.arange(n) / float(n - 1)
        betacdf = beta.cdf
        #
        for (i, p) in enumerate(prob):
            _w = betacdf(vv, (n + 1) * p, (n + 1) * (1 - p))
            w = _w[1:] - _w[:-1]
            mx_ = np.fromiter(
                [np.dot(w, xsorted[np.r_[range(0, k), range(k + 1, n)].astype(int_)]) for k in range(n)], dtype=float_
            )
            mx_var = np.array(mx_.var(), copy=False, ndmin=1) * n / float(n - 1)
            hdsd[i] = float(n - 1) * np.sqrt(np.diag(mx_var).diagonal() / float(n))
        return hdsd

    # Initialization & checks ---------
    data = ma.array(data, copy=False, dtype=float_)
    p = np.array(prob, copy=False, ndmin=1)
    # Computes quantiles along axis (or globally)
    if axis is None:
        result = _hdsd_1D(data, p)
    else:
        assert data.ndim <= 2, "Array should be 2D at most !"
        result = ma.apply_along_axis(_hdsd_1D, axis, data, p)
    #
    return ma.fix_invalid(result, copy=False).ravel()

    def test_kstwosamp(self):
        x = [
            [nan, nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1],
            [4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3],
            [3, 2, 5, 6, 18, 4, 9, 1, 1, nan, 1, 1, nan],
            [nan, 6, 11, 4, 17, nan, 6, 1, 1, 2, 5, 1, 1],
        ]
        x = ma.fix_invalid(x).T
        (winter, spring, summer, fall) = x.T

        assert_almost_equal(np.round(mstats.ks_twosamp(winter, spring), 4), (0.1818, 0.9892))
        assert_almost_equal(np.round(mstats.ks_twosamp(winter, spring, "g"), 4), (0.1469, 0.7734))
        assert_almost_equal(np.round(mstats.ks_twosamp(winter, spring, "l"), 4), (0.1818, 0.6744))

def chi2(data=None, back_calc=None, errors=None, params=None):
    """Function to calculate the chi-squared value."""

    # Calculate the chi-squared statistic.
    if raise_warnings:
        try:
            t_chi2 = sum((1.0 / errors * (data - back_calc))**2)
        except RuntimeWarning:
            # Handle if algorithm takes wrong step.
            #print "Oppps. np=%i, sim=%i, R2=%3.2f, I0=%3.2f" % (np_i, sim_j, params[0], params[1])
            t_chi2 =  1e100
    else:
        t_chi2 = sum((1.0 / errors * (data - back_calc))**2)
        if 0:
            fix_invalid(t_chi2, copy=False, fill_value=1e100)
            t_chi2 = nan_to_num( t_chi2 )
        if not isfinite(t_chi2):
            t_chi2_2 = nan_to_num( t_chi2 )
            #print "Oppps. np=%i, sim=%i, R2=%3.2f, I0=%3.2f %s %s" % (np_i, sim_j, params[0], params[1], t_chi2, t_chi2_2)
            t_chi2 = t_chi2_2

    return t_chi2

 def test_kstwosamp(self):
     "Tests the Kolmogorov-Smirnov 2 samples test"
     x = [[nan,nan,  4,  2, 16, 26,  5,  1,  5,  1,  2,  3,  1],
          [  4,  3,  5,  3,  2,  7,  3,  1,  1,  2,  3,  5,  3],
          [  3,  2,  5,  6, 18,  4,  9,  1,  1,nan,  1,  1,nan],
          [nan,  6, 11,  4, 17,nan,  6,  1,  1,  2,  5,  1,  1]]
     x = ma.fix_invalid(x).T
     (winter,spring,summer,fall) = x.T
     #
     assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring),4),
                         (0.1818,0.9892))
     assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring,'g'),4),
                         (0.1469,0.7734))
     assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring,'l'),4),
                         (0.1818,0.6744))

    def test_spearmanr(self):
        # Tests some computations of Spearman's rho
        (x, y) = ([5.05,6.75,3.21,2.66],[1.65,2.64,2.64,6.95])
        assert_almost_equal(mstats.spearmanr(x,y)[0], -0.6324555)
        (x, y) = ([5.05,6.75,3.21,2.66,np.nan],[1.65,2.64,2.64,6.95,np.nan])
        (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y))
        assert_almost_equal(mstats.spearmanr(x,y)[0], -0.6324555)

        x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1,
              1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7]
        y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6,
              0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4]
        assert_almost_equal(mstats.spearmanr(x,y)[0], 0.6887299)
        x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1,
              1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7, np.nan]
        y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6,
              0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4, np.nan]
        (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y))
        assert_almost_equal(mstats.spearmanr(x,y)[0], 0.6887299)

        # test for namedtuple attributes
        res = mstats.spearmanr(x, y)
        attributes = ('correlation', 'pvalue')
        check_named_results(res, attributes, ma=True)

def densitystep(S, T, P):
    """
    """
    assert S.shape == T.shape
    assert S.shape == P.shape
    try:
        import gsw
        rho0 = gsw.pot_rho_t_exact(S, T, P, 0)
        assert S.ndim == 1, "Not able to densitystep an array ndim > 1"
        ds = ma.concatenate([ma.masked_all(1),
                np.sign(np.diff(P))*np.diff(rho0)])
        return ma.fix_invalid(ds)

    except ImportError:
        print("Package gsw is required and is not available.")

 def test_avf_masked(self):
     presidents = ma.fix_invalid(self.presidents)
     # periodogram : True
     avfp = avf(presidents)
     assert_almost_equal(avfp[:21].round(2),
                         [241.74,185.75,159.63,116.92, 95.91, 60.36, 45.69,
                           34.97, 31.74, 10.91,  7.48,  1.32, 11.70,  7.71,
                           13.57,  4.16, -1.05, -9.76,-11.24,-15.67,-12.32])
     # preiodogram : False
     avfp = avf(presidents,0)
     pz = presidents.anom()
     mz = (~pz.mask).astype(int)
     assert_almost_equal(avfp[:21],
                         np.r_[[pz.var()],
                               [(pz[k:]*pz[:-k]).sum()/(mz[k:]*mz[:-k]).sum()
                                for k in range(1,21)]])

def estimate_anomaly(features, params, method='produtorium'):
    """ Estimate probability from PDF defined by params

        The output is the natural logarithm of the estimated probability.

        params are the parameters that define the PDF for each feature
          in features. This function estimate the combined probability of
          each row in features as the produtorium between the probabilities
          of the different features on the same row.

        ATENTION!! I should think more about what would I like from this
          function. What should happens in case of a masked feature? And
          if all features for one measurement are masked? Right now it
          simply don't add for the estimate, so that all features masked
          would lead to an expectation of 100% it's good.
    """
    assert hasattr(params, 'keys')
    assert hasattr(features, 'keys')
    for k in params.keys():
        assert k in features.keys(), "features doesn't have: %s" % k

    prob = ma.masked_all(len(features[features.keys()[0]]))

    for t in params.keys():
        param = params[t]['param']
        valid = ~ma.fix_invalid(features[t]).mask

        tmp = exponweib.sf(np.asanyarray(features[t]),
                *param[:-2], loc=param[-2], scale=param[-1])
        # Arbitrary solution. No value can have a probability of 0.
        tmp[tmp == 0] = 1e-15
        p = ma.log(tmp)

        # Update prob if new value is valid and prob is masked
        ind = prob.mask & valid
        prob[ind] = p[ind]

        # If both are valid, operate as choosed method.
        ind = ~prob.mask & valid
        if method == 'produtorium':
            prob[ind] = prob[ind] + p[ind]
        elif method == 'min':
            prob[ind] = min(prob[ind], p[ind])
        else:
            return

    return prob

 def test_friedmanchisq(self):
     # No missing values
     args = ([9.0,9.5,5.0,7.5,9.5,7.5,8.0,7.0,8.5,6.0],
             [7.0,6.5,7.0,7.5,5.0,8.0,6.0,6.5,7.0,7.0],
             [6.0,8.0,4.0,6.0,7.0,6.5,6.0,4.0,6.5,3.0])
     result = mstats.friedmanchisquare(*args)
     assert_almost_equal(result[0], 10.4737, 4)
     assert_almost_equal(result[1], 0.005317, 6)
     # Missing values
     x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1],
          [4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3],
          [3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan],
          [nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]]
     x = ma.fix_invalid(x)
     result = mstats.friedmanchisquare(*x)
     assert_almost_equal(result[0], 2.0156, 4)
     assert_almost_equal(result[1], 0.5692, 4)