1 | |
---|
2 | |
---|
3 | import time |
---|
4 | import numpy |
---|
5 | import copy |
---|
6 | import math |
---|
7 | import sys |
---|
8 | import wx |
---|
9 | |
---|
10 | from data_util.calcthread import CalcThread |
---|
11 | |
---|
12 | class CalcChisqr1D(CalcThread): |
---|
13 | """ |
---|
14 | Compute chisqr |
---|
15 | """ |
---|
16 | def __init__(self, x, y,dy, model, |
---|
17 | smearer=None, |
---|
18 | qmin=None, |
---|
19 | qmax=None, |
---|
20 | completefn = None, |
---|
21 | updatefn = None, |
---|
22 | yieldtime = 0.01, |
---|
23 | worktime = 0.01 |
---|
24 | ): |
---|
25 | CalcThread.__init__(self,completefn, |
---|
26 | updatefn, |
---|
27 | yieldtime, |
---|
28 | worktime) |
---|
29 | self.smearer =smearer |
---|
30 | self.y = numpy.array(y) |
---|
31 | self.x = numpy.array(x) |
---|
32 | self.dy= copy.deepcopy(dy) |
---|
33 | self.model = model |
---|
34 | self.qmin = qmin |
---|
35 | self.qmax = qmax |
---|
36 | self.smearer = smearer |
---|
37 | self.starttime = 0 |
---|
38 | |
---|
39 | def isquit(self): |
---|
40 | """ |
---|
41 | @raise KeyboardInterrupt: when the thread is interrupted |
---|
42 | """ |
---|
43 | try: |
---|
44 | CalcThread.isquit(self) |
---|
45 | except KeyboardInterrupt: |
---|
46 | raise KeyboardInterrupt |
---|
47 | |
---|
48 | |
---|
49 | def compute(self): |
---|
50 | """ |
---|
51 | Compute the data given a model function |
---|
52 | """ |
---|
53 | self.starttime = time.time() |
---|
54 | |
---|
55 | x,y = [numpy.asarray(v) for v in (self.x,self.y)] |
---|
56 | if self.dy==None or self.dy==[]: |
---|
57 | self.dy= numpy.zeros(len(self.y)) |
---|
58 | self.dy[self.dy==0]=1 |
---|
59 | |
---|
60 | if self.qmin==0.0 and not numpy.isfinite(self.y[self.qmin]): |
---|
61 | self.qmin = min(self.x[sel.x!=0]) |
---|
62 | elif self.qmin==None: |
---|
63 | self.qmin= min(self.x) |
---|
64 | |
---|
65 | if self.qmax==None: |
---|
66 | self.qmax= max(self.x) |
---|
67 | |
---|
68 | fx = numpy.zeros(len(self.x)) |
---|
69 | |
---|
70 | output= None |
---|
71 | res=[] |
---|
72 | try: |
---|
73 | |
---|
74 | for i_x in range(len(self.x)): |
---|
75 | |
---|
76 | # Check whether we need to bail out |
---|
77 | self.isquit() |
---|
78 | |
---|
79 | fx[i_x]=self.model.run(self.x[i_x]) |
---|
80 | |
---|
81 | if self.smearer!=None: |
---|
82 | fx= self.smearer(fx) |
---|
83 | |
---|
84 | for i_y in range(len(fx)): |
---|
85 | # Check whether we need to bail out |
---|
86 | self.isquit() |
---|
87 | |
---|
88 | temp=(self.y[i_y] - fx[i_y])/self.dy[i_y] |
---|
89 | res.append(temp*temp) |
---|
90 | #sum of residuals |
---|
91 | sum=0 |
---|
92 | for item in res: |
---|
93 | # Check whether we need to bail out |
---|
94 | self.isquit() |
---|
95 | if numpy.isfinite(item): |
---|
96 | sum +=item |
---|
97 | output = sum/ len(res) |
---|
98 | elapsed = time.time()-self.starttime |
---|
99 | self.complete(output= output, elapsed=elapsed) |
---|
100 | |
---|
101 | except KeyboardInterrupt: |
---|
102 | # Thread was interrupted, just proceed and re-raise. |
---|
103 | # Real code should not print, but this is an example... |
---|
104 | raise |
---|
105 | except: |
---|
106 | raise |
---|
107 | |
---|
108 | class CalcChisqr2D(CalcThread): |
---|
109 | """ |
---|
110 | Compute chisqr |
---|
111 | """ |
---|
112 | |
---|
113 | def __init__(self, x_bins, y_bins,data,err_data, model, |
---|
114 | qmin, |
---|
115 | qmax, |
---|
116 | completefn = None, |
---|
117 | updatefn = None, |
---|
118 | yieldtime = 0.01, |
---|
119 | worktime = 0.01 |
---|
120 | ): |
---|
121 | CalcThread.__init__(self,completefn, |
---|
122 | updatefn, |
---|
123 | yieldtime, |
---|
124 | worktime) |
---|
125 | |
---|
126 | self.y_bins = y_bins |
---|
127 | self.x_bins = x_bins |
---|
128 | self.data= data |
---|
129 | self.err_data= copy.deepcopy(err_data) |
---|
130 | self.model = model |
---|
131 | self.qmin = qmin |
---|
132 | self.qmax = qmax |
---|
133 | |
---|
134 | self.starttime = 0 |
---|
135 | |
---|
136 | def isquit(self): |
---|
137 | """ |
---|
138 | @raise KeyboardInterrupt: when the thread is interrupted |
---|
139 | """ |
---|
140 | try: |
---|
141 | CalcThread.isquit(self) |
---|
142 | except KeyboardInterrupt: |
---|
143 | raise KeyboardInterrupt |
---|
144 | |
---|
145 | |
---|
146 | def compute(self): |
---|
147 | """ |
---|
148 | Compute the data given a model function |
---|
149 | """ |
---|
150 | self.starttime = time.time() |
---|
151 | if self.model ==None: |
---|
152 | return |
---|
153 | if self.data==None: |
---|
154 | return |
---|
155 | if self.err_data==None or self.err_data==[]: |
---|
156 | self.err_data= numpy.zeros(len(self.x_bins),len(self.y_bins)) |
---|
157 | |
---|
158 | self.err_data[self.err_data==0]=1 |
---|
159 | |
---|
160 | output= None |
---|
161 | res=[] |
---|
162 | try: |
---|
163 | |
---|
164 | for i in range(len(self.x_bins)): |
---|
165 | # Check whether we need to bail out |
---|
166 | self.isquit() |
---|
167 | for j in range(len(self.y_bins)): |
---|
168 | #Check the range containing data between self.qmin_x and self.qmax_x |
---|
169 | value = math.pow(self.x_bins[i],2)+ math.pow(self.y_bins[j],2) |
---|
170 | if value >= math.pow(self.qmin,2) and value <= math.pow(self.qmax,2): |
---|
171 | |
---|
172 | temp = [self.x_bins[i],self.y_bins[j]] |
---|
173 | error= self.err_data[j][i] |
---|
174 | chisqrji = (self.data[j][i]- self.model.runXY(temp ))/error |
---|
175 | #Vector containing residuals |
---|
176 | res.append( math.pow(chisqrji,2) ) |
---|
177 | |
---|
178 | sum=0 |
---|
179 | for item in res: |
---|
180 | # Check whether we need to bail out |
---|
181 | self.isquit() |
---|
182 | if numpy.isfinite(item): |
---|
183 | sum +=item |
---|
184 | output = sum/ len(res) |
---|
185 | elapsed = time.time()-self.starttime |
---|
186 | self.complete(output= output, elapsed=elapsed) |
---|
187 | |
---|
188 | except KeyboardInterrupt: |
---|
189 | # Thread was interrupted, just proceed and re-raise. |
---|
190 | # Real code should not print, but this is an example... |
---|
191 | raise |
---|