Python birth_death_chain.BirthDeathChain类代码示例

    def setUp(self):
        """Store state of the rng"""
        self.state = np.random.mtrand.get_state()

        """Reseed the rng to enforce 'deterministic' behavior"""
        np.random.mtrand.seed(42)
    
        """Meta-stable birth-death chain"""
        b = 2
        q = np.zeros(7)
        p = np.zeros(7)
        q[1:] = 0.5
        p[0:-1] = 0.5
        q[2] = 1.0-10**(-b)
        q[4] = 10**(-b)
        p[2] = 10**(-b)
        p[4] = 1.0-10**(-b)

        bdc = BirthDeathChain(q, p)
        P = bdc.transition_matrix()
        self.dtrajs = [generate_traj(P, 10000, start=0)]
        self.tau = 1

        """Estimate MSM"""
        self.C_MSM = cmatrix(self.dtrajs, self.tau, sliding=True).toarray()
        self.lcc_MSM = largest_connected_set(self.C_MSM)
        self.Ccc_MSM = connected_cmatrix(self.C_MSM, lcc=self.lcc_MSM)
        self.P_MSM = tmatrix(self.Ccc_MSM, reversible=True)
        self.mu_MSM = statdist(self.P_MSM)
        self.k = 3
        self.ts = timescales(self.P_MSM, k=self.k+1, tau=self.tau)[1:]

class TestCommittorDense(unittest.TestCase):
    def setUp(self):
        p = np.zeros(10)
        q = np.zeros(10)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[4] = 0.01
        q[6] = 0.1

        self.bdc = BirthDeathChain(q, p)

    def tearDown(self):
        pass

    def test_forward_comittor(self):
        P = self.bdc.transition_matrix()
        un = committor(P, [0, 1], [8, 9], forward=True)
        u = self.bdc.committor_forward(1, 8)
        assert_allclose(un, u)

    def test_backward_comittor(self):
        P = self.bdc.transition_matrix()
        un = committor(P, [0, 1], [8, 9], forward=False)
        u = self.bdc.committor_backward(1, 8)
        assert_allclose(un, u)

class TestCommittorSparse(unittest.TestCase):
    def setUp(self):
        p = np.zeros(100)
        q = np.zeros(100)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[49] = 0.01
        q[51] = 0.1

        self.bdc = BirthDeathChain(q, p)

    def tearDown(self):
        pass

    def test_forward_comittor(self):
        P = self.bdc.transition_matrix_sparse()
        un = committor(P, range(10), range(90, 100), forward=True)
        u = self.bdc.committor_forward(9, 90)
        assert_allclose(un, u)

    def test_backward_comittor(self):
        P = self.bdc.transition_matrix_sparse()
        un = committor(P, range(10), range(90, 100), forward=False)
        u = self.bdc.committor_backward(9, 90)
        assert_allclose(un, u)

class TestRelaxationDense(unittest.TestCase):
    def setUp(self):
        p=np.zeros(10)
        q=np.zeros(10)
        p[0:-1]=0.5
        q[1:]=0.5
        p[4]=0.01
        q[6]=0.1

        self.bdc=BirthDeathChain(q, p)
        
        self.mu = self.bdc.stationary_distribution()
        self.T = self.bdc.transition_matrix()

        """Test matrix-vector product against spectral decomposition"""        
        R, D, L=rdl_decomposition(self.T)
        self.L=L
        self.R=R
        self.ts=timescales(self.T)
        self.times=np.array([1, 5, 10, 20, 100])

        ev=np.diagonal(D)
        self.ev_t=ev[np.newaxis,:]**self.times[:,np.newaxis]

        self.k=4

        """Observable"""
        obs1 = np.zeros(10)
        obs1[0] = 1
        obs1[1] = 1
        self.obs=obs1

        """Initial distribution"""
        w0=np.zeros(10)
        w0[0:4]=0.25
        self.p0=w0     

    def test_relaxation(self):        
        """k=None"""
        relax_amp=np.dot(self.p0, self.R)*np.dot(self.L, self.obs)
        relax=np.dot(self.ev_t, relax_amp)        
        relaxn=relaxation(self.T, self.p0, self.obs, times=self.times)        
        self.assertTrue(np.allclose(relaxn, relax))        
        
        """k=4"""
        k=self.k
        relax_amp=np.dot(self.p0, self.R[:,0:k])*np.dot(self.L[0:k,:], self.obs)
        relax=np.dot(self.ev_t[:,0:k], relax_amp)        
        relaxn=relaxation(self.T, self.p0, self.obs, k=k, times=self.times) 
        self.assertTrue(np.allclose(relaxn, relax))

class TestCorrelationSparse(unittest.TestCase):
    def setUp(self):
        self.k=4

        p=np.zeros(10)
        q=np.zeros(10)
        p[0:-1]=0.5
        q[1:]=0.5
        p[4]=0.01
        q[6]=0.1

        self.bdc=BirthDeathChain(q, p)
        
        self.mu = self.bdc.stationary_distribution()
        self.T = self.bdc.transition_matrix_sparse()
        R, D, L=rdl_decomposition(self.T, k=self.k)
        self.L=L
        self.R=R
        self.ts=timescales(self.T, k=self.k)
        self.times=np.array([1, 5, 10, 20, 100])

        ev=np.diagonal(D)
        self.ev_t=ev[np.newaxis,:]**self.times[:,np.newaxis]

        obs1 = np.zeros(10)
        obs1[0] = 1
        obs1[1] = 1
        obs2 = np.zeros(10)
        obs2[8] = 1
        obs2[9] = 1

        self.obs1=obs1
        self.obs2=obs2
        self.one_vec=np.ones(10)

    def test_correlation(self):
        """Auto-correlation"""
        acorr_amp=np.dot(self.mu*self.obs1, self.R)*np.dot(self.L, self.obs1)
        acorr=np.dot(self.ev_t, acorr_amp)
        acorrn=correlation(self.T, self.obs1, k=self.k, times=self.times)
        self.assertTrue(np.allclose(acorrn, acorr))

  
        """Cross-correlation"""
        corr_amp=np.dot(self.mu*self.obs1, self.R)*np.dot(self.L, self.obs2)
        corr=np.dot(self.ev_t, corr_amp)    
        corrn=correlation(self.T, self.obs1, obs2=self.obs2, k=self.k, times=self.times)
        self.assertTrue(np.allclose(corrn, corr))    

    def setUp(self):
        self.k=4

        p=np.zeros(10)
        q=np.zeros(10)
        p[0:-1]=0.5
        q[1:]=0.5
        p[4]=0.01
        q[6]=0.1

        self.bdc=BirthDeathChain(q, p)
        
        self.mu = self.bdc.stationary_distribution()
        self.T = self.bdc.transition_matrix_sparse()

        """Test matrix-vector product against spectral decomposition"""        
        R, D, L=rdl_decomposition(self.T, k=self.k)
        self.L=L
        self.R=R
        self.ts=timescales(self.T, k=self.k)
        self.times=np.array([1, 5, 10, 20, 100])

        ev=np.diagonal(D)
        self.ev_t=ev[np.newaxis,:]**self.times[:,np.newaxis]

        """Observable"""
        obs1 = np.zeros(10)
        obs1[0] = 1
        obs1[1] = 1
        self.obs=obs1

        """Initial distribution"""
        w0=np.zeros(10)
        w0[0:4]=0.25
        self.p0=w0     

    def setUp(self):
        p = np.zeros(100)
        q = np.zeros(100)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[49] = 0.01
        q[51] = 0.1

        self.bdc = BirthDeathChain(q, p)

    def setUp(self):
        p = np.zeros(10)
        q = np.zeros(10)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[4] = 0.01
        q[6] = 0.1

        self.bdc = BirthDeathChain(q, p)
        self.T = self.bdc.transition_matrix()
        self.mu = self.bdc.stationary_distribution()

        self.A = create_rate_matrix()

class TestAssessmentDense(unittest.TestCase):
    def setUp(self):
        p = np.zeros(10)
        q = np.zeros(10)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[4] = 0.01
        q[6] = 0.1

        self.bdc = BirthDeathChain(q, p)
        self.T = self.bdc.transition_matrix()
        self.mu = self.bdc.stationary_distribution()

        self.A = create_rate_matrix()

    def test_IsRateMatrix(self):
        self.assert_(is_rate_matrix(self.A), 'A should be a rate matrix')

        # manipulate matrix so it isn't a rate matrix any more
        self.A[0][0] = 3
        self.assertFalse(is_rate_matrix(self.A), 'matrix is not a rate matrix')

    def test_IsReversible(self):
        # create a reversible matrix
        self.assertTrue(is_reversible(self.T, self.mu),
                        "T should be reversible")

    def test_is_transition_matrix(self):
        self.assertTrue(is_transition_matrix(self.T))

        """Larger test-case to prevent too restrictive tolerance settings"""
        X = np.random.random((2000, 2000))
        Tlarge = X / X.sum(axis=1)[:, np.newaxis]
        self.assertTrue(is_transition_matrix(Tlarge))

    def test_is_connected(self):
        self.assertTrue(is_connected(self.T))
        self.assertTrue(is_connected(self.T, directed=False))

class TestExpectation(unittest.TestCase):
    def setUp(self):
        p=np.zeros(10)
        q=np.zeros(10)
        p[0:-1]=0.5
        q[1:]=0.5
        p[4]=0.01
        q[6]=0.1

        self.bdc=BirthDeathChain(q, p)
        
        self.mu = self.bdc.stationary_distribution()
        self.T = self.bdc.transition_matrix()

        obs1 = np.zeros(10)
        obs1[0] = 1
        obs1[1] = 1

        self.obs1=obs1

    def test_expectation(self):
        exp=np.dot(self.mu, self.obs1)
        expn=expectation(self.T, self.obs1)
        self.assertTrue(np.allclose(exp, expn)) 

    def setUp(self):
        self.dim=100
        self.k=10
        
        """Set up meta-stable birth-death chain"""
        p=np.zeros(self.dim)
        p[0:-1]=0.5
        
        q=np.zeros(self.dim)
        q[1:]=0.5

        p[self.dim/2-1]=0.001
        q[self.dim/2+1]=0.001
        
        self.bdc=BirthDeathChain(q, p)

class TestAssessmentDense(unittest.TestCase):
    
    def setUp(self):
        p = np.zeros(10)
        q = np.zeros(10)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[4] = 0.01
        q[6] = 0.1

        self.bdc = BirthDeathChain(q, p)
        self.T = self.bdc.transition_matrix()
        self.mu = self.bdc.stationary_distribution()

        self.A = create_rate_matrix()

    def test_IsRateMatrix(self):
        self.assert_(is_rate_matrix(self.A), \
                     'A should be a rate matrix')
        
        # manipulate matrix so it isn't a rate matrix any more
        self.A[0][0] = 3
        self.assertFalse(is_rate_matrix(self.A), \
                        'matrix is not a rate matrix')
    
    def test_IsReversible(self):
        # create a reversible matrix
        self.assertTrue(is_reversible(self.T, self.mu),
                        "T should be reversible")

    def test_is_transition_matrix(self):
        self.assertTrue(is_transition_matrix(self.T))

    def test_is_connected(self):
        self.assertTrue(is_connected(self.T))
        self.assertTrue(is_connected(self.T, directed=False))

class TestFingerprintSparse(unittest.TestCase):
    def setUp(self):
        self.k=4

        p=np.zeros(10)
        q=np.zeros(10)
        p[0:-1]=0.5
        q[1:]=0.5
        p[4]=0.01
        q[6]=0.1

        self.bdc=BirthDeathChain(q, p)
        
        self.mu = self.bdc.stationary_distribution()
        self.T = self.bdc.transition_matrix_sparse()
        R, D, L=rdl_decomposition(self.T, k=self.k)
        self.L=L
        self.R=R
        self.ts=timescales(self.T, k=self.k)
        self.times=np.array([1, 5, 10, 20])

        ev=np.diagonal(D)
        self.ev_t=ev[np.newaxis,:]**self.times[:,np.newaxis]

        self.tau=7.5

        """Observables"""
        obs1 = np.zeros(10)
        obs1[0] = 1
        obs1[1] = 1
        obs2 = np.zeros(10)
        obs2[8] = 1
        obs2[9] = 1

        self.obs1=obs1
        self.obs2=obs2

        """Initial vector for relaxation"""
        w0=np.zeros(10)
        w0[0:4]=0.25
        self.p0=w0     

    def test_fingerprint_correlation(self):
        """Autocorrelation"""

        """k=4, tau=1"""
        k=self.k
        acorr_amp=np.dot(self.mu*self.obs1, self.R)*np.dot(self.L,self.obs1)
        tsn, acorr_ampn=fingerprint_correlation(self.T, self.obs1, k=k)
        self.assertTrue(np.allclose(tsn, self.ts))
        self.assertTrue(np.allclose(acorr_ampn, acorr_amp))

        """k=4, tau=7.5"""
        tau=self.tau
        tsn, acorr_ampn=fingerprint_correlation(self.T, self.obs1, k=k, tau=tau)
        self.assertTrue(np.allclose(tsn, tau*self.ts))
        self.assertTrue(np.allclose(acorr_ampn, acorr_amp))

        """Cross-correlation"""

        """k=4, tau=1"""
        k=self.k
        corr_amp=np.dot(self.mu*self.obs1, self.R)*np.dot(self.L,self.obs2)
        tsn, corr_ampn=fingerprint_correlation(self.T, self.obs1, obs2=self.obs2, k=k)
        self.assertTrue(np.allclose(tsn, self.ts))
        self.assertTrue(np.allclose(corr_ampn, corr_amp))

        """k=4, tau=7.5"""
        tau=self.tau
        tsn, corr_ampn=fingerprint_correlation(self.T, self.obs1, obs2=self.obs2, k=k, tau=tau)
        self.assertTrue(np.allclose(tsn, tau*self.ts))
        self.assertTrue(np.allclose(corr_ampn, corr_amp))

    def test_fingerprint_relaxation(self):
        one_vec=np.ones(self.T.shape[0])

        relax_amp=np.dot(self.p0, self.R)*np.dot(self.L, self.obs1)
        tsn, relax_ampn=fingerprint_relaxation(self.T, self.p0, self.obs1, k=self.k)        
        self.assertTrue(np.allclose(tsn, self.ts))
        self.assertTrue(np.allclose(relax_ampn, relax_amp))

class TestDecompositionDense(unittest.TestCase):
    def setUp(self):
        self.dim = 100
        self.k = 10

        """Set up meta-stable birth-death chain"""
        p = np.zeros(self.dim)
        p[0:-1] = 0.5

        q = np.zeros(self.dim)
        q[1:] = 0.5

        p[self.dim / 2 - 1] = 0.001
        q[self.dim / 2 + 1] = 0.001

        self.bdc = BirthDeathChain(q, p)

    def test_statdist(self):
        P = self.bdc.transition_matrix()
        mu = self.bdc.stationary_distribution()
        mun = stationary_distribution(P)
        assert_allclose(mu, mun)

    def test_eigenvalues(self):
        P = self.bdc.transition_matrix()
        ev = eigvals(P)
        """Sort with decreasing magnitude"""
        ev = ev[np.argsort(np.abs(ev))[::-1]]

        """k=None"""
        evn = eigenvalues(P)
        assert_allclose(ev, evn)

        """k is not None"""
        evn = eigenvalues(P, k=self.k)
        assert_allclose(ev[0:self.k], evn)

    def test_eigenvectors(self):
        P = self.bdc.transition_matrix()

        # k==None
        ev = eigvals(P)
        ev = ev[np.argsort(np.abs(ev))[::-1]]
        Dn = np.diag(ev)

        # right eigenvectors
        Rn = eigenvectors(P)
        assert_allclose(np.dot(P,Rn),np.dot(Rn,Dn))
        # left eigenvectors
        Ln = eigenvectors(P, right=False)
        assert_allclose(np.dot(Ln.T,P),np.dot(Dn,Ln.T))
        # orthogonality
        Xn = np.dot(Ln.T, Rn)
        di = np.diag_indices(Xn.shape[0])
        Xn[di] = 0.0
        assert_allclose(Xn,0)

        # k!=None
        Dnk = Dn[:,0:self.k][0:self.k,:]
        # right eigenvectors
        Rn = eigenvectors(P, k=self.k)
        assert_allclose(np.dot(P,Rn),np.dot(Rn,Dnk))
        # left eigenvectors
        Ln = eigenvectors(P, right=False, k=self.k)
        assert_allclose(np.dot(Ln.T,P),np.dot(Dnk,Ln.T))
        # orthogonality
        Xn = np.dot(Ln.T, Rn)
        di = np.diag_indices(self.k)
        Xn[di] = 0.0
        assert_allclose(Xn,0)

    def test_rdl_decomposition(self):
        P = self.bdc.transition_matrix()
        mu = self.bdc.stationary_distribution()

        """Non-reversible"""

        """k=None"""
        Rn, Dn, Ln = rdl_decomposition(P)
        Xn = np.dot(Ln, Rn)
        """Right-eigenvectors"""
        assert_allclose(np.dot(P, Rn), np.dot(Rn, Dn))
        """Left-eigenvectors"""
        assert_allclose(np.dot(Ln, P), np.dot(Dn, Ln))
        """Orthonormality"""
        assert_allclose(Xn, np.eye(self.dim))
        """Probability vector"""
        assert_allclose(np.sum(Ln[0, :]), 1.0)

        """k is not None"""
        Rn, Dn, Ln = rdl_decomposition(P, k=self.k)
        Xn = np.dot(Ln, Rn)
        """Right-eigenvectors"""
        assert_allclose(np.dot(P, Rn), np.dot(Rn, Dn))
        """Left-eigenvectors"""
        assert_allclose(np.dot(Ln, P), np.dot(Dn, Ln))
        """Orthonormality"""
        assert_allclose(Xn, np.eye(self.k))
        """Probability vector"""
        assert_allclose(np.sum(Ln[0, :]), 1.0)
#.........这里部分代码省略.........

class TestCorrelationDense(unittest.TestCase):
    def setUp(self):
        p = np.zeros(10)
        q = np.zeros(10)
        p[0:-1] = 0.5
        q[1:] = 0.5
        p[4] = 0.01
        q[6] = 0.1

        self.bdc = BirthDeathChain(q, p)

        self.mu = self.bdc.stationary_distribution()
        self.T = self.bdc.transition_matrix()
        R, D, L = rdl_decomposition(self.T, norm='reversible')
        self.L = L
        self.R = R
        self.ts = timescales(self.T)
        self.times = np.array([1, 5, 10, 20, 100])

        ev = np.diagonal(D)
        self.ev_t = ev[np.newaxis, :] ** self.times[:, np.newaxis]

        self.k = 4

        obs1 = np.zeros(10)
        obs1[0] = 1
        obs1[1] = 1
        obs2 = np.zeros(10)
        obs2[8] = 1
        obs2[9] = 1

        self.obs1 = obs1
        self.obs2 = obs2
        self.one_vec = np.ones(10)

    def test_correlation(self):
        """Auto-correlation"""

        """k=None"""
        acorr_amp = np.dot(self.mu * self.obs1, self.R) * np.dot(self.L, self.obs1)
        acorr = np.dot(self.ev_t, acorr_amp)
        acorrn = correlation(self.T, self.obs1, times=self.times)
        assert_allclose(acorrn, acorr)

        """k=4"""
        k = self.k
        acorr_amp = np.dot(self.mu * self.obs1, self.R[:, 0:k]) * np.dot(self.L[0:k, :], self.obs1)
        acorr = np.dot(self.ev_t[:, 0:k], acorr_amp)
        acorrn = correlation(self.T, self.obs1, times=self.times, k=k)
        assert_allclose(acorrn, acorr)

        """Cross-correlation"""

        """k=None"""
        corr_amp = np.dot(self.mu * self.obs1, self.R) * np.dot(self.L, self.obs2)
        corr = np.dot(self.ev_t, corr_amp)
        corrn = correlation(self.T, self.obs1, obs2=self.obs2, times=self.times)
        assert_allclose(corrn, corr)

        """k=4"""
        k = self.k
        corr_amp = np.dot(self.mu * self.obs1, self.R[:, 0:k]) * np.dot(self.L[0:k, :], self.obs2)
        corr = np.dot(self.ev_t[:, 0:k], corr_amp)
        corrn = correlation(self.T, self.obs1, obs2=self.obs2, times=self.times, k=k)
        assert_allclose(corrn, corr)