qubx.model

1 """Hidden Markov model and math. 2 3 Copyright 2008-2015 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 qubx.fast.model 23 import qubx.table 24 from qubx.model_constraints import * 25 from qubx.util_types import * 26 from qubx.accept import * 27 from qubx.tree import node_data_or_def 28 from qubx.undo import UndoStack 29 from numpy import * 30 import scipy 31 import scipy.linalg 32 import qubx.tree 33 import collections 34 import itertools 35 import json 36 import os 37 import re 38 import sys 39 40 MIN_TOTAL_RATE = 1e-9 41

42 -class QubModel(object):

43 """ 44 Hidden Markov model, e.g. of one ion channel. 45 46 @ivar undoStack: L{qubx.undo.UndoStack} for all model edits 47 @ivar states: L{qubx.table.SimpleTable} of Class, Pr, x, y, Label 48 @ivar rates: L{qubx.table.SimpleTable} of From, To, k0, k1, k2, Ligand, Voltage, Pressure 49 @ivar classes: L{qubx.table.SimpleTable} of Amp, Std, Cond, CondStd, vRev 50 @ivar constraints_kin: L{qubx.table.SimpleTable} of Type, States 51 @ivar channelCount: number of identical mechanisms 52 @ivar IeqFv: whether to use Cond instead of Amp 53 @ivar vRev: reversal potential in millivolts 54 @ivar k0format: format expr for all k0s 55 @ivar k1format: format expr for all k1s and k2s 56 @ivar changed: (bool) if different from on disk 57 @ivar path: folder/file name 58 @ivar variables: names of active variables e.g. Ligand, Voltage 59 @ivar rate_connecting: dict: (a,b) -> index of rate connecting from state a to state b 60 @ivar should_enforce_constraints: generally True; temporarily False to set many rates at once 61 @ivar balance_loops: True to maintain detailed balance 62 63 @ivar OnSetChannelCount: L{WeakEvent}C{(QubModel, channelCount)} 64 @ivar OnSetIeqFv: L{WeakEvent}C{(QubModel, IeqFv)} 65 @ivar OnSetVRev: L{WeakEvent}C{(QubModel, vRev)} 66 @ivar OnSetK0Format: L{WeakEvent}C{(QubModel, k0format)} 67 @ivar OnSetK1Format: L{WeakEvent}C{(QubModel, k1format)} 68 @ivar OnChangeChanged: L{WeakEvent}C{(QubModel, is_changed)} 69 @ivar OnChangePath: L{WeakEvent}C{(QubModel, path or "")} 70 @ivar OnChangeField: L{WeakEvent}C{(QubModel, 'Name' or 'In Folder' or 'Channel Count' or 'vRev' or 'k*format')} for ObjectTable 71 @ivar OnChangeVariables: L{WeakEvent}C{(QubModel, [var_name])} 72 @ivar OnChangeBalanceLoops: L{WeakEvent}C{(QubModel, balance_loops)} 73 """

74 - def __init__(self):

75 self.ref = Reffer() 76 self.undoStack = UndoStack() 77 self.undoStack.OnChange += self.ref(self._onUndoChange) 78 self.states = qubx.table.SimpleTable('States', global_name='QubX.Models.file.states') 79 self.states.add_field('Class', 0, acceptIntGreaterThanOrEqualTo(0), str, '') 80 self.states.add_field('Pr', 0.0, acceptFloatGreaterThanOrEqualTo(0.0), '%.2f', '') 81 self.states.add_field('x', 50.0, acceptFloatBetween(0.0, 100.0), '%.2f', '%') 82 self.states.add_field('y', 50.0, acceptFloatBetween(0.0, 100.0), '%.2f', '%') 83 self.states.add_field('Label', '', acceptString, str, '') 84 self.states_undo = qubx.table.TableUndo(self.states, self.undoStack) 85 self.states.OnRemoving += self.ref(self._onRemovingState), 'Model.onRemovingState' 86 self.states.OnSet += self.ref(self._onSetState), 'Model.onSetState' 87 self.rates = qubx.table.SimpleTable('Rates', auto_add_fields=True, global_name='QubX.Models.file.rates') 88 self.rates.add_field('From', 0, acceptNothing, self.formatStateAsLbl, '') 89 self.rates.add_field('To', 0, acceptNothing, self.formatStateAsLbl, '') 90 self.rates.add_field('k0', 100.0, acceptFloatGreaterThan(0.0), '%.4g', '') 91 self.rates.add_field('k0 std', 0.0, acceptFloat, '%.3g', '') 92 self.rates.add_field('k1', 0.0, acceptFloat, '%.4g', '') 93 self.rates.add_field('k1 std', 0.0, acceptFloat, '%.3g', '') 94 self.rates.add_field('k2', 0.0, acceptFloat, '%.4g', '') 95 self.rates.add_field('k2 std', 0.0, acceptFloat, '%.3g', '') 96 self.rates.add_field('Ligand', '', acceptString, str, '') 97 self.rates.add_field('Voltage', '', acceptString, str, '') 98 self.rates.add_field('Pressure', '', acceptString, str, '') 99 self.rates.add_field('K', 0.0, acceptNothing, '%.4g', '') 100 self.rates.add_field('Format', '', acceptString, str, '') 101 self.rates_undo = qubx.table.TableUndo(self.rates, self.undoStack) 102 self.rates.OnInsert += self.ref(self._onInsertRate), 'Model.onInsertRate' 103 self.rates.OnRemoving += self.ref(self._onRemovingRate), 'Model.onRemovingRate' 104 self.rates.OnSet += self.ref(self._onSetRate), 'Model.onSetRate' 105 self.classes = qubx.table.SimpleTable('Classes', global_name='QubX.Models.file.classes') 106 self.classes.add_field('Amp', 0.0, acceptFloat, '%.4g', 'pA') 107 self.classes.add_field('Std', 0.1, acceptFloatGreaterThanOrEqualTo(0.0), '%.4g', 'pA') 108 self.classes.add_field('Cond', 0.0, acceptFloat, '%.4g', 'pS') 109 self.classes.add_field('CondStd', 0.1, acceptFloatGreaterThanOrEqualTo(0.0), '%.4g', 'pS') 110 self.classes.add_field('vRev', UNSET_VALUE, acceptFloatOrUnset, formatFloatOrUnset('%.4g'), 'mV') 111 self.classes_undo = qubx.table.TableUndo(self.classes, self.undoStack) 112 self.constraints_kin = ConstraintsKinTable('Kinetic Constraints', 'QubX.Models.file.constraints_kin', self.rates) 113 self.constraints_kin_undo = qubx.table.TableUndo(self.constraints_kin, self.undoStack) 114 self.constraints_kin.OnEdit += self.ref(self._onEditConstraint) 115 # maybe on-bogus-constraint-states-set we nullify and/or complain? 116 # self.classes.add_field('Color', ...) ? 117 # some events e.g. to create a freshly referenced class 118 self._channelCount = 1 119 self._IeqFv = False 120 self._vRev = 0.0 121 self._k0format = '' 122 self._k1format = '' 123 self._changed = False 124 self._path = "Untitled" 125 self._variables = [] 126 self._varcount = {} 127 self._balance_loops = False 128 self.OnSetChannelCount = WeakEvent() # (Model, channelCount) 129 self.OnSetIeqFv = WeakEvent() 130 self.OnSetVRev = WeakEvent() # (Model, vRev) 131 self.OnSetK0Format = WeakEvent() # (Model, k0format) 132 self.OnSetK1Format = WeakEvent() # (Model, k1format) 133 self.OnChangeChanged = WeakEvent() # (Model, is_changed) 134 self.OnChangePath = WeakEvent() # (Model, path or "") 135 self.OnChangeField = WeakEvent() # (Model, 'Name' or 'In Folder' or 'Channel Count' or 'vRev' or 'k*format') for ObjectTable 136 self.OnChangeVariables = WeakEvent() # (Model, variables) 137 self.OnChangeBalanceLoops = WeakEvent() 138 self.rate_connecting = collections.defaultdict() # (a, b) : {} 139 self.constraints_kin.rate_connecting = self.rate_connecting 140 self.should_enforce_constraints = True 141 # not (yet) supported, but carried along: 142 self.constraints_ampvar = \ 143 self.extra_kin_pars = \ 144 self.nAr = \ 145 self.ars = \ 146 self.properties_node = qubx.tree.NullNode() 147 self.__setting_rate = False 148 self.__reading = False 149 self.group = 1

150 channelCount = property(lambda self: self._channelCount, lambda self, x: self.set_channelCount(x)) 151 IeqFv = property(lambda self: self._IeqFv, lambda self, x: self.set_IeqFv(x)) 152 vRev = property(lambda self: self._vRev, lambda self, x: self.set_vRev(x)) 153 k0format = property(lambda self: self._k0format, lambda self, x: self.set_k0format(x)) 154 k1format = property(lambda self: self._k1format, lambda self, x: self.set_k1format(x)) 155 changed = property(lambda self: self._changed, lambda self, x: self.set_changed(x)) 156 path = property(lambda self: self._path, lambda self, x: self.set_path(x)) 157 variables = property(lambda self: self._variables) 158 balance_loops = property(lambda self: self._balance_loops, lambda self, x: self.set_balance_loops(x)) 159

160 - def dispose(self):

161 for tbl in (self.states, self.classes, self.rates, self.constraints_kin): 162 tbl.dispose() 163 self.ref.clear() 164 self.undoStack.clear()

165 - def __str__(self):

166 return '\n'.join(['Channel count: %i' % self.channelCount, 167 'I = f(V): %s' % repr(self.IeqFv), 168 'vRev: %.3g' % self.vRev] + 169 ['-------> %12s <-----------------------------------------------------\n%s' % (lbl, tbl) 170 for lbl, tbl in [('Classes', self.classes), 171 ('States', self.states), 172 ('Rates', self.rates), 173 ('Constraints', self.constraints_kin)] 174 if tbl.size])

175

176 - def set_channelCount(self, x, seal=True):

177 if self._channelCount == x: return 178 prev = self._channelCount 179 self._channelCount = x 180 self.OnSetChannelCount(self, x) 181 self.OnChangeField(self, 'Channel Count') 182 self.undoStack.push_undo(lambda: self.set_channelCount(prev), 183 lambda: self.set_channelCount(x)) 184 if seal: 185 self.undoStack.seal_undo('set channel count')

186 - def set_IeqFv(self, x):

187 x = bool(x) 188 if self._IeqFv == x: return 189 self._IeqFv = x 190 self.OnSetIeqFv(self, x) 191 if x: 192 self._onSetRate(-1, 'Voltage', 'Voltage', '', False) 193 else: 194 self._onSetRate(-1, 'Voltage', '', 'Voltage', False) 195 self.undoStack.push_undo(lambda: self.set_IeqFv(not x), 196 lambda: self.set_IeqFv(x)) 197 self.undoStack.seal_undo('set I=f(V)')

198 - def set_vRev(self, x):

199 if self._vRev == x: return 200 prev = self._vRev 201 self._vRev = x 202 self.OnSetVRev(self, x) 203 self.OnChangeField(self, 'vRev') 204 self.undoStack.push_undo(lambda: self.set_vRev(prev), 205 lambda: self.set_vRev(x)) 206 self.undoStack.seal_undo('set vRev')

207 - def set_k0format(self, x):

208 if self._k0format == x: return 209 prev = self._k0format 210 self._k0format = x 211 self.OnSetK0Format(self, x) 212 self.OnChangeField(self, 'k0 Format') 213 self.undoStack.push_undo(lambda: self.set_k0format(prev), 214 lambda: self.set_k0format(x)) 215 self.undoStack.seal_undo('set k0format')

216 - def set_k1format(self, x):

217 if self._k1format == x: return 218 prev = self._k1format 219 self._k1format = x 220 self.OnSetK1Format(self, x) 221 self.OnChangeField(self, 'k1 Format') 222 self.undoStack.push_undo(lambda: self.set_k1format(prev), 223 lambda: self.set_k1format(x)) 224 self.undoStack.seal_undo('set k1format')

225 - def set_changed(self, x):

226 if self._changed == x: return 227 self._changed = x 228 self.OnChangeChanged(self, x)

229 - def _onUndoChange(self, undoStack):

230 self.changed = True

231 - def set_path(self, x):

232 if self._path == x: return 233 self._path = x 234 self.OnChangePath(self, x) 235 self.OnChangeField(self, 'Name') 236 self.OnChangeField(self, 'In Folder')

237 - def set_balance_loops(self, x):

238 if self._balance_loops == x: return 239 self._balance_loops = x 240 self.constraints_kin.balance_loops = x 241 self.enforce_constraints() 242 self.OnChangeBalanceLoops(self, x)

243 - def acceptStateIndex(self, x):

244 # TODO: use? 245 i = int(x) 246 if not (0 <= i < self.states.size): 247 raise ValueError('State index %s is out of bounds' % x) 248 return i

249 - def formatStateAsLbl(self, i):

250 return self.states.get(i, 'Label') or str(i)

251 - def _gen_connecting(self):

252 self.rate_connecting.clear() 253 rr = self.rates 254 for i in xrange(rr.size): 255 self.rate_connecting[(rr.get(i, 'From'), rr.get(i, 'To'))] = i

256 - def _onEditConstraint(self, need_enforced):

257 if need_enforced and not self.__reading: 258 self.enforce_constraints()

259 - def enforce_constraints(self):

260 """Modifies rates as necessary to satisfy constraints_kin such as detailed (im)balance.""" 261 if not self._balance_loops: 262 return 263 # derive constrained params 264 wanted_names = list(set([self.rates.get(r, 'Ligand') for r in xrange(self.rates.size)] + 265 [self.rates.get(r, 'Voltage') for r in xrange(self.rates.size)] + 266 [self.rates.get(r, 'Pressure') for r in xrange(self.rates.size)])) 267 K0, K1, K2, L, V, P = ModelToRateMatricesP(self, wanted_names) 268 K = K012toK(K0, K1, K2) 269 A, B = self.constraints_kin.get_matrices_p(K0, K1, K2, L, V, P) 270 Kr, AinR, BinR, Kr_indices = qubx.fast.model.EliminateFixedParams(K, A, B) 271 A, B = qubx.fast.model.ReduceConstraints(AinR, BinR) # reduce_constraints(A, B) 272 Ac, Bc, Aci, pars = qubx.fast.model.SetupLinearConstraints(A, B, Kr) # linear_constraints(A, B, K) 273 # derive full params 274 Kr = parsToK(pars, Ac, Aci, Bc) 275 qubx.fast.model.Kr_to_K(Kr, K, Kr_indices) 276 # params to K0, K1 277 K0, K1, K2 = KtoK012(K, K0) 278 # set all rates 279 self.should_enforce_constraints = False 280 for j in xrange(self.rates.size): 281 x,y = [self.rates.get(j, f) for f in ['From', 'To']] 282 self.rates.set(j, 'k0', K0[x,y]) 283 self.rates.set(j, 'k1', K1[x,y]) 284 self.rates.set(j, 'k2', K2[x,y]) 285 self.should_enforce_constraints = True

286 - def get_stimulus(self):

287 stimulus = {} 288 rr = self.rates 289 for r in xrange(rr.size): 290 l = rr.get(r, 'Ligand') 291 if l: 292 stimulus[l] = 1.0 293 v = rr.get(r, 'Voltage') 294 if v: 295 stimulus[v] = 0.0 296 p = rr.get(r, 'Pressure') 297 if p: 298 stimulus[p] = 0.0 299 if self.IeqFv: 300 stimulus['Voltage'] = 0.0 301 return stimulus

302 - def _onSetRate(self, i, field, val, prev, undoing):

303 a, b = [self.rates.get(i, f) for f in ['From', 'To']] 304 if field in ['Ligand', 'Voltage', 'Pressure']: 305 if val == prev: 306 return 307 if prev: 308 ix, ct = self._varcount[prev] 309 if ct > 1: 310 self._varcount[prev][1] -= 1 311 else: 312 del self._variables[ix] 313 del self._varcount[prev] 314 if val: 315 if val in self._varcount: 316 self._varcount[val][1] += 1 317 else: 318 ix = len(self._variables) 319 self._variables.append(val) 320 self._varcount[val] = [ix, 1] 321 self.OnChangeVariables(self, self._variables) 322 elif (not undoing) and self.should_enforce_constraints and (self.constraints_kin.size or (self.balance_loops and ((a, b) in self.constraints_kin.loop_rates))) and (field in ['k0', 'k1', 'k2']): 323 # get model matrices 324 wanted_names = list(set([self.rates.get(r, 'Ligand') for r in xrange(self.rates.size)] + 325 [self.rates.get(r, 'Voltage') for r in xrange(self.rates.size)] + 326 [self.rates.get(r, 'Pressure') for r in xrange(self.rates.size)])) 327 K0, K1, K2, L, V, P = ModelToRateMatricesP(self, wanted_names) 328 K = K012toK(K0, K1, K2) 329 A, B = self.constraints_kin.get_matrices_p(K0, K1, K2, L, V, P) 330 rows, sums = [A], [B] 331 constraints = [] 332 # pre-scale any scaled rates/exps 333 def get_scale_class(ctyp): 334 classes = {} 335 for i in xrange(self.constraints_kin.size): 336 if self.constraints_kin[i, 'Type'] == ctyp: 337 p, q, r, s = self.constraints_kin[i, 'States'] 338 cpq = classes[(p,q)] if (p,q) in classes else None 339 crs = classes[(r,s)] if (r,s) in classes else None 340 if cpq and crs: 341 if not (cpq is crs): 342 for x in crs: 343 cpq[x] = 1 344 classes[(r,s)] = cpq 345 elif cpq: 346 classes[(r,s)] = cpq 347 elif crs: 348 classes[(p,q)] = cpq = crs 349 else: 350 classes[(r,s)] = classes[(p,q)] = cpq = {} 351 cpq[(p,q)] = cpq[(r,s)] = 1 352 if (a, b) in classes: 353 return classes[(a, b)].keys() 354 return []

355 def do_scales(ctyp, mat): 356 factor = val/prev 357 for p,q in get_scale_class(ctyp): 358 if (p,q) != (a,b): 359 mat[p,q] *= factor

360 361 #scale0_sts = [self.constraints_kin.get(i, 'States') for i in xrange(self.constraints_kin.size) 362 # if ('ScaleRate' == self.constraints_kin.get(i, 'Type'))] 363 #scale1_sts = [self.constraints_kin.get(i, 'States') for i in xrange(self.constraints_kin.size) 364 # if ('ScaleExp' == self.constraints_kin.get(i, 'Type'))] 365 #def do_scales(scale_sts, mat): 366 # for sts in scale_sts: 367 # if a == sts[0] and b == sts[1]: 368 # x, y = sts[2:4] 369 # elif a == sts[2] and b == sts[3]: 370 # x, y = sts[0:2] 371 # else: 372 # continue 373 # mat[x,y] *= val/prev 374 if field == 'k0': 375 do_scales('ScaleRate', K0) 376 constraints.append( ('FixRate', (a, b)) ) 377 elif field == 'k1': 378 do_scales('ScaleExp', K1) 379 constraints.append( ('FixExp', (a, b)) ) 380 elif field == 'k2': 381 do_scales('ScalePress', K2) 382 constraints.append( ('FixPress', (a, b)) ) 383 if self._balance_loops: 384 # derive constrained params 385 if constraints: 386 Aextra, Bextra = UserConstraints(K0, K1, K2, L, V, P, constraints) 387 A, B = vstack([A, Aextra]), vstack([B, Bextra]) 388 Kr, AinR, BinR, Kr_indices = qubx.fast.model.EliminateFixedParams(K, A, B) 389 A, B = qubx.fast.model.ReduceConstraints(AinR, BinR) # reduce_constraints(A, B) 390 Ac, Bc, Aci, pars = qubx.fast.model.SetupLinearConstraints(AinR, BinR, Kr) # linear_constraints(A, B, K) 391 # derive full params 392 Kr = parsToK(pars, Ac, Aci, Bc) 393 K = zeros(shape=K.shape, dtype='float64') 394 qubx.fast.model.Kr_to_K(Kr, K, Kr_indices) 395 # params to K0, K1 396 K0, K1, K2 = KtoK012(K, K0) 397 # set all rates 398 try: 399 self.should_enforce_constraints = False 400 for j in xrange(self.rates.size): 401 x,y = [self.rates.get(j, f) for f in ['From', 'To']] 402 self.rates.set(j, 'k0', K0[x,y]) 403 self.rates.set(j, 'k1', K1[x,y]) 404 self.rates.set(j, 'k2', K2[x,y]) 405 finally: 406 self.should_enforce_constraints = True

407 - def _onInsertRate(self, i, undoing):

408 self._gen_connecting() 409 add_stim = False 410 for name in [self.rates.get(i, x) for x in ['Ligand', 'Voltage', 'Pressure']]: 411 if name: 412 if name in self._varcount: 413 self._varcount[name][1] += 1 414 else: 415 ix = len(self._variables) 416 self._variables.append(name) 417 self._varcount[name] = [ix, 1] 418 add_stim = True 419 if add_stim: 420 self.OnChangeVariables(self, self._variables)

421 - def _onRemovingRate(self, i, undoing):

422 qubx.pyenv.env.call_later(self._gen_connecting) 423 del_stim = False 424 for name in [self.rates.get(i, x) for x in ['Ligand', 'Voltage', 'Pressure']]: 425 if name: 426 ix, ct = self._varcount[name] 427 if ct > 1: 428 self._varcount[name][1] -= 1 429 else: 430 del self._variables[ix] 431 del self._varcount[name] 432 del_stim = True 433 for jx in xrange(ix, len(self._variables)): 434 nm = self._variables[jx] 435 kx, count = self._varcount[nm] 436 self._varcount[nm] = [jx, count] 437 if del_stim: 438 self.OnChangeVariables(self, self._variables)

439 - def _onRemovingState(self, i, undoing):

440 r = self.rates 441 for j in reversed(xrange(self.rates.size)): 442 f, t = r.get(j, 'From'), r.get(j, 'To') 443 if (i == f) or (i == t): 444 if j%2 == 0: 445 r.remove(2*(j/2)+1) 446 r.remove(2*(j/2)) 447 else: 448 if r.get(j, 'From') > i: 449 r.set(j, 'From', r.get(j, 'From') - 1) 450 else: 451 r.set(j, 'From', r.get(j, 'From')) 452 if r.get(j, 'To') > i: 453 r.set(j, 'To', r.get(j, 'To') - 1) 454 else: 455 r.set(j, 'To', r.get(j, 'To')) 456 cc = self.constraints_kin 457 for j in xrange(cc.size): 458 sts = map(dec_if_gt(i), cc.get(j, 'States')) 459 cc.set(j, 'States', sts) 460 self._gen_connecting()

461 - def _onSetState(self, i, field, val, prev, undoing):

462 if field == 'Label': 463 for r in xrange(self.rates.size): 464 for rf in ('From', 'To'): 465 if self.rates.get(r, rf) == i: 466 self.rates.set(r, rf, i) 467 elif field == 'Class': 468 while val >= self.classes.size: 469 std = self.classes[self.classes.size-1, 'Std'] if self.classes.size else 0.1 470 cond_std = self.classes[self.classes.size-1, 'CondStd'] if self.classes.size else 0.1 471 self.classes.append({'Amp' : float(self.classes.size), 'Std' : std, 472 'Cond' : float(10*self.classes.size), 'CondStd' : cond_std})

473 - def saveAs(self, spath=None, progressf=lambda pct: True):

474 """Saves the model at as_path.""" 475 path = spath or self.path 476 if not path: 477 raise IOError("Can't save: no file name") 478 ext = os.path.splitext(path)[1].lower() 479 if not (ext in ['.qmf', '.qmj']): 480 ext = '.qmf' 481 path = '%s%s' % (path, ext) 482 if ext == '.qmf': 483 f = self.as_tree() 484 # TODO: save alongside, check, replace/move 485 if not f.saveAs(path, as_copy=True): 486 raise IOError("Error saving model %s" % path) 487 elif ext == '.qmj': 488 open(path, 'w').write(repr(self.as_json())) 489 self.path = path 490 self.changed = False

491 - def revert_to_saved(self):

492 """Re-reads the model from path.""" 493 if not self.path: 494 raise IOError("Can't revert model that was never saved.") 495 saved = OpenQubModel(self.path) 496 self.clear() 497 self.channelCount = saved.channelCount 498 self.IeqFv = saved.IeqFv 499 self.vRev = saved.vRev 500 for x in saved.classes.all_rows(): 501 self.classes.append(x.fields) 502 for x in saved.states.all_rows(): 503 self.states.append(x.fields) 504 for x in saved.rates.all_rows(): 505 self.rates.append(x.fields) 506 for x in saved.constraints_kin.all_rows(): 507 self.constraints_kin.append(x.fields) 508 self.undoStack.clear() 509 self.changed = False

510 - def pre_close(self):

511 pass

512 - def clear(self):

513 """Removes all states, rates and constraints.""" 514 self.constraints_kin.clear() 515 self.rates.clear() 516 self.states.clear() 517 self.classes.clear() 518 self.undoStack.seal_undo('clear')

519 - def balance_all_loops(self):

520 self.balance_loops = True

521 - def as_json(self):

522 json = JsonAnon(states = self.states.to_json(), 523 classes = self.classes.to_json(), 524 rates = self.rates.to_json(), 525 constraints_kin = self.constraints_kin.to_json(), 526 channelCount = self.channelCount, 527 IeqFv = self.IeqFv, 528 vRev = self.vRev, 529 balanceLoops = self.balance_loops) 530 return json

531 - def from_json(self, json):

532 self.__reading = True 533 if 'classes' in json: 534 self.classes.from_json(json['classes']) 535 if 'states' in json: 536 self.states.from_json(json['states']) 537 if 'rates' in json: 538 self.rates.from_json(json['rates']) 539 if 'constraints_kin' in json: 540 self.constraints_kin.from_json(json['constraints_kin']) 541 self.channelCount = json['channelCount'] if ('channelCount' in json) else 1 542 self.IeqFv = json['IeqFv'] if ('IeqFv' in json) else False 543 self.vRev = json['vRev'] if ('vRev' in json) else 0 544 self.balance_loops = json['balanceLoops'] if ('balanceLoops' in json) else False 545 self.__reading = False

546 - def as_tree(self):

547 """Returns the model as a U{QMF<http://www.qub.buffalo.edu/qubdoc/files/qmf.html>} tree. 548 549 @return: L{qubx.tree.Node_numpy} 550 """ 551 f = qubx.tree.Node('ModelFile') 552 states = f['States'] 553 rates = f['Rates'] 554 constraints = f['Constraints'] 555 for i in xrange(self.states.size): 556 s = states.append('State') 557 s['x'].data = self.states.get(i, 'x') 558 s['y'].data = self.states.get(i, 'y') 559 s['Class'].data = self.states.get(i, 'Class') 560 s['Label'].data = self.states.get(i, 'Label') 561 s['Pr'].data = self.states.get(i, 'Pr') 562 s['Gr'].data = self.states.get(i, 'Group') 563 for j in xrange(1, self.rates.size, 2): 564 i = j - 1 565 r = rates.append('Rate') 566 r['States'].data = (self.rates.get(i, 'From'), self.rates.get(i, 'To')) 567 r['k0'].data = (self.rates.get(i, 'k0'), self.rates.get(j, 'k0')) 568 r['k1'].data = (self.rates.get(i, 'k1'), self.rates.get(j, 'k1')) 569 r['k2'].data = (self.rates.get(i, 'k2'), self.rates.get(j, 'k2')) 570 r['P'].data = (self.rates.get(i, 'Ligand') and 1 or 0, self.rates.get(j, 'Ligand') and 1 or 0) 571 r['Q'].data = (self.rates.get(i, 'Voltage') and 1 or 0, self.rates.get(j, 'Voltage') and 1 or 0) 572 r['Pressure'].data = (self.rates.get(i, 'Pressure') and 1 or 0, self.rates.get(j, 'Pressure') and 1 or 0) 573 r['PNames'].append('PName').data = self.rates.get(i, 'Ligand') 574 r['PNames'].append('PName').data = self.rates.get(j, 'Ligand') 575 r['QNames'].append('QName').data = self.rates.get(i, 'Voltage') 576 r['QNames'].append('QName').data = self.rates.get(j, 'Voltage') 577 r['PressureNames'].append('PressureName').data = self.rates.get(i, 'Pressure') 578 r['PressureNames'].append('PressureName').data = self.rates.get(j, 'Pressure') 579 r['RateFormats'].append('RateFormat').data = self.rates.get(i, 'Format') 580 r['RateFormats'].append('RateFormat').data = self.rates.get(j, 'Format') 581 for i in xrange(self.constraints_kin.size): 582 constraints.append(self.constraints_kin.get_for_qmf(i)) 583 f['ChannelCount'].data = self.channelCount 584 f['VRev'].data = self.vRev 585 f['BalanceLoops'].data = self.balance_loops 586 f['Amps'].data = [self.classes.get(i, 'Amp') for i in xrange(self.classes.size)] 587 f['Stds'].data = [self.classes.get(i, 'Std') for i in xrange(self.classes.size)] 588 f['Conds'].data = [self.classes.get(i, 'Cond') for i in xrange(self.classes.size)] 589 f['CondStds'].data = [self.classes.get(i, 'CondStd') for i in xrange(self.classes.size)] 590 append_clone_if(f, self.constraints_ampvar) 591 append_clone_if(f, self.extra_kin_pars) 592 if not append_clone_if(f, self.nAr): 593 f['NAr'].data = self.classes.size * [0] 594 append_clone_if(f, self.ars) 595 append_clone_if(f, self.properties_node) 596 f['Properties']['IeqFv'].data = int(self._IeqFv) 597 f['Properties']['PreFormat'].data = self._k0format 598 f['Properties']['ExpFormat'].data = self._k1format 599 return f

600 - def from_tree(self, f):

601 """Reads in the model from f, a {qubx.tree.Node_numpy} formatted as U{QMF<http://www.qub.buffalo.edu/qubdoc/files/qmf.html>}.""" 602 self.__reading = True 603 self.clear() 604 605 amps, stds = f['Amps'], f['Stds'] 606 conds, condstds = f['Conds'], f['CondStds'] 607 vrevs = f['VRevs'] 608 nclass = min(node.data.count for node in [amps, stds]) 609 for i in xrange(nclass): 610 self.classes.append({'Amp' :float(node_data_or_def(amps, float(i), i)), 611 'Std' :float(node_data_or_def(stds, 0.15, i)), 612 'Cond' :float(node_data_or_def(conds, 10.0*i, i)), 613 'CondStd' :float(node_data_or_def(condstds, 1.0*i, i)), 614 'vRev' :float(node_data_or_def(vrevs, UNSET_VALUE, i))}) 615 for state in qubx.tree.children(f['States'], 'State'): 616 self.states.append({'Class':node_data_or_def(state.find('Class'), 0), 617 'x' :node_data_or_def(state.find('x'), 50.0), 618 'y' :node_data_or_def(state.find('y'), 50.0), 619 'Pr' :node_data_or_def(state.find('Pr'), 0.0), 620 'Label':node_data_or_def(state.find('Label'), ""), 621 'Group':node_data_or_def(state.find('Gr'), 0)}) 622 for rate in qubx.tree.children(f['Rates'], 'Rate'): 623 states, k0, k1, k2, p, q, press, pnames, qnames, pressnames, fmts = [rate.find(nm) for nm in ('States', 'k0', 'k1', 'k2', 'P', 'Q', 'Pressure', 'PNames', 'QNames', 'PressureNames', 'RateFormats')] 624 pname = pnames.find('PName') 625 qname = qnames.find('QName') 626 pressname = pressnames.find('PressureName') 627 fmt = fmts.find('RateFormat') 628 for i in (0,1): 629 j = (i+1)%2 630 lig = node_data_or_def(p, 0, i) and str(pname.data) or "" 631 volt = node_data_or_def(q, 0, i) and str(qname.data) or "" 632 pressure = node_data_or_def(press, 0, i) and str(pressname.data) or "" 633 self.rates.append({'From' : node_data_or_def(states, i, i), 634 'To' : node_data_or_def(states, j, j), 635 'k0' : float(node_data_or_def(k0, 10.0, i)), 636 'k1' : float(node_data_or_def(k1, 0.0, i)), 637 'k2' : float(node_data_or_def(k2, 0.0, i)), 638 'Ligand' : lig, 639 'Voltage' : volt, 640 'Pressure' : pressure, 641 'Format' : str(fmt.data)}) 642 pname = pname.nextSameName() 643 qname = qname.nextSameName() 644 pressname = pressname.nextSameName() 645 fmt = fmt.nextSameName() 646 for cns in qubx.tree.children(f.find('Constraints')): 647 if cns.name == 'LoopBal': 648 self.balance_loops = True 649 elif cns.name == 'LoopImbal': 650 pass 651 else: 652 self.constraints_kin.append_from_qmf(cns) 653 self.channelCount = node_data_or_def(f.find('ChannelCount'), 1) 654 self.IeqFv = node_data_or_def(f.find('Properties').find('IeqFv'), self._IeqFv) 655 self.k0format = node_data_or_def(f.find('Properties').find('PreFormat'), self._k0format) 656 self.k1format = node_data_or_def(f.find('Properties').find('ExpFormat'), self._k1format) 657 self.vRev = node_data_or_def(f.find('VRev'), self._vRev) 658 self.balance_loops = node_data_or_def(f.find('BalanceLoops'), False) 659 660 # setting aside: 661 self.constraints_ampvar = f.find('ConstraintsAmpVar') 662 self.extra_kin_pars = f.find('ExtraKineticPars') 663 self.nAr = f.find('NAr') 664 self.ars = f.find('Ars') 665 self.properties_node = f.find('Properties') 666 667 # discarding: 668 # Constraints:*:HasValue, Value, Coeffs 669 670 self.__reading = False

671 - def from_mdl(self, f):

672 self.__reading = True 673 self.clear() 674 675 self.channelCount = int(re.search(r"([0-9]+).*hannels", f.next()).group(1)) 676 nclass = int(re.search(r"([0-9]+).*lass", f.next()).group(1)) 677 f.next() #header 678 for i in xrange(nclass): 679 match = re.search(r"(\S+)\s+(\S+)\s+(\S+)", f.next()) 680 self.classes.append({'Amp' :float(match.group(1)), 681 'Std' :float(match.group(2)), 682 'Cond' :10.0*i, 683 'CondStd' :1.0*i, 684 'Vrev' :UNSET_VALUE}) 685 nstate = int(re.search(r"([0-9]+).*tates", f.next()).group(1)) 686 f.next() #header 687 for i in xrange(nstate): 688 match = re.search(r"(\S+)\s+(\S+)\s+(\S+)\s+(\S+)", f.next()) 689 self.states.append({'Class' :int(match.group(2)), 690 'x' :float(match.group(3)), 691 'y' :float(match.group(4)), 692 'Pr' :float(match.group(1)), 693 'Label' :'', 694 'Group' :0}) 695 nrate = int(re.search(r"([0-9]+).*onnections", f.next()).group(1)) 696 for i in xrange(nrate): 697 match = re.search(r"(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)", f.next()) 698 a, b = int(match.group(1)) - 1, int(match.group(2)) - 1 699 self.rates.append({'From' :a, 700 'To' :b, 701 'k0' :float(match.group(3)), 702 'k1' :float(match.group(4)), 703 'Ligand' :int(match.group(5)) and "Ligand" or "", 704 'Voltage' :float(match.group(4)) and "Voltage" or ""}) 705 self.rates.append({'From' :b, 706 'To' :a, 707 'k0' :float(match.group(6)), 708 'k1' :float(match.group(7)), 709 'Ligand' :int(match.group(8)) and "Ligand" or "", 710 'Voltage' :float(match.group(7)) and "Voltage" or ""}) 711 ncns = int(re.search(r"([0-9]+).*onstrain", f.next()).group(1)) 712 f.next() 713 f.next() 714 for i in xrange(ncns): 715 ints = [int(tok) for tok in re.split(r"\s+", f.next())] 716 if ints[0] == 2: # LoopBal 717 self.balance_loops = True 718 else: 719 self.constraints_kin.append({'Type' : {0 : 'FixRate', 720 1 : 'ScaleRate', 721 # 2 : 'LoopBal', 722 3 : 'FixExp', 723 4 : 'ScaleExp'}[ints[0]], 724 'States' : [intt - 1 for intt in ints[1:]]}) 725 726 # ignoring Ar coeffs. 727 728 self.__reading = False

729 730

731 -class ConstraintsKinTable(qubx.table.SimpleTable):

732 - def __init__(self, label, global_name, rates):

733 qubx.table.SimpleTable.__init__(self, label, auto_add_fields=True, global_name=global_name) 734 self.max_fields = 2048 735 self.add_field('Type', 'Generalized', acceptNothing, str, '', ['FixRate', 'FixExp', 'ScaleRate', 'ScaleExp', 'Generalized']) 736 self.add_field('States', (), acceptIntList(acceptIntGreaterThanOrEqualTo(0)), 737 lambda ll: ' '.join(str(l) for l in ll), '') 738 self.add_field('Sum', 0.0, acceptEval(), '%.3g', '') 739 self.OnEdit = WeakEvent() # () 740 self.__in_edit = False 741 self.__v_sens = self.__p_sens = 0 742 self.__ref = Reffer() 743 self.OnAfterInsert += self.__ref(self.__onAfterInsert) 744 self.OnRemoved += self.__ref(self.__onRemoved) 745 self.OnSet += self.__ref(self.__onSet) 746 self.rates = rates 747 self.rates.OnInsert += self.__ref(self.__onInsertRate) 748 self.rates.OnRemoving += self.__ref(self.__onRemovingRate) 749 self.rates.OnRemoved += self.__ref(self.__onRemovedRate) 750 self.rates.OnSet += self.__ref(self.__onSetRate) 751 self.__balance_loops = False 752 for i in xrange(rates.size): 753 self.__onInsertRate(i, undoing=False)

754 - def dispose(self):

755 qubx.table.SimpleTable.dispose(self) 756 self.__ref.clear()

757 - def set_balance_loops(self, x):

758 if self.__balance_loops == x: return 759 self.__balance_loops = x 760 if x: 761 self.build_loop_constraints()

762 balance_loops = property(lambda self: self.__balance_loops, lambda self, x: self.set_balance_loops(x))

763 - def build_loop_constraints(self):

764 """Finds the fundamental cycles in the model and constrains them all to detailed balance.""" 765 self.loop_constraints = [] 766 # find edges 767 edges = [tuple([self.rates.get(i, field) for field in ('From', 'To')]) 768 for i in xrange(self.rates.size)] 769 edges = [(fr, to) for fr, to in edges if fr < to] 770 self.loop_constraints = [('LoopBal', loop) for loop in qubx.fast.model.find_loops(edges)] 771 loop_rates = [] 772 for name, states in self.loop_constraints: 773 ns = len(states) 774 rates = [] 775 for i,s in enumerate(states): 776 s2 = states[(i+1)%ns] 777 rates.extend([(s, s2), (s2, s)]) 778 loop_rates.extend(rates) 779 self.loop_rates = set(loop_rates)

780 - def append_from_qmf(self, cns):

781 typ = cns.name 782 states = cns.data[:] 783 rhs = cns.find('Value').data and cns['Value'].data[0] or 0.0 784 coeffs = cns.find('Coeffs').data and cns.find('Coeffs').data[:] or [] 785 entry = {'Type' : typ, 'States' : states, 'Sum' : rhs} 786 if coeffs: 787 try: 788 for r in xrange(self.rates.size): 789 sa, sb = self.rates.get(r, 'From'), self.rates.get(r, 'To') 790 entry[fieldname_k(0, sa, sb)] = coeffs[2*r] 791 if coeffs[2*r+1]: 792 entry[fieldname_k(1, sa, sb)] = coeffs[2*r+1] 793 except: 794 traceback.print_exc() 795 #pass 796 self.append(entry)

797 - def get_for_qmf(self, i):

798 row = self.get_row(i) 799 cns = qubx.tree.Node(row.Type) 800 cns.data = row.States 801 if row.Type == 'Generalized': 802 coeffs = [] 803 for r in xrange(self.rates.size): 804 sa, sb = self.rates.get(r, 'From'), self.rates.get(r, 'To') 805 coeffs.append(row.__dict__[fieldname_k(0, sa, sb)]) 806 coeffs.append(self.__v_sens and row.__dict__[fieldname_k(1, sa, sb)] or 0.0) 807 cns['HasValue'].data = 1 808 cns['Value'].data = row.Sum 809 cns['Coeffs'].data = coeffs 810 return cns

811 - def get_matrices(self, K0, K1, L, V, ix=None, auto=True):

812 Nstate = K0.shape[0] 813 K2 = zeros(dtype='float64', shape=(Nstate,Nstate)) 814 P = zeros(dtype='int32', shape=(Nstate,Nstate)) 815 return self.get_matrices_p(K0, K1, K2, L, V, P, ix, auto)

816 - def get_matrices_p(self, K0, K1, K2, L, V, P, ix=None, auto=True):

817 Nstate = K0.shape[0] 818 if not (ix is None): 819 ixs = K_index(Nstate) 820 ixo = dict([(ab, i) for i, ab in enumerate(ixs)]) 821 N = len(ixs) 822 A = zeros(dtype='float64', shape=(1,3*N)) 823 for a in xrange(Nstate): 824 for b in xrange(Nstate): 825 if a != b: 826 fieldname = fieldname_k(0, a, b) 827 if fieldname in self.fields: 828 Ai = ixo[(a, b)] 829 A[0,Ai] = self.get(ix, fieldname) 830 if self.__v_sens: 831 A[0,N+Ai] = self.get(ix, fieldname_k(1, a, b)) 832 if self.__p_sens: 833 A[0,2*N+Ai] = self.get(ix, fieldname_k(2, a, b)) 834 return A, self.get(ix, 'Sum') 835 else: 836 rows, sums = [], [] 837 if auto: 838 Aauto, Bauto = AutoConstraints(K0, K1, K2, L, V, P) 839 if Aauto.shape[0]: 840 rows, sums = [Aauto], [Bauto] 841 for i in xrange(self.size): 842 row, rhs = self.get_matrices_p(K0, K1, K2, L, V, P, i) 843 rows.append(row) 844 sums.append(rhs) 845 if self.balance_loops: 846 Aloop, Bloop = UserConstraints(K0, K1, K2, L, V, P, self.loop_constraints) 847 if Aloop.shape[0]: 848 rows.append(Aloop) 849 sums.append(Bloop) 850 if not rows: 851 return None, None 852 return vstack(rows), vstack(sums)

853 - def list_involving(self, sa, sb, field):

854 """Returns a list of indices of constraints which include the rate from sa to sb. 855 856 @param field: either 'k0' or 'k1' 857 """ 858 cns = [] 859 for i in xrange(self.size): 860 typ = self.get(i, 'Type') 861 sts = self.get(i, 'States') 862 if typ == 'FixRate': 863 if (field == 'k0') and (sa == sts[0]) and (sb == sts[1]): 864 cns.append(i) 865 continue 866 elif typ == 'FixExp': 867 if (field == 'k1') and (sa == sts[0]) and (sb == sts[1]): 868 cns.append(i) 869 continue 870 elif typ == 'FixPress': 871 if (field == 'k2') and (sa == sts[0]) and (sb == sts[1]): 872 cns.append(i) 873 continue 874 elif typ == 'ScaleRate': 875 if (field == 'k0'): 876 if ((sa == sts[0]) and (sb == sts[1])) or ((sa == sts[2]) and (sb == sts[3])): 877 cns.append(i) 878 continue 879 elif typ == 'ScaleExp': 880 if (field == 'k1'): 881 if ((sa == sts[0]) and (sb == sts[1])) or ((sa == sts[2]) and (sb == sts[3])): 882 cns.append(i) 883 continue 884 elif typ == 'ScalePress': 885 if (field == 'k2'): 886 if ((sa == sts[0]) and (sb == sts[1])) or ((sa == sts[2]) and (sb == sts[3])): 887 cns.append(i) 888 continue 889 elif typ == 'Generalized': 890 ki = {'k0':0, 'k1':1, 'k2':2}[field] 891 if self.get(i, fieldname_k(ki, sa, sb)): 892 cns.append(i) 893 continue 894 return cns

895 - def __update_v_sens(self, delta):

896 v_prev, self.__v_sens = self.__v_sens, self.__v_sens+delta 897 if self.__v_sens and not v_prev: 898 for r in xrange(self.rates.size): 899 sa, sb = self.rates.get(r, 'From'), self.rates.get(r, 'To') 900 self.add_field(fieldname_k1(sa, sb), self.auto_default, self.auto_accept, self.auto_format, '') 901 elif v_prev and not self.__v_sens: 902 for r in xrange(self.rates.size): 903 sa, sb = self.rates.get(r, 'From'), self.rates.get(r, 'To') 904 self.remove_field(fieldname_k1(sa, sb))

905 - def __update_p_sens(self, delta):

906 p_prev, self.__p_sens = self.__p_sens, self.__p_sens+delta 907 if self.__p_sens and not p_prev: 908 for r in xrange(self.rates.size): 909 sa, sb = self.rates.get(r, 'From'), self.rates.get(r, 'To') 910 self.add_field(fieldname_k2(sa, sb), self.auto_default, self.auto_accept, self.auto_format, '') 911 elif p_prev and not self.__p_sens: 912 for r in xrange(self.rates.size): 913 sa, sb = self.rates.get(r, 'From'), self.rates.get(r, 'To') 914 self.remove_field(fieldname_k2(sa, sb))

915 - def __update_coeffs(self, i):

916 coeffs = {} 917 def sts_at_size(ns=None): 918 sts = list(self.get(i, 'States')) 919 if (not (ns is None)) and (ns != len(sts)): 920 while len(sts) < ns: 921 sts.append(len(sts)) 922 if len(sts) > ns: 923 del sts[ns:] 924 self.set(i, 'States', tuple(sts)) 925 return sts

926 typ = self.get(i, 'Type') 927 rr = self.rates 928 rate_connecting = self.rate_connecting 929 try: 930 if typ == 'FixRate': 931 sa, sb = sts_at_size(2) 932 self.set(i, 'Sum', log(rr.get(rate_connecting[(sa, sb)], 'k0'))) 933 coeffs[(0, sa, sb)] = 1.0 934 elif typ == 'FixExp': 935 sa, sb = sts_at_size(2) 936 self.set(i, 'Sum', rr.get(rate_connecting[(sa, sb)], 'k1')) 937 coeffs[(1, sa, sb)] = 1.0 938 elif typ == 'FixPress': 939 sa, sb = sts_at_size(2) 940 self.set(i, 'Sum', rr.get(rate_connecting[(sa, sb)], 'k2')) 941 coeffs[(2, sa, sb)] = 1.0 942 elif typ == 'ScaleRate': 943 sts = sts_at_size(4) 944 self.set(i, 'Sum', (log(rr.get(rate_connecting[(sts[0], sts[1])], 'k0')) 945 - log(rr.get(rate_connecting[(sts[2], sts[3])], 'k0')))) 946 coeffs[(0, sts[0], sts[1])] = 1.0 947 coeffs[(0, sts[2], sts[3])] = -1.0 948 elif typ == 'ScaleExp': 949 sts = sts_at_size(4) 950 self.set(i, 'Sum', 0.0) 951 coeffs[(1, sts[0], sts[1])] = 1.0 952 coeffs[(1, sts[2], sts[3])] = - (rr.get(rate_connecting[(sts[0], sts[1])], 'k1') 953 / rr.get(rate_connecting[(sts[2], sts[3])], 'k1')) 954 elif typ == 'ScalePress': 955 sts = sts_at_size(4) 956 self.set(i, 'Sum', 0.0) 957 coeffs[(2, sts[0], sts[1])] = 1.0 958 coeffs[(2, sts[2], sts[3])] = - (rr.get(rate_connecting[(sts[0], sts[1])], 'k2') 959 / rr.get(rate_connecting[(sts[2], sts[3])], 'k2')) 960 except: 961 traceback.print_exc() 962 #pass 963 if coeffs: 964 for field in self.fields: 965 if not (field in set(('Index', 'Type', 'States', 'Sum', 'Group'))): 966 if fieldinfo_k(field) in coeffs: 967 val = coeffs[fieldinfo_k(field)] 968 if abs(val) < 1e-12: 969 val = 0.0 970 self.set(i, field, val) 971 elif self.get(i, field): 972 self.set(i, field, 0.0)

973 - def __onAfterInsert(self, i, undoing):

974 if self.__in_edit: return 975 self.__in_edit = True 976 typ = self.get(i, 'Type') 977 if typ != 'Generalized': 978 self.__update_coeffs(i) 979 self.__in_edit = False 980 self.OnEdit(not (typ in ('FixRate', 'FixExp', 'ScaleRate', 'ScaleExp')))

981 - def __onRemoved(self, index, undoing):

982 self.OnEdit(False)

983 - def __onSet(self, i, field, val, prev, undoing):

984 if self.__in_edit: return 985 self.__in_edit = True 986 if self.get(i, 'Type') != 'Generalized': 987 if field in ('Type', 'States'): 988 self.__update_coeffs(i) 989 elif not (field in ('Index', 'Group')): 990 self.set(i, 'Type', 'Generalized') 991 self.__in_edit = False

992 # self.OnEdit(False?) prevent nonsense enforcement of half-entered constraint

993 - def __onInsertRate(self, i, undoing):

994 if self.__in_edit: return 995 self.__in_edit = True 996 if self.rates.get(i, 'Voltage'): 997 self.__update_v_sens(1) 998 if self.rates.get(i, 'Pressure'): 999 self.__update_p_sens(1) 1000 sa, sb = self.rates.get(i, 'From'), self.rates.get(i, 'To') 1001 self.add_field(fieldname_k0(sa, sb), self.auto_default, self.auto_accept, self.auto_format, '') 1002 if self.__v_sens: 1003 self.add_field(fieldname_k1(sa, sb), self.auto_default, self.auto_accept, self.auto_format, '') 1004 if self.__p_sens: 1005 self.add_field(fieldname_k2(sa, sb), self.auto_default, self.auto_accept, self.auto_format, '') 1006 self.__in_edit = False 1007 if self.__balance_loops: 1008 self.build_loop_constraints()

1009 - def __onRemovingRate(self, i, undoing):

1010 if self.__in_edit: return 1011 self.__in_edit = True 1012 changed = False 1013 r = self.rates 1014 if r.get(i, 'Voltage'): 1015 self.__update_v_sens(-1) 1016 if r.get(i, 'Pressure'): 1017 self.__update_p_sens(-1) 1018 sa, sb = r.get(i, 'From'), r.get(i, 'To') 1019 self.remove_field(fieldname_k0(sa, sb)) 1020 if self.__v_sens: 1021 self.remove_field(fieldname_k1(sa, sb)) 1022 if self.__p_sens: 1023 self.remove_field(fieldname_k2(sa, sb)) 1024 # remove invalid constraints 1025 for j in reversed(xrange(self.size)): 1026 typ = self.get(j, 'Type') 1027 sts = self.get(j, 'States') 1028 if (typ == 'FixRate') or (typ == 'FixExp'): 1029 if (sa == sts[0]) and (sb == sts[1]): 1030 self.remove(j) 1031 changed = True 1032 continue 1033 elif (typ == 'ScaleRate') or (typ == 'ScaleExp'): 1034 if ((sa == sts[0]) and (sb == sts[1])) or ((sa == sts[2]) and (sb == sts[3])): 1035 self.remove(j) 1036 changed = True 1037 continue 1038 self.__in_edit = False 1039 if changed: 1040 self.OnEdit(False)

1041 - def __onRemovedRate(self, i, undoing):

1042 if self.__balance_loops: 1043 self.build_loop_constraints()

1044 - def __onSetRate(self, i, field, val, prev, undoing):

1045 if self.__in_edit: return 1046 self.__in_edit = True 1047 if field == 'Voltage': 1048 if val and not prev: 1049 self.__update_v_sens(1) 1050 elif prev and not val: 1051 self.__update_v_sens(-1) 1052 if field == 'Pressure': 1053 if val and not prev: 1054 self.__update_p_sens(1) 1055 elif prev and not val: 1056 self.__update_p_sens(-1) 1057 if field in ('k0', 'k1', 'k2'): # may alter sum, coeffs: 1058 for j in xrange(self.size): 1059 if self.get(j, 'Type') != 'Generalized': 1060 self.__update_coeffs(j) 1061 self.__in_edit = False

1062

1063 -def fieldname_k0_f(sa, sb):

1064 return 'ln(k0),%i,%i'%(sa, sb)

1065 fieldname_k0 = memoize(fieldname_k0_f) 1066

1067 -def fieldname_k1_f(sa, sb):

1068 return 'k1,%i,%i'%(sa, sb)

1069 fieldname_k1 = memoize(fieldname_k1_f) 1070

1071 -def fieldname_k2_f(sa, sb):

1072 return 'k2,%i,%i'%(sa, sb)

1073 fieldname_k2 = memoize(fieldname_k2_f) 1074

1075 -def fieldname_k_f(ki, sa, sb):

1076 if ki == 0: 1077 return fieldname_k0_f(sa, sb) 1078 elif ki == 1: 1079 return fieldname_k1_f(sa, sb) 1080 else: 1081 return fieldname_k2_f(sa, sb)

1082 fieldname_k = memoize(fieldname_k_f) 1083

1084 -def fieldinfo_k_f(s):

1085 fields = re.search(r"(.*)k([012])[^,]*,([^,]+),(.*)", s) 1086 return fields and (int(fields.group(2)), int(fields.group(3)), int(fields.group(4))) or (0,0,0)

1087 fieldinfo_k = memoize(fieldinfo_k_f) 1088

1089 -def append_clone_if(parent, node):

1090 """if node: parent.appendClone(node)""" 1091 if not node.isNull: 1092 parent.appendClone(node) 1093 return True 1094 return False

1095

1096 -def dec_if_gt(n):

1097 """if x > n: return x-1; else: return x.""" 1098 def dec(x): 1099 if x > n: 1100 return x-1 1101 return x

1102 return dec 1103

1104 -def NewQubModel():

1105 """Returns a new L{QubModel} with two states.""" 1106 m = QubModel() 1107 m.classes.append({'Amp' : 0.0, 'Std' : 0.1, 'Cond' : 0.0, 'CondStd' : 0.0}) 1108 m.classes.append({'Amp' : 1.0, 'Std' : 0.15, 'Cond' : 10.0, 'CondStd' : 1.0}) 1109 m.states.append({'Class' : 0, 'x' : 30.0, 'y' : 50.0}) 1110 m.states.append({'Class' : 1, 'x' : 70.0, 'y' : 50.0}) 1111 m.rates.append({'From' : 0, 'To' : 1, 'k0' : 10.0}) 1112 m.rates.append({'From' : 1, 'To' : 0, 'k0' : 10.0}) 1113 m.undoStack.clear() 1114 m.changed = False 1115 return m

1116

1117 -def OpenQubModel(path):

1118 """Reads a U{QMF<http://www.qub.buffalo.edu/qubdoc/files/qmf.html>} file from disk and returns a L{QubModel}.""" 1119 m = QubModel() 1120 ext = os.path.splitext(path)[1].lower() 1121 if ext == '.qmf': 1122 f = qubx.tree.Open(path, True) 1123 f.close() 1124 if f.isNull: 1125 raise IOError("Can't open %s" % path) 1126 if f.name != 'ModelFile': 1127 raise IOError("%s is not a qub model file." % path) 1128 m.from_tree(f) 1129 elif ext == '.qmj': 1130 m.from_json(json.loads(open(path, 'r').read())) 1131 else: 1132 try: 1133 m.from_mdl(open(path, 'r')) 1134 except: 1135 raise IOError("%s is not a qub model file." % path) 1136 m.undoStack.clear() 1137 m.changed = False 1138 m.path = path 1139 return m

1140 1141

1142 -class QubModels(qubx.table.ObjectTable):

1143 """Table-like collection of open model files."""

1144 - def __init__(self):

1145 qubx.table.ObjectTable.__init__(self, ' Models', global_name='QubX.Models.table') 1146 self.add_field('Name', "", acceptNothing, '%s', '', get=self.get_name, set=None) 1147 self.add_field('Channel Count', 1, acceptIntGreaterThan(0), '%s', '', get=self.get_channelCount, set=self.set_channelCount) 1148 self.add_field('vRev', 0.0, acceptFloat, '%.3g', 'mV', get=self.get_vRev, set=self.set_vRev) 1149 self.add_field('In Folder', '', acceptNothing, '%s', '', get=self.get_folder, set=None) 1150 self.add_field('k0 Format', '', str, str, '', get=self.get_k0format, set=self.set_k0format) 1151 self.add_field('k1 Format', '', str, str, '', get=self.get_k1format, set=self.set_k1format)

1152 - def get_name(self, obj):

1153 return os.path.split(obj.path)[1]

1154 - def get_folder(self, obj):

1155 return os.path.split(obj.path)[0]

1156 - def get_channelCount(self, obj):

1157 return obj.channelCount

1158 - def set_channelCount(self, obj, val):

1159 obj.channelCount = val

1160 - def get_vRev(self, obj):

1161 return obj.vRev

1162 - def set_vRev(self, obj, val):

1163 obj.vRev = val

1164 - def get_k0format(self, obj):

1165 return obj.k0format

1166 - def set_k0format(self, obj, val):

1167 obj.k0format = val

1168 - def get_k1format(self, obj):

1169 return obj.k1format

1170 - def set_k1format(self, obj, val):

1171 obj.k1format = val

1172 1173

1174 -def ModelToClazz(model):

1175 """Returns the array([class of state s, for s in (0..nstate-1)]).""" 1176 return array([model.states.get(s, 'Class') for s in xrange(model.states.size)])

1177 1178 1179

1180 -def ModelToRateMatrices(model, signal_names=[]):

1181 K0, K1, K2, L, V, P = ModelToRateMatricesP(model, signal_names) 1182 return K0, K1, L, V

1183

1184 -def ModelToRateMatricesP(model, signal_names=[]):

1185 """Returns the (nstate x nstate) matrices (K0, K1, K2, Ligand, Voltage, Pressure).""" 1186 signal_ix = dict([(nm, ix) for ix, nm in enumerate(signal_names)]) 1187 nstate = model.states.size 1188 K0 = matrix(zeros(shape=(nstate, nstate))) 1189 K1 = matrix(zeros(shape=(nstate, nstate))) 1190 K2 = matrix(zeros(shape=(nstate, nstate))) 1191 L = matrix(zeros(shape=(nstate, nstate), dtype='int32')) 1192 V = matrix(zeros(shape=(nstate, nstate), dtype='int32')) 1193 P = matrix(zeros(shape=(nstate, nstate), dtype='int32')) 1194 rr = model.rates 1195 for r in xrange(rr.size): 1196 a = rr.get(r, 'From') 1197 b = rr.get(r, 'To') 1198 K0[a,b] = rr.get(r, 'k0') 1199 K1[a,b] = rr.get(r, 'k1') 1200 K2[a,b] = rr.get(r, 'k2') 1201 nm = rr.get(r, 'Ligand') 1202 try: 1203 if nm: L[a,b] = signal_ix[nm] or 0 1204 except KeyError: 1205 pass # print 'Ignoring %s-dependence on %i->%i: no such signal.' % (nm, a, b) 1206 nm = rr.get(r, 'Voltage') 1207 try: 1208 if nm: V[a,b] = signal_ix[nm] or 0 1209 except KeyError: 1210 pass # print 'Ignoring %s-dependence on %i->%i: no such signal.' % (nm, a, b) 1211 nm = rr.get(r, 'Pressure') 1212 try: 1213 if nm: P[a,b] = signal_ix[nm] or 0 1214 except KeyError: 1215 pass # print 'Ignoring %s-dependence on %i->%i: no such signal.' % (nm, a, b) 1216 return K0, K1, K2, L, V, P

1217

1218 -def K01ToModel(K0, K1, model):

1219 """Copies rate constants from matrices K0 and K1 into model.rates.""" 1220 model.should_enforce_constraints = False 1221 rr = model.rates 1222 for r in xrange(rr.size): 1223 a = rr.get(r, 'From') 1224 b = rr.get(r, 'To') 1225 rr.set(r, 'k0', K0[a,b]) 1226 rr.set(r, 'k1', K1[a,b]) 1227 model.should_enforce_constraints = True

1228

1229 -def K012ToModel(K0, K1, K2, model):

1230 """Copies rate constants from matrices K0 and K1 into model.rates.""" 1231 model.should_enforce_constraints = False 1232 rr = model.rates 1233 for r in xrange(rr.size): 1234 a = rr.get(r, 'From') 1235 b = rr.get(r, 'To') 1236 rr.set(r, 'k0', float(K0[a,b])) 1237 rr.set(r, 'k1', float(K1[a,b])) 1238 rr.set(r, 'k2', float(K2[a,b])) 1239 model.should_enforce_constraints = True

1240

1241 -def ModelToFast(model, signal_names=[], on_status=sys.stdout.write):

1242 """Returns a new L{qubx.fast.model.Model}.""" 1243 clazz = ModelToClazz(model) 1244 Nclass = int(max(clazz) + 1) if clazz.shape[0] else 0 1245 fast = qubx.fast.model.Model(model.states.size, Nclass, clazz, on_status) 1246 if model.states: 1247 fast.P0[:] = ModelToP0(model) 1248 fast.IeqFv = model.IeqFv 1249 fast.Nchannel = model.channelCount 1250 fast.Vrev = model.vRev 1251 if model.IeqFv: 1252 for i in xrange(Nclass): 1253 fast.amp[i] = model.classes.get(i, 'Cond') 1254 fast.var[i] = model.classes.get(i, 'CondStd')**2 1255 else: 1256 for i in xrange(Nclass): 1257 fast.amp[i] = model.classes.get(i, 'Amp') 1258 fast.var[i] = model.classes.get(i, 'Std')**2 1259 if model.classes: 1260 fast.amp0 = model.classes.get(0, 'Amp') 1261 fast.var0 = model.classes.get(0, 'Std')**2 1262 if model.rates: 1263 fast.K0[:], fast.K1[:], fast.K2[:], fast.L[:], fast.V[:], fast.P[:] = ModelToRateMatricesP(model, signal_names) 1264 return fast

1265 1266 1267

1268 -def RandomWeighted(weights):

1269 """Returns a random index into weights (number between 0 and len(weights)-1), with likelihood proportional to weights[i].""" 1270 r = random.uniform() # in [0, 1) 1271 tot = weights[0] 1272 i = 0 1273 last = len(weights)-1 1274 while (tot <= r) and (i < last): 1275 i += 1 1276 tot += weights[i] 1277 return i

1278

1279 -def RandomLT(high):

1280 """Returns a uniform random number between 0 and high.""" 1281 return random.uniform(high=high)

1282

1283 -class RandomNormalsSmallN(object):

1284 """Generates batches of gaussian data; mean 0, std 1."""

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

1286 self.count = buffer_size 1287 self.new_buffer()

1288 - def new_buffer(self):

1289 self.norms = random.randn(self.count) 1290 self.i = 0

1291 - def next(self, N):

1292 """Returns an array or slice containing N gaussian random numbers.""" 1293 if N > self.count: 1294 return random.randn(N) 1295 if self.i + N > self.count: 1296 self.new_buffer() 1297 last, self.i = self.i, self.i + N 1298 return self.norms[last:self.i]

1299 1300 1301

1302 -def GetNewState(q, state):

1303 """Returns the next state, given the state it's coming from.""" 1304 nstate = q.shape[0] 1305 rand = RandomLT( -q[state, state] ) 1306 tot = 0.0 1307 for s in xrange(nstate): 1308 if s != state: 1309 tot += q[state, s] 1310 if tot > rand: 1311 return s 1312 return state

1313

1314 -def GenerateEvents(model):

1315 """Yields (state, duration) of a Markov process, starting from fixed probabilities (model.states:Pr).""" 1316 p = ModelToP0(model) 1317 q = ModelToQ(model) 1318 s = RandomWeighted(p) 1319 for x in GenerateEventsX(q, s): 1320 yield x

1321

1322 -def GenerateEventsX(q, start_state):

1323 """Yields (state, duration) of a Markov process, starting from a particular state.""" 1324 s = start_state 1325 while True: 1326 dur = random.exponential()/(- q[s, s]) 1327 yield s, dur 1328 s = GetNewState(q, s)

1329

1330 -def GenerateSampleBunchesX(events, sampling, amps, stds):

1331 """Yields (sample_array, state_index), one event at a time from a hidden Markov process. 1332 1333 @param amps: mean current of each state 1334 @param stds: standard deviation of current in each state 1335 """ 1336 sam = 0.0 1337 norms = RandomNormalsSmallN() 1338 for state, dur in events: 1339 sam += dur/sampling 1340 count = int(floor(sam)) 1341 sam -= count 1342 samples = norms.next(count) * stds[state] 1343 samples += amps[state] 1344 yield samples, state

1345

1346 -def GenerateSampleBunches(model, sampling, events):

1347 """Yields (sample_array, state_index), one event at a time, from a hidden Markov process.""" 1348 sam = 0.0 1349 amp = [model.classes.get(model.states.get(i, 'Class'), 'Amp') for i in xrange(model.states.size)] 1350 std = [model.classes.get(model.states.get(i, 'Class'), 'Std') for i in xrange(model.states.size)] 1351 norms = RandomNormalsSmallN() 1352 for state, dur in events: 1353 sam += dur/sampling 1354 count = int(floor(sam)) 1355 sam -= count 1356 samples = norms.next(count) * std[state] 1357 samples += amp[state] 1358 yield samples

1359

1360 -def GenerateSamples(model, sampling, events):

1361 """Yields one sample at a time from a hidden Markov process.""" 1362 sam = 0.0 1363 cls_amp = [model.classes.get(i, 'Amp') for i in xrange(model.classes.size)] 1364 cls_std = [model.classes.get(i, 'Std') for i in xrange(model.classes.size)] 1365 # minus baseline and background noise, for channel addition: 1366 cls_amp = [a - cls_amp[0] for a in cls_amp] 1367 cls_std = [s - cls_std[0] for s in cls_std] 1368 amp = [cls_amp[model.states.get(i, 'Class')] for i in xrange(model.states.size)] 1369 std = [cls_std[model.states.get(i, 'Class')] for i in xrange(model.states.size)] 1370 norms = RandomNormalsSmallN() 1371 for state, dur in events: 1372 sam += dur/sampling 1373 count = int(floor(sam)) 1374 sam -= count 1375 if state: 1376 for value in (amp[state] + (norms.next(count) * std[state])): 1377 yield value 1378 else: 1379 for i in xrange(count): 1380 yield 0.0

1381 1382 1383 1384

1385 -def ModelToP0(model):

1386 """Returns the array[nstate] of fixed starting probabilities.""" 1387 p = zeros(shape=(model.states.size,)) 1388 for s in xrange(model.states.size): 1389 p[s] = model.states.get(s, 'Pr') 1390 tot = sum(p) 1391 if tot: 1392 p *= 1/tot 1393 else: 1394 for s in xrange(model.states.size): 1395 p[s] = 1.0/model.states.size 1396 return p

1397 1398

1399 -def ModelToQ(model, stimulus=None):

1400 """Returns the matrix Q of combined rate constants.""" 1401 nstate = model.states.size 1402 q = matrix(zeros(shape=(nstate, nstate))) 1403 rr = model.rates 1404 for r in xrange(rr.size): 1405 rate = rr.get(r, 'k0') 1406 if stimulus: 1407 lig, volt, press = rr.get(r, 'Ligand'), rr.get(r, 'Voltage'), rr.get(r, 'Pressure') 1408 if lig and stimulus.has_key(lig): 1409 rate *= stimulus[lig] 1410 kExp = 0.0 1411 if volt and stimulus.has_key(volt): 1412 kExp += stimulus[volt] * rr.get(r, 'k1') 1413 if press and stimulus.has_key(press): 1414 kExp += stimulus[press] * rr.get(r, 'k2') 1415 rate *= exp(kExp) 1416 if not rate: 1417 rate = MIN_TOTAL_RATE 1418 q[rr.get(r, 'From'), rr.get(r, 'To')] = rate 1419 for s in xrange(nstate): 1420 q[s, s] = - sum(q[s]) 1421 return q

1422 1423

1424 -def QtoA(Q, dt):

1425 """Returns the matrix A = exp( dt * Q ).""" 1426 #Q = matrix(Q, copy=False) 1427 return scipy.linalg.matfuncs.expm( dt * Q )

1428

1429 -def QtoPe(Q):

1430 """Returns the array[nstate] of equilibrium probabilities.""" 1431 # S = [Q | 1] 1432 S = matrix(ones(shape=(Q.shape[0], Q.shape[1]+1))) 1433 S[:,:-1] = Q 1434 u = matrix(ones(shape=(1,Q.shape[0]))) 1435 p = u * (S*S.T).I 1436 return array(p).reshape((Q.shape[0],))

1437 1438 1439 # begin_html 1440 # "Qe" aka \({}^eQ\) is the "apparent" rate constant matrix, given that events 1441 # with duration <= tdead are not recorded. MIL computes the probability of 1442 # staying in class <b>a</b> for duration <b>t</b> using the submatrix \({}^eQ_{aa}\) 1443 # \[A(a, t) = e^{{}^eQ_{aa} t'}\] 1444 # end_html 1445

1446 -def QtoQe(Q, clazz, td):

1447 """Returns the apparent rate constant matrix, given that events with duration <= td are not recorded.""" 1448 # z = a-bar (not a) 1449 expm = scipy.linalg.matfuncs.expm 1450 eQ = matrix(zeros(shape=Q.shape)) 1451 for a in set(clazz): 1452 # partitions: 1453 A = clazz == a 1454 Z = clazz != a 1455 Qaa = Q[ix_(A, A)] 1456 Qaz = Q[ix_(A, Z)] 1457 Qza = Q[ix_(Z, A)] 1458 Qzz = Q[ix_(Z, Z)] 1459 Iz = identity(Qzz.shape[0]) 1460 # staying and returning: 1461 missed_returning = ((Qaz * (Iz - expm(td*Qzz))) * Qzz.I) * Qza 1462 eQ[ix_(A, A)] = Qaa - missed_returning 1463 1464 for b in (set(clazz) - set([a])): 1465 # partitions 1466 B = clazz == b 1467 C = (clazz != a) * (clazz != b) 1468 Qab = Q[ix_(A, B)] 1469 Qbb = Q[ix_(B, B)] 1470 if any(C): 1471 Qac = Q[ix_(A, C)] 1472 Qcc = Q[ix_(C, C)] 1473 Qcb = Q[ix_(C, B)] 1474 Ic = identity(Qcc.shape[0]) 1475 # switching class: 1476 escaped_indirect = ((Qac * (Ic - expm(td * Qcc))) * Qcc.I) * Qcb 1477 else: 1478 escaped_indirect = matrix(zeros(shape=Qab.shape)) 1479 # as published: 1480 # eQ[ix_(A,B)] = expm(td * missed_returning) * (Qab - escaped_indirect) 1481 # as coded, since before 2000: 1482 eQ[ix_(A,B)] = (Qab - escaped_indirect) * expm(td * Qbb) 1483 return eQ

1484 1485

1486 -def MeanFirstPassage(A, Peq, dt):

1487 """ 1488 Returns the matrix M of mean first passage times M[i,j] between states i and j. 1489 1490 >>> from http://www.scopenvironment.org/downloadpubs/scope34/contents.html 1491 1492 Z = inv(I - A + col(1)*row(pEq)) 1493 FL[i][j] = (Z[j][j] - Z[i][j]) / pEq[j] 1494 1495 @param A: sampled transition probability matrix A = expm(dt * Q) 1496 @param Peq: equilibrium probability vector 1497 @param dt: sampling interval in seconds 1498 """ 1499 n = A.shape[0] 1500 I = identity(n) 1501 U = matrix(ones(shape=(n,1))) 1502 P = matrix(Peq).reshape(1,n) 1503 Z = (identity(n) - A + U*P).I 1504 FL = zeros(shape=(n,n)) 1505 for i in xrange(n): 1506 for j in xrange(n): 1507 FL[i][j] = dt * (Z[j,j] - Z[i,j]) / P[0,j] 1508 return FL

1509

Source Code for Module qubx.model