def change_pmin_by_innovation(self, x, f):
        Lx, _ = self._gp_innovation_local(x)
        dMdb = Lx
        dVdb = -Lx.dot(Lx.T)
        stoch_changes = dMdb.dot(self.W)
        Mb_new = self.Mb[:, None] + stoch_changes

        Vb_new = self.Vb + dVdb

        Vb_new[np.diag_indices(Vb_new.shape[0])] = np.clip(Vb_new[np.diag_indices(Vb_new.shape[0])], np.finfo(Vb_new.dtype).eps, np.inf)

        Vb_new[np.where((Vb_new < np.finfo(Vb_new.dtype).eps) & (Vb_new > -np.finfo(Vb_new.dtype).eps))] = 0
        try:
            cVb_new = np.linalg.cholesky(Vb_new)
        except np.linalg.LinAlgError:
            try:
                cVb_new = np.linalg.cholesky(Vb_new + 1e-10 * np.eye(Vb_new.shape[0]))
            except np.linalg.LinAlgError:
                try:
                    cVb_new = np.linalg.cholesky(Vb_new + 1e-6 * np.eye(Vb_new.shape[0]))
                except np.linalg.LinAlgError:
                    cVb_new = np.linalg.cholesky(Vb_new + 1e-3 * np.eye(Vb_new.shape[0]))
        f_new = np.dot(cVb_new, self.F.T)
        f_new = f_new[:, :, None]
        Mb_new = Mb_new[:, None, :]
        f_new = Mb_new + f_new
        return self.calc_pmin(f_new)

    def test_eigenvectors(self):
        P = self.bdc.transition_matrix()

        # k==None
        ev = eigvals(P)
        ev = ev[np.argsort(np.abs(ev))[::-1]]
        Dn = np.diag(ev)

        # right eigenvectors
        Rn = eigenvectors(P)
        assert_allclose(np.dot(P,Rn),np.dot(Rn,Dn))
        # left eigenvectors
        Ln = eigenvectors(P, right=False)
        assert_allclose(np.dot(Ln.T,P),np.dot(Dn,Ln.T))
        # orthogonality
        Xn = np.dot(Ln.T, Rn)
        di = np.diag_indices(Xn.shape[0])
        Xn[di] = 0.0
        assert_allclose(Xn,0)

        # k!=None
        Dnk = Dn[:,0:self.k][0:self.k,:]
        # right eigenvectors
        Rn = eigenvectors(P, k=self.k)
        assert_allclose(np.dot(P,Rn),np.dot(Rn,Dnk))
        # left eigenvectors
        Ln = eigenvectors(P, right=False, k=self.k)
        assert_allclose(np.dot(Ln.T,P),np.dot(Dnk,Ln.T))
        # orthogonality
        Xn = np.dot(Ln.T, Rn)
        di = np.diag_indices(self.k)
        Xn[di] = 0.0
        assert_allclose(Xn,0)

def _make_rdm1(mycc, d1, with_frozen=True, ao_repr=False):
    '''dm1[p,q] = <q_alpha^\dagger p_alpha> + <q_beta^\dagger p_beta>

    The convention of 1-pdm is based on McWeeney's book, Eq (5.4.20).
    The contraction between 1-particle Hamiltonian and rdm1 is
    E = einsum('pq,qp', h1, rdm1)
    '''
    doo, dov, dvo, dvv = d1
    nocc, nvir = dov.shape
    nmo = nocc + nvir
    dm1 = numpy.empty((nmo,nmo), dtype=doo.dtype)
    dm1[:nocc,:nocc] = doo + doo.conj().T
    dm1[:nocc,nocc:] = dov + dvo.conj().T
    dm1[nocc:,:nocc] = dm1[:nocc,nocc:].conj().T
    dm1[nocc:,nocc:] = dvv + dvv.conj().T
    dm1[numpy.diag_indices(nocc)] += 2

    if with_frozen and not (mycc.frozen is 0 or mycc.frozen is None):
        nmo = mycc.mo_occ.size
        nocc = numpy.count_nonzero(mycc.mo_occ > 0)
        rdm1 = numpy.zeros((nmo,nmo), dtype=dm1.dtype)
        rdm1[numpy.diag_indices(nocc)] = 2
        moidx = numpy.where(mycc.get_frozen_mask())[0]
        rdm1[moidx[:,None],moidx] = dm1
        dm1 = rdm1

    if ao_repr:
        mo = mycc.mo_coeff
        dm1 = lib.einsum('pi,ij,qj->pq', mo, dm1, mo.conj())
    return dm1

    def predict_wishard_embedding(self, Xnew, kern=None, mean=True, covariance=True):
        """
        Predict the wishard embedding G of the GP. This is the density of the
        input of the GP defined by the probabilistic function mapping f.
        G = J_mean.T*J_mean + output_dim*J_cov.

        :param array-like Xnew: The points at which to evaluate the magnification.
        :param :py:class:`~GPy.kern.Kern` kern: The kernel to use for the magnification.

        Supplying only a part of the learning kernel gives insights into the density
        of the specific kernel part of the input function. E.g. one can see how dense the
        linear part of a kernel is compared to the non-linear part etc.
        """
        if kern is None:
            kern = self.kern

        mu_jac, var_jac = self.predict_jacobian(Xnew, kern, full_cov=False)
        mumuT = np.einsum('iqd,ipd->iqp', mu_jac, mu_jac)
        Sigma = np.zeros(mumuT.shape)
        if var_jac.ndim == 3:
            Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = var_jac.sum(-1)
        else:
            Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = self.output_dim*var_jac
        G = 0.
        if mean:
            G += mumuT
        if covariance:
            G += Sigma
        return G

def _gamma1_intermediates(mycc, t1, t2, l1, l2, eris=None):
    doo, dov, dvo, dvv = gccsd_rdm._gamma1_intermediates(mycc, t1, t2, l1, l2)

    if eris is None: eris = mycc.ao2mo()

    nocc, nvir = t1.shape
    bcei = numpy.asarray(eris.ovvv).conj().transpose(3,2,1,0)
    majk = numpy.asarray(eris.ooov).conj().transpose(2,3,0,1)
    bcjk = numpy.asarray(eris.oovv).conj().transpose(2,3,0,1)

    mo_e = eris.mo_energy
    eia = mo_e[:nocc,None] - mo_e[nocc:]
    d3 = lib.direct_sum('ia+jb+kc->ijkabc', eia, eia, eia)

    t3c =(numpy.einsum('jkae,bcei->ijkabc', t2, bcei)
        - numpy.einsum('imbc,majk->ijkabc', t2, majk))
    t3c = t3c - t3c.transpose(0,1,2,4,3,5) - t3c.transpose(0,1,2,5,4,3)
    t3c = t3c - t3c.transpose(1,0,2,3,4,5) - t3c.transpose(2,1,0,3,4,5)
    t3c /= d3

    t3d = numpy.einsum('ia,bcjk->ijkabc', t1, bcjk)
    t3d += numpy.einsum('ai,jkbc->ijkabc', eris.fock[nocc:,:nocc], t2)
    t3d = t3d - t3d.transpose(0,1,2,4,3,5) - t3d.transpose(0,1,2,5,4,3)
    t3d = t3d - t3d.transpose(1,0,2,3,4,5) - t3d.transpose(2,1,0,3,4,5)
    t3d /= d3

    goo = numpy.einsum('iklabc,jklabc->ij', (t3c+t3d).conj(), t3c) * (1./12)
    gvv = numpy.einsum('ijkacd,ijkbcd->ab', t3c+t3d, t3c.conj()) * (1./12)
    doo[numpy.diag_indices(nocc)] -= goo.diagonal()
    dvv[numpy.diag_indices(nvir)] += gvv.diagonal()
    dvo += numpy.einsum('ijab,ijkabc->ck', t2.conj(), t3c) * (1./4)

    return doo, dov, dvo, dvv

def test_diag_indices():
    di = diag_indices(4)
    a = array([[1, 2, 3, 4],
               [5, 6, 7, 8],
               [9, 10, 11, 12],
               [13, 14, 15, 16]])
    a[di] = 100
    yield (assert_array_equal, a,
           array([[100,   2,   3,   4],
                  [  5, 100,   7,   8],
                  [  9,  10, 100,  12],
                  [ 13,  14,  15, 100]]))

    # Now, we create indices to manipulate a 3-d array:
    d3 = diag_indices(2, 3)

    # And use it to set the diagonal of a zeros array to 1:
    a = zeros((2, 2, 2),int)
    a[d3] = 1
    yield (assert_array_equal, a,
           array([[[1, 0],
                   [0, 0]],

                  [[0, 0],
                   [0, 1]]]) )

def test_naf_layer_full():
    batch_size = 2
    for nb_actions in (1, 3):
        # Construct single model with NAF as the only layer, hence it is fully deterministic
        # since no weights are used, which would be randomly initialized.
        L_flat_input = Input(shape=((nb_actions * nb_actions + nb_actions) // 2,))
        mu_input = Input(shape=(nb_actions,))
        action_input = Input(shape=(nb_actions,))
        x = NAFLayer(nb_actions, mode='full')([L_flat_input, mu_input, action_input])
        model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
        model.compile(loss='mse', optimizer='sgd')
        
        # Create random test data.
        L_flat = np.random.random((batch_size, (nb_actions * nb_actions + nb_actions) // 2)).astype('float32')
        mu = np.random.random((batch_size, nb_actions)).astype('float32')
        action = np.random.random((batch_size, nb_actions)).astype('float32')

        # Perform reference computations in numpy since these are much easier to verify.
        L = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
        LT = np.copy(L)
        for l, l_T, l_flat in zip(L, LT, L_flat):
            l[np.tril_indices(nb_actions)] = l_flat
            l[np.diag_indices(nb_actions)] = np.exp(l[np.diag_indices(nb_actions)])
            l_T[:, :] = l.T
        P = np.array([np.dot(l, l_T) for l, l_T in zip(L, LT)]).astype('float32')
        A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
        A_ref *= -.5

        # Finally, compute the output of the net, which should be identical to the previously
        # computed reference.
        A_net = model.predict([L_flat, mu, action]).flatten()
        assert_allclose(A_net, A_ref, rtol=1e-5)

def _laplacian_dense(csgraph, normed='geometric', symmetrize=True, scaling_epps=0., renormalization_exponent=1, return_diag=False, return_lapsym = False):
    n_nodes = csgraph.shape[0]
    if symmetrize:
        lap = (csgraph + csgraph.T)/2.
    else:
        lap = csgraph.copy()
    degrees = np.asarray(lap.sum(axis=1)).squeeze()
    di = np.diag_indices( lap.shape[0] )  # diagonal indices

    if normed == 'symmetricnormalized':
        w = np.sqrt(degrees)
        w_zeros = (w == 0)
        w[w_zeros] = 1
        lap /= w
        lap /= w[:, np.newaxis]
        di = np.diag_indices( lap.shape[0] )
        lap[di] -= (1 - w_zeros).astype(lap.dtype)
    if normed == 'geometric':
        w = degrees.copy()     # normalize once symmetrically by d
        w_zeros = (w == 0)
        w[w_zeros] = 1
        lap /= w
        lap /= w[:, np.newaxis]
        w = np.asarray(lap.sum(axis=1)).squeeze() #normalize again asymmetricall
        if return_lapsym:
            lapsym = lap.copy()
        lap /= w[:, np.newaxis]
        lap[di] -= (1 - w_zeros).astype(lap.dtype)
    if normed == 'renormalized':
        w = degrees**renormalization_exponent;
        # same as 'geometric' from here on
        w_zeros = (w == 0)
        w[w_zeros] = 1
        lap /= w
        lap /= w[:, np.newaxis]
        w = np.asarray(lap.sum(axis=1)).squeeze() #normalize again asymmetricall
        if return_lapsym:
            lapsym = lap.copy()
        lap /= w[:, np.newaxis]
        lap[di] -= (1 - w_zeros).astype(lap.dtype)
    if normed == 'unnormalized':
        dum = lap[di]-degrees[np.newaxis,:]
        lap[di] = dum[0,:]
    if normed == 'randomwalk':
        lap /= degrees[:,np.newaxis]
        lap -= np.eye(lap.shape[0])

    if scaling_epps > 0.:
        lap *= 4/(scaling_epps**2)

    if return_diag:
        diag = np.array( lap[di] )
        if return_lapsym:
            return lap, diag, lapsym, w
        else: 
            return lap, diag
    elif return_lapsym:
        return lap, lapsym, w
    else:
        return lap

    def kNN_graph(self, k, metric, mutual=False):
#        self.latex = []
        nn = NearestNeighbors(k, algorithm="brute", metric=metric, n_jobs=-1).fit(self.X)
        UAM = nn.kneighbors_graph(self.X).toarray() #unweighted adjacency matrix
        m = UAM.shape[0]
        self.W = np.zeros((m, m)) #(weighted) adjecancy matrix
        self.D = np.zeros((m, m)) #degree matrix
        if mutual == False:
            if self.full_calculated:
                indices = np.where(UAM == 1)
                self.W[indices] = self.full_W[indices]
                self.D[np.diag_indices(m)] = np.sum(self.W, 1)
            else:
                for i in range(m):
                    for j in range(m):
                        if UAM[i,j] == 1:
                            sim = self.s(self.X[i], self.X[j], self.d)
                            self.W[i,j] = sim
                            self.D[i,i] += sim
        else:
            if self.full_calculated:
                indices = np.where(np.logical_and(UAM == 1, UAM.T == 1).astype(int) == 1)
                self.W[indices] = self.full_W[indices]
                self.D[np.diag_indices(m)] = np.sum(self.W != 0, 1)
            else:
                for i in range(m):
                    for j in range(m):
                        if UAM[i,j] == 1 and UAM[j,i] == 1:
                            sim = self.s(self.X[i], self.X[j], self.d)
                            self.W[i,j] = sim
                            self.D[i,i] += sim
        self.W = np.nan_to_num(self.W)
        self.graph = "kNN graph, k = " + str(k) + ", mutual:" + str(mutual)

    def predict(self, X_test, **kwargs):

        # Normalize input data to 0 mean and unit std
        X_ = (X_test - self.X_mean) /  self.X_std

        # Get features from the net

        layers = lasagne.layers.get_all_layers(self.network)
        theta = lasagne.layers.get_output(layers[:-1], X_)[-1].eval()

        # Marginalise predictions over hyperparameters of the BLR
        mu = np.zeros([len(self.models), X_test.shape[0]])
        var = np.zeros([len(self.models), X_test.shape[0]])

        for i, m in enumerate(self.models):
            mu[i], var[i] = m.predict(theta)

        # See the algorithm runtime prediction paper by Hutter et al
        # for the derivation of the total variance
        m = np.array([[mu.mean()]])
        v = np.mean(mu ** 2 + var) - m ** 2

        # Clip negative variances and set them to the smallest
        # positive float value
        if v.shape[0] == 1:
            v = np.clip(v, np.finfo(v.dtype).eps, np.inf)
        else:
            v[np.diag_indices(v.shape[0])] = \
                    np.clip(v[np.diag_indices(v.shape[0])],
                            np.finfo(v.dtype).eps, np.inf)
            v[np.where((v < np.finfo(v.dtype).eps) & (v > -np.finfo(v.dtype).eps))] = 0


        return m, v

def test_calculate2():
    # Two mutations in one cluster, two in second
    c = np.zeros((4,2))
    c[0:2,0] = 1
    c[2:4,1] = 1
    c = np.dot(c,c.T)

    # Identical
    assert round(calculate2(c,c), 2) == 1.00

    #Inverted
    c2 = np.abs(c-1)
    c2[np.diag_indices(4)] = 1
    assert round(calculate2(c,c2), 2) == -0.68

    # Differences first 3 SSMs in first cluster, 4th ssm in second cluster
    c3 = np.zeros((4,2))
    c3[0:3,0] = 1
    c3[3:4,1] = 1
    c3 = np.dot(c3,c3.T)
    assert round(calculate2(c,c3), 2) == -0.20

    # Metric doesn't count the diagnonal
    c4 = c+0
    c4[np.diag_indices(4)] = 0
    assert round(calculate2(c,c4), 2) == 0.74

    c2[np.diag_indices(4)] = 0
    assert round(calculate2(c,c2), 2) == -0.23

def compute_genetic_distance(X, **kwargs):
    """Given genotype matrix X, returns pairwise genetic distance between individuals
    using the estimator described in Theorem 1.

    Args:
        X: n * p matrix of 0/1/2/nan, n is #individuals, p is #SNPs
    """

    n, p = X.shape
    missing = np.isnan(X)
    col_sums = np.nansum(X, axis=0)
    col_counts = np.sum(~missing, axis=0)
    mu_hat = col_sums / 2. / col_counts     # p dimensional

    eta0_hat = np.nansum(X**2 - X, axis=0) / 2. / col_counts - mu_hat**2

    X_tmp = X/2.
    X_tmp[missing] = 0
    non_missing = np.array(~missing, dtype=float)

    X_shifted = X_tmp - mu_hat
    gdm_squared = 2. * np.mean(eta0_hat) - 2. * np.dot(X_shifted, X_shifted.T) / np.dot(non_missing, non_missing.T)
    gdm_squared[np.diag_indices(n)] = 0.

    if len(gdm_squared[gdm_squared < 0]) > 0:
        # shift all entries by the smallest amount that makes them non-negative
        shift = - np.min(gdm_squared[gdm_squared < 0])
        gdm_squared += shift
        gdm_squared[np.diag_indices(n)] = 0.

    gdm = np.sqrt(np.maximum(gdm_squared, 0.))

    return gdm

 def derivative_of_marginalLikelihood(self, hyp):
     """
     This method calculates the derivative of marginal likelihood.
     
     You may refer to this code to see what methods in numpy is useful
     
     while you are implementing other functions.
     
     Do not modify this method.
     """
     trainingData = self.trainingData
     trainingLabels = self.trainingLabels
     c = len(self.legalLabels)
     [n,d] = self.trainingShape
     hyp = np.reshape(hyp, self.hypSize)
     
     [mode,_] = self.findMode(trainingData, trainingLabels, hyp)
     [t,_] = self.trainingLabels2t(trainingLabels)
       
     Ks = self.calculateCovariance(trainingData, hyp)
     [_,[E, M, R, b, totpi, K],_] = self.calculateIntermediateValues(t, mode, Ks)
       
     MRE = np.linalg.solve(M,R.T.dot(E))
     MMRE = np.linalg.solve(M.T,MRE)
     KWinvinv = E-E.dot(R.dot(MMRE))
       
     KinvWinv = K-K.dot(KWinvinv.dot(K))
     partitioned_KinvWinv = np.transpose(np.array(np.split(np.array(np.split(KinvWinv, c)),c,2)),[2,3,1,0])
       
     s2 = np.zeros([n,c])
     for i in range(n):
         pi_n = softmax(np.reshape(mode,[c,n])[:,i:i+1].T).T
         pipj = pi_n.dot(pi_n.T)
         pi_3d = np.zeros([c,c,c])
         pi_3d[np.diag_indices(c,3)] = pi_n.ravel()
         pipjpk = np.tensordot(pi_n,np.reshape(pipj,(1,c,c)),(1,0))
         pipj_3d = np.zeros([c,c,c])
         pipj_3d[np.diag_indices(c)] = pipj
         W_3d = pi_3d + 2 * pipjpk - pipj_3d - np.transpose(pipj_3d,[2,1,0]) - np.transpose(pipj_3d,[1,2,0])
         s2[i,:] = -0.5*np.trace(partitioned_KinvWinv[i,i].dot(W_3d))
           
     b_rs = np.reshape(b, [c,n])
     dZ = np.zeros(hyp.shape)
     for j in range(2):
         cs = []
         zeroCs = [np.zeros([n,n]) for i in range(c)]
         for i in range(c):
             C = self.covARD(hyp[i,:],trainingData,None,j)
             dZ[i,j] = 0.5*b_rs[i,:].T.dot(C.dot(b_rs[i,:]))
             zeroCs[i] = C
             cs.append(self.block_diag(zeroCs))
             zeroCs[i] = np.zeros([n,n])
             
         for i in range(c):
             dd = cs[i].dot(t-totpi)
             s3 = dd - K.dot(KWinvinv.dot(dd))
             dZ[i,j] +=  - 0.5 * np.trace(KWinvinv.dot(cs[i])) + s2.T.ravel().dot(s3) # 
               
     return -dZ.ravel()

def main():
    order = 4
    cov, inv_cov, inv_cov_est, t_direct, t_patches = get_mats(order)
    header = '=== inv_scale_length:{:.1f}: ==='
    header = header.format(inv_scale_length)
    print(header, file=logf)
    print('Time to invert directly:   {:.4f} s'.format(t_direct), file=logf)
    maybe_zero = inv_cov_est - inv_cov_est.T
    for i in range(npix):
        for j in range(npix):
            print(maybe_zero[i][j])
    return 0
    #patch_identity.shape = shape

    t2 = time()
    print('Time to invert by patches: {:.4f} s'.format(t2-t1), file=logf)
    I = np.dot(inv_cov, cov)
    I_est = np.dot(inv_cov_est, cov)

    t3 = time()
    print('Time to calculate checks:  {:.4f} s'.format(t3-t2), file=logf)

    I_diag = I[np.diag_indices(I.shape[0])]
    I_diag_dev = np.max(np.abs(I_diag-1.))
    I[np.diag_indices(I.shape[0])] = 0.
    I_off_diag_dev = np.max(np.abs(I))

    print('C^-1 C:', file=logf)
    print('  * max. abs. dev. along diagonal: {:.3g}'.format(I_diag_dev), file=logf)
    print('  * max. abs. dev. off diagonal: {:.3g}'.format(I_off_diag_dev), file=logf)

    I_diag = I_est[np.diag_indices(I_est.shape[0])]
    I_diag_dev = np.max(np.abs(I_diag-1.))
    I_est[np.diag_indices(I_est.shape[0])] = 0.
    I_off_diag_dev = np.max(np.abs(I_est))

    print('(C^-1)_{est} C:', file=logf)
    print('  * max. abs. dev. along diagonal: {:.3g}'.format(I_diag_dev), file=logf)
    print('  * max. abs. dev. off diagonal: {:.3g}'.format(I_off_diag_dev), file=logf)

    #plt.imsave('dist.png', dist, cmap=colors.inferno_r)
    row2fig('dist.png', dist)
    plt.imsave('dist_matrix.png', dist_mat, cmap=colors.inferno_r)
    plt.imsave('cov.png', cov, cmap=colors.inferno_r)
    vmax = np.max(np.abs(inv_cov))
    plt.imsave('inv_cov.png', inv_cov, cmap='coolwarm_r', vmin=-vmax, vmax=vmax)
    vmax = np.max(np.abs(inv_cov_est))
    plt.imsave('inv_cov_est.png', inv_cov_est, cmap='coolwarm_r', vmin=-vmax, vmax=vmax)
    row2fig('inv_cov_est_row.png', inv_cov_est[1387])
    with open('inv_cov_est_diag.txt', 'w') as icer:
        for i in range(npix):
            print(np.max(inv_cov_est[i]), file=icer)
    with open('inv_cov_diag.txt', 'w') as icer:
        for i in range(npix):
            print(np.max(inv_cov[i]), file=icer)
    #plt.imsave('patches.png', patch_identity, cmap=colors.inferno_r)
    row2fig('patches.png', patch_identity)

    return 0

def beam_search(dec,state,y,data,beam_width,mydict_inv):  
    beam_width=beam_width
    xp=cuda.cupy
    batchsize=data.shape[0]
    vocab_size=len(mydict_inv)
    topk=20
    route = np.zeros((batchsize,beam_width,50)).astype(np.int32)
    
    for j in range(50):
        if j == 0:
            y = Variable(xp.array(np.argmax(y.data.get(), axis=1)).astype(xp.int32))
            state,y = dec(y, state, train=False)
            h=state['h1'].data
            c=state['c1'].data
            h=xp.tile(h.reshape(batchsize,1,-1), (1,beam_width,1))
            c=xp.tile(c.reshape(batchsize,1,-1), (1,beam_width,1))
            ptr=F.log_softmax(y).data.get()
            pred_total_city = np.argsort(ptr)[:,::-1][:,:beam_width]
            pred_total_score = np.sort(ptr)[:,::-1][:,:beam_width]
            route[:,:,j] = pred_total_city
            pred_total_city=pred_total_city.reshape(batchsize,beam_width,1)
        else:
            pred_next_score=np.zeros((batchsize,beam_width,topk))
            pred_next_city=np.zeros((batchsize,beam_width,topk)).astype(np.int32)
            score2idx=np.zeros((batchsize,beam_width,topk)).astype(np.int32)
            for b in range(beam_width):
                state={'c1':Variable(c[:,b,:]), 'h1':Variable(h[:,b,:])}
                cur_city = xp.array([pred_total_city[i,b,j-1] for i in range(batchsize)]).astype(xp.int32)
                state,y = dec(cur_city,state, train=False)
                h[:,b,:]=state['h1'].data
                c[:,b,:]=state['c1'].data
                ptr=F.log_softmax(y).data.get()
                pred_next_score[:,b,:]=np.sort(ptr, axis=1)[:,::-1][:,:topk]
                pred_next_city[:,b,:]=np.argsort(ptr, axis=1)[:,::-1][:,:topk]

            h=F.stack([h for i in range(topk)], axis=2).data
            c=F.stack([c for i in range(topk)], axis=2).data
            
            pred_total_city = np.tile(route[:,:,:j],(1,1,topk)).reshape(batchsize,beam_width,topk,j)
            pred_next_city = pred_next_city.reshape(batchsize,beam_width,topk,1)
            pred_total_city = np.concatenate((pred_total_city,pred_next_city),axis=3)

            pred_total_score = np.tile(pred_total_score.reshape(batchsize,beam_width,1),(1,1,topk)).reshape(batchsize,beam_width,topk,1)
            pred_next_score = pred_next_score.reshape(batchsize,beam_width,topk,1)
            pred_total_score += pred_next_score

            idx = pred_total_score.reshape(batchsize,beam_width * topk).argsort(axis=1)[:,::-1][:,:beam_width]

            pred_total_city = pred_total_city[:,idx//topk, np.mod(idx,topk), :][np.diag_indices(batchsize,ndim=2)].reshape(batchsize,beam_width,j+1)
            pred_total_score = pred_total_score[:,idx//topk, np.mod(idx,topk), :][np.diag_indices(batchsize,ndim=2)].reshape(batchsize,beam_width,1)
            h = h[:,idx//topk, np.mod(idx,topk), :][np.diag_indices(batchsize,ndim=2)].reshape(batchsize,beam_width,-1)
            c = c[:,idx//topk, np.mod(idx,topk), :][np.diag_indices(batchsize,ndim=2)].reshape(batchsize,beam_width,-1)

            route[:,:,:j+1] =pred_total_city
            if (pred_total_city[:,:,j] == 15).all():
                break


    return route[:,0,:j+1].tolist()

 def convert(self, network):
     """ Generates an n-dimensional connectivity matrix. """
     if not isinstance(network, NeuralNetwork):
         network = NeuralNetwork(network)
     
     os = np.atleast_1d(self.substrate_shape)
     # Unpack the genotype
     w, f = network.cm.copy(), network.node_types[:]
     
     # Create substrate
     if len(os) == 1:
         cm = np.mgrid[-1:1:os[0]*1j,-1:1:os[0]*1j].transpose((1,2,0))
     elif len(os) == 2:
         cm = np.mgrid[-1:1:os[0]*1j,-1:1:os[1]*1j,-1:1:os[0]*1j,-1:1:os[1]*1j].transpose(1,2,3,4,0)
     else:
         raise NotImplementedError("3+D substrates not supported yet.")
     # Insert a bias
     cm = np.insert(cm, 0, 1.0, -1)
     # Check if the genotype has enough weights
     if w.shape[0] < cm.shape[-1]:
         raise Exception("Genotype weight matrix is too small (%s)" % (w.shape,) )
     # Append zeros
     n_elems = len(f)
     nvals = np.zeros(cm.shape[:-1] + (n_elems - cm.shape[-1],))
     cm = np.concatenate((cm, nvals), -1)
     shape = cm.shape
     
     # Fix the input elements
     frozen = len(os) * 2 + 1
     w[:frozen] = 0.0
     w[np.diag_indices(frozen, 2)] = 1.0
     f[:frozen] = [lambda x: x] * frozen
     w[np.diag_indices(n_elems)] = (1 - self.recursion) * w[np.diag_indices(n_elems)] + self.recursion
     
     # Compute the reaction
     self._steps = []
     laplacian = np.empty_like(cm[..., frozen:])
     kernel = self.diffusion * np.array([1.,2.,1.])
     for _ in range(self.reaction_steps):
         cm = np.dot(w, cm.reshape((-1, n_elems)).T)
         cm = np.clip(cm, self.cm_range[0], self.cm_range[1])
         for el in xrange(cm.shape[0]):
             cm[el,:] = f[el](cm[el,:])
         cm = cm.T.reshape(shape)            
         # apply diffusion
         laplacian[:] = 0.0
         for ax in xrange(cm.ndim - 1):
             laplacian += scipy.ndimage.filters.convolve1d(cm[..., frozen:], kernel, axis=ax, mode='constant')
         cm[..., frozen:] += laplacian
         self._steps.append(cm[...,-1])
         
     # Return the values of the last element (indicating connectivity strength)
     output = cm[..., -1]
     # Build a network object
     net = NeuralNetwork().from_matrix(output)
     if self.sandwich:
         net.make_sandwich()
     return net

def update_amps(cc, t1, t2, eris):
    assert(isinstance(eris, _PhysicistsERIs))
    nocc, nvir = t1.shape
    fock = eris.fock

    fov = fock[:nocc,nocc:]
    mo_e_o = eris.mo_energy[:nocc]
    mo_e_v = eris.mo_energy[nocc:] + cc.level_shift

    tau = imd.make_tau(t2, t1, t1)

    Fvv = imd.cc_Fvv(t1, t2, eris)
    Foo = imd.cc_Foo(t1, t2, eris)
    Fov = imd.cc_Fov(t1, t2, eris)
    Woooo = imd.cc_Woooo(t1, t2, eris)
    Wvvvv = imd.cc_Wvvvv(t1, t2, eris)
    Wovvo = imd.cc_Wovvo(t1, t2, eris)

    # Move energy terms to the other side
    Fvv[np.diag_indices(nvir)] -= mo_e_v
    Foo[np.diag_indices(nocc)] -= mo_e_o

    # T1 equation
    t1new  =  einsum('ie,ae->ia', t1, Fvv)
    t1new += -einsum('ma,mi->ia', t1, Foo)
    t1new +=  einsum('imae,me->ia', t2, Fov)
    t1new += -einsum('nf,naif->ia', t1, eris.ovov)
    t1new += -0.5*einsum('imef,maef->ia', t2, eris.ovvv)
    t1new += -0.5*einsum('mnae,mnie->ia', t2, eris.ooov)
    t1new += fov.conj()

    # T2 equation
    Ftmp = Fvv - 0.5*einsum('mb,me->be', t1, Fov)
    tmp = einsum('ijae,be->ijab', t2, Ftmp)
    t2new = tmp - tmp.transpose(0,1,3,2)
    Ftmp = Foo + 0.5*einsum('je,me->mj', t1, Fov)
    tmp = einsum('imab,mj->ijab', t2, Ftmp)
    t2new -= tmp - tmp.transpose(1,0,2,3)
    t2new += np.asarray(eris.oovv).conj()
    t2new += 0.5*einsum('mnab,mnij->ijab', tau, Woooo)
    t2new += 0.5*einsum('ijef,abef->ijab', tau, Wvvvv)
    tmp = einsum('imae,mbej->ijab', t2, Wovvo)
    tmp -= -einsum('ie,ma,mbje->ijab', t1, t1, eris.ovov)
    tmp = tmp - tmp.transpose(1,0,2,3)
    tmp = tmp - tmp.transpose(0,1,3,2)
    t2new += tmp
    tmp = einsum('ie,jeba->ijab', t1, np.array(eris.ovvv).conj())
    t2new += (tmp - tmp.transpose(1,0,2,3))
    tmp = einsum('ma,ijmb->ijab', t1, np.asarray(eris.ooov).conj())
    t2new -= (tmp - tmp.transpose(0,1,3,2))

    eia = mo_e_o[:,None] - mo_e_v
    eijab = lib.direct_sum('ia,jb->ijab', eia, eia)
    t1new /= eia
    t2new /= eijab

    return t1new, t2new

def test_general_metric(seed=1234, N=2, ndim=3):
    np.random.seed(seed)

    _general_metric(np.eye(ndim), N=N, ndim=ndim)

    L = np.random.randn(ndim, ndim)
    L[np.diag_indices(ndim)] = np.exp(L[np.diag_indices(ndim)])
    L[np.triu_indices(ndim, 1)] = 0.0
    metric = np.dot(L, L.T)
    _general_metric(metric, N=N, ndim=ndim)

def parse_dataFrame(df, df_name):
    """
    removes overlapping cells by x, y, and source image. 
    """

    x = np.array(df['x'])
    y = np.array(df['y'])
    i = np.array(df['image'])

    # create linkage table 
    dx = x[:, np.newaxis] - x
    dy = y[:, np.newaxis] - y
    di = (i[:, np.newaxis] != i)*1e9

    # set diagonal of linkage table to large number
    dx[np.diag_indices(len(dx))] = 1e9
    dy[np.diag_indices(len(dy))] = 1e9
    di[np.diag_indices(len(di))] = 1e9

    # get absolute values
    dx = np.abs(dx)
    dy = np.abs(dy)

    # sum vector of x, y, and image
    d = dx + dy + di

    # set variable = to min of 0 for vector summation
    b = d.min(0)

    # create array with zeros of length b 
    exclude = np.zeros(len(b))

    # iterate through index of length b 
    for i in range(len(b)):
        # if True/anything present at exclude[i]; skip
        if exclude[i]:
            continue
            # add True to index position of exclude if d[i] < 30     
        exclude[(d[i] < 20).nonzero()] = True
        
     
    d = {}

    # create overlap frame
    df_overlap = df[exclude == True]
    df_overlap = df_overlap.reset_index(drop=True)
    

    # create non_overlap frame
    df_nonoverlap = df[exclude == False]
    df_nonoverlap = df_nonoverlap.reset_index(drop=True)
    d['{}_non_overlap'.format(df_name)] = df_nonoverlap

    return d

        def on_stiffness(self, factor):

            diag3 = np.diag_indices( 3 )
            diag6 = np.diag_indices( 6 )            

            # diagonal
            self.stiffness[ diag6 ] = - factor * self.k

            # off-diagonal
            self.stiffness[ 3:, :3 ][diag3] = factor * self.k
            self.stiffness[ :3, 3: ][diag3] = factor * self.k