Python numpy.diag_indices函数代码示例

    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 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 _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 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()