calculate the sensitivity matrix for index of the forward committor from A to B given transition matrix T.
    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    A : array like
        List of integer state labels for set A
    B : array like
        List of integer state labels for set B
    index : entry of the committor for which the sensitivity is to be computed

    Returns
    -------
    x : ndarray, shape=(n, n)
        Sensitivity matrix for entry index around transition matrix T. Reversibility is not assumed.

def forward_committor_sensitivity(T, A, B, index):
    """
    calculate the sensitivity matrix for index of the forward committor from A to B given transition matrix T.
    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    A : array like
        List of integer state labels for set A
    B : array like
        List of integer state labels for set B
    index : entry of the committor for which the sensitivity is to be computed

    Returns
    -------
    x : ndarray, shape=(n, n)
        Sensitivity matrix for entry index around transition matrix T. Reversibility is not assumed.
    """

    n = len(T)
    set_X = numpy.arange(n)  # set(range(n))
    set_A = numpy.unique(A)  # set(A)
    set_B = numpy.unique(B)  # set(B)
    set_AB = numpy.union1d(set_A, set_B)  # set_A | set_B
    notAB = numpy.setdiff1d(set_X, set_AB, True)  # list(set_X - set_AB)
    m = len(notAB)

    K = T - numpy.diag(numpy.ones(n))

    U = K[numpy.ix_(notAB.tolist(), notAB.tolist())]

    v = numpy.zeros(m)

    # for i in xrange(0, m):
    #   for k in xrange(0, len(set_B)):
    #       v[i] = v[i] - K[notAB[i], B[k]]
    v[:] = v[:] - K[notAB[:], B[:]]

    qI = numpy.linalg.solve(U, v)

    q_forward = numpy.zeros(n)
    #q_forward[set_A] = 0 # double assignment.
    q_forward[set_B] = 1
    #for i in range(len(notAB)):
    q_forward[notAB[:]] = qI[:]

    target = numpy.eye(1, n, index)
    target = target[0, notAB]

    UinvVec = numpy.linalg.solve(U.T, target)
    Siab = numpy.zeros((n, n))

    for i in range(m):
        Siab[notAB[i]] = - UinvVec[i] * q_forward

    return Siab