qubx_testing.qubx

1 """Base classes for QUB Express plugins, and sample implementations. 2 3 This unit is not required to make a plugin (see qubx.pyenv), 4 but is intended to simplify and explain the process. 5 6 To make a plugin: 7 * make a folder in qub-express's Plugins folder (e.g. "my_plugin") 8 * make a subfolder with a unique name for your plugin's package (e.g. "my_plugin/my_uniquestr_package"); 9 the my_plugin folder will be part of sys.path; if there is any other my_uniquestr_package in sys.path, it 10 may be loaded instead. 11 * copy this file into my_uniquestr_package; by putting it there, you're guaranteed to import this copy, 12 and not one from some other plugin. 13 * make a python module, e.g. "my_uniquestr_package/my_module.py", which starts with "from . import qubx_plugin", 14 and defines a subclass of one of the plugin base classes. 15 * make a python module, "my_uniquestr_package/__init__.py" along the lines of: 16 __all__ = ['qubx_plugin', 'my_module'] 17 or if you want, leave it blank, but make sure it exists. 18 * make these text files in the my_plugin folder: 19 - qubx_plugin_about.txt - information for the Admin:Plugins panel 20 - qubx_plugin_imports.txt - newline-separated list of modules, in the order they should be re-loaded 21 - qubx_plugin_init.py - python script executed on load; 22 typically creates and adds one or more "Face" to the tool/about panels. 23 - qubx_plugin_fini.py - python script executed on unload; typically removes the "Face"(s) 24 - qubx_plugin_name.txt - name of plugin, as shown in Admin:Plugins panel 25 26 For examples, see qubx_sample_plugin, qubx_macroscopic, and qubx_single_molecule. 27 """ 28 29 import gobject 30 import gtk 31 import pango 32 import traceback 33 34 import qubx.notebook 35 import qubx.optimize 36 import qubx.optimizeGTK 37 import qubx.pyenv 38 import qubx.trialset 39 from qubx.accept import * 40 from qubx.faces import Face 41 from qubx.GTK import Requestable, pack_item, pack_space, pack_label, pack_scrolled, pack_check, pack_button 42 from qubx.model import QubModel, ModelToQ, ModelToP0, ModelToClazz, QtoPe, K012ToModel, ModelToRateMatricesP 43 from qubx.task import Robot, Tasks 44 from qubx.util_types import Anon, Product, Reffer, WeakEvent 45 46 from itertools import izip, count 47 from numpy import array, matrix, zeros, any, sum 48 import scipy.linalg.matfuncs 49 expm = scipy.linalg.matfuncs.expm 50

51 -class AboutPlugin(Face):

52 - def __init__(self, label, global_name=""):

53 super(AboutPlugin, self).__init__(label, global_name=global_name)

54

55 -class ToolPlugin(Face):

56 """Base class for top-pane plugin interface with optional worker thread (robot). 57 58 @ivar robot: L{qubx.task.Robot} managing a background work thread 59 """ 60 61 # these names are visible in the Explore dialog (Admin:Scripts): 62 __explore_featured = ['robot', 'on_robot_begin', 'on_robot_end', 'on_exception', 'on_interrupt'] 63

64 - def __init__(self, label, global_name="", has_robot=True):

65 """ 66 @param label: name on this tab in the gui 67 @param has_robot: if True, starts a background thread called self.robot 68 """ 69 super(ToolPlugin, self).__init__(label, global_name=global_name) 70 self.__ref = Reffer() 71 self.set_size_request(200, 50) 72 73 if has_robot: 74 self.robot = Robot(label, self.__ref(self.__on_robot_begin), self.__ref(self.__on_robot_end)) 75 self.robot.OnException += self.__ref(self.on_exception) 76 self.robot.OnInterrupt += self.__ref(self.on_interrupt) 77 QubX = qubx.pyenv.env.globals['QubX'] 78 QubX.OnQuit += self.__ref(self.__qubx_finalize) 79 else: 80 self.robot = None

81 - def __qubx_finalize(self):

82 self.robot.stop()

83 - def __on_robot_begin(self):

84 Tasks.add_task(self.robot) 85 self.on_robot_begin()

86 - def __on_robot_end(self):

87 Tasks.remove_task(self.robot) 88 self.on_robot_end()

89 - def on_robot_begin(self):

90 """Override this method to respond when the robot starts working.""" 91 pass

92 - def on_robot_end(self):

93 """Override this method to respond when the robot finishes working.""" 94 pass

95 - def on_exception(self, task, typ, val, tb):

96 """Override this method to do something with exceptions thrown on robot thread; default: print to stdout.""" 97 traceback.print_exception(typ, val, tb)

98 - def on_interrupt(self, task, cancel):

99 """Override this method to customize the effect of the Stop button (e.g. set a flag); default: KeyboardInterrupt raised in robot thread.""" 100 # if you don't call cancel, qubx will throw a KeyboardInterrupt into the robot thread. 101 pass

102

103 -class OptToolPlugin(ToolPlugin):

104 """Base class for optimization plugins. 105 106 Override at minimum: robot_calculate_score, robot_get_pars, robot_set_pars. 107 Optionally, override robot_optimize to replace the built-in optimizers. 108 109 They are called robot_* because they run on the worker (robot) thread. When called, 110 the latest model has been copied into self.robot.model, a L{QubModel}. 111 112 Add any additional widgets to self.top_panel and self.right_panel, or rearrange layout as needed. 113 114 model_prep has .source (the QubModel in GUI), .model (private copy for robot computation), 115 and .set_rates2(K0, K1, K2) [faster than making new model_prep] 116 117 @ivar right_panel: L{gtk.VBox} along right edge, with optimize controls already at the bottom 118 @ivar top_panel: L{gtk.HBox} along top edge, with score display already at far right 119 @ivar txtStatus: L{gtk.TextView} with text output; convenience methods self.append_output(s), self.clear_output() 120 @ivar model_prep: L{ModelPrep} with source model and private copy for robot use 121 122 """ 123 124 __explore_featured = ['right_panel', 'top_panel', 'txtStatus', 'model_prep', 'grads', 'hessian', 125 'OnLL', 'OnIteration', 'OnOptimized', 'model_tracker', 'data', 'data_prep', 126 'btnGetScore', 'txtScore', 'lblScore', 'btnAddToTrials', 'txtStatus', 127 'btnOpt', 'optimizing', 'score', 'model', 'get_model_prep', 'get_data_prep', 'request_data_prep', 128 'append_output', 'clear_output', 'add_to_trials', 'make_trial', 'get_ndata', 'add_trial_fields', 129 'calculate_score', 'optimize', 'robot_setup_data', 'update_score', 'robot_calculate_score', 130 'robot_get_pars', 'robot_set_pars', 'robot_optimize'] 131

132 - def __init__(self, label, global_name="", opt_name=None, opt_options={}, have_grads=False):

133 """ 134 @param label: tab label, becomes self.face_name 135 @param opt_name: optimizer preferences selector, defaults to label 136 @param opt_options: {option_label : default} to override specific defaults of the built-in optimizers (see qubx.optimize) 137 @param have_grads: True if your robot_calculate_score can also calculate dScore/dPar[i] 138 """ 139 super(OptToolPlugin, self).__init__(label, global_name=global_name) 140 self.__ref = Reffer() 141 self.__score = 0.0 142 self.__optimizing = False 143 self.have_grads = have_grads 144 self.opt_name = opt_name or label 145 qubx.optimize.setup_opt(self.opt_name, opt_options) 146 147 self.grads = [] 148 self.hessian = None 149 self.OnLL = WeakEvent() 150 self.OnIteration = WeakEvent() 151 self.OnOptimized = WeakEvent() 152 153 self.model_tracker = ModelTracker() 154 self.model_tracker.OnSwitch += self.__ref(self.__onSwitchModel) 155 self.model_tracker.OnChange += self.__ref(self.__onSetModel) 156 self.model_tracker.OnSetRate += self.__ref(self.__onSetRate) 157 self.model_prep = Product(self.get_model_prep) 158 self.model_prep.OnInvalidate += self.__ref(self.__onInvalidateModel) 159 self.rates = Product(self.__get_rates) 160 self.rates.OnInvalidate += self.__ref(self.__onInvalidateModel) 161 self.__model = None 162 self.data = OptToolData(self.get_data_prep) 163 self.data.OnInvalidate += self.__ref(self.__onInvalidateData) 164 self.data_prep = None 165 166 self.top_panel = pack_item(gtk.HBox(), self) 167 v = pack_item(gtk.VBox(), self.top_panel, at_end=True) 168 vh = pack_item(gtk.HBox(), v) 169 self.btnGetScore = pack_button('Get Score', vh, at_end=True, on_click=self.__onClickGetScore) 170 self.txtScore = pack_item(gtk.Entry(), vh, at_end=True) 171 self.txtScore.set_width_chars(10) 172 self.lblScore = pack_label('Score:', vh, at_end=True) 173 self.btnAddToTrials = pack_button('Add to Trials', v, expand=True, on_click=self.__onClickAddToTrials) 174 pack_space(self.top_panel, x=20, at_end=True) 175 176 box = pack_item(gtk.HBox(), self, expand=True) 177 self.txtStatus = gtk.TextView() 178 self.txtStatus.set_editable(False) 179 self.txtStatus.modify_font(pango.FontDescription('monospace 10')) 180 pack_scrolled(self.txtStatus, box, expand=True) 181 182 self.right_panel = pack_item(gtk.VBox(), box) 183 h = pack_item(gtk.HBox(True), self.right_panel, at_end=True) 184 self.btnOpt = pack_button('Optimize', h, on_click=self.__onClickOptimize) 185 186 self.optOptions = qubx.optimizeGTK.OptOptionsFace(label, self.opt_name, 'QubX.About.%s'%label) 187 # set self.optOptions to None if you're not using qubx.optimize 188 189 try: 190 qubx.global_namespace.QubX.Modeling.Utils.add_optimizer(label, self.get_model_prep, self.get_data_prep, 191 self.calculate_score, self.optimize, 192 self.OnLL, self.OnIteration, self.OnOptimized) 193 except: 194 traceback.print_exc() 195 196 self.__onSwitchModel(qubx.global_namespace.QubX.Models.file)

197 optimizing = property(lambda self: self.__optimizing)

198 - def set_score(self, x):

199 gobject.idle_add(self.txtScore.set_text, "%.5g" % x) 200 self.__score = x

201 score = property(lambda self: self.__score, lambda self, x: self.set_score(x), 202 "quantity to maximize, e.g. log-likelihoo; auto-displayed at top-right")

203 - def onShow(self, showing):

204 # if using built-in optimizers, show their options in the About panel 205 if not self.optOptions: 206 return 207 QubX = qubx.pyenv.env.globals['QubX'] 208 if showing: 209 QubX.About.append_face(self.optOptions) 210 QubX.About.show_face(self.optOptions.face_name) 211 else: 212 QubX.About.remove_face(QubX.About.index(self.optOptions.face_name))

213 - def __onSwitchModel(self, model):

214 if self.__optimizing: # put it off until optimization finishes 215 self.robot.do(gobject.idle_add, self.__onSwitchModel, model) 216 else: 217 self.__model = model 218 self.__onSetModel()

219 model = property(lambda self: self.__model)

220 - def __onInvalidateModel(self, product):

221 pass

222 #if not self.__optimizing: 223 # self.robot.do(gobject.idle_add, self.invalidate_controls)

224 - def __onSetModel(self, *args):

225 self.model_prep.invalidate() 226 self.rates.invalidate()

227 - def __onSetRate(self, r, field, val, prev, undoing):

228 self.rates.invalidate()

229 - def get_model_prep(self, model=None):

230 # override to customize the ModelPrep 231 m = model or self.__model 232 if not m: 233 return None 234 return ModelPrep(m)

235 - def __onInvalidateData(self, product):

236 self.model_prep.invalidate() 237 self.rates.invalidate()

238 - def get_data_prep(self, segs, model_prep, wait=True, receiver=lambda dp: None):

239 # leave this one alone; override the next one to customize DataPrep 240 if wait: 241 return qubx.pyenv.call_async_wait(self.request_data_prep, segs, model_prep, last_data_prep=self.data_prep) 242 else: 243 self.request_data_prep(receiver, segs, model_prep, last_data_prep=self.data_prep)

244 - def request_data_prep(self, receiver, segs=None, model_prep=None, last_data_prep=None):

245 # override to customize the DataPrep 246 data_prep = Anon(segs=segs, model_prep=model_prep) 247 receiver(data_prep)

248 - def __get_rates(self):

249 if self.__model == None: 250 return None, None 251 signal_names = qubx.global_namespace.QubX.Modeling.Stimulus.active_names 252 K0, K1, K2, L, V, P = ModelToRateMatricesP(self.__model, signal_names) 253 return K0, K1, K2

254 - def on_robot_begin(self):

255 def do_on_idle(): 256 self.btnGetScore.set_sensitive(False) 257 self.btnOpt.set_sensitive(False)

258 gobject.idle_add(do_on_idle)

259 - def on_robot_end(self):

260 def do_on_idle(): 261 self.btnGetScore.set_sensitive(True) 262 self.btnOpt.set_sensitive(True)

263 gobject.idle_add(do_on_idle)

264 - def __onClickGetScore(self, btn):

265 qubx.pyenv.env.scriptable_if_matching('%s.calculate_score()' % self.global_name, 266 [(self.global_name, self)]) 267 self.calculate_score(wait=False, receiver=self.set_score)

268 - def __onClickAddToTrials(self, btn):

269 self.add_to_trials(*self.make_trial())

270 - def __onClickOptimize(self, btn):

271 qubx.pyenv.env.scriptable_if_matching('%s.optimize(on_trial=%s.add_to_trials)' % (self.global_name, self.global_name), 272 [(self.global_name, self)]) 273 self.__mtree_pre_opt = self.__model.as_tree() 274 self.optimize(wait=False, on_trial=self.add_to_trials)

275 - def append_output(self, s):

276 """Appends a string to the output panel.""" 277 buf = self.txtStatus.get_buffer() 278 buf.insert(buf.get_end_iter(), s) 279 self.txtStatus.scroll_to_mark(buf.get_insert(), 0)

280 - def clear_output(self):

281 """Clears the text output panel.""" 282 self.txtStatus.set_text("")

283 - def on_exception(self, task, typ, val, tb):

284 self.append_output('\n'.join(traceback.format_exception(typ, val, tb)))

285 - def add_to_trials(self, mtree, row):

286 qubx.pyenv.env.OnScriptable('%s.add_to_trials(*%s.make_trial())' % (self.global_name, self.global_name)) 287 Trials = qubx.global_namespace.QubX.Trials 288 ts = Trials.get_trialset(self.face_name, in_qsf=True) 289 ts.append_trial(mtree, row) 290 Trials.trialset = ts

291 - def make_trial(self, model_prep=None, data_prep=None, fields={}, initial_tree=None):

292 try: 293 QubX = qubx.global_namespace.QubX 294 if model_prep: 295 model = model_prep.source or model_prep.model 296 else: 297 model = self.QubX.Models.file 298 if data_prep and data_prep.segs: 299 data = data_prep.segs[0].file 300 else: 301 data = self.QubX.Data.file 302 mtree = model.as_tree() 303 nParam = fields['nParam'] if ('nParam' in fields) else 0 304 if 'nData' in fields: 305 nData = fields['nData'] 306 else: 307 nData = self.get_ndata(data_prep) 308 row = qubx.trialset.make_trial_row(model.path, data.path, self.score, nParam, nData=nData) 309 row['nData'] = nData 310 if initial_tree: 311 self.add_trial_fields(row, initial_tree, 'Initial ', model_prep, data_prep) 312 self.add_trial_fields(row, mtree, "", model_prep, data_prep) 313 for k,v in fields.iteritems(): 314 row[k] = v 315 return mtree, row 316 except: 317 traceback.print_exc()

318 - def get_ndata(self, data_prep=None):

319 """Override to tell how many datapoints/events used, for Trials.""" 320 return 0

321 - def add_trial_fields(self, row, mtree, prefix='', model_prep=None, data_prep=None):

322 """Override to pick info from the model.""" 323 pass # e.g. qubx.trialset.add_trial_rates(row, mtree, prefix)

324 - def calculate_score(self, model_prep=None, data_prep=None, wait=True, receiver=lambda score: None, on_trial=None, **kw):

325 """Starts the robot calculating the score. 326 @param receiver: func(score) called when score is finished calculating 327 """ 328 if wait: 329 qubx.pyenv.call_async_wait(self.__calculate_score, model_prep, data_prep, on_trial=on_trial, **kw) 330 receiver(self.__score) 331 return self.__score 332 else: 333 self.__calculate_score(receiver, model_prep, data_prep, on_trial=on_trial, **kw)

334 - def __calculate_score(self, receiver=lambda score: None, model_prep=None, data_prep=None, on_trial=None, **kw):

335 QubX = qubx.pyenv.env.globals['QubX'] 336 self.__set_keys(**kw) 337 mp = model_prep 338 if mp is None: 339 mp = self.model_prep.val 340 mp.set_rates2(*self.rates.val) 341 self.robot.do(self.__robot_calculate_score, receiver, mp, data_prep, on_trial)

342 - def optimize(self, model_prep=None, data_prep=None, wait=True, receiver=None, ext_return=False, on_trial=None, **kw):

343 """Starts the robot optimizing. 344 @param receiver: func() called when optimization is finished 345 """ 346 if wait: 347 qubx.pyenv.call_async_wait(self.__optimize, model_prep, data_prep, on_trial=on_trial, **kw) 348 if ext_return: 349 if receiver: 350 receiver(self.score, self.iterations, self.grads, self.hessian) 351 return self.score, self.iterations, self.grads, self.hessian 352 else: 353 if receiver: 354 receiver(self.score) 355 return self.score 356 else: 357 self.__optimize(receiver, model_prep, data_prep, ext_return, on_trial=on_trial, **kw)

358 - def __optimize(self, receiver=None, model_prep=None, data_prep=None, ext_return=False, on_trial=None, **kw):

359 QubX = qubx.pyenv.env.globals['QubX'] 360 self.__set_keys(**kw) 361 mp = model_prep 362 if mp is None: 363 mp = self.model_prep.val 364 if mp.source: 365 mp.source.undoStack.enable(False) 366 mp.set_rates2(*self.rates.val) 367 self.__mtree_pre_opt = mp.model.as_tree() 368 QubX.Trials.get_trialset("%s Iterations"%self.face_name, in_qsf=True).clear() 369 self.__optimizing = True 370 self.robot.do(self.__robot_optimize, receiver, mp, data_prep, ext_return, on_trial)

371 - def __set_keys(self, **kw):

372 for k, v in kw.iteritems(): 373 try: 374 setattr(self, k, v) 375 except: 376 print 'MacRates optimizer has no option named %s' % k

377 - def robot_setup_data(self, model_prep, data_prep):

378 self.robot.model_prep = model_prep 379 if data_prep: 380 self.data_prep = data_prep 381 else: 382 self.data.model_prep = model_prep 383 self.data_prep = self.robot.gui_call_recv(self.data.request, self.robot)[0]

384 - def __robot_calculate_score(self, on_finish=lambda: None, model_prep=None, data_prep=None, on_trial=None):

385 """Calculates score on worker thread, then passes it back to GUI thread.""" 386 self.robot_setup_data(model_prep, data_prep) 387 try: 388 score = self.robot_calculate_score() 389 except: 390 score = 0.0 391 self.robot.send_exception() 392 if on_trial: 393 mtree, row = self.make_trial(model_prep, data_prep) 394 row['LL'] = score 395 gobject.idle_add(on_trial, mtree, row) 396 gobject.idle_add(self.update_score, score, lambda: on_finish(score)) 397 gobject.idle_add(self.OnLL, score)

398 - def __robot_optimize(self, on_finish=None, model_prep=None, data_prep=None, ext_return=False, on_trial=None):

399 """Optimizes on worker thread, then notifies GUI thread.""" 400 self.robot_setup_data(model_prep, data_prep) 401 mtree_initial = model_prep.model.as_tree() 402 try: 403 score, pars, iterations, grads, Hessian = self.robot_optimize() 404 self.iterations = iterations 405 self.grads = grads 406 self.hessian = Hessian 407 except: 408 self.robot.send_exception() 409 score = 0.0 410 self.score = score 411 self.__optimizing = False 412 if on_trial: 413 row = {'nParam' : len(self.grads), 414 'Iterations' : self.iterations} 415 for i, g in enumerate(self.grads): 416 row['Gradient %i' % i] = g 417 gobject.idle_add(on_trial, *self.make_trial(model_prep, data_prep, row, mtree_initial)) 418 if ext_return: 419 gobject.idle_add(self.update_score, score, lambda: on_finish and on_finish((score, iterations, grads, Hessian))) 420 else: 421 gobject.idle_add(self.update_score, score, lambda: on_finish and on_finish(score)) 422 gobject.idle_add(self.__onOptimized, iterations, score, grads, Hessian)

423 - def __onOptimized(self, iterations, score, grads, Hessian):

424 model_prep = self.robot.model_prep 425 if model_prep.source: 426 model_prep.source.undoStack.enable() 427 model_prep.source.undoStack.push_undo(bind(model_prep.source.from_tree, self.__mtree_pre_opt), 428 bind(model_prep.source.from_tree, model_prep.source.as_tree())) 429 model_prep.source.undoStack.seal_undo('set rates') 430 ts = qubx.global_namespace.QubX.Trials.get_trialset("%s Iterations"%self.face_name, in_qsf=True) 431 for mtree, row in self.iter_trials: 432 ts.append_trial(mtree, row) 433 self.OnOptimized(iterations, score, grads, Hessian)

434 - def update_score(self, score, call_me=lambda: None):

435 """Updates score and txtScore then calls call_me(). Intended as GUI callback from robot thread; e.g. 436 >>> gobject.idle_add(self.update_score, score, on_finish) 437 """ 438 self.score = score 439 call_me()

440 - def robot_calculate_score(self, grads=None):

441 """Calculates and returns the score, on worker thread; override this method. 442 self.model_prep.model has just been updated. 443 @param grads: either None (no gradient requested) or an array to fill with dScore/dPar[i] 444 """ 445 if grads != None: 446 grads[0] = 0.0 447 grads[1] = 0.0 448 return 1.0

449 - def robot_get_pars(self):

450 """Returns the list of initial parameter values, on worker thread; override this method. 451 Derive parameter values from self.robot.model_prep.model and/or from your own widgets.""" 452 return [1.0, -1.0]

453 - def robot_set_pars(self, pars):

454 """Updates self.robot.model_prep.model and/or other state from new param values, on worker thread; override this method.""" 455 pass

456 - def robot_optimize(self):

457 """Finds improved param values, starting from robot_get_pars(); override if you don't want the built-in optimizers. 458 If you override, also do self.optOptions = None. 459 self.robot.model_prep.model has just been updated. 460 """ 461 pars = self.robot_get_pars() 462 self.robot.iter = 0 463 self.iter_trials = [] 464 self.robot.maxiter = qubx.optimize.get_options(self.opt_name)['max iterations'] 465 if self.optOptions: 466 optimize = qubx.optimize.get_opts(self.opt_name, 467 on_start=lambda rec, opts: gobject.idle_add(self.optOptions.optimizers.onStartOpt, rec, opts), 468 on_end=lambda rec, opts: gobject.idle_add(self.optOptions.optimizers.onEndOpt, rec, opts)) 469 else: 470 optimize = qubx.optimize.get_opts(self.opt_name) 471 score, pars, iterations, grads, Hessian = optimize(array(pars), self.__robot_func, self.__robot_iter, self.have_grads and self.__robot_grad_func or None) 472 # important: use idle_add to pass GUI-changing stuff to the main thread: 473 gobject.idle_add(self.append_output, "Optimization finished after %d iterations.\nScore = %.9g\n" % (iterations, -score)) 474 # (on linux, it can be any thread, as long as the main thread is paused [with robot.main_hold:], 475 # but windows will crash unless all GUI API calls come from one thread.) 476 477 return -score, pars, iterations, grads, Hessian

478 # callbacks for built-in optimizers:

479 - def __robot_func(self, pars):

480 self.robot_set_pars(pars) 481 return - self.robot_calculate_score()

482 - def __robot_grad_func(self, pars, grads):

483 self.robot_set_pars(pars) 484 return - self.robot_calculate_score(grads)

485 - def __robot_iter(self, pars, grads, iterations, score):

486 self.robot.iter = iterations 487 self.robot.progress = self.robot.iter * 100.0 / self.robot.maxiter 488 gobject.idle_add(self.update_score, -score) 489 gobject.idle_add(self.OnIteration, self.robot.iter, score, grads) 490 self.robot_on_iteration(pars, grads, iterations, score) 491 mtree, row = self.make_trial(model_prep=self.robot.model_prep, data_prep=self.data_prep, fields={"Iteration": iterations}) 492 self.iter_trials.append((mtree, row))

493 - def robot_on_iteration(self, pars, grads, iterations, score):

494 """Responds to a completed iteration; override this function to send model params etc back to gui.""" 495 pass

496 497

498 -class ModelPrep(Anon):

499 """Represents a model ready for computation with QubX.Modeling.Utils."""

500 - def __init__(self, model=None, modeltree=None): # one or the other, not both; modeltree is more efficient

501 Anon.__init__(self, source=model, modeltree=modeltree or model.as_tree(), model=QubModel()) 502 self.model.from_tree(self.modeltree) 503 self.K0, self.K1, self.K2, L, V, P = ModelToRateMatricesP(self.model)

504 - def set_rates2(self, K0, K1, K2):

505 self.K0, self.K1, self.K2 = K0, K1, K2 506 K012ToModel(K0, K1, K2, self.model)

507 508

509 -class OptToolData(Requestable):

510 """Caches last-used data until changed."""

511 - def __init__(self, get_data_prep):

512 Requestable.__init__(self) 513 self.get_data_prep = get_data_prep 514 self.__ref = Reffer() 515 self.__use_group = None 516 self.__data = None 517 QubX = qubx.global_namespace.QubX 518 QubX.DataSource.OnChangeIdealization += self.__ref(self.__onChangeIdl) 519 QubX.Data.table.OnInsert += self.__ref(self.__onChangeFileList) 520 QubX.Data.table.OnRemoved += self.__ref(self.__onChangeFileList) 521 QubX.Data.table.OnSet += self.__ref(self.__onSetFileList) 522 QubX.Data.OnSwitch += self.__ref(self.__onSwitchDataFile) 523 self.__datafile = None 524 self.__onSwitchDataFile(QubX.Data, QubX.Data.file) 525 QubX.Modeling.Stimulus.OnInsert += self.__ref(self.__onChangeStimulus) 526 QubX.Modeling.Stimulus.OnChange += self.__ref(self.__onChangeStimulus) 527 QubX.Modeling.Stimulus.OnRemoving += self.__ref(self.__onChangeStimulus) 528 self.invalidate()

529 - def set_use_group(self, x):

530 if self.__use_group != x: 531 self.__use_group = x 532 self.invalidate()

533 use_group = property(lambda self: self.__use_group, lambda self, x: self.set_use_group(x))

534 - def __onChangeIdl(self, *args):

535 self.invalidate()

536 - def __onChangeFileList(self, *args):

537 if self.__use_group != None: 538 self.invalidate()

539 - def __onSetFileList(self, i, field_name, val, prev, undoing):

540 if (self.__use_group != None) and (field_name == 'Group'): 541 self.invalidate()

542 - def __onChangeStimulus(self, *args):

543 self.invalidate()

544 - def __onSwitchDataFile(self, Data, file):

545 if file == self.__datafile: return 546 if self.__datafile: 547 self.__datafile.constants.OnSet -= self.__ref(self.__onSetConstant) 548 self.__datafile = file 549 if self.__datafile: 550 self.__datafile.constants.OnSet += self.__ref(self.__onSetConstant)

551 - def __onSetConstant(self, r, field, val, prev, undoing):

552 if self.__datafile.constants[r,'Name'] == 'Tdead (ms)': 553 self.invalidate()

554 - def get(self, receiver):

555 QubX = qubx.global_namespace.QubX 556 if self.__use_group != None: 557 view_segm = QubX.DataSource.get_segmentation_files(group=self.__use_group) 558 else: 559 view_segm = [(QubX.Data.view, QubX.DataSource.get_segmentation())] 560 segs = [] 561 for view, segm in view_segm: 562 segs.extend(segm) 563 self.get_data_prep(segs, self.model_prep, wait=False, receiver=receiver)

564 565

566 -class AboutPluginExample(AboutPlugin):

567 __explore_featured = ['lblModelStatus', 'modelTracker']

568 - def __init__(self, name='Example', global_name='QubX.About.Example'):

569 super(AboutPluginExample, self).__init__(name, global_name=global_name) 570 571 # qubx uses WeakEvent for synchronous notification of listeners. 572 # WeakEvent uses weak references, so it's not a big deal if a listener never unsubscribes, 573 # but this means the references must be held somewhere else; 574 # it requires special handling for lambda expressions, as well as bound methods, 575 # which are ordinarily created and destroyed as needed. Thus the following weak 576 # reference becomes invalid by the end of the line: 577 # QubX.OnSomeEvent += self.__onSomeEvent 578 # We keep a permanent reference to event handlers in a Reffer object. 579 # self.__ref(x) == x 580 # and the act of calling self.__ref(x) adds it to the Reffer's internal dict (as a strong reference). 581 # Thus the idiom: 582 # QubX.OnSomeEvent += self.__ref(self.__onSomeEvent) 583 self.__ref = Reffer() 584 585 # add components to self (it is a subclass of gtk.VBox): 586 self.lblModelStatus = pack_label('', self, expand=True) 587 self.lblDataStatus = pack_label('', self, expand=True) 588 # they will show updated info about data and model, as they change: 589 self.modelTracker = ModelTracker() # defined below 590 self.modelTracker.OnSwitch += self.__ref(self.__onSwitchModel) 591 self.modelTracker.OnChange += self.__ref(self.__onChangeModel) 592 self.modelTracker.OnSetRate += self.__ref(self.__onSetModelRate) 593 self.__onChangeModel(self.modelTracker.model) # initialize 594 # QubX is the main window. DataSource represents the active portion of data 595 # (frontmost file, screen or whole file, which signal) 596 QubX = qubx.pyenv.env.globals['QubX'] 597 self.DataSource = QubX.DataSource 598 QubX.DataSource.OnChangeSamples += self.__ref(self.__onChangeSamples) 599 QubX.DataSource.OnChangeIdealization += self.__ref(self.__onChangeIdealization) 600 self.__onChangeData() # initialize

601 - def __onSwitchModel(self, models, model):

602 self.__onChangeModel(model)

603 - def __onSetModelRate(self, r, field, val, prev, undoing):

604 self.__onChangeModel(self.__model)

605 - def __onChangeModel(self, model):

606 self.lblModelStatus.set_text("%d states" % model.states.size)

607 - def __onChangeSamples(self):

608 self.__onChangeData()

609 - def __onChangeIdealization(self):

610 self.__onChangeData()

611 - def __onChangeData(self):

612 segs = self.DataSource.get_segmentation() 613 nseg = len(segs) 614 npoint = 0 615 nideal = 0 616 for seg in segs: 617 for chunk in seg.chunks: 618 if chunk.included: 619 npoint += chunk.n 620 ff, ll, cc = chunk.get_idealization() 621 nideal += len(ff) 622 self.lblDataStatus.set_text("%d points in %d segments\n(%d events)" % (npoint, nseg, nideal))

623 624

625 -class ToolPluginExample(ToolPlugin):

626 - def __init__(self, name='Test', global_name='QubX.Modeling.Test'):

627 super(ToolPluginExample, self).__init__(name, global_name=global_name) 628 # unused: self.__ref = Reffer() 629 630 self.lblMsg = pack_label('This plugin finds the max and min value in the Data Source.', self, expand=True) 631 self.panButtons = pack_item(gtk.HBox(), self) 632 self.btnRun = pack_button('Run', self.panButtons, at_end=True, on_click=self.__onClickRun)

633 - def on_robot_begin(self):

634 gobject.idle_add(self.btnRun.set_sensitive, False)

635 - def on_robot_end(self):

636 gobject.idle_add(self.btnRun.set_sensitive, True)

637 - def __onClickRun(self, btn):

638 # calls self.__robot_compute() on the robot (worker) thread: 639 self.robot.do(self.__robot_compute)

640 - def __robot_compute(self):

641 QubX = qubx.pyenv.env.globals['QubX'] 642 with self.robot.main_hold: # pause the main thread to read/write volatile info: 643 units = QubX.Data.file.signals.get(QubX.DataSource.signal, 'Units') 644 lo = 1e20 645 hi = -lo 646 segments = QubX.DataSource.get_segmentation() 647 for s, seg in enumerate(segments): 648 # A segment is divided in chunks, alternating included/excluded. 649 # The chunks don't have sampled data yet. 650 # gen_samples yields sampled chunks, and also subdivides any chunks bigger than maxlen. 651 # It also takes care of pausing the main thread to read data. 652 for chunk in QubX.DataSource.gen_samples(seg.chunks, self.robot.main_hold, maxlen=(1<<20)): 653 if chunk.included: 654 lo = min(lo, min(chunk.samples)) 655 hi = max(hi, max(chunk.samples)) 656 self.robot.progress = s * 100.0 / len(segments) 657 # for short chunks, gen_samples is equivalent to: 658 # for chunk in seg.chunks: 659 # if chunk.included: 660 # with self.robot.main_hold: 661 # datachunk = chunk.get_samples() 662 # lo = min(lo, min(datachunk.samples)) 663 # hi = max(hi, max(datachunk.samples)) 664 665 # gobject.idle_add(a, b, c, ...) calls a(b, c, ...) on the main thread, sometime later. 666 # It's good for displaying intermediate and final results, messages, etc. which must 667 # be done on the main thread for Win32 compatibility. 668 gobject.idle_add(self.__show_results, 'Lo: %.6g %s\nHi: %.6g %s' % (lo, units, hi, units))

669 - def __show_results(self, results):

670 self.lblMsg.set_text(results)

671 672

673 -class OptToolPluginExample(OptToolPlugin):

674 """Sample L{OptToolPlugin} which optimizes rate constants given idealized single-molecule data. 675 The algorithm is a simplification of MIL (Qin et al). 676 """ 677 __explore_featured = ['chkPeq']

678 - def __init__(self, label='Example', global_name='QubX.Modeling.Example'):

679 super(OptToolPluginExample, self).__init__(label, global_name=global_name, have_grads=False) # it will auto-numerical-differentiate when required 680 self.__ref = Reffer() 681 682 # add a checkbox: equilibrium or fixed entry probability 683 self.chkPeq = pack_check('start at equilibrium (else use States:P0)', self.top_panel, active=True) 684 685 self.append_output('This plugin demonstrates log-likelihood optimization of rate constants from idealized data, in just 85 lines.\nThe algorithm is a simplification of MIL (Qin et al, 1996).\n')

686

687 - def robot_get_pars(self):

688 # this plugin doesn't handle stimulus-dependent rates, just k0 689 model = self.robot.model_prep.model 690 for i in xrange(model.rates.size): 691 if (model.rates.get(i, 'k1') or model.rates.get(i, 'k2') 692 or model.rates.get(i, 'Ligand') 693 or model.rates.get(i, 'Voltage') or model.rates.get(i, 'Pressure')): 694 gobject.idle_add(self.append_output, "Warning: Ignoring Ligand- and/or Voltage-dependence.\n") 695 # the params are log(k0), to constrain k0 > 0 696 return [log(model.rates.get(i, 'k0')) for i in xrange(model.rates.size)]

697 - def robot_set_pars(self, pars):

698 model = self.robot.model_prep.model 699 k0 = [exp(pars[i]) for i in xrange(model.rates.size)] 700 # update rates in the worker-thread's model: 701 for i in xrange(model.rates.size): 702 model.rates.set(i, 'k0', k0[i]) 703 # and in the visible one, on the main thread: 704 gobject.idle_add(self.__update_gui_rates, k0)

705 - def __update_gui_rates(self, k0):

706 if self.robot.model_prep.source: 707 model = self.robot.model_prep.source 708 for i in xrange(model.rates.size): 709 model.rates.set(i, 'k0', k0[i])

710

711 - def request_data_prep(self, receiver, segs=None, model_prep=None, last_data_prep=None):

712 # this one runs in the gobject thread, as needed (cached by OptToolPlugin) 713 idl = [] 714 nevent = 0 715 dt = segs[0].sampling if segs else 1e-4 716 for seg in segs: 717 # first, last, class, duration -- we need only class and duration arrays 718 # i.e. a segment is a sequence of dwells (events): 719 # - in a certain class (color; collection of indistinguishable states) 720 # - for some number of seconds (by assumption, durations are even multiples of dt) 721 ff, ll, cc, dd = seg.get_idealization(mark_excluded=True, get_durations=True) 722 # "excluded" and un-idealized portions will have class < 0 723 dd *= 1e-3 # convert from ms to sec 724 idl.append( (cc, dd) ) 725 if any(cc >= 0): # don't count segments with un-idealized or all-excluded data 726 nevent += len(cc) 727 receiver(Anon(segs=segs, model_prep=model_prep, idl=idl, nevent=nevent, dt=dt))

728

729 - def robot_calculate_score(self, grads=None):

730 QubX = qubx.pyenv.env.globals['QubX'] 731 idl = self.data_prep.idl 732 nevent = self.data_prep.nevent 733 dt = self.data_prep.dt 734 if nevent == 0: 735 gobject.idle_add(self.append_output, "No idealized data available (check the data source, or idealize).\n") 736 return 0.0 737 # read widgets while the main thread is paused: 738 with self.robot.main_hold: 739 usePeq = self.chkPeq.get_active() 740 # get rate matrix (Q) and state colors (class) 741 model = self.robot.model_prep.model 742 Q = ModelToQ(model) 743 N = Q.shape[0] 744 class_of_state = ModelToClazz(model) 745 # get equilibrium or fixed entry prob 746 if usePeq: 747 P0 = QtoPe(Q) 748 else: 749 P0 = ModelToP0(model) 750 751 # math subroutines for transition probability and log likelihood accumulation: 752 753 # Prob(state b at time t+dt |given| state a at time t) = A[a,b] = e^(Q*dt)[a,b] 754 A_base = expm(Q*dt) 755 # except we zero out columns which are ruled out by data (wrong class [color]), in this work matrix: 756 A_work = matrix(zeros(shape=(N,N), dtype='float64')) 757 # cache exponentiated rate matrices A(dt*nsample), since there will be a limited number of multiples of dt 758 A_cache = {} 759 def get_A(dur, dest_class): 760 """Returns A'^(dur/dt), where A' = e^(Q*dt), except with A'[r,c]==0.0 if c!=dest_class; 761 caches all answers in A_cache and returns answers from cache when possible.""" 762 if not ((dur, dest_class) in A_cache): 763 # start with one-sample transition prob. matrix: 764 A_work[:] = A_base[:] 765 # if data specifies destination class, then zero out other columns 766 # (once per sample, it is observed in dest_class, so prob(going anywhere else in one sample) is zero) 767 if dest_class >= 0: 768 A_work[:,class_of_state != c] = 0.0 769 # exponentiate to repeat by number of samples 770 A_cache[(dur, dest_class)] = A_work ** int(round(dur/dt)) 771 return A_cache[(dur, dest_class)]

772 773 # accumulating likelihood quickly results in underflow; following Qin et al (1996), 774 # we accumulate -log(scale) at each dwell, and rescale Pt 775 def Scale(Pt): 776 """Returns scale=1.0/sum(Pt); multiplies each component of Pt by scale.""" 777 mag = sum(Pt) 778 if mag == 0.0: 779 raise Exception('Impossible!') 780 scale = 1.0 / mag 781 Pt *= scale 782 return scale

783 784 # main loop: 785 LL = 0.0 786 for classes, durations in idl: # each segment 787 for i, c, d in izip(count(), classes, durations): # each dwell: index, class, duration 788 # transition probability into class c: 789 if i == 0: # first dwell of segment: Pt[i] = P0[i], if state i is in class c 790 Pt = matrix(zeros(shape=(1,N), dtype='float64')) 791 for i in xrange(N): 792 Pt[0,i] = (class_of_state[i] == c) and P0[i] or 0.0 793 else: # apply transition matrix for one sample duration (rough approximation) 794 transition_dur = min(dt, durations[i-1]) 795 Pt *= get_A(transition_dur, c) 796 # transition probability for dwell duration (minus transition of dt): 797 dwell_dur = d - dt 798 if dwell_dur > 0.0: 799 Pt *= get_A(dwell_dur, c) 800 # rescale Pt toward 1.0, and accumulate log likelihood: 801 LL -= log(Scale(Pt)) 802 return LL 803 804 805

806 -class ModelTracker(object):

807 """Follows the frontmost L{qubx.model.QubModel} and L{qubx.modelGTK.ModelView}, 808 with events when anything changes. 809 810 @ivar OnSwitch: L{WeakEvent}(QubX.Models, model) when a new file is brought to front 811 @ivar OnChange: L{WeakEvent}(model) when something in the model changes (except a rate constant) 812 @ivar OnSetRate: L{WeakEvent}(rate_index, field_name, value, prev, is_undoing) 813 """

814 - def __init__(self):

815 self.__ref = Reffer() 816 QubX = qubx.pyenv.env.globals['QubX'] 817 self.Models = QubX.Models 818 self.Models.OnSwitch += self.__ref(self.__onSwitchModel) 819 self.__model = self.__view = None 820 self.OnSwitch = WeakEvent() 821 self.OnChange = WeakEvent() 822 self.OnSetRate = WeakEvent() 823 self.__onSwitchModel(self.Models, self.Models.file)

824 - def __onSwitchModel(self, models, model):

825 if model == self.__model: return 826 if self.__model: 827 for event in (self.__model.OnSetChannelCount, 828 self.__model.OnSetIeqFv, 829 self.__model.OnSetVRev, 830 self.__model.states.OnSet, 831 self.__model.states.OnInsert, 832 self.__model.states.OnRemoved, 833 self.__model.classes.OnSet, 834 self.__model.rates.OnInsert, 835 self.__model.rates.OnRemoved, 836 self.__model.OnChangeBalanceLoops, 837 self.__model.constraints_kin.OnEdit): 838 event -= self.__ref(self.__onChangeModel) 839 self.__model.rates.OnSet -= self.__ref(self.__onSetRate) 840 self.__model = model 841 self.__view = models.view 842 if self.__model: 843 for event in (self.__model.OnSetChannelCount, 844 self.__model.OnSetIeqFv, 845 self.__model.OnSetVRev, 846 self.__model.states.OnSet, 847 self.__model.states.OnInsert, 848 self.__model.states.OnRemoved, 849 self.__model.classes.OnSet, 850 self.__model.rates.OnInsert, 851 self.__model.rates.OnRemoved, 852 self.__model.OnChangeBalanceLoops, 853 self.__model.constraints_kin.OnEdit): 854 event += self.__ref(self.__onChangeModel) 855 self.__model.rates.OnSet += self.__ref(self.__onSetRate) 856 self.OnSwitch(model)

857 model = property(lambda self: self.__model) 858 view = property(lambda self: self.__view)

859 - def __onChangeModel(self, *args):

860 self.OnChange(self.__model)

861 - def __onSetRate(self, r, field, val, prev, undoing):

862 self.OnSetRate(r, field, val, prev, undoing)

863

Source Code for Module qubx_testing.qubx_plugin