Python MSMLib.permute_mat方法代码示例

# 需要导入模块: from msmbuilder import MSMLib [as 别名]
# 或者: from msmbuilder.MSMLib import permute_mat [as 别名]
def calculate_fraction_visits(tprob, waypoint, source, sink, return_cond_Q=False):
    """
    Calculate the fraction of times a walker on `tprob` going from `sources`
    to `sinks` will travel through the set of states `waypoints` en route.

    Computes the conditional committors q^{ABC^+} and uses them to find the
    fraction of paths mentioned above. The conditional committors can be

    Note that in the notation of Dickson et. al. this computes h_c(A,B), with
        sources   = A
        sinks     = B
        waypoint  = C

    Parameters
    ----------
    tprob : matrix
        The transition probability matrix
    waypoint : int
        The index of the intermediate state
    sources : nd_array, int or int
        The indices of the source state(s)
    sinks : nd_array, int or int
        The indices of the sink state(s)
    return_cond_Q : bool
        Whether or not to return the conditional committors

    Returns
    -------
    fraction_paths : float
        The fraction of times a walker going from `sources` -> `sinks` stops
        by `waypoints` on its way.
    cond_Q : nd_array, float (optional)
        Optionally returned (`return_cond_Q`)

    See Also
    --------
    calculate_hub_score : function
        Compute the 'hub score', the weighted fraction of visits for an
        entire network.
    calculate_all_hub_scores : function
        Wrapper to compute all the hub scores in a network.

    Notes
    -----
    Employs dense linear algebra,
      memory use scales as N^2
      cycle use scales as N^3

    References
    ----------
    ..[1] Dickson & Brooks (2012), J. Chem. Theory Comput.,
          Article ASAP DOI: 10.1021/ct300537s
    """

    # do some typechecking - we need to be sure that the lumped sources are in
    # the second to last row, and the lumped sinks are in the last row
    # check `tprob`
    msm_analysis.check_transition(tprob)
    if type(tprob) != np.ndarray:
        try:
            tprob = tprob.todense()
        except AttributeError as e:
            raise TypeError('Argument `tprob` must be convertable to a dense'
                            'numpy array. \n%s' % e)

    # typecheck
    for data in [source, sink, waypoint]:
        if type(data) == int:
            pass
        elif hasattr(data, 'len'):
            if len(data) == 1:
                data = data[0]
        else:
            raise TypeError('Arguments source/sink/waypoint must be an int')

    if (source == waypoint) or (sink == waypoint) or (sink == source):
        raise ValueError('source, sink, waypoint must all be disjoint!')

    N = tprob.shape[0]
    Q = calculate_committors([source], [sink], tprob)

    # permute the transition matrix into cannonical form - send waypoint the the
    # last row, and source + sink to the end after that
    Bsink_indices = [source, sink, waypoint]
    perm = np.arange(N)
    perm = np.delete(perm, Bsink_indices)
    perm = np.append(perm, Bsink_indices)
    T = MSMLib.permute_mat(tprob, perm)

    # extract P, R
    n = N - len(Bsink_indices)
    P = T[:n, :n]
    R = T[:n, n:]

    # calculate the conditional committors ( B = N*R ), B[i,j] is the prob
    # state i ends in j, where j runs over the source + sink + waypoint
    # (waypoint is position -1)
    B = np.dot(np.linalg.inv(np.eye(n) - P), R)
    # Not sure if this is sparse or not...

#.........这里部分代码省略.........

# 需要导入模块: from msmbuilder import MSMLib [as 别名]
# 或者: from msmbuilder.MSMLib import permute_mat [as 别名]
    def evaluate_partition_functions(self):
        """
        Computes the partition function of the cut-based free energy profile based
        on transition network and reaction coordinate.

        Employs an MSM to do so, thus taking input in the form of a discrete state
        space and the observed dynamics in that space (assignments). Further, one
        can provide a "reaction coordiante" who's values will be employed to
        find the min-cut/max-flow coordinate describing dynamics on that space.

        Generously contributed by Sergei Krivov
        Optimizations and modifications due to TJL

        Stores
        ------
        zc : nd_array, float
            The cut-based free energy profile along the reaction coordinate. The
            negative log of this function gives the free energy profile along the
            coordinate.
        zh : nd_array, float
            The histogram-based free energy profile.

        See Also
        --------
        optimize_cut_based_coordinate : function
            Optimize a flexible reaction coordinate to maximize the free energy
            barrier along that coordinate.
        """

        if TIME:
            starttime = time.clock()
        self._check_coordinate_values()

        # permute the counts matrix to order it with the rxn coordinate
        state_order_along_RC = np.argsort(self.reaction_coordinate_values)
        perm_counts = MSMLib.permute_mat(self.counts, state_order_along_RC)

        # set up the three variables for weave -- need to be locals
        data = perm_counts.data
        indices = perm_counts.indices
        indptr = perm_counts.indptr

        zc = np.zeros(self.N)
        N = self.N

        import scipy.weave
        scipy.weave.inline(r"""
        Py_BEGIN_ALLOW_THREADS
        int i, j, k;
        double nij, incr;

        #pragma omp parallel for private(nij, incr, j, k) shared(N, indptr, indices, data, zc)
        for (i = 0; i < N; i++) {
            for (k = indptr[i]; k < indptr[i+1]; k++) {
                j = indices[k];
                if (i == j) { continue; }

                nij = data[k];

                // iterate through all values in the matrix in csr format
                // i, j are the row and column indices
                // nij is the entry

                incr = j < i ? nij : -nij;

                #pragma omp critical(zc_update)
                {
                    zc[j] += incr;
                    zc[i] -= incr;
                }
            }
        }
        Py_END_ALLOW_THREADS
        """, ['N', 'data', 'indices', 'indptr', 'zc'], extra_link_args=['-lgomp'],
            extra_compile_args=["-O3", "-fopenmp"])

        zc /= 2.0  # we overcounted in the above - fix that

        # put stuff back in the original order
        inv_ordering = np.argsort(state_order_along_RC)
        zc = zc[inv_ordering]

        # calculate the histogram-based partition function
        zh = np.array(self.counts.sum(axis=0)).flatten()

        self.zc = zc
        self.zh = zh

        if TIME:
            endtime = time.clock()
            print("Time spent in zc:", endtime - starttime)

        return