-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathEq.R
155 lines (98 loc) · 3.84 KB
/
Eq.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
library(pracma)
# Summary: this method calculates the theoretical equilibria of a set of transition PROBABILITY matrices
#
# Parameters:
# transitions: an SxS matrix M of floats between 0 and 1 where M[i,j] is the
# probabilty of transitioning from state j to state i over
# the course of the experiment, OR an array of T such matrices
#
# Returns: An Sx1 matrix of the theoretical equilibiria of the transitions
calcEQs = function(transitions) {
dimen = dim(transitions)
mats = 1
if(length(dimen) == 3) {mats = dimen[3]}
print("dfsa")
print(transitions)
print((dim(transitions)))
print(mats)
eqs = matrix(nrow =(dim(transitions)[1]),ncol = mats)
modtrans = array(transitions,c(dim(transitions)[1:2],mats))
for(k in 1:mats){
cleanup = modtrans[,,k]
for(j in 1:ncol(cleanup)){
cleanup[,j] = cleanup[,j]/sum(cleanup[,j])
}
eigs = eigen(cleanup)
eqs[,k] = (1/sum(eigs$vector[,abs(eigs$values-1)<1e-6]))*eigs$vector[,abs(eigs$values-1)<1e-6]
}
return(apply(eqs,1,mean))
}
# Summary: this method calculates the theoretical equilibria of a set of transition RATE matrices
#
# Parameters:
# transitions: an SxS matrix M of floats where M[i,j] is the
# rate constant of the reaction from state j to state i
# in units 1/sec, OR an array of T such matrices
#
# Returns: An Sx1 matrix of the theoretical equilibiria of the transition rates
theoreticalratesEq = function(rates){
dimen = dim(rates)
mats = 1
if(length(dimen) == 3) {mats = dimen[3]}
eqs = matrix(nrow =dimen[1],ncol = mats)
modtrans = array(rates,c(dimen[1:2],mats))
for(k in 1:mats){
cleanup = modtrans[,,k]
for(j in 1:ncol(cleanup)){
cleanup[j,j] = -sum(cleanup[,j])
}
cleanup = cbind(cleanup, rep(0, dimen[1]))
cleanup[sample.int(dimen[1],1),] = rep(1,dimen[2]+1)
eqs[,k] = rref(cleanup)[,dimen[2]+1]
print(rref(cleanup))
}
return(list(dat=eqs, eq=rowMeans(eqs)))
}
# Summary: this method calculates the chemical rates of transition if the
# experimental probabilities of transition over the experimental
# time frame truly correspond to single transitions over that time
# frame
#
# Parameters:
# transitions: an array of T SxS matrices M of floats between 0 and 1 where M[i,j] is the
# probabilty of transitioning from state j to state i over
# the course of the experiment
#
# Returns: An SxSxT array R, where R[i,j,k] = the 1/sec rate constant of going from state j to state i
naiveRates = function(transitions){
tau = .04
naives = array(dim=dim(transitions))
for(k in 1:dim(transitions)[3]){
for(j in 1:dim(transitions)[2]){
bigL = -log(transitions[j,j,k])/tau
naives[,j,k] = bigL*transitions[,j,k]/(1-exp(-bigL*tau))
naives[j,j,k] = 0
}
}
return(naives)
}
# Summary: this method calculates the theoretical experimental probabilities if each
# transition reaction can occur at most once over the course of the experiment
#
# Parameters:
# transitions: an array of T SxS matrices M of floats where M[i,j] is the
# 1/sec rate constant of going from state j to state i
# the course of the experiment
#
# Returns: An SxSxT array P, where P[i,j,k] = the probability of a single transition
# from state j to state i from the kth set of rate constants
naiveTrans = function(rates){
tau = .04
naives = array(dim=dim(rates))
for(j in 1:dim(rates)[2]){
bigL = sum(rates[,j])
naives[,j] = (rates[,j]/bigL)*(1-exp(-bigL*tau))
naives[j,j] = exp(-bigL*tau)
}
return(naives)
}