-
Notifications
You must be signed in to change notification settings - Fork 1
/
stochastic_trait_matrix.py
48 lines (37 loc) · 1.48 KB
/
stochastic_trait_matrix.py
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
import numpy as np
import sys
import utils
def CreateRandomQ(S, U, num_non_cumu_traits=0):
# Initial all the elements of the matrix to random values
A = np.ones([int(S/2), U]) * 2
B = np.ones([(S - int(S/2)), U]) * 10
C = np.concatenate((A, B), axis=0)
Q = np.random.rand(S, U) * C
# Binary representation for non-cumulative traits
for k in range(num_non_cumu_traits):
Q.T[k] = np.random.randint(2, size=S)
# make the matrix a little sparse.
num = np.random.randint(U/2, (U - num_non_cumu_traits) * S + 1)
for n in range(num):
i = np.random.randint(0, S)
j = np.random.randint(num_non_cumu_traits, U)
Q[i, j] = 0
return Q.astype(np.float)
def CreateRankedQ(S, U, num_non_cumu_traits=0):
# Guarantees that Q has maximum rank (== U).
assert U <= S
Q = CreateRandomQ(S, U, num_non_cumu_traits)
while np.linalg.matrix_rank(Q) != U:
Q = CreateRandomQ(S, U, num_non_cumu_traits)
# Define the variance of each element of Q
var_Q = np.random.rand(Q.shape[0], Q.shape[1])
return Q, var_Q
if __name__ == '__main__':
num_species = 20
num_traits = 10
sys.stdout.write('Generating random matrix with maximum ranks...\t')
sys.stdout.flush()
for i in range(num_traits):
Q = CreateRankedQ(num_species, num_traits)
assert np.linalg.matrix_rank(Q) == num_traits
sys.stdout.write(utils.Highlight('[DONE]\n', utils.GREEN, bold=True))