qubx.fast.model

1 """Compiled routines for graph theory and linear constraints. 2 3 Copyright 2008-2014 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 from qubx.fast.fast_utils import * 23 import scipy.linalg 24 import sys 25 26 qubfast.qub_findloops.argtypes = (c_int, c_int, c_int_p, c_int_p, c_int_p) 27 qubfast.qub_findloops.restype = c_int 28

29 -def find_loops(edges):

30 """Returns a list of the fundamental cycles in a bidirectional graph. 31 @param edges: a list of pairs (From, To) of connected state indices; should have both (From, To) and (To, From). 32 @return: a list of cycles; each cycle is a list of state indices 33 34 See David Colquhoun, Kathryn Dowsland, Marco Beato and Andrew Plested. (2004). How to impose microscopic reversibility in complex reaction mechanisms. (with appendix by Kathryn. A. Dowsland and Frank G. Ball) Biophysical Journal, 86, 3510-3518.U{Available here<http://www.ucl.ac.uk/Pharmacology/dc-bits/dcpubs.html>}. 35 """ 36 Nedge = len(edges) 37 if Nedge == 0: 38 return [] 39 edges_arg = numpy.zeros(dtype='int32', shape=(Nedge*2,)) 40 Nvertex = max( max(edge) for edge in edges ) 41 for i, edge in enumerate(edges): 42 edges_arg[2*i:2*(i+1)] = edge 43 Nstate = numpy.zeros(dtype='int32', shape=(Nedge,)) 44 loops = numpy.zeros(dtype='int32', shape=(Nvertex*Nedge,)) 45 Nloop = qubfast.qub_findloops(Nvertex, Nedge, edges_arg.ctypes.data_as(c_int_p), 46 Nstate.ctypes.data_as(c_int_p), loops.ctypes.data_as(c_int_p)) 47 loop_lst = [ [x for x in loops[i*Nvertex:i*Nvertex+Nstate[i]]] 48 for i in xrange(Nloop) ] 49 return loop_lst

50 51 52 53 #TODO: qubx.fast.*constraint* 54 55 56 qubfast.IdentifyFixedParams_p.argtypes = (c_int, c_int, c_double_p, c_double_p, c_double_p, c_int_p, c_int_p, c_int_p, c_int_p) 57 qubfast.IdentifyFixedParams_p.restype = None 58 qubfast.EliminateFixedParams_p.argtypes = (c_int, c_int, c_double_p, c_double_p, c_double_p, c_int, c_int_p, c_int, c_int_p, c_double_p, c_double_p, c_double_p) 59 qubfast.EliminateFixedParams_p.restype = None 60 qubfast.Kr_to_K.argtypes = (c_double_p, c_int, c_double_p, c_int_p) 61 qubfast.Kr_to_K.restype = None 62

63 -def EliminateFixedParams(K, Ain, Bin):

64 """Returns K, Ain, Bin, kR_indices with all fixed rates (and corresponding constraint rows) removed.""" 65 Nc = Ain.shape[0] 66 Nk = Ain.shape[1] 67 Kr_indices = numpy.zeros(shape=(max(1, Nk),), dtype='int32') 68 constraint_indices = numpy.zeros(shape=(max(1, Nc),), dtype='int32') 69 NkR = c_int() 70 NcnsR = c_int() 71 qubfast.IdentifyFixedParams_p(Nk, Nc, K.ctypes.data_as(c_double_p), Ain.ctypes.data_as(c_double_p), Bin.ctypes.data_as(c_double_p), 72 byref(NkR), Kr_indices.ctypes.data_as(c_int_p), byref(NcnsR), constraint_indices.ctypes.data_as(c_int_p)) 73 NkR = NkR.value 74 NcnsR = NcnsR.value 75 Kr = numpy.matrix(numpy.zeros(shape=(NkR,1), dtype='float64')) 76 AinR = numpy.matrix(numpy.zeros(shape=(NcnsR, NkR), dtype='float64')) 77 BinR = numpy.matrix(numpy.zeros(shape=(NcnsR, 1), dtype='float64')) 78 qubfast.EliminateFixedParams_p(Nk, Nc, K.ctypes.data_as(c_double_p), Ain.ctypes.data_as(c_double_p), Bin.ctypes.data_as(c_double_p), 79 NkR, Kr_indices.ctypes.data_as(c_int_p), NcnsR, constraint_indices.ctypes.data_as(c_int_p), 80 Kr.ctypes.data_as(c_double_p), AinR.ctypes.data_as(c_double_p), BinR.ctypes.data_as(c_double_p)) 81 return Kr, AinR, BinR, Kr_indices

82

83 -def Kr_to_K(Kr, K, Kr_indices):

84 qubfast.Kr_to_K(Kr.ctypes.data_as(c_double_p), Kr.shape[0], K.ctypes.data_as(c_double_p), Kr_indices.ctypes.data_as(c_int_p))

85 86 qubfast.ReduceConstraints_p.argtypes = (c_int, c_int, c_double_p, c_double_p) 87 qubfast.ReduceConstraints_p.restype = c_int 88 qubfast.SetupLinearConstraints_p.argtypes = (c_int, c_double_p, c_double_p, 89 c_int, c_double_p, c_double_p, 90 c_double_p, c_double_p, c_double_p, 91 ReportFunc, c_void_p) 92 qubfast.SetupLinearConstraints_p.restype = c_int 93

94 -def ReduceConstraints(Ain, Bin):

95 Ncns = qubfast.ReduceConstraints_p(Ain.shape[1], Ain.shape[0], Ain.ctypes.data_as(c_double_p), Bin.ctypes.data_as(c_double_p)) 96 return numpy.matrix(Ain[:Ncns,:]), numpy.matrix(Bin[:Ncns])

97

98 -def SetupLinearConstraints(Ain, Bin, K):

99 Nc = Ain.shape[0] 100 Nk = Ain.shape[1] 101 Np = Nk - Nc 102 if Np <= 0: 103 print 'number of constraints >= number of parameters:' 104 print Ain 105 print Bin 106 print K 107 P = numpy.matrix(numpy.zeros(shape=(Np,1), dtype='float64')) 108 Acns = numpy.matrix(numpy.zeros(shape=(Nk,Np), dtype='float64')) 109 Bcns = numpy.matrix(numpy.zeros(shape=(Nk,1), dtype='float64')) 110 Ainv = numpy.matrix(numpy.zeros(shape=(Np,Nk), dtype='float64')) 111 def do_report(msg, obj): 112 print msg 113 return 0

114 qubfast.SetupLinearConstraints_p(Nc, Ain.ctypes.data_as(c_double_p), Bin.ctypes.data_as(c_double_p), 115 Nk, K.ctypes.data_as(c_double_p), P.ctypes.data_as(c_double_p), 116 Acns.ctypes.data_as(c_double_p), Bcns.ctypes.data_as(c_double_p), 117 Ainv.ctypes.data_as(c_double_p), ReportFunc(do_report), None) 118 return Acns, Bcns, Ainv, P 119 120 121 # scipy.linalg.eig is better than qub's old eigen 122 EigenFunc = CFUNCTYPE(c_int, c_int, c_double_p, c_double_p, c_double_p, c_double_p, c_double_p) 123 # monkeypatch to allow None as null pointer:

124 -def ef_from_param(cls, obj):

125 if obj is None: 126 return None # return a NULL pointer 127 from ctypes import _CFuncPtr 128 return _CFuncPtr.from_param(obj)

129 EigenFunc.from_param = classmethod(ef_from_param) 130 # could instead wrap each None in a cast: 131 # nullptr = cast(None, EigenFunc) 132

133 -def Eigen_Scipy(n, a_, wr_, wi_, v_, v_inv_):

134 135 #traceback.print_stack() 136 #sys.stdout.flush() 137 138 try: 139 a = numpy.frombuffer(pybuf(a_, n*n*sizeof(c_double)), dtype='float64', count=n*n).reshape((n,n)) 140 v = numpy.frombuffer(pybuf(v_, n*n*sizeof(c_double)), dtype='float64', count=n*n).reshape((n,n)) 141 v_inv = numpy.frombuffer(pybuf(v_inv_, n*n*sizeof(c_double)), dtype='float64', count=n*n).reshape((n,n)) 142 eig, vecs = scipy.linalg.eig(a) 143 for i in xrange(n): 144 wr_[i] = numpy.real(eig[i]) 145 wi_[i] = numpy.imag(eig[i]) 146 for j in xrange(n): 147 v[i,j] = numpy.real(vecs[i,j]) 148 v_inv[:,:] = numpy.real(numpy.matrix(vecs).I) 149 #print 150 #sys.stdout.flush() 151 return 0 152 except: 153 traceback.print_exc() 154 sys.stdout.flush() 155 return -1

156 Eigen_Scipy_Ptr = EigenFunc(Eigen_Scipy) 157

158 -class ModelStruct(Structure):

159 _pack_ = 1 160 _fields_ = [('version', c_int), 161 ('storage', c_void_p), 162 ('Nstate', c_int), 163 ('Nclass', c_int), 164 ('clazz', c_int_p), 165 ('usePeq', c_int), 166 ('IeqFv', c_int), 167 ('Nchannel', c_int), 168 ('iVoltage', c_int), 169 ('Gleak', c_double), 170 ('Vleak', c_double), 171 ('Vrev', c_double), 172 ('amp0', c_double), 173 ('var0', c_double), 174 ('amp', c_double_p), 175 ('var', c_double_p), 176 ('P0', c_double_p), 177 ('K0', c_double_pp), 178 ('K1', c_double_pp), 179 ('L', c_int_pp), 180 ('V', c_int_pp), 181 ('stop_flag', c_int), 182 ('progress', c_int), 183 ('d_Nchannel', c_double), 184 ('d_amp0', c_double), 185 ('d_var0', c_double), 186 ('d_amp', c_double_p), 187 ('d_var', c_double_p), 188 ('d_K0', c_double_pp), 189 ('d_K1', c_double_pp), 190 ('K2', c_double_pp), 191 ('d_K2', c_double_pp), 192 ('P', c_int_pp)]

193 PModel = POINTER(ModelStruct) 194 195 qubfast.qubx_model_create.argtypes = (c_int, c_int, c_int_p, ReportFunc, c_void_p) 196 qubfast.qubx_model_create.restype = PModel 197 qubfast.qubx_model_free.argtypes = (PModel,) 198 qubfast.qubx_model_free.restype = None 199

200 -class Model(object):

201 - def __init__(self, Nstate, Nclass, class_of_state, on_status):

202 self.on_status = on_status 203 self.report_cb = ReportFunc(self.do_report) 204 clazz = numpy.array(class_of_state, dtype='int32') 205 self.obj = qubfast.qubx_model_create(Nstate, Nclass, cdata(clazz, c_int_p), self.report_cb, None) 206 if Nstate: 207 self.clazz_ = numpy.frombuffer(pybuf(self.obj[0].clazz, Nstate*sizeof(c_int)), dtype='int32', count=Nstate) 208 if Nclass: 209 self.amp_ = numpy.frombuffer(pybuf(self.obj[0].amp, Nclass*sizeof(c_double)), dtype='float64', count=Nclass) 210 self.var_ = numpy.frombuffer(pybuf(self.obj[0].var, Nclass*sizeof(c_double)), dtype='float64', count=Nclass) 211 self.amp_err_ = numpy.frombuffer(pybuf(self.obj[0].d_amp, Nclass*sizeof(c_double)), dtype='float64', count=Nclass) 212 self.var_err_ = numpy.frombuffer(pybuf(self.obj[0].d_var, Nclass*sizeof(c_double)), dtype='float64', count=Nclass) 213 if Nstate: 214 self.P0_ = numpy.frombuffer(pybuf(self.obj[0].P0, Nstate*sizeof(c_double)), dtype='float64', count=Nstate) 215 self.K0_ = numpy.frombuffer(pybuf(self.obj[0].K0[0], Nstate*Nstate*sizeof(c_double)), dtype='float64', count=Nstate*Nstate).reshape((Nstate,Nstate)) 216 self.K1_ = numpy.frombuffer(pybuf(self.obj[0].K1[0], Nstate*Nstate*sizeof(c_double)), dtype='float64', count=Nstate*Nstate).reshape((Nstate,Nstate)) 217 self.K2_ = numpy.frombuffer(pybuf(self.obj[0].K2[0], Nstate*Nstate*sizeof(c_double)), dtype='float64', count=Nstate*Nstate).reshape((Nstate,Nstate)) 218 self.K0_err_ = numpy.frombuffer(pybuf(self.obj[0].d_K0[0], Nstate*Nstate*sizeof(c_double)), dtype='float64', count=Nstate*Nstate).reshape((Nstate,Nstate)) 219 self.K1_err_ = numpy.frombuffer(pybuf(self.obj[0].d_K1[0], Nstate*Nstate*sizeof(c_double)), dtype='float64', count=Nstate*Nstate).reshape((Nstate,Nstate)) 220 self.K2_err_ = numpy.frombuffer(pybuf(self.obj[0].d_K2[0], Nstate*Nstate*sizeof(c_double)), dtype='float64', count=Nstate*Nstate).reshape((Nstate,Nstate)) 221 self.L_ = numpy.frombuffer(pybuf(self.obj[0].L[0], Nstate*Nstate*sizeof(c_int)), dtype='int32', count=Nstate*Nstate).reshape((Nstate,Nstate)) 222 self.V_ = numpy.frombuffer(pybuf(self.obj[0].V[0], Nstate*Nstate*sizeof(c_int)), dtype='int32', count=Nstate*Nstate).reshape((Nstate,Nstate)) 223 self.P_ = numpy.frombuffer(pybuf(self.obj[0].P[0], Nstate*Nstate*sizeof(c_int)), dtype='int32', count=Nstate*Nstate).reshape((Nstate,Nstate))

224

225 - def __del__(self):

226 qubfast.qubx_model_free(self.obj)

227 - def dispose(self):

228 """Prevents this object from landing in gc.garbage by breaking a circular reference involving self.report_cb -- must call when you're done using it.""" 229 del self.report_cb # makes circular ref

230 - def do_report(self, msg, obj):

231 self.on_status(msg) 232 return 0

233 version = property(lambda self: self.obj[0].version) 234 Nstate = property(lambda self: self.obj[0].Nstate, doc="state count") 235 Nclass = property(lambda self: self.obj[0].Nclass, doc="class (color) count") 236 clazz = property(lambda self: self.clazz_, doc="class index of each state")

237 - def set_usePeq(self, x): self.obj[0].usePeq = x

238 usePeq = property(lambda self: self.obj[0].usePeq, lambda self, x: self.set_usePeq(x), 239 doc="True to start each segment at equilibrium; False to use normalized States:Pr")

240 - def set_IeqFv(self, x): self.obj[0].IeqFv = x

241 IeqFv = property(lambda self: self.obj[0].IeqFv, lambda self, x: self.set_IeqFv(x), 242 doc="True to model current as a function of voltage, using Cond and CondStd")

243 - def set_Nchannel(self, x): self.obj[0].Nchannel = x

244 Nchannel = property(lambda self: self.obj[0].Nchannel, lambda self, x: self.set_Nchannel(x), 245 doc="number of molecules") 246 Nchannel_err = property(lambda self: self.obj[0].d_Nchannel, doc="estimated error in Nchannel (after some optimizers)")

247 - def set_iVoltage(self, x): self.obj[0].iVoltage = x

248 iVoltage = property(lambda self: self.obj[0].iVoltage, lambda self, x: self.set_iVoltage(x), 249 doc="index of voltage signal (for IeqFv)")

250 - def set_Gleak(self, x): self.obj[0].Gleak = x

251 Gleak = property(lambda self: self.obj[0].Gleak, lambda self, x: self.set_Gleak(x), 252 doc="leakage conductance (pS)")

253 - def set_Vleak(self, x): self.obj[0].Vleak = x

254 Vleak = property(lambda self: self.obj[0].Vleak, lambda self, x: self.set_Vleak(x), 255 doc="leakage voltage (mV)")

256 - def set_Vrev(self, x): self.obj[0].Vrev = x

257 Vrev = property(lambda self: self.obj[0].Vrev, lambda self, x: self.set_Vrev(x), 258 doc="reversal potential (mV)")

259 - def set_amp0(self, x): self.obj[0].amp0 = x

260 amp0 = property(lambda self: self.obj[0].amp0, lambda self, x: self.set_amp0(x), 261 doc="background current (pA) if IeqFv") 262 amp0_err = property(lambda self: self.obj[0].d_amp0, doc="estimated error in amp0 (after some optimizers)")

263 - def set_var0(self, x): self.obj[0].var0 = x

264 var0 = property(lambda self: self.obj[0].var0, lambda self, x: self.set_var0(x), 265 doc="background current variance if IeqFv") 266 var0_err = property(lambda self: self.obj[0].d_var0, doc="estimated error in var0 (after some optimizers)") 267 amp = property(lambda self: self.amp_, doc="current amplitude (pA) or conductance (pS -- if IeqFv) of each class") 268 var = property(lambda self: self.var_, doc="standard deviation of amp (pA or pS) of each class") 269 amp_err = property(lambda self: self.amp_err_, doc="estimated error in each class's amp (after some optimizers)") 270 var_err = property(lambda self: self.var_err_, doc="estimated error in each class's var (after some optimizers)") 271 P0 = property(lambda self: self.P0_, doc="entry probability of each state, if not usePeq") 272 K0 = property(lambda self: self.K0_, doc="matrix[a,b] of pre-exponential rate constants between states") 273 K1 = property(lambda self: self.K1_, doc="matrix[a,b] of exponential v-sensitive rate constants") 274 K2 = property(lambda self: self.K2_, doc="matrix[a,b] of exponential pressure-sensitive rate constants") 275 K0_err = property(lambda self: self.K0_err_, doc="matrix[a,b] of estimated error in K0 (after some optimizers)") 276 K1_err = property(lambda self: self.K1_err_, doc="matrix[a,b] of estimated error in K1 (after some optimizers)") 277 K2_err = property(lambda self: self.K2_err_, doc="matrix[a,b] of estimated error in K2 (after some optimizers)") 278 L = property(lambda self: self.L_, doc="matrix[a,b] of Ligand signal index") 279 V = property(lambda self: self.V_, doc="matrix[a,b] of Voltage signal index") 280 P = property(lambda self: self.P_, doc="matrix[a,b] of Pressure signal index")

281 - def set_stop_flag(self, x): self.obj[0].stop_flag = x

282 stop_flag = property(lambda self: self.obj[0].stop_flag, lambda self, x: self.set_stop_flag(x), 283 doc="advisory to an algorithm to stop and return")

284 - def set_progress(self, x): self.obj[0].progress = x

285 progress = property(lambda self: self.obj[0].progress, lambda self, x: self.set_progress(x), 286 doc="available for algorithms to report percent complete")

287 288

289 -class MultiModelStruct(Structure):

290 _fields_ = [('version', c_int), 291 ('storage', c_void_p), 292 ('single_model', PModel), 293 ('multi_model', PModel), 294 ('nmutistate', c_int), 295 ('nmutipath', c_int), 296 ('nmuticlass', c_int), 297 ('nclass_order', c_int), 298 ('npath', c_int), 299 ('class_order', c_int_p), 300 ('cmutistate', c_int_p), 301 ('path', c_int_p), 302 ('mutistate', c_int_p), 303 ('mutipath', c_int_p)]

304 PMultiModel = POINTER(MultiModelStruct) 305 306 qubfast.qubx_multi_model_create.argtypes = (PModel,) 307 qubfast.qubx_multi_model_create.restype = PMultiModel 308 qubfast.qubx_multi_model_free.argtypes = (PMultiModel,) 309 qubfast.qubx_multi_model_free.restype = None 310 qubfast.qubx_multi_model_update.argtypes = (PMultiModel,) 311 qubfast.qubx_multi_model_update.restype = None 312

313 -class MultiModel(object):

314 - def __init__(self, model):

315 self.__single_model = model # keep it alive for the duration 316 self.obj = qubfast.qubx_multi_model_create(model.obj)

317 - def __del__(self):

318 qubfast.qubx_multi_model_free(self.obj)

319 - def update(self):

320 """Re-derives multi_model rates and amps from single_model.""" 321 qubfast.qubx_multi_model_update(self.obj)

322 version = property(lambda self: self.obj[0].version) 323 single_model = property(lambda self: self.obj[0].single_model) 324 multi_model = property(lambda self: self.obj[0].multi_model) 325 nmutistate = property(lambda self: self.obj[0].nmutistate) 326 nmutipath = property(lambda self: self.obj[0].nmutipath) 327 nmuticlass = property(lambda self: self.obj[0].nmuticlass) 328 nclass_order = property(lambda self: self.obj[0].nclass_order) 329 npath = property(lambda self: self.obj[0].npath) 330 class_order = property(lambda self: self.obj[0].class_order) 331 cmutistate = property(lambda self: self.obj[0].cmutistate) 332 path = property(lambda self: self.obj[0].path) 333 mutistate = property(lambda self: self.obj[0].mutistate) 334 mutipath = property(lambda self: self.obj[0].mutipath)

335 336

337 -class ConstrainedStruct(Structure):

338 _fields_ = [('version', c_int), 339 ('storage', c_void_p), 340 ('Npar', c_int), 341 ('NparR', c_int), 342 ('NfPar', c_int), 343 ('Ncns', c_int), 344 ('NcnsR', c_int), 345 ('Ain', c_double_pp), 346 ('Bin', c_double_p), 347 ('pars', c_double_p), 348 ('fPars', c_double_p), 349 ('Acns', c_double_pp), 350 ('Ascal', c_double_pp), 351 ('Bcns', c_double_p), 352 ('Bscal', c_double_p)]

353 PConstrained = POINTER(ConstrainedStruct) 354 355 qubfast.qubx_constrained_create.argtypes = (c_int,) 356 qubfast.qubx_constrained_create.restype = PConstrained 357 qubfast.qubx_constrained_free.argtypes = (PConstrained,) 358 qubfast.qubx_constrained_free.restype = None 359 qubfast.qubx_constrained_setup.argtypes = (PConstrained,) 360 qubfast.qubx_constrained_setup.restype = c_int 361 qubfast.qubx_constrained_free_to_pars.argtypes = (PConstrained,) 362 qubfast.qubx_constrained_free_to_pars.restype = None 363 qubfast.qubx_constrained_pars_to_free.argtypes = (PConstrained,) 364 qubfast.qubx_constrained_pars_to_free.restype = None 365 qubfast.qubx_constrained_calc_std.argtypes = (PConstrained, c_double_p, c_double_p) 366 qubfast.qubx_constrained_calc_std.restype = None 367

368 -class ConstrainedBase(object):

369 - def __init__(self, obj):

370 self.obj = obj 371 self.Npar_ = obj[0].Npar 372 self.Ain_ = numpy.frombuffer(pybuf(self.obj[0].Ain[0], self.Npar_*self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_*self.Npar_).reshape((self.Npar_,self.Npar_)) 373 self.Bin_ = numpy.frombuffer(pybuf(self.obj[0].Bin, self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_) 374 self.pars_ = numpy.frombuffer(pybuf(self.obj[0].pars, self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_) 375 self.fPars_ = numpy.frombuffer(pybuf(self.obj[0].fPars, self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_)

376 Npar = property(lambda self: self.Npar_) 377 NparR = property(lambda self: self.obj[0].NparR) 378 NfPar = property(lambda self: self.obj[0].NfPar)

379 - def set_Ncns(self, x):

380 self.obj[0].Ncns = x

381 Ncns = property(lambda self: self.obj[0].Ncns, lambda self, x: self.set_Ncns(x)) 382 NcnsR = property(lambda self: self.obj[0].NcnsR) 383 Ain = property(lambda self: self.Ain_) 384 Bin = property(lambda self: self.Bin_) 385 pars = property(lambda self: self.pars_) 386 fPars = property(lambda self: self.fPars_[:self.NfPar]) 387 Acns = property(lambda self: self.Acns_) 388 Ascal = property(lambda self: self._Ascal_) 389 Bcns = property(lambda self: self.Bcns_) 390 Bscal = property(lambda self: self.Bscal_)

391 - def setup(self):

392 result = qubfast.qubx_constrained_setup(self.obj) 393 self.Acns_ = numpy.frombuffer(pybuf(self.obj[0].Acns[0], self.Npar_*self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_*self.Npar_).reshape((self.Npar_,self.Npar_)) 394 self.Ascal_ = numpy.frombuffer(pybuf(self.obj[0].Ascal[0], self.Npar_*self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_*self.Npar_).reshape((self.Npar_,self.Npar_)) 395 self.Bcns_ = numpy.frombuffer(pybuf(self.obj[0].Bin, self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_) 396 self.Bscal_ = numpy.frombuffer(pybuf(self.obj[0].Bin, self.Npar_*sizeof(c_double)), dtype='float64', count=self.Npar_) 397 return result

398 - def free_to_pars(self):

399 qubfast.qubx_constrained_free_to_pars(self.obj)

400 - def pars_to_free(self):

401 qubfast.qubx_constrained_pars_to_free(self.obj)

402 - def calc_std(self, InvHessian):

403 hess = numpy.array(InvHessian, dtype='float64', copy=True) 404 dPar = numpy.zeros(shape=(self.Npar,), dtype='float64') 405 qubfast.qubx_constrained_calc_std(self.obj, cdata(hess, c_double_p), cdata(dPar, c_double_p)) 406 return dPar

407

408 -class Constrained(ConstrainedBase):

409 - def __init__(self, Npar):

410 ConstrainedBase.__init__(self, qubfast.qubx_constrained_create(Npar))

411 - def __del__(self):

412 qubfast.qubx_constrained_free(self.obj)

413

414 -class StimAmpsStruct(Structure):

415 _fields_ = [('version', c_int), 416 ('storage', c_void_p), 417 ('model', PModel), 418 ('Nstim', c_int), 419 ('Nsig', c_int), 420 ('stimclasses', c_double_p), 421 ('sc_frac', c_double_p), 422 ('cns', PConstrained), 423 ('sc_amp', c_double_pp), 424 ('sc_var', c_double_pp), 425 ('sc_std', c_double_pp)]

426 PStimAmps = POINTER(StimAmpsStruct) 427 428 qubfast.qubx_stim_amps_create.argtypes = (PModel, c_int, c_int, c_double_p, c_double_p) 429 qubfast.qubx_stim_amps_create.restype = PStimAmps 430 qubfast.qubx_stim_amps_free.argtypes = (PStimAmps,) 431 qubfast.qubx_stim_amps_free.restype = None 432 qubfast.qubx_stim_amps_setup_scs.argtypes = (PStimAmps,) 433 qubfast.qubx_stim_amps_setup_scs.restype = None 434 qubfast.qubx_stim_amps_pars_to_model.argtypes = (PStimAmps,) 435 qubfast.qubx_stim_amps_pars_to_model.restype = None 436 qubfast.qubx_stim_amps_model_to_pars.argtypes = (PStimAmps,) 437 qubfast.qubx_stim_amps_model_to_pars.restype = None 438 qubfast.qubx_stim_amps_fix_amp0.argtypes = (PStimAmps,) 439 qubfast.qubx_stim_amps_fix_amp0.restype = None 440 qubfast.qubx_stim_amps_fix_var0.argtypes = (PStimAmps,) 441 qubfast.qubx_stim_amps_fix_var0.restype = None 442 qubfast.qubx_stim_amps_fix_amp.argtypes = (PStimAmps,c_int) 443 qubfast.qubx_stim_amps_fix_amp.restype = None 444 qubfast.qubx_stim_amps_fix_var.argtypes = (PStimAmps,c_int) 445 qubfast.qubx_stim_amps_fix_var.restype = None 446

447 -class StimAmps(object):

448 - def __init__(self, model, Nstim, Nsig, stimclasses, sc_frac):

449 model_obj = model.obj if isinstance(model, Model) else model 450 self.obj = qubfast.qubx_stim_amps_create(model_obj, Nstim, Nsig, stimclasses, sc_frac) 451 self.Nstim_ = Nstim 452 self.Nsig_ = Nsig 453 self.stimclasses_ = numpy.frombuffer(pybuf(self.obj[0].stimclasses, Nstim*Nsig*sizeof(c_double)), dtype='float64', count=Nstim*Nsig).reshape((Nstim,Nsig)) 454 self.sc_frac_ = numpy.frombuffer(pybuf(self.obj[0].sc_frac, Nstim*sizeof(c_double)), dtype='float64', count=Nstim) 455 self.cns_ = ConstrainedBase(self.obj[0].cns) 456 self.sc_amp_ = numpy.frombuffer(pybuf(self.obj[0].sc_amp, Nstim*model_obj[0].Nclass*sizeof(c_double)), dtype='float64', count=Nstim*model_obj[0].Nclass).reshape((Nstim, model_obj[0].Nclass)) 457 self.sc_var_ = numpy.frombuffer(pybuf(self.obj[0].sc_var, Nstim*model_obj[0].Nclass*sizeof(c_double)), dtype='float64', count=Nstim*model_obj[0].Nclass).reshape((Nstim, model_obj[0].Nclass)) 458 self.sc_std_ = numpy.frombuffer(pybuf(self.obj[0].sc_std, Nstim*model_obj[0].Nclass*sizeof(c_double)), dtype='float64', count=Nstim*model_obj[0].Nclass).reshape((Nstim, model_obj[0].Nclass))

459 - def __del__(self):

460 qubfast.qubx_stim_amps_free(self.obj)

461 Nstim = property(lambda self: self.Nstim_) 462 Nsig = property(lambda self: self.Nsig_) 463 stimclasses = property(lambda self: self.stimclasses_) 464 sc_frac = property(lambda self: self.sc_frac_) 465 cns = property(lambda self: self.cns_) 466 sc_amp = property(lambda self: self.sc_amp_) 467 sc_var = property(lambda self: self.sc_var_) 468 sc_std = property(lambda self: self.sc_std_)

469 - def setup_scs(self):

470 qubfast.qubx_stim_amps_setup_scs(self.obj)

471 - def pars_to_model(self):

472 qubfast.qubx_stim_amps_pars_to_model(self.obj)

473 - def model_to_pars(self):

474 qubfast.qubx_stim_amps_model_to_pars(self.obj)

475 - def fix_amp0(self):

476 qubfast.qubx_stim_amps_fix_amp0(self.obj)

477 - def fix_var0(self):

478 qubfast.qubx_stim_amps_fix_var0(self.obj)

479 - def fix_amp(self, cls):

480 qubfast.qubx_stim_amps_fix_amp(self.obj, cls)

481 - def fix_var(self, cls):

482 qubfast.qubx_stim_amps_fix_var(self.obj, cls)

483 484 485

486 -class StimRatesStruct(Structure):

487 _fields_ = [('version', c_int), 488 ('storage', c_void_p), 489 ('model', PModel), 490 ('custom_eigen', EigenFunc), 491 ('Nstim', c_int), 492 ('Nsig', c_int), 493 ('stimclasses', c_double_p), 494 ('sc_frac', c_double_p), 495 ('cns', PConstrained), 496 ('QQ', c_double_ppp), 497 ('AA', c_double_ppp), 498 ('QQ_eig', c_double_pp), 499 ('QQ_imag', c_double_pp), 500 ('QQ_evec', c_double_ppp), 501 ('QQ_einv', c_double_ppp), 502 ('A', c_double_pp)]

503 PStimRates = POINTER(StimRatesStruct) 504 505 qubfast.qubx_stim_rates_create.argtypes = (PModel, c_int, c_int, c_double_p, c_double_p, c_double, EigenFunc) 506 qubfast.qubx_stim_rates_create.restype = PStimRates 507 qubfast.qubx_stim_rates_free.argtypes = (PStimRates,) 508 qubfast.qubx_stim_rates_free.restype = None 509 qubfast.qubx_stim_rates_setup_QQ.argtypes = (PStimRates,) 510 qubfast.qubx_stim_rates_setup_QQ.restype = None 511 qubfast.qubx_stim_rates_setup_AA.argtypes = (PStimRates,) 512 qubfast.qubx_stim_rates_setup_AA.restype = None 513 qubfast.qubx_stim_rates_setup_A.argtypes = (PStimRates, c_int, c_double) 514 qubfast.qubx_stim_rates_setup_A.restype = None 515 qubfast.qubx_stim_rates_setup_constraints.argtypes = (PStimRates,) 516 qubfast.qubx_stim_rates_setup_constraints.restype = None 517 qubfast.qubx_stim_rates_pars_to_model.argtypes = (PStimRates,) 518 qubfast.qubx_stim_rates_pars_to_model.restype = None 519 qubfast.qubx_stim_rates_model_to_pars.argtypes = (PStimRates,) 520 qubfast.qubx_stim_rates_model_to_pars.restype = None 521

522 -class StimRates(object):

523 - def __init__(self, model, Nstim, Nsig, stimclasses, sc_frac, sampling, custom_eigen=Eigen_Scipy_Ptr):

524 model_obj = model.obj if isinstance(model, Model) else model 525 self.obj = qubfast.qubx_stim_rates_create(model_obj, Nstim, Nsig, stimclasses, sc_frac, sampling, custom_eigen) 526 self.Nstim_ = Nstim 527 self.Nsig_ = Nsig 528 self.stimclasses_ = numpy.frombuffer(pybuf(self.obj[0].stimclasses, Nstim*Nsig*sizeof(c_double)), dtype='float64', count=Nstim*Nsig).reshape((Nstim, Nsig)) 529 self.sc_frac_ = numpy.frombuffer(pybuf(self.obj[0].sc_frac, Nstim*sizeof(c_double)), dtype='float64', count=Nstim) 530 self.cns_ = ConstrainedBase(self.obj[0].cns) 531 Ns = model_obj[0].Nstate 532 self.QQ_ = [numpy.frombuffer(pybuf(self.obj[0].QQ[sc][0], Ns*Ns*sizeof(c_double)), dtype='float64', count=Ns*Ns).reshape((Ns,Ns)) for sc in xrange(Nstim)] 533 self.AA_ = [numpy.frombuffer(pybuf(self.obj[0].AA[sc][0], Ns*Ns*sizeof(c_double)), dtype='float64', count=Ns*Ns).reshape((Ns,Ns)) for sc in xrange(Nstim)] 534 self.QQ_eig_ = [numpy.frombuffer(pybuf(self.obj[0].QQ_eig[sc], Ns*sizeof(c_double)), dtype='float64', count=Ns) for sc in xrange(Nstim)] 535 self.QQ_imag_ = [numpy.frombuffer(pybuf(self.obj[0].QQ_imag[sc], Ns*sizeof(c_double)), dtype='float64', count=Ns) for sc in xrange(Nstim)] 536 self.QQ_evec_ = [numpy.frombuffer(pybuf(self.obj[0].QQ_evec[sc][0], Ns*Ns*sizeof(c_double)), dtype='float64', count=Ns*Ns).reshape((Ns,Ns)) for sc in xrange(Nstim)] 537 self.QQ_einv_ = [numpy.frombuffer(pybuf(self.obj[0].QQ_einv[sc][0], Ns*Ns*sizeof(c_double)), dtype='float64', count=Ns*Ns).reshape((Ns,Ns)) for sc in xrange(Nstim)] 538 self.AA_ = numpy.frombuffer(pybuf(self.obj[0].A[0], Ns*Ns*sizeof(c_double)), dtype='float64', count=Ns*Ns).reshape((Ns,Ns))

539 - def __del__(self):

540 qubfast.qubx_stim_rates_free(self.obj)

541 Nstim = property(lambda self: self.Nstim_) 542 Nsig = property(lambda self: self.Nsig_) 543 stimclasses = property(lambda self: self.stimclasses_) 544 sc_frac = property(lambda self: self.sc_frac_) 545 cns = property(lambda self: self.cns_) 546 QQ = property(lambda self: self.QQ_) 547 AA = property(lambda self: self.AA_) 548 QQ_eig = property(lambda self: self.QQ_eig_) 549 QQ_imag = property(lambda self: self.QQ_imag_) 550 QQ_evec = property(lambda self: self.QQ_evec_) 551 QQ_einv = property(lambda self: self.QQ_einv_) 552 A = property(lambda self: self.A_)

553 - def setup_QQ(self):

554 qubfast.qubx_stim_rates_setup_QQ(self.obj)

555 - def setup_AA(self):

556 qubfast.qubx_stim_rates_setup_AA(self.obj)

557 - def setup_A(self, sc, dt):

558 qubfast.qubx_stim_rates_setup_A(self.obj, sc, dt)

559 - def setup_constraints(self):

560 qubfast.qubx_stim_rates_setup_constraints(self.obj)

561 - def pars_to_model(self):

562 qubfast.qubx_stim_rates_pars_to_model(self.obj)

563 - def model_to_pars(self):

564 qubfast.qubx_stim_rates_model_to_pars(self.obj)

565

Source Code for Module qubx.fast.model