我正在写一个monte carlo模拟,它计算使n顶点图连接所需的最小有向边数的期望值。它通过从一个全0邻接矩阵开始,添加一个有向边,然后测试该矩阵是否表示一个连通图这个过程是循环的,直到一个连通图被构造出来,然后它所需要的迭代次数就是一次试验的样本量。该程序对于小图似乎是精确的,但是一旦图超过~10个顶点,就越来越清楚地表明,它在一个连通图已经构造完成之后,仍然在添加顶点。从本质上说,该算法在大型图上并没有足够早地停止。
看起来可能的罪魁祸首是不相关的功能,但我不完全确定。任何建议都将不胜感激。

import numpy as np
import math
import random

# Randomly adds an edge to the graph by
# choosing a random 0 from the adjecency
# matrix and changing it to a 1.
# @param mat - the list type matrix
# @param k   - the dimension of the matrix
def addEdge(mat, k):

    flag = False

    while flag == False:

        colNum = random.randint(0, k-1)
        rowNum = random.randint(0, k-1)

        if mat[rowNum][colNum] == 0:
            mat[rowNum][colNum] = 1
            flag = True

# Runs singleTrial multiple times and finds
# the average of their sample sizes
def getExpectedValue(size, trials):

    sampleSum = 0.0

    flag = True

    for i in range(trials):
        sample = singleTrial(size)
        sampleSum += sample

    expectedValue = float(sampleSum/trials)

    return expectedValue


# Adds edges to an initially edgeless
# graph UNTIL the graph becomes connected.
def singleTrial(size):

    # Create all zero matrix
    mat = np.random.randint(0,1,size=(size,size))

    sample = 0

    flag = True

    while flag:
        # Uncomment this code to see each matrix that is
        # generated in a single trial. Upon viewing it
        # is clear that this while loop does not terminate
        # early enough.
        #print mat
        #print "\n"
        addEdge(mat, size)
        sample += 1
        if isConnected(mat, size):
            print mat
            flag = False

    print sample
    return sample


# Checks if a given matrix is connected by
# calculating the sum of the number of 1 step
# paths though the number of k step paths
# for every pair of verticies. If this number
# is zero for any pair, then the graph is
# not connected.
def isConnected(A, k):

    B = A

    for i in range(2, k-1):
        B = B + np.linalg.matrix_power(A, i)
    if np.prod(B) != 0:
        return True
    else:
        return False


if __name__ == '__main__' :


    # Perform 1 trial on an 11 vertex graph
    print getExpectedValue(11, 1)

最佳答案

我解决了这个问题。我使用函数numpy.prod(mat)检查矩阵是否包含零,因为它返回矩阵中所有项的乘积。我没有考虑溢出的可能性。我假设,如果numpy.prod的值变得太大(这是肯定的),它将重置为零。这会造成许多错误的否定。

10-07 16:36