qubx.maths

1 """Mathematical utilities. 2 3 Copyright 2008-2012 Research Foundation State University of New York 4 This file is part of QUB Express. 5 6 QUB Express is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation, either version 3 of the License, or 9 (at your option) any later version. 10 11 QUB Express is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License, 17 named LICENSE.txt, in the QUB Express program directory. If not, see 18 <http://www.gnu.org/licenses/>. 19 20 """ 21 22 import sys 23 24 from numpy import * 25 import scipy 26 import scipy.integrate 27 import scipy.linalg 28 import scipy.optimize 29 import scipy.special 30 erf = scipy.special.erf 31 32 import qubx.accept 33 34 """ 35 Here is a technique to estimate mean and std.dev of an unbounded sequence:: 36 The basic equation you need for a single pass update is: 37 new S2 = old S2 + (x - old_M).(x - new_M) 38 where S2 = sum of squares of deviations from the mean M. 39 40 It is implemented as follows: 41 Start: n = 0, M = 0, S2 = 0 42 Then for each observation x: 43 dev = x - M 44 n = n + 1 45 M = M + dev/n 46 S2 = S2 + dev.(x - M) 47 48 When you want to calculate the sample std. devn. : 49 std.dev. = sqrt( S2 / (n-1) ) 50 -- 51 Alan Miller, Retired Scientist (Statistician) 52 CSIRO Mathematical & Information Sciences 53 a...@vic.cmis.csiro.au, mille...@ozemail.com.au 54 http://www.ozemail.com.au/~milleraj 55 """ 56

57 -class RunningMean(object):

58 """ 59 Accumulates the approximate mean and standard deviation of a partial dataset. 60 61 Example: 62 63 >>> accumulator = RunningMean() 64 >>> for x in dataset: 65 >>> accumulator.add(x) 66 >>> print "after %d points: %f +/- %f" % (accumulator.n, accumulator.mean, accumulator.stdDev) 67 """

68 - def __init__(self):

69 self.n = 0 70 self.mean = 0 71 self.s2 = 0

72 - def add(self, x):

73 dev = x - self.mean 74 self.n += 1 75 self.mean += dev/self.n 76 self.s2 += dev*(x - self.mean)

77 - def get_std(self):

78 return (self.n > 1) and sqrt(self.s2 / (self.n - 1)) or 0.0

79 stdDev = property(get_std)

80 81 82 83 # random utils: 84

85 -def RandomWeighted(weights):

86 """Returns a random index into weights, with likelihood proportional to the elements of weights.""" 87 r = random.uniform() # in [0, 1) 88 tot = weights[0] 89 i = 0 90 last = len(weights)-1 91 while (tot <= r) and (i < last): 92 i += 1 93 tot += weights[i] 94 return i

95

96 -def RandomLT(high):

97 """Returns a uniform random number between 0 and high.""" 98 return random.uniform(high=high)

99

100 -class RandomNormalsSmallN(object):

101 """Generates batches of gaussian random numbers, mean 0, std 1."""

102 - def __init__(self, buffer_size=32768):

103 self.count = buffer_size 104 self.new_buffer()

105 - def new_buffer(self):

106 self.norms = random.randn(self.count) 107 self.i = 0

108 - def next(self, N):

109 """Returns an array of N Gaussian random numbers.""" 110 if N > self.count: 111 return random.randn(N) 112 if self.i + N > self.count: 113 self.new_buffer() 114 last, self.i = self.i, self.i + N 115 return self.norms[last:self.i]

116 117 118

119 -class SampledGaussian(object):

120 """ 121 """

122 - def __init__(self, mean=0.0, std=1.0, resolution=1e-3, cutoff_std=6.0, max_for_complete=256):

123 """ 124 """ 125 self.mean = mean 126 self.std = std 127 self.resolution = float(resolution) 128 self.cutoff_std = cutoff_std 129 self.max_for_complete = max_for_complete 130 self.s2s2 = std * sqrt(2) 131 self.N = int(ceil(2 * std * cutoff_std / resolution)) 132 if self.N > max_for_complete: 133 self.N = 1 134 self.G = self.G_centered 135 self.r2 = resolution / 2.0 136 else: 137 self.G = self.G_sampled 138 self.A = self.make_A() 139 self.A_out = self.G_area(mean+(cutoff_std+1)*std, mean+(cutoff_std+1)*std + resolution) 140 self.y2i_m = self.N / (2 * std * cutoff_std) 141 self.y2i_b = self.N/2 - mean*self.y2i_m

142 - def G_area(self, a, b):

143 mean, s2s2 = self.mean, self.s2s2 144 return 0.5 * (erf((b-mean) / s2s2) - erf((a-mean) / s2s2))

145 - def G_centered(self, x):

146 r2 = self.r2 147 return self.G_area(x - r2, x + r2)

148 - def make_A(self):

149 mean, std, s2s2, cutoff_std, resolution, N = self.mean, self.std, self.s2s2, self.cutoff_std, self.resolution, self.N 150 def G_area2(a, b, cdf_a=None): 151 if cdf_a is None: 152 cdfa = erf((a-mean) / s2s2) 153 else: 154 cdfa = cdf_a 155 cdf_b = erf((b-mean) / s2s2) 156 return (0.5 * (cdf_b - cdfa), cdf_b)

157 A = zeros(shape=(N,), dtype=float32) 158 cdf_a = None 159 b = mean - std * cutoff_std 160 for i in xrange(N): 161 a, b = b, b + resolution 162 A[i], cdf_a = G_area2(a, b, cdf_a) 163 return A

164 - def G_sampled(self, x):

165 i = int(floor(self.y2i_m * x + self.y2i_b)) 166 if 0 <= i < self.N: 167 return self.A[i] 168 else: 169 return self.A_out

170 - def __call__(self, x):

171 """ 172 """ 173 return self.G(x)

174 175 176 177 (METHOD_ODEINT, METHOD_VODE) = (0, 1) 178

179 -class DrivenODESystem(object):

180 - def __init__(self, expr, method=METHOD_ODEINT, locals=None):

181 """Initializes a differential integrator. 182 183 @param expr: ODE system suitable for L{qubx.accept.acceptODEs}(static=['x'], custom=True) 184 """ 185 self.ode = qubx.accept.acceptODEs(custom=True, locals=locals)(expr) 186 self.method = method 187 self.diffeq_names = self.ode.ynames 188 self.input_names = self.ode.args[1:]

189 - def __call__(self, xx, *var0, **signals):

190 """Integrates the ODEs. 191 192 @param xx: numpy.array of x (aka t) values 193 @param var0: one positional argument per differential equation: float initial value 194 @param signals: name -> numpy.array or float for all names in self.input_names 195 @returns: numpy.array, len(xx) rows, len(self.diffeq_names) columns 196 """ 197 if len(xx) == 0: 198 return array([]) 199 if len(var0) != len(self.diffeq_names): 200 raise Exception("Number of equations doesn't match number of initial values.") 201 for name in self.input_names: 202 if not (name in signals): 203 raise Exception("Undefined signal: %s; please provide a numpy.array or constant as a keyword arg." % name) 204 205 xsig = [signals[name] for name in self.input_names] 206 isig = [] 207 for i in xrange(len(xsig)): 208 try: 209 xsig[i][0] 210 isig.append(i) 211 except TypeError: 212 pass 213 214 if self.method == METHOD_ODEINT: 215 return self.solve_odeint(xx, var0, isig, xsig) 216 elif self.method == METHOD_VODE: 217 return self.solve_vode(xx, var0, isig, xsig) 218 else: 219 raise Exception('Unknown method: %s' % self.method)

220 - def solve_odeint(self, xx, var0, driven_ixs, inputs):

221 # the integrator generates yy,x at sub-samples of the data. 222 # if user equations include series other than x, 223 # we need to get the integer+fractional index of x 224 # and interpolate (linearly) the relevant series(s) 225 frac_pos = FracPos(xx) 226 ny = len(self.diffeq_names) 227 dy = self.ode.dy 228 args = [0.0] * (len(self.ode.args) + ny) 229 args[1:len(inputs)+1] = inputs # apply constants; arrays will be replaced with subsamples later 230 def dyy(yy, x): 231 args[0] = x 232 if driven_ixs: 233 frac_pos.find(x) 234 for i in driven_ixs: 235 args[i+1] = frac_pos.subsample(inputs[i]) 236 args[-ny:] = yy 237 return array( dy(*args) )

238 return scipy.integrate.odeint(dyy, var0, xx) #, h0=0.1*(xx[1]-xx[0]))

239 - def solve_vode(self, xx, var0, driven_ixs, inputs):

240 frac_pos = FracPos(xx) 241 ny = len(self.diffeq_names) 242 dy = self.ode.dy 243 args = [0.0] * (len(self.ode.args) + ny) 244 args[1:len(inputs)+1] = inputs 245 def dyy(x, yy): 246 args[0] = x 247 if driven_ixs: 248 frac_pos.find(x) 249 for i in driven_ixs: 250 args[i+1] = frac_pos.subsample(inputs[i]) 251 args[-ny:] = yy 252 return array( dy(*args) )

253 ode = scipy.integrate.ode(dyy) 254 # choosing 'vode' generates an annoying console message; we hide it 255 stdout, sys.stdout = sys.stdout, lambda x: x 256 ode.set_integrator('vode') #, first_step=0.1*(xx[1]-xx[0])) 257 sys.stdout = stdout 258 ode.set_initial_value(var0, xx[0]) 259 ff = zeros(shape=(len(xx), ny), dtype=xx.dtype) 260 ff[0,:] = var0 261 for i in xrange(1, len(xx)): 262 x = xx[i] 263 ff[i,:] = ode.integrate(x) 264 if not ode.successful(): 265 raise Exception("VODE error after %d of %d points" % (i+1, len(xx))) 266 return ff 267 268

269 -class FracPos(object):

270 """Keeps a pointer in an array of data, for random but generally small-offset access by ODE solver. 271 272 @ivar data: list sorted increasing 273 @ivar n: len(data) 274 @ivar ix: current position 275 @ivar frac_ix: fractional position 276 """

277 - def __init__(self, data):

278 # data is sorted increasing 279 self.data = data 280 self.n = len(data) 281 self.ix = 0 282 self.frac_ix = 0.0 283 if len(data) < 2: 284 self.find = lambda x: None 285 self.update = None 286 self.hw = 1.0 287 self.mid = data and data[0] or 0.0 288 else: 289 self.update()

290 - def update(self):

291 ix = self.ix 292 data = self.data 293 if ix == (self.n - 1): 294 self.hw = (data[ix] - data[ix-1]) / 2 295 else: 296 self.hw = (data[ix+1] - data[ix]) / 2 297 self.mid = data[ix] + self.hw

298 #print " at %i (%.3g) / %i\t\tmid: %.3g\thw: %.3g" % (ix, data[ix], self.n, self.mid, self.hw)

299 - def find(self, x):

300 """Moves ix and frac_ix.""" 301 d = x - self.mid; ad = abs(d); # print " dx: %.3g" % d 302 # binary search while far 303 if (-d) > (5*self.hw): 304 step = - max(1, int(ceil((self.ix+1)/2.0))) 305 elif d > (5*self.hw): 306 step = max(1, int(ceil((self.n - self.ix + 1)/2.0))) 307 while ad > (5*self.hw): 308 oldix = self.ix 309 self.ix = min(self.n-1, max(0, self.ix + step)) 310 #print "bsearch %d+%d:\tx = %.3g\td = %.3g\thw = %.3g" % (oldix, step, x, d, self.hw) 311 self.update() 312 d = x - self.mid; ad = abs(d); # print " dx: %.3g" % d 313 if step == 1: break 314 step = max(1, int(ceil(step / 2.0))) 315 if d*step < 0: 316 step = - step 317 # fine adjustment 318 if not (-self.hw <= d < self.hw): 319 if d > 0: 320 while (self.ix < (self.n-1)) and (self.data[self.ix+1] < x): 321 self.ix += 1 322 elif d < 0: 323 while (self.ix > 0) and (self.data[self.ix] > x): 324 self.ix -= 1 325 self.update() 326 d = x - self.mid; ad = abs(d); # print " dx: %.3g" % d 327 if (-1.0001*self.hw) <= d <= (1.0001*self.hw): 328 self.frac_ix = .5 + d/(2*self.hw) 329 else: 330 self.frac_ix = 0.0

331 #print "d: %.3g\thw: %.3g\tx: %.3g\tdata: %.3g" % (d, self.hw, x, self.data[self.ix]) 332 #print " -- > %d, %.3f" % (self.ix, self.frac_ix)

333 - def subsample(self, dd):

334 """Returns the interpolated value from array dd at the fractional index.""" 335 ix, fi = self.ix, self.frac_ix 336 if (ix < (self.n-1)) and fi: # the bastard pushes beyond xx[n-1] 337 return dd[ix] + fi*(dd[ix+1]-dd[ix]) 338 else: # but we treat it like xx[n-1] 339 return dd[ix]

340 341 342 343 344 ## pos' = vel; vel' = (-2*damping*natfreq*vel) - natfreq**2 * (pos - inp) 345

Source Code for Module qubx.maths