Python MatrixUtil.assert_symmetric方法代码示例

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def get_endpoint_conditioned_expected_occupancy(R, v, a, b, T):
    """
    Holmes and Rubin 2002.
    @param R: rate matrix
    @param v: stationary distribution
    @param a: integer state index of initial state
    @param b: integer state index of final state
    @param T: elapsed time
    @return: endpoint conditioned expected amount of time spent in each state
    """
    n = len(v)
    psi = np.sqrt(v)
    S = (R.T * psi).T / psi
    MatrixUtil.assert_symmetric(S)
    w, U = scipy.linalg.eigh(S)
    if not np.allclose(np.dot(U, U.T), np.eye(n)):
        raise Exception('U should be orthogonal')
    P = scipy.linalg.expm(T*R)
    # the Mab is Holmes and Rubin 2002 notation
    Mab = (psi[b] / psi[a]) * np.sum(U[a] * U[b] * np.exp(T*w))
    if not np.allclose(P[a,b], Mab):
        raise Exception('not close: %s %s' % (P[a,b], Mab))
    coeff = (psi[b] / psi[a]) / Mab
    K = _holmes_rubin_2002_kernel(w, T)
    occupancy = coeff * np.array([
        _holmes_rubin_2002_summation(U, a, b, i, K) for i in range(n)])
    if not np.allclose(T, np.sum(occupancy)):
        raise Exception(
                'the expectected occupancy times should add up '
                'to the total time')
    return occupancy

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def fiedler_cut_valuator(A):
    """
    @param A: a symmetric weighted adjacency matrix
    @return: a valuation of the indices of the adjacency matrix
    """
    MatrixUtil.assert_symmetric(A)
    MatrixUtil.assert_nonnegative(A)
    MatrixUtil.assert_hollow(A)
    nverts = A.shape[0]
    L = np.diag(np.sum(A, axis=1)) - A
    w, v = scipy.linalg.eigh(L)
    return v[:, 1]

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def check_generic_cut(valuator, extendor, A):
    """
    The input matrix is expected to have a certain block structure.
    In particular, the leaf vertices are expected to
    precede the points of articulation.
    Because the tree is expected to be an unrooted binary tree,
    the relative number of leaves and points of articulation
    is determined by the size of the adjacency matrix.
    @param A: adjacency matrix of an unrooted edge-weighted binary tree
    @return: True if a counterexample is found
    """
    MatrixUtil.assert_symmetric(A)
    MatrixUtil.assert_nonnegative(A)
    MatrixUtil.assert_hollow(A)
    nverts = A.shape[0]
    if nverts < 4:
        raise Exception('expected at least four vertices')
    if nverts % 2 != 0:
        raise Exception('expected an even number of vertices')
    ntips = nverts / 2 + 1
    narts = nverts / 2 - 1
    # get the schur complement laplacian and its associated adjacency matrix
    L = np.diag(np.sum(A, axis=1)) - A
    L_tips = L[:ntips, :ntips] - ndot(
            L[:ntips, -narts:],
            scipy.linalg.inv(L[-narts:, -narts:]),
            L[-narts:, :ntips],
            )
    A_tips = np.diag(np.diag(L_tips)) - L_tips
    tip_valuations = valuator(A_tips)
    #tip_valuations -= np.mean(tip_valuations)
    #tip_valuations /= np.linalg.norm(tip_valuations)
    art_valuations = extendor(A, tip_valuations)
    #ntip_pos = sum(1 for v in tip_valuations if v > 0)
    #ntip_neg = sum(1 for v in tip_valuations if v < 0)
    #nart_pos = sum(1 for v in art_valuations if v > 0)
    #nart_neg = sum(1 for v in art_valuations if v < 0)
    #print ((ntip_pos, ntip_neg), (nart_pos, nart_neg))
    valuations = np.concatenate((tip_valuations, art_valuations))
    ncrossings = 0
    for i in range(nverts):
        for j in range(i+1, nverts):
            if valuations[i] * valuations[j] * A[i, j] < 0:
                ncrossings += 1
    if ncrossings != 1:
        # found a counterexample!
        print ncrossings
        print A
        return True

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def min_cut_valuator(A):
    """
    @param A: a symmetric weighted adjacency matrix
    @return: a valuation of the indices of the adjacency matrix
    """
    MatrixUtil.assert_symmetric(A)
    MatrixUtil.assert_nonnegative(A)
    MatrixUtil.assert_hollow(A)
    nverts = A.shape[0]
    v_pos = StoerWagner.stoer_wagner_min_cut(A)
    v_neg = set(range(nverts)) - v_pos
    valuations = np.zeros(nverts, dtype=float)
    for v in v_pos:
        valuations[v] = 1
    for v in v_neg:
        valuations[v] = -1
    return valuations

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def stoer_wagner_min_cut(weight_matrix):
    """
    The input matrix is assumed to be a numpy ndarray with float dtype.
    @param weight_matrix: non-negative symmetric weighted adjacency matrix
    @return: the set of indices belonging to one of the two clusters
    """
    w = weight_matrix.copy()
    n = w.shape[0]
    MatrixUtil.assert_symmetric(w)
    MatrixUtil.assert_nonnegative(w)
    MatrixUtil.assert_hollow(w)
    # no cut has been observed so far
    min_cut = None
    min_cut_weight = None
    remaining_indices = set(range(n))
    index_sets = [set([i]) for i in range(n)]
    # reduce the number of remaining indices by one each iteration
    while len(remaining_indices) > 1:
        # initialize the borg
        assimilated_indices = [0]
        unassimilated_indices = remaining_indices - set(assimilated_indices)
        weight_index_pairs = [(w[i, 0], i) for i in unassimilated_indices]
        while len(assimilated_indices) < len(remaining_indices):
            max_weight, best_index = max(weight_index_pairs)
            assimilated_indices.append(best_index)
            unassimilated_indices.remove(best_index)
            next_weight_index_pairs = []
            for weight, index in weight_index_pairs:
                if index is not best_index:
                    next_weight = weight + w[index, best_index]
                    next_weight_index_pairs.append((next_weight, index))
            weight_index_pairs = next_weight_index_pairs
        # the cut between the last two assimilated indices is a possible min cut
        cut = set(index_sets[assimilated_indices[-1]])
        s, t = list(sorted(assimilated_indices[-2:]))
        cut_weight = max_weight
        if min_cut is None or cut_weight < min_cut_weight:
            min_cut, min_cut_weight = cut, cut_weight
        # combine the last two assimilated indices
        for i in remaining_indices:
            w[i, s] += w[i, t]
            w[s, i] += w[t, i]
        index_sets[s].update(index_sets[t])
        remaining_indices.remove(t)
    return min_cut

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def blinkize(pre_Q, blink_birth, blink_death):
    """
    @param pre_Q: a pre-rate matrix
    @param blink_birth: a rate
    @param blink_death: a rate
    @return: a new 'blink-ized' pre-rate matrix
    """
    n = pre_Q.shape[0]
    
    # pre-compute some properties of the new rate matrix
    expo_states = get_expanded_states(n)
    adj_toggle_on = get_toggle_on_adjacency(expo_states)
    adj_toggle_off = get_toggle_off_adjacency(expo_states)
    adj_observable = get_observable_adjacency(expo_states)

    # check some invariants
    if np.any(adj_toggle_on * adj_toggle_off):
        raise Exception
    adj_toggle_either = adj_toggle_on + adj_toggle_off
    MatrixUtil.assert_symmetric(adj_toggle_either)
    MatrixUtil.assert_symmetric(adj_observable)
    if np.any(adj_toggle_either * adj_observable):
        raise Exception

    # initialize the new rate matrix
    nblinks = len(expo_states)
    pre_blink = np.zeros((nblinks, nblinks), dtype=float)
    
    # add the rates that control the toggling between the hidden states
    pre_blink += adj_toggle_off * blink_death
    pre_blink += adj_toggle_on * blink_birth

    # add the rates that control the observable state change
    for i in range(nblinks):
        for j in range(nblinks):
            perm_i, mask_i = expo_states[i]
            perm_j, mask_j = expo_states[j]
            if adj_observable[i, j]:
                a = perm_i[0]
                b = perm_j[0]
                pre_blink[i, j] = pre_Q[a, b]

    # return the newly constructed pre-rate matrix
    return pre_blink

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def harmonic_extension(A, tip_valuations):
    """
    This extension uses the 'Laplacian matrix' and the 'accompanying matrix'.
    @param A: weighted undirected adjacency matrix with tips first
    @param valuations: valuations associated with the first part of A
    @return: extension of valuations to the remaining vertices
    """
    MatrixUtil.assert_symmetric(A)
    MatrixUtil.assert_nonnegative(A)
    MatrixUtil.assert_hollow(A)
    nverts = A.shape[0]
    ntips = tip_valuations.shape[0]
    narts = nverts - ntips
    L = np.diag(np.sum(A, axis=1)) - A
    acc = -np.dot(L[:ntips, -narts:], scipy.linalg.inv(L[-narts:, -narts:]))
    if acc.shape != (ntips, narts):
        raise Exception('unexpected shape of accompanying matrix')
    MatrixUtil.assert_positive(acc)
    return np.dot(tip_valuations, acc)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_symmetric [as 别名]
def combinatorial_extension(A, tip_valuations):
    """
    This extension uses {-1, +1} valuations to minimize zero-crossings.
    It currently uses brute force.
    A more sophisticated approach would minimize not the number of
    zero-crossings, but would instead minimize the number of
    connected components induced by the valuation sign cut.
    @param A: weighted undirected adjacency matrix with tips first
    @param valuations: valuations associated with the first part of A
    @return: extension of valuations to the remaining vertices
    """
    MatrixUtil.assert_symmetric(A)
    MatrixUtil.assert_nonnegative(A)
    MatrixUtil.assert_hollow(A)
    nverts = A.shape[0]
    ntips = tip_valuations.shape[0]
    narts = nverts - ntips
    # get some arbitrarily directed edges from the adjacency matrix
    edges = []
    for i in range(nverts):
        for j in range(i+1, nverts):
            if A[i, j] > 0:
                edges.append((i, j))
    # use brute force to get the best sign valuation of internal vertices
    best_art_vals = None
    best_ncrossings = None
    for art_vals in itertools.product((-1, 1), repeat=narts):
        vfull = np.concatenate((tip_valuations, art_vals))
        ncrossings = 0
        for i, j in edges:
            if A[i, j] * vfull[i] * vfull[j] < 0:
                ncrossings += 1
        if best_art_vals is None or ncrossings < best_ncrossings:
            best_art_vals = art_vals
            best_ncrossings = ncrossings
    return np.array(best_art_vals, dtype=float)