-
Notifications
You must be signed in to change notification settings - Fork 0
/
tabulate_like.py
executable file
·121 lines (99 loc) · 3.12 KB
/
tabulate_like.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
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
#!/usr/bin/env python
import toy
import time, glob
import cPickle
from scipy import *
from math import *
import sys
class Gridder:
def __init__ (self,N):
self.T=toy.ToyGenerator(0,0)
self.shearlist=[(0.0,0.0),(0.05, 0.00), (0.1, 0.0), (0.05,0.05),
(0.1,0.05) ]
self.grid=[]
self.N=N
self.Ns=len(self.shearlist)
for i in range(self.Ns):
self.grid.append((zeros((N,N)),zeros((N,N))))
def add_samples(self,N):
for i in range(N):
lst=self.T.generate_multishear(self.shearlist)
for e12,grd in zip(lst, self.grid):
if e12!=None:
e1,e2=e12
i1=int((1+e1)/2.*self.N)
i2=int((1+e2)/2.*self.N)
grd[0][i1,i2]+=1
esq=e1**2+e2**2
theta=atan2(e2,e1)+pi
i1=int(esq*self.N)
i2=int(theta/(2*pi)*self.N)
try:
grd[1][i1,i2]+=1
except:
print i1,i2,e1,e2,esq
def dump_to_file(self,fn):
cPickle.dump (self,open(fn,'w'),-1)
def load_many(self,root):
fnl=glob.glob(root+"_*.pickle")
fnl.sort()
for i,fn in enumerate(fnl):
print "Loading ",fn
if i==0:
cself=cPickle.load(open(fn))
# convert tuples to list:
ng=[]
for g in cself.grid:
ng.append(list(g))
cself.grid=ng
print len(ng)
else:
tmp=cPickle.load(open(fn))
for grd,grd2 in zip(cself.grid,tmp.grid):
grd[0]+=grd2[0]
grd[1]+=grd2[1]
print cself.grid[0][1]
self.grid=cself.grid
self.N=cself.N
self.Ns=cself.Ns
self.shearlist=cself.shearlist
class TabLike:
def __init__(self, N,P, L_i, L_ij, Fisher, E3, M13I):
self.N=N
self.P=P
self.L_i=L_i
self.L_ij=L_ij
self.Fisher=Fisher
self.E3=E3
self.M13I=M13I
## we assume Fisher is diagonal here
self.IFisher=array([[1/Fisher[0,0], 0],[0,1/Fisher[1,1]]])
def P(self,e1m,e2m):
i1=int((1+e1m)/2.*self.N)
i2=int((1+e2m)/2.*self.N)
P=self.P[i,i2]
def FD(self,e1m,e2m):
i1=int((1+e1m)/2.*self.N)
i2=int((1+e2m)/2.*self.N)
return array([self.L_i[0][i1,i2], self.L_i[1][i1,i2]])
def E3D(self,e1m,e2m):
i1=int((1+e1m)/2.*self.N)
i2=int((1+e2m)/2.*self.N)
return self.E3[i1,i2]
def SD(self,e1m,e2m):
i1=int((1+e1m)/2.*self.N)
i2=int((1+e2m)/2.*self.N)
i00=self.L_ij[0][0][i1,i2]
i10=self.L_ij[1][0][i1,i2]
i11=self.L_ij[1][1][i1,i2]
return array([[i00,i10],[i10,i11]])
def main():
G=Gridder(200)
Ni=200
N=400000
for i in range(Ni):
print i,'/',Ni
G.add_samples(N)
G.dump_to_file('grids4/grid_%s.pickle'%sys.argv[1])
if __name__=="__main__":
main()