Python MatrixUtil.m_to_mathematica_string方法代码示例

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import m_to_mathematica_string [as 别名]
def main(args):
    alpha = args.alpha
    N = args.N
    k = 3
    print 'alpha:', alpha
    print 'N:', N
    print 'k:', k
    print
    M = np.array(list(multinomstate.gen_states(N, k)), dtype=int)
    T = multinomstate.get_inverse_map(M)
    R_mut = wrightcore.create_mutation_abc(M, T)
    R_drift = wrightcore.create_moran_drift_rate_k3(M, T)
    Q = alpha * R_mut + R_drift
    # pick out the correct eigenvector
    W, V = scipy.linalg.eig(Q.T)
    w, v = min(zip(np.abs(W), V.T))
    print 'rate matrix:'
    print Q
    print
    print 'transpose of rate matrix:'
    print Q.T
    print
    print 'eigendecomposition of transpose of rate matrix as integers:'
    print scipy.linalg.eig(Q.T)
    print
    print 'transpose of rate matrix in mathematica notation:'
    print MatrixUtil.m_to_mathematica_string(Q.T.astype(int))
    print
    print 'abs eigenvector corresponding to smallest abs eigenvalue:'
    print np.abs(v)
    print

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import m_to_mathematica_string [as 别名]
def main(args):
    alpha = args.alpha
    N = args.N
    k = 4
    print 'alpha:', alpha
    print 'N:', N
    print 'k:', k
    print
    print 'defining the state vectors...'
    M = np.array(list(gen_states_for_induction(N)), dtype=int)
    #M = np.array(list(multinomstate.gen_states(N, k)), dtype=int)
    print 'M.shape:', M.shape
    print 'M:'
    print M
    print
    if args.dense:
        print 'defining the state vector inverse map...'
        T = multinomstate.get_inverse_map(M)
        print 'T.shape:', T.shape
        print 'constructing dense mutation rate matrix...'
        R_mut = wrightcore.create_mutation(M, T)
        print 'constructing dense drift rate matrix...'
        R_drift = wrightcore.create_moran_drift_rate_k4(M, T)
        Q = alpha * R_mut + R_drift
        # pick out the correct eigenvector
        print 'taking eigendecomposition of dense rate matrix...'
        W, V = scipy.linalg.eig(Q.T)
        w, v = min(zip(np.abs(W), V.T))
        if args.eigvals:
            print 'eigenvalues:'
            print W
            # get integer approximations of eigenvalues
            d = collections.defaultdict(int)
            for raw_eigval in W:
                int_eigval = int(np.round(raw_eigval.real))
                d[int_eigval] += 1
            arr = []
            for int_eigval in reversed(sorted(d)):
                s = '%d^%d' % (-int_eigval, d[int_eigval])
                arr.append(s)
            print ' '.join(arr)
        else:
            print 'rate matrix:'
            print Q
            print
            print 'transpose of rate matrix:'
            print Q.T
            print
            print 'eigendecomposition of transpose of rate matrix as integers:'
            print (W, V)
            print
            print 'rate matrix in mathematica notation:'
            print MatrixUtil.m_to_mathematica_string(Q.astype(int))
            print
            print 'abs eigenvector corresponding to smallest abs eigenvalue:'
            print np.abs(v)
            print
    if args.sparse or args.shift_invert:
        print 'defining the state vector inverse dict...'
        T = multinomstate.get_inverse_dict(M)
        print 'sys.getsizeof(T):', sys.getsizeof(T)
        print 'constructing sparse coo mutation+drift rate matrix...'
        R_coo = create_coo_moran(M, T, alpha)
        print 'converting to sparse csr transpose rate matrix...'
        RT_csr = scipy.sparse.csr_matrix(R_coo.T)
        if args.shift_invert:
            print 'compute an eigenpair using shift-invert mode...'
            W, V = scipy.sparse.linalg.eigs(RT_csr, k=1, sigma=1)
        else:
            print 'compute an eigenpair using "small magnitude" mode...'
            W, V = scipy.sparse.linalg.eigs(RT_csr, k=1, which='SM')
        #print 'dense form of sparsely constructed matrix:'
        #print RT_csr.todense()
        #print
        print 'sparse eigenvalues:'
        print W
        print
        print 'sparse stationary distribution eigenvector:'
        print V[:, 0]
        print
        v = abs(V[:, 0])
        v /= np.sum(v)
        autosave_filename = 'full-moran-autosave.txt'
        print 'writing the stationary distn to', autosave_filename, '...'
        with open(autosave_filename, 'w') as fout:
            for p, (X, Y, Z, W) in zip(v, M):
                print >> fout, X, Y, Z, W, p

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import m_to_mathematica_string [as 别名]
def main(args):
    alpha = args.alpha
    beta = args.beta
    if alpha is None:
        alpha = beta
    if beta is None:
        beta = alpha
    if alpha == math.floor(alpha) and beta == math.floor(beta):
        alpha = int(alpha)
        beta = int(beta)
    N = args.N
    print 'alpha:', alpha
    print 'beta:', alpha
    print 'N:', N
    print
    pre_Q = np.zeros((N+1, N+1), dtype=type(alpha))
    print 'dtype of pre-rate-matrix:'
    print pre_Q.dtype
    print
    # Construct a tridiagonal rate matrix.
    # This rate matrix is scaled by pop size so its entries become huge,
    # but this does not affect the equilibrium distribution.
    for i in range(N+1):
        if i > 0:
            # add the drift rate
            pre_Q[i, i-1] += (i*(N-i))
            # add the mutation rate
            pre_Q[i, i-1] += beta * i
        if i < N:
            # add the drift rate
            pre_Q[i, i+1] += ((N-i)*i)
            # add the mutation rate
            pre_Q[i, i+1] += alpha * (N-i)
    Q = pre_Q - np.diag(np.sum(pre_Q, axis=1))
    # define the domain
    #x = np.linspace(0, 1, N+2)
    #x = 0.5 * (x[:-1] + x[1:])
    #x = np.linspace(0, 1, N+1)
    #
    #x = np.linspace(0, 1, N+3)[1:-1]
    #y = scipy.stats.beta.pdf(x, alpha, alpha)
    #y /= scipy.linalg.norm(y)
    # pick out the correct eigenvector
    W, V = scipy.linalg.eig(Q.T)
    w, v = min(zip(np.abs(W), V.T))
    #
    print 'rate matrix:'
    print Q
    print
    print 'transpose of rate matrix:'
    print Q.T
    print
    print 'transpose of rate matrix in mathematica notation:'
    print MatrixUtil.m_to_mathematica_string(Q.T)
    print
    print 'abs eigenvector corresponding to smallest abs eigenvalue:'
    print np.abs(v)
    print
    v_exact = exact_scaled_distn_asymmetric(alpha, beta, N)
    v_exact_normalized = v_exact / scipy.linalg.norm(v_exact)
    print 'exact unnormalized solution:'
    print v_exact
    print
    print 'exact normalized solution:'
    print v_exact_normalized
    print
    v_bernstein = get_bernstein_approximation(alpha, beta, N)
    v_bernstein_normalized = v_bernstein / scipy.linalg.norm(v_bernstein)
    print 'bernstein unnormalized solution:'
    print v_bernstein
    print
    print 'bernstein normalized solution:'
    print v_bernstein_normalized
    print
    v_bernstein = x_get_bernstein_approximation(alpha, beta, N)
    v_bernstein_normalized = v_bernstein / scipy.linalg.norm(v_bernstein)
    print 'bernstein unnormalized solution:'
    print v_bernstein
    print
    print 'bernstein normalized solution:'
    print v_bernstein_normalized
    print