def compute_nmi(self, z1, z2):
        # compute normalized mutual information
        
        n = np.size(z1)
        k1 = np.size(np.unique(z1))
        k2 = np.size(np.unique(z2))
        
        nk1 = np.zeros((k1,1))
        nk2 = np.zeros((k2,1))

        for kk in range(k1):
            nk1[kk] = np.sum(z1==kk)
        for kk in range(k2):
            nk2[kk] = np.sum(z2==kk)
            
        pk1 = nk1/float(np.sum(nk1))
        pk2 = nk2/float(np.sum(nk2))
        
        nk12 = np.zeros((k1,k2))
        for ii in range(k1):
            for jj in range(k2):
                nk12[ii,jj] = np.sum((z1==ii)*(z2==jj))
        pk12 = nk12/float(n)        
        
        Hx = -np.sum(pk1 * np.log(pk1 + np.finfo(float).eps))
        Hy = -np.sum(pk2 * np.log(pk2 + np.finfo(float).eps))
        
        Hxy = -np.sum(pk12 * np.log(pk12 + np.finfo(float).eps))
        
        MI = Hx + Hy - Hxy;
        nmi = MI/float(0.5*(Hx+Hy))
        
        return nmi

def recognition(obs, oracle, order=1, smooth=False):
    hmm_tensor = extract_hmm_tensor(oracle, max_lrs=order, smooth=smooth)

    cluster_means = np.array([np.median(oracle.f_array.data[np.array(c), :].T, axis=1)
                              for c in oracle.latent])

    cluster_means += np.finfo('float').eps
    cluster_means = (cluster_means.T / np.sum(cluster_means, axis=1)).T

    a = hmm_tensor[-1]
    a += np.finfo('float').eps
    a += 1.0
    divider = np.sum(a, axis=1)
    a = np.divide(a.T, divider).T
    log_a = np.log(a)
    # hist = np.array([len(c) for c in oracle.latent])/float(oracle.n_states-1)

    v = np.zeros((len(obs), len(oracle.latent)))
    p = np.zeros(v.shape)
    v[0] = np.log(np.dot(cluster_means, obs[0])) + np.log(1.0/len(oracle.latent))
    # v[0] = np.log(np.dot(cluster_means, obs[0])) + np.log(hist)
    # p[0] = np.arange(len(oracle.latent))
    for t in range(1, len(obs)):
        s = v[t-1]+log_a.T
        v[t] = np.max(s, axis=1)+np.log(np.dot(cluster_means, obs[t]))
        p[t-1] = np.argmax(s, axis=1)

    return v, p

def _gpinv(p, k, sigma):
    """Inverse Generalized Pareto distribution function"""
    x = np.full_like(p, np.nan)
    if sigma <= 0:
        return x
    ok = (p > 0) & (p < 1)
    if np.all(ok):
        if np.abs(k) < np.finfo(float).eps:
            x = - np.log1p(-p)
        else:
            x = np.expm1(-k * np.log1p(-p)) / k
        x *= sigma
    else:
        if np.abs(k) < np.finfo(float).eps:
            x[ok] = - np.log1p(-p[ok])
        else:
            x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
        x *= sigma
        x[p == 0] = 0
        if k >= 0:
            x[p == 1] = np.inf
        else:
            x[p == 1] = - sigma / k

    return x

def quatAverage(q_in, qsym):
    """
    """
    assert q_in.ndim == 2, 'input must be 2-s hstacked quats'

    # renormalize
    q_in = unitVector(q_in)

    # check to see num of quats is > 1
    if q_in.shape[1] < 3:
        if q_in.shape[1] == 1:
            q_bar = q_in
        else:
            ma, mq = misorientation(q_in[:, 0].reshape(4, 1),
                                    q_in[:, 1].reshape(4, 1), (qsym,))
            q_bar = quatProduct(q_in[:, 0].reshape(4, 1),
                                quatOfExpMap(0.5*ma*unitVector(mq[1:].reshape(3, 1))))
    else:
        # use first quat as initial guess
        phi = 2. * arccos(q_in[0, 0])
        if phi <= finfo(float).eps:
            x0 = zeros(3)
        else:
            n = unitVector(q_in[1:, 0].reshape(3, 1))
            x0 = phi*n.flatten()
        results = leastsq(quatAverage_obj, x0, args=(q_in, qsym))
        phi = sqrt(sum(results[0]*results[0]))
        if phi <= finfo(float).eps:
            q_bar = c_[1., 0., 0., 0.].T
        else:
            n     = results[0] / phi
            q_bar = hstack([cos(0.5*phi), sin(0.5*phi)*n]).reshape(4, 1)
    return q_bar

 def test_nan(self):
     # Test that nan is 'far' from small, tiny, inf, max and min
     for dt in [np.float32, np.float64]:
         if dt == np.float32:
             maxulp = 1e6
         else:
             maxulp = 1e12
         inf = np.array([np.inf]).astype(dt)
         nan = np.array([np.nan]).astype(dt)
         big = np.array([np.finfo(dt).max])
         tiny = np.array([np.finfo(dt).tiny])
         zero = np.array([np.PZERO]).astype(dt)
         nzero = np.array([np.NZERO]).astype(dt)
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, inf,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, big,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, tiny,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, zero,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, nzero,
                                                            maxulp=maxulp))

def test_ulp():
    assert_equal(ulp(), np.finfo(np.float64).eps)
    assert_equal(ulp(1.0), np.finfo(np.float64).eps)
    assert_equal(ulp(np.float32(1.0)), np.finfo(np.float32).eps)
    assert_equal(ulp(np.float32(1.999)), np.finfo(np.float32).eps)
    # Integers always return 1
    assert_equal(ulp(1), 1)
    assert_equal(ulp(2**63-1), 1)
    # negative / positive same
    assert_equal(ulp(-1), 1)
    assert_equal(ulp(7.999), ulp(4.0))
    assert_equal(ulp(-7.999), ulp(4.0))
    assert_equal(ulp(np.float64(2**54-2)), 2)
    assert_equal(ulp(np.float64(2**54)), 4)
    assert_equal(ulp(np.float64(2**54)), 4)
    # Infs, NaNs return nan
    assert_true(np.isnan(ulp(np.inf)))
    assert_true(np.isnan(ulp(-np.inf)))
    assert_true(np.isnan(ulp(np.nan)))
    # 0 gives subnormal smallest
    subn64 = np.float64(2**(-1022-52))
    subn32 = np.float32(2**(-126-23))
    assert_equal(ulp(0.0), subn64)
    assert_equal(ulp(np.float64(0)), subn64)
    assert_equal(ulp(np.float32(0)), subn32)
    # as do multiples of subnormal smallest
    assert_equal(ulp(subn64 * np.float64(2**52)), subn64)
    assert_equal(ulp(subn64 * np.float64(2**53)), subn64*2)
    assert_equal(ulp(subn32 * np.float32(2**23)), subn32)
    assert_equal(ulp(subn32 * np.float32(2**24)), subn32*2)

    def _hessian_main(self, params):
        params_infl = params[:self.k_inflate]
        params_main = params[self.k_inflate:]

        y = self.endog
        w = self.model_infl.predict(params_infl)
        w = np.clip(w, np.finfo(float).eps, 1 - np.finfo(float).eps)
        score = self.score(params)
        zero_idx = np.nonzero(y == 0)[0]
        nonzero_idx = np.nonzero(y)[0]

        mu = self.model_main.predict(params_main)

        hess_arr = np.zeros((self.k_exog, self.k_exog))

        coeff = (1 + w[zero_idx] * (np.exp(mu[zero_idx]) - 1))

        #d2l/dp2
        for i in range(self.k_exog):
            for j in range(i, -1, -1):
                hess_arr[i, j] = ((
                    self.exog[zero_idx, i] * self.exog[zero_idx, j] *
                    mu[zero_idx] * (w[zero_idx] - 1) * (1 / coeff -
                    w[zero_idx] * mu[zero_idx] * np.exp(mu[zero_idx]) /
                    coeff**2)).sum() - (mu[nonzero_idx] * self.exog[nonzero_idx, i] *
                    self.exog[nonzero_idx, j]).sum())

        return hess_arr

def system_diagnostic():
    """ return various helpful/informative information about the
    current system. For instance versions of python & available packages.

    Mostly undocumented function...
    """

    # There is probably a more clever way to do the following via introspection?

    import platform
    import os
    import poppy
    import numpy
    from .version import version
    try:
        import ttk
        ttk_version = ttk.__version__
    except ImportError:
        ttk_version = 'not found'

    try:
        import wx
        wx_version = wx.__version__
    except ImportError:
        wx_version = 'not found'

    try:
        import pyfftw
        pyfftw_version = pyfftw.version
    except ImportError:
        pyfftw_version = 'not found'

    try:
        import pysynphot
        pysynphot_version = pysynphot.__version__
    except ImportError:
        pysynphot_version = 'not found'


    try:
        import astropy
        astropy_version = astropy.__version__
    except ImportError:
        astropy_version = 'not found'

    result = DIAGNOSTIC_REPORT.format(
        os=platform.platform(),
        numpy=numpy.__version__,
        python=sys.version.replace("\n", " "),
        poppy=poppy.__version__,
        webbpsf=version,
        tkinter=ttk_version,
        wxpython=wx_version,
        pyfftw=pyfftw_version,
        pysyn=pysynphot_version,
        astropy=astropy_version,
        finfo_float=numpy.finfo(numpy.float),
        finfo_complex=numpy.finfo(numpy.complex),
    )
    return result

def get_matrix_2d_ragged(workspace, distribution, histogram2D=False, transpose=False):
    num_hist = workspace.getNumberHistograms()
    delta = numpy.finfo(numpy.float64).max
    min_value = numpy.finfo(numpy.float64).max
    max_value = numpy.finfo(numpy.float64).min
    for i in range(num_hist):
        xtmp = workspace.readX(i)
        if workspace.isHistogramData():
            #input x is edges
            xtmp = mantid.plots.helperfunctions.points_from_boundaries(xtmp)
        else:
            #input x is centers
            pass
        min_value = min(min_value, xtmp.min())
        max_value = max(max_value, xtmp.max())
        diff = xtmp[1:] - xtmp[:-1]
        delta = min(delta, diff.min())
    num_edges = int(numpy.ceil((max_value - min_value)/delta)) + 1
    x_centers = numpy.linspace(min_value, max_value, num=num_edges)
    y = mantid.plots.helperfunctions.boundaries_from_points(workspace.getAxis(1).extractValues())
    z = numpy.empty([num_hist, num_edges], dtype=numpy.float64)
    for i in range(num_hist):
        centers, ztmp, _, _ = mantid.plots.helperfunctions.get_spectrum(workspace, i, distribution=distribution, withDy=False, withDx=False)
        f = interp1d(centers, ztmp, kind='nearest', bounds_error=False, fill_value=numpy.nan)
        z[i] = f(x_centers)
    if histogram2D:
        x = mantid.plots.helperfunctions.boundaries_from_points(x_centers)
    else:
        x = x_centers
    if transpose:
        return y.T,x.T,z.T
    else:
        return x,y,z

    def __init__(self):
        self.gravity = 9.8
        self.masscart = 1.0
        self.masspole = 0.1
        self.total_mass = (self.masspole + self.masscart)
        self.length = 0.5 # actually half the pole's length
        self.polemass_length = (self.masspole * self.length)
        self.force_mag = 10.0
        self.tau = 0.02  # seconds between state updates

        # Angle at which to fail the episode
        self.theta_threshold_radians = 12 * 2 * math.pi / 360
        self.x_threshold = 2.4

        # Angle limit set to 2 * theta_threshold_radians so failing observation is still within bounds
        high = np.array([
            self.x_threshold * 2,
            np.finfo(np.float32).max,
            self.theta_threshold_radians * 2,
            np.finfo(np.float32).max])

        self.action_space = spaces.Discrete(2)
        self.observation_space = spaces.Box(-high, high)

        self.seed()
        self.viewer = None
        self.state = None

        self.steps_beyond_done = None

 def almost_equal(self, other,
                  rtol=2**4 * np.finfo(np.zeros(1).dtype).eps,
                  atol=0.):
     # assert atol == 0., 'Not supported'
     assert type(other) == type(self)
     return self._impl.almost_equal(other._impl,
                                    rtol if rtol is not None else 2**4 *np.finfo(np.zeros(1.).dtype).eps)

    def get_cost(self, weights, data, display=False):
        # print 'getting cost...'
        N = float(len(data))
        reg = (self.lmbda / (2.0 * N)) * np.sum(weights ** 2)
        # reg = (self.lmbda / self.N) * np.sum(np.abs(weights))
        # self.set_network_weights(weights)
        layers = self.convert_weights_to_layers(weights)
        cost = 0.0
        for d, l in data[:]:
            z = d
            for idx, layer in enumerate(layers):
                if idx == len(layers) - 1:
                    # this is a output layer
                    prediction = layer.get_z(z)
                    prediction[prediction >= 1.0] = 1.0 - np.finfo(float).eps  # to avoid nan showing up
                    prediction[prediction <= 0.0] = 0.0 + np.finfo(float).eps
                    l1p = -l * np.log(prediction)
                    l2p = -(1.0 - l) * np.log((1.0 - prediction))
                    lcost = np.sum(l1p + l2p)
                    cost += lcost * (1.0 / float(N))

                else:

                    z = layer.get_z(z)
        if display:
            print 'cost', cost + reg
        return cost + reg

  def _testZeroDensity(self, alpha):
    """Zero isn't in the support of the gamma distribution.

    But quantized floating point math has its limits.
    TODO(bjp): Implement log-gamma sampler for small-shape distributions.

    Args:
      alpha: float shape value to test
    """
    try:
      from scipy import stats  # pylint: disable=g-import-not-at-top
    except ImportError as e:
      tf_logging.warn("Cannot test zero density proportions: %s" % e)
      return
    allowable_zeros = {
        dtypes.float16: stats.gamma(alpha).cdf(np.finfo(np.float16).tiny),
        dtypes.float32: stats.gamma(alpha).cdf(np.finfo(np.float32).tiny),
        dtypes.float64: stats.gamma(alpha).cdf(np.finfo(np.float64).tiny)
    }
    failures = []
    for use_gpu in [False, True]:
      for dt in dtypes.float16, dtypes.float32, dtypes.float64:
        sampler = self._Sampler(
            10000, alpha, 1.0, dt, use_gpu=use_gpu, seed=12345)
        x = sampler()
        allowable = allowable_zeros[dt] * x.size
        allowable = allowable * 2 if allowable < 10 else allowable * 1.05
        if np.sum(x <= 0) > allowable:
          failures += [(use_gpu, dt)]
      self.assertEqual([], failures)

    def test_exp(self):
        source = np.array([0, 0, 0, 1.0])
        result = quaternion.exp(source)
        expected = np.array([0, 0, 0, np.exp(1)])
        np.testing.assert_almost_equal(result, expected)

        source = quaternion.create_from_eulers([np.pi, 0, 0])
        result = quaternion.exp(source)
        expected = np.array([0.84147098, 0, 0, 0.54030231])
        np.testing.assert_almost_equal(result, expected)

        # Tests from the boost::math::quaternion
        source = np.array([4 * np.arctan(1), 0, 0, 0])
        result = quaternion.exp(source) + [0, 0, 0, 1.0]
        result = np.linalg.norm(result)
        expected = 2 * np.finfo(result.dtype).eps
        np.testing.assert_almost_equal(result, expected)

        source = np.array([0, 4 * np.arctan(1), 0, 0])
        result = quaternion.exp(source) + [0, 0, 0, 1.0]
        result = np.linalg.norm(result)
        expected = 2 * np.finfo(result.dtype).eps
        np.testing.assert_almost_equal(result, expected)

        source = np.array([0, 0, 4 * np.arctan(1), 0])
        result = quaternion.exp(source) + [0, 0, 0, 1.0]
        result = np.linalg.norm(result)
        expected = 2 * np.finfo(result.dtype).eps
        np.testing.assert_almost_equal(result, expected)

def elop(X, Y, op):
    """
    Compute element-wise operation of matrix :param:`X` and matrix :param:`Y`.
    
    :param X: First input matrix.
    :type X: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil, dia or :class:`numpy.matrix`
    :param Y: Second input matrix.
    :type Y: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil, dia or :class:`numpy.matrix`
    :param op: Operation to be performed. 
    :type op: `func` 
    """
    try:
        zp1 = op(0, 1) if sp.isspmatrix(X) else op(1, 0)
        zp2 = op(0, 0) 
        zp = zp1 != 0 or zp2 != 0
    except:
        zp = 0
    if sp.isspmatrix(X) or sp.isspmatrix(Y):
        return _op_spmatrix(X, Y, op) if not zp else _op_matrix(X, Y, op)
    else:
        try:
            X[X == 0] = np.finfo(X.dtype).eps
            Y[Y == 0] = np.finfo(Y.dtype).eps
        except ValueError:
            return op(np.asmatrix(X), np.asmatrix(Y))
        return op(np.asmatrix(X), np.asmatrix(Y))

    def _wave_fit(self, tseries, candidate_start_points, candidate_peak_volumes):
        '''
        This method has the fitting startegy to minimize the cost.
        ''' 
        
        period = self.period
        threshold = self.threshold
        residual_metric = self.residual_metric
        score_func = self.score_func
        fit_audience = self.fit_audience

        curr_score = np.finfo('d').max
        best_score = np.finfo('d').max
        best_params = None

        params = None
        for i in xrange(len(candidate_start_points)):
            sp = candidate_start_points[i]
            pv = candidate_peak_volumes[i]
            
            params = _fit_one(tseries, period, residual_metric, sp, pv, \
                    fit_audience, params, candidate_start_points[0])
            model = phoenix_r_with_period(params, tseries.shape[0])
            num_params = 5 * (i + 1) + 2
            curr_score = self.score_func(model, tseries, num_params, params)

            if (curr_score <= best_score):
                best_params = params
                best_score = curr_score
            else:
                increased_score = (curr_score - best_score) / best_score
                if increased_score > threshold:
                    break

        return best_score, best_params

            def __init__(self, n, n_folds=3, shuffle=False, random_state=None):
                if abs(n - int(n)) >= np.finfo('f').eps:
                    raise ValueError("n must be an integer")
                self.n = int(n)

                if abs(n_folds - int(n_folds)) >= np.finfo('f').eps:
                    raise ValueError("n_folds must be an integer")
                self.n_folds = n_folds = int(n_folds)

                if n_folds <= 1:
                    raise ValueError(
                        "k-fold cross validation requires at least one"
                        " train / test split by setting n_folds=2 or more,"
                        " got n_folds={0}.".format(n_folds))
                if n_folds > self.n:
                    raise ValueError(
                        ("Cannot have number of folds n_folds={0} greater"
                         " than the number of samples: {1}.").format(n_folds, n))

                if not isinstance(shuffle, bool):
                    raise TypeError("shuffle must be True or False;"
                                    " got {0}".format(shuffle))
                self.shuffle = shuffle
                self.random_state = random_state

                self.idxs = np.arange(n)
                if shuffle:
                    rng = check_random_state(self.random_state)
                    rng.shuffle(self.idxs)

def fbank(signal, samplerate=16000, winlen=0.025, winstep=0.01,
          nfilt=26, nfft=512, lowfreq=0, highfreq=None, preemph=0.97):
    """Compute Mel-filterbank energy features from an audio signal.
    :param signal: the audio signal from which to compute features. Should be an N*1 array
    :param samplerate: the samplerate of the signal we are working with.
    :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds)
    :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds)
    :param nfilt: the number of filters in the filterbank, default 26.
    :param nfft: the FFT size. Default is 512.
    :param lowfreq: lowest band edge of mel filters. In Hz, default is 0.
    :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2
    :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97.
    :returns: 2 values. The first is a numpy array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. The
        second return value is the energy in each frame (total energy, unwindowed)
    """
    highfreq = highfreq or samplerate / 2
    signal = sigproc.preemphasis(signal, preemph)
    frames = sigproc.framesig(signal, winlen * samplerate, winstep * samplerate)
    pspec = sigproc.powspec(frames, nfft)
    energy = numpy.sum(pspec, 1)  # this stores the total energy in each frame
    energy = numpy.where(energy == 0, numpy.finfo(float).eps, energy)  # if energy is zero, we get problems with log

    fb = get_filterbanks(nfilt, nfft, samplerate, lowfreq, highfreq)
    feat = numpy.dot(pspec, fb.T)  # compute the filterbank energies
    feat = numpy.where(feat == 0, numpy.finfo(float).eps, feat)  # if feat is zero, we get problems with log

    return feat, energy

    def _accumulate_sufficient_statistics(self, stats, obs, framelogprob,
                                          posteriors, fwdlattice, bwdlattice,
                                          params):
        super(GMMHMM, self)._accumulate_sufficient_statistics(
            stats, obs, framelogprob, posteriors, fwdlattice, bwdlattice,
            params)

        for state, g in enumerate(self.gmms):
            lgmm_posteriors = np.log(g.eval(obs)[1] + np.finfo(np.float).eps)
            lgmm_posteriors += np.log(posteriors[:, state][:, np.newaxis]
                                      + np.finfo(np.float).eps)
            gmm_posteriors = np.exp(lgmm_posteriors)
            tmp_gmm = GMM(g.n_components, cvtype=g.cvtype)
            tmp_gmm.n_features = g.n_features
            tmp_gmm.covars = _distribute_covar_matrix_to_match_cvtype(
                                np.eye(g.n_features), g.cvtype, g.n_components)
            norm = tmp_gmm._do_mstep(obs, gmm_posteriors, params)

            if np.any(np.isnan(tmp_gmm.covars)):
                raise ValueError

            stats['norm'][state] += norm
            if 'm' in params:
                stats['means'][state] += tmp_gmm.means * norm[:, np.newaxis]
            if 'c' in params:
                if tmp_gmm.cvtype == 'tied':
                    stats['covars'][state] += tmp_gmm._covars * norm.sum()
                else:
                    cvnorm = np.copy(norm)
                    shape = np.ones(tmp_gmm._covars.ndim)
                    shape[0] = np.shape(tmp_gmm._covars)[0]
                    cvnorm.shape = shape
                    stats['covars'][state] += tmp_gmm._covars * cvnorm

    def test_complex128_fail(self):
        nulp = 5
        x = np.linspace(-20, 20, 50, dtype=np.float64)
        x = 10**x
        x = np.r_[-x, x]
        xi = x + x*1j

        eps = np.finfo(x.dtype).eps
        y = x + x*eps*nulp*2.
        self.assertRaises(AssertionError, assert_array_almost_equal_nulp,
                          xi, x + y*1j, nulp)
        self.assertRaises(AssertionError, assert_array_almost_equal_nulp,
                          xi, y + x*1j, nulp)
        # The test condition needs to be at least a factor of sqrt(2) smaller
        # because the real and imaginary parts both change
        y = x + x*eps*nulp
        self.assertRaises(AssertionError, assert_array_almost_equal_nulp,
                          xi, y + y*1j, nulp)

        epsneg = np.finfo(x.dtype).epsneg
        y = x - x*epsneg*nulp*2.
        self.assertRaises(AssertionError, assert_array_almost_equal_nulp,
                          xi, x + y*1j, nulp)
        self.assertRaises(AssertionError, assert_array_almost_equal_nulp,
                          xi, y + x*1j, nulp)
        y = x - x*epsneg*nulp
        self.assertRaises(AssertionError, assert_array_almost_equal_nulp,
                          xi, y + y*1j, nulp)