qubx.modeling

1 """QubX.Modeling.Utils panel; hosts a registry of generic model/data-based optimizers, and sub-panels which can use them. 2 3 Optimizer methods should be called from the gobject (gui) thread. 4 To call from a qubx.task.Robot thread, use robot.gui_call_recv with wait=False. 5 6 Copyright 2003-2014 Research Foundation State University of New York 7 This file is part of QUB Express. 8 9 QUB Express is free software; you can redistribute it and/or modify 10 it under the terms of the GNU General Public License as published by 11 the Free Software Foundation, either version 3 of the License, or 12 (at your option) any later version. 13 14 QUB Express is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 GNU General Public License for more details. 18 19 You should have received a copy of the GNU General Public License, 20 named LICENSE.txt, in the QUB Express program directory. If not, see 21 <http://www.gnu.org/licenses/>. 22 23 """ 24 25 import collections 26 import gobject 27 import gtk 28 from gtk import gdk 29 import numpy 30 import random 31 import traceback 32 from itertools import izip 33 from math import * 34 35 import qubx.faces 36 import qubx.global_namespace 37 import qubx.GTK 38 import qubx.settings 39 from qubx.settings import Propertied, Property 40 from qubx.accept import * 41 from qubx.GTK import pack_item, pack_space, pack_hsep, pack_vsep, pack_label, pack_button, pack_check, pack_radio, pack_scrolled, build_menuitem 42 from qubx.trialset import make_trial_row, add_trial_rates 43 from qubx.util_types import *

44 45 -class Optimizer(Anon):

46 """Presents a consistent interface on one model/data-based optimization routine, e.g. MIL, MSL, MacRates. 47 Created by QubX.Modeling.Utils.add_optimizer, listed in QubX.Modeling.Utils.optimizers, assignable as QubX.Modeling.Utils.optimizer. 48 Used by sub-panels of QubX.Modeling.Utils."""

49 - def __init__(self, name, get_model_prep, get_data_prep, get_ll, optimize, OnLL, OnIteration, OnOptimized):

50 """ 51 52 @param name: string to identify the algorithm 53 @param get_model_prep: func(model=None, modeltree=None) -> model_prep 54 model: qubx.model.QubModel 55 modeltree: qubx.tree.Node containing QMF (only if no model provided) 56 model_prep: object derived from model, ready for computation, with fields: 57 .source: the input QubModel, if provided 58 .model: a private QubModel copy 59 .modeltree: QMF tree, even if not provided 60 since rates change so much more often, my algorithms re-use the rest of the prep when possible, by defining these fields: 61 .K0, .K1, .K2: numpy.matrix of rate constants 62 .set_rates(K0, K1): updates .K0, .K1, and .model 63 .set_rates2(K0, K1, K2): preferred newer form for pressure sensitivity 64 @param get_data_prep: func(segs, model_prep, wait=True, receiver=lambda dp: None) -> data_prep 65 segs: list of qubx.data_types.SourceSeg 66 model_prep: result of get_model_prep, since signal layout may depend on model stimuli 67 wait: False to return immediately 68 receiver: if wait, this function is called with data_prep when it's ready (on gobject thread) 69 data_prep: object with prepared data, ready for computation, with fields: 70 .segs: input the list of SourceSeg 71 .model_prep: input prepared model 72 and whatever else the algorithm needs, e.g. sampled or idealized data 73 @param get_ll: func(model_prep=None, data_prep=None, wait=True, receiver=lambda LL: None, on_trial, **kw) -> LL 74 evaluates the score for given model and data; fills them in from onscreen sources if not provided. 75 wait: False to return immediately 76 receiver: if wait == False: func(LL) called when computation is done 77 on_trial: func(model_tree, row_dict) with completed trial info 78 **kw: other keywords to apply to properties before running; e.g. use_Peq=True 79 @param optimize: func(model_prep=None, data_prep=None, wait=True, receiver=lambda LL: None, ext_return=False, on_trial, **kw) -> LL 80 iterates, adjusting parameters to improve LL for given model and data; fills them in from onscreen sources if not provided. 81 wait: False to return immediately 82 receiver: if wait == False: func(result--see ext_return) called when computation is done 83 ext_return: True to return the tuple (LL, iterations, grads, hessian); False to return LL 84 on_trial: func(model_tree, row_dict) with completed trial info 85 **kw: other keywords to apply to properties before running; e.g. use_Peq=True 86 @param OnLL: L{WeakEvent}(LL) when get_ll completes 87 @param OnIteration: L{WeakEvent}(iterations, LL, grads) each iteration during opt 88 @param OnOptimized: L{WeakEvent}(iterations, LL, grads, hessian) when optimize completes 89 """ 90 Anon.__init__(self, name=name, get_model_prep=get_model_prep, get_data_prep=get_data_prep, get_ll=get_ll, optimize=optimize, 91 OnLL=OnLL, OnIteration=OnIteration, OnOptimized=OnOptimized, LL=0.0, iterations=0, grads=[], hessian=None) 92 self.__ref = Reffer() 93 OnLL += self.__ref(self.__onLL) 94 OnIteration += self.__ref(self.__onIteration) 95 OnOptimized += self.__ref(self.__onOptimized)

96 - def __onLL(self, LL):

97 self.LL = LL

98 - def __onIteration(self, iterations, LL, grads):

99 self.iterations = iterations 100 self.LL = LL 101 self.grads = grads

102 - def __onOptimized(self, iterations, LL, grads, hessian):

103 self.__onIteration(iterations, LL, grads) 104 self.hessian = hessian

105

106 107 -class ModelingUtilsFace(qubx.faces.ToolsFace):

108 """Hosts a registry of generic model/data-based optimizers, and sub-panels which can use them. 109 110 @ivar optimizers: list of L{Optimizer} 111 """ 112 __explore_featured = ['optimizers', 'optimizer', 'optimizer_name', 'add_optimizer']

113 - def __init__(self, name='Utils', global_name="QubX.Modeling.Utils"):

114 qubx.faces.ToolsFace.__init__(self, name, pos=qubx.faces.POS_DROPDOWN, dropdown_label="Utility:", global_name=global_name) 115 self.__ref = Reffer() 116 self.__prefs = qubx.settings.SettingsMgr['qubx.modeling_utils'].active 117 pack_space(self.mnu_line, x=10, expand=True) 118 pack_label('Optimizer:', self.mnu_line) 119 self.mnuOptimizer = pack_item(qubx.GTK.StaticComboList([]), self.mnu_line, expand=True) 120 self.mnuOptimizer.OnChange += self.__ref(self.__onChangeOptimizer) 121 self.__optimizer = None 122 self.optimizers = [] 123 self.append_face(ChainFace()) 124 self.append_face(ClassifyFace()) 125 self.append_face(ModelSearchFace()) 126 self.append_face(PerturbFace()) 127 self.append_face(StarFace())

128 optimizer = property(lambda self: self.__optimizer, lambda self, x: self.set_optimizer(x), doc='the L{Optimizer} chosen in the gui') 129 optimizer_name = property(lambda self: self.__optimizer.name, lambda self, x: self.set_optimizer_name(x), doc='name of the chosen L{Optimizer} (set this to choose by name)')

130 - def add_optimizer(self, name, get_model_prep, get_data_prep, get_ll, optimize, OnLL, OnIteration, OnOptimized):

131 """Registers a model/data-based optimization method for sub-panel use. 132 133 @param name: string to identify the algorithm 134 @param get_model_prep: func(model=None, modeltree=None) -> model_prep 135 model: qubx.model.QubModel 136 modeltree: qubx.tree.Node containing QMF (only if no model provided) 137 model_prep: object derived from model, ready for computation, with fields: 138 .source: the input QubModel, if provided 139 .model: a private QubModel copy 140 .modeltree: QMF tree, even if not provided 141 since rates change so much more often, my algorithms re-use the rest of the prep when possible, by defining these fields: 142 .K0, .K1, .K2: numpy.matrix of rate constants 143 .set_rates(K0, K1): updates .K0, .K1, and .model 144 .set_rates2(K0, K1, K2): preferred newer form for pressure sensitivity 145 @param get_data_prep: func(segs, model_prep, wait=True, receiver=lambda dp: None) -> data_prep 146 segs: list of qubx.data_types.SourceSeg 147 model_prep: result of get_model_prep, since signal layout may depend on model stimuli 148 wait: False to return immediately 149 receiver: if wait, this function is called with data_prep when it's ready (on gobject thread) 150 data_prep: object with prepared data, ready for computation, with fields: 151 .segs: input the list of SourceSeg 152 .model_prep: input prepared model 153 and whatever else the algorithm needs, e.g. sampled or idealized data 154 @param get_ll: func(model_prep=None, data_prep=None, wait=True, receiver=lambda LL: None, **kw) -> LL 155 evaluates the score for given model and data; fills them in from onscreen sources if not provided. 156 wait: False to return immediately 157 receiver: if wait == False: func(LL) called when computation is done 158 **kw: other keywords to apply to properties before running; e.g. use_Peq=True 159 @param optimize: func(model_prep=None, data_prep=None, wait=True, receiver=lambda LL: None, ext_return=False, **kw) -> LL 160 iterates, adjusting parameters to improve LL for given model and data; fills them in from onscreen sources if not provided. 161 wait: False to return immediately 162 receiver: if wait == False: func(result--see ext_return) called when computation is done 163 ext_return: True to return the tuple (LL, iterations, grads, hessian); False to return LL 164 **kw: other keywords to apply to properties before running; e.g. use_Peq=True 165 @param OnLL: L{WeakEvent}(LL) when get_ll completes 166 @param OnIteration: L{WeakEvent}(iterations, LL, grads) each iteration during opt 167 @param OnOptimized: L{WeakEvent}(iterations, LL, grads, hessian) when optimize completes 168 """ 169 opt = Optimizer(name, get_model_prep, get_data_prep, get_ll, optimize, OnLL, OnIteration, OnOptimized) 170 for i, existing in enumerate(self.optimizers): 171 if existing.name == name: 172 self.optimizers[i] = opt 173 break 174 else: 175 self.optimizers.append(opt) 176 self.mnuOptimizer.add_choice(name) 177 if name == str(self.__prefs['optimizer_name'].data): 178 self.optimizer = opt

179 - def __onChangeOptimizer(self, mnu, name):

180 qubx.pyenv.env.OnScriptable('%s.optimizer_name = %s' % (self.global_name, repr(name))) 181 self.optimizer_name = name

182 - def set_optimizer(self, opt):

183 self.mnuOptimizer.active_i = self.optimizers.index(opt) 184 self.__optimizer = opt

185 - def set_optimizer_name(self, name):

186 for i, opt in enumerate(self.optimizers): 187 if opt.name == name: 188 self.__optimizer = opt 189 self.mnuOptimizer.active_i = i 190 self.__prefs['optimizer_name'].data = name 191 return 192 raise ValueError('No optimizer named %s' % repr(name))

193

194 195 -def task_exception_to_console(task, typ, val, tb):

196 traceback.print_exception(typ, val, tb)

197 198 199 200 @Propertied(Property('optimize', False, 'False to evaluate LL without optimizing'), 201 Property('model_group', 1, 'Group index of models to consider'))

202 -class ClassifyFace(qubx.faces.Face):

203 __explore_featured = ['run']

204 - def __init__(self):

205 qubx.faces.Face.__init__(self, 'Classify', 'QubX.Modeling.Utils.Classify') 206 self.propertied_connect_settings('Modeling.Utils.Classify') 207 line = pack_item(gtk.HBox(), self) 208 chkEval = pack_radio('evaluate or', line) 209 chkOpt = pack_radio('optimize', line, group=chkEval) 210 self.propertied_connect_radios('optimize', [(False, chkEval), (True, chkOpt)]) 211 pack_label('the LL of each selection in the active List', line) 212 line = pack_item(gtk.HBox(), self) 213 pack_label('for each model in group', line) 214 txtGroup = pack_item(qubx.GTK.NumEntry(self.model_group, acceptIntGreaterThanOrEqualTo(0), width_chars=4), line) 215 self.propertied_connect_NumEntry('model_group', txtGroup) 216 line = pack_item(gtk.HBox(), self) 217 pack_label("Each selection's Group will be set to the Index of the highest-scoring model.", line) 218 line = pack_item(gtk.HBox(), self) 219 pack_space(line, expand=True) 220 pack_space(line, expand=True) 221 self.btnRun = pack_button('Run', line, self.__onClickRun) 222 pack_space(line, expand=True)

223 - def __onClickRun(self, btn):

224 qubx.pyenv.env.OnScriptable('%s.run()' % self.global_name) 225 self.run(wait=False)

226 - def run(self, receiver=lambda: None, wait=True):

227 if wait: 228 return qubx.pyenv.call_async_wait(self.run, wait=False) 229 self.btnRun.set_sensitive(False) 230 task = ClassifyTask(qubx.global_namespace.QubX.Data.file.list, self.optimize, self.model_group, receiver) 231 task.OnException += task_exception_to_console 232 task.start()

233

234 -class ClassifyTask(qubx.task.Task):

235 - def __init__(self, lst, optimize=False, model_group=0, receiver=lambda: None):

236 qubx.task.Task.__init__(self, 'Classify') 237 self.__ref_onInterrupt = self.__onInterrupt 238 self.OnInterrupt += self.__ref_onInterrupt 239 self.__stop_flag = False 240 QubX = qubx.global_namespace.QubX 241 QubX.Trials.get_trialset('Classify').clear() 242 QubX.Trials.trialset_by_name = 'Classify' 243 self.lst = lst 244 self.receiver = receiver 245 self.optimizer = QubX.Modeling.Utils.optimizer 246 self.optimize = optimize 247 self.models_prep = [] 248 for i, view in enumerate(QubX.Models.views): 249 if QubX.Models.table[i, 'Group'] == model_group: 250 model_prep = self.optimizer.get_model_prep(view.file) 251 model_prep.Index = i 252 self.models_prep.append(model_prep) 253 self.segs = [] 254 for i in xrange(len(lst)): 255 self.segs.append(QubX.Data.file.get_segmentation_indexed(lst[i, 'From'], lst[i, 'To'], signal=QubX.DataSource.signal)[0])

256 - def __onInterrupt(self, task, cancel):

257 cancel() 258 self.__stop_flag = True

259 - def get_score(self, model_prep, data_prep, **kw):

260 def on_trial(model_tree, row): 261 row['Label'] = data_prep.segs[0].Label 262 qubx.global_namespace.QubX.Trials.get_trialset('Classify').append_trial(model_tree, row)

263 if self.optimize: 264 bk = model_prep.source.as_tree() 265 result = self.optimizer.optimize(model_prep, data_prep, on_trial=on_trial, **kw) 266 model_prep.source.from_tree(bk) 267 else: 268 result = self.optimizer.get_ll(model_prep, data_prep, on_trial=on_trial, **kw) 269 return result

270 - def run(self):

271 qubx.task.Tasks.add_task(self) 272 self.status = 'Running...' 273 try: 274 scores = [] 275 for s, seg in enumerate(self.segs): 276 seg.Label = 'List[%i]' % s 277 progress = self.progress = s * 100.0 / len(self.segs) 278 row = [] 279 scores.append(row) 280 for m, model_prep in enumerate(self.models_prep): 281 if self.__stop_flag: 282 raise KeyboardInterrupt 283 self.progress = progress + m * 100.0 / (len(self.segs) * len(self.models_prep)) 284 data_prep = self.gui_call_recv(self.optimizer.get_data_prep, [seg], model_prep, wait=False, receiver=self)[0] 285 ll = self.gui_call_recv(self.get_score, model_prep, data_prep, wait=False, receiver=self)[0] 286 row.append(ll) 287 groups = [numpy.argmax(row) for row in scores] if (len(scores) and len(scores[0])) else [] 288 gobject.idle_add(self.__set_groups, groups, scores) 289 finally: 290 qubx.global_namespace.QubX.Modeling.Utils.Classify.btnRun.set_sensitive(True) 291 gobject.idle_add(qubx.global_namespace.QubX.Modeling.Utils.request_show) 292 gobject.idle_add(self.receiver) 293 qubx.task.Tasks.remove_task(self)

294 - def __set_groups(self, groups, scores):

295 l = self.lst 296 mdl_lbls = ['%s %s' % (self.optimizer.name, model_prep.source and os.path.split(model_prep.source.path)[1] or str(model_prep.Index)) 297 for model_prep in self.models_prep] 298 for i, g in enumerate(groups): 299 l[i, 'Group'] = g 300 for j, score in enumerate(scores[i]): 301 l[i, mdl_lbls[j]] = score

302 303 304 305 @Propertied(Property('delta_ll', 1.0, 'minimum improvement in LL, or else it stops adding states'), 306 Property('hub_state', 0, 'Index of hub state (connected to added states)'))

307 -class StarFace(qubx.faces.Face):

308 __explore_featured = ['run']

309 - def __init__(self):

310 qubx.faces.Face.__init__(self, 'Star', 'QubX.Modeling.Utils.Star') 311 self.propertied_connect_settings('Modeling.Utils.Star') 312 line = pack_item(gtk.HBox(), self) 313 pack_label("Repeatedly adds states, connected from hub state, until log-likelihood improves by less than delta LL", line) 314 line = pack_item(gtk.HBox(), self) 315 pack_label('Hub state:', line) 316 txtHub = pack_item(qubx.GTK.NumEntry(self.hub_state, acceptIntGreaterThanOrEqualTo(0), width_chars=3), line) 317 self.propertied_connect_NumEntry('hub_state', txtHub) 318 pack_label('Delta LL:', line) 319 txtDeltaLL = pack_item(qubx.GTK.NumEntry(self.delta_ll, acceptFloatNonzero, '%.3g', width_chars=5), line) 320 self.propertied_connect_NumEntry('delta_ll', txtDeltaLL) 321 pack_space(line, x=20) 322 self.btnRun = pack_button('Run', line, self.__onClickRun)

323 - def __onClickRun(self, btn):

324 qubx.pyenv.env.OnScriptable('%s.run()' % self.global_name) 325 self.run(wait=False)

326 - def run(self, receiver=lambda: None, wait=True):

327 if wait: 328 return qubx.pyenv.call_async_wait(self.run, wait=False) 329 self.btnRun.set_sensitive(False) 330 task = StyledBuildTask('Star', AddStarState(self.hub_state), self.delta_ll, self, receiver) 331 task.OnException += task_exception_to_console 332 task.start()

333

334 @Propertied(Property('delta_ll', 1.0, 'minimum improvement in LL, or else it stops adding states')) 335 -class ChainFace(qubx.faces.Face):

336 __explore_featured = ['run']

337 - def __init__(self):

338 qubx.faces.Face.__init__(self, 'Chain', 'QubX.Modeling.Utils.Chain') 339 self.propertied_connect_settings('Modeling.Utils.Chain') 340 line = pack_item(gtk.HBox(), self) 341 pack_label("Repeatedly adds linear states until log-likelihood improves by less than delta LL", line) 342 line = pack_item(gtk.HBox(), self) 343 pack_label('Delta LL:', line) 344 txtDeltaLL = pack_item(qubx.GTK.NumEntry(self.delta_ll, acceptFloatNonzero, '%.3g', width_chars=5), line) 345 self.propertied_connect_NumEntry('delta_ll', txtDeltaLL) 346 pack_space(line, x=20) 347 self.btnRun = pack_button('Run', line, self.__onClickRun)

348 - def __onClickRun(self, btn):

349 qubx.pyenv.env.OnScriptable('%s.run()' % self.global_name) 350 self.run(wait=False)

351 - def run(self, receiver=lambda: None, wait=True):

352 if wait: 353 return qubx.pyenv.call_async_wait(self.run, wait=False) 354 self.btnRun.set_sensitive(False) 355 task = StyledBuildTask('Chain', AddChainState, self.delta_ll, self, receiver) 356 task.OnException += task_exception_to_console 357 task.start()

358

359 -class StyledBuildTask(qubx.task.Task):

360 - def __init__(self, name, add_state, delta_ll, face, receiver=lambda: None):

361 qubx.task.Task.__init__(self, name) 362 self.add_state = add_state 363 self.delta_ll = delta_ll 364 self.face = face 365 self.receiver = receiver 366 self.__ref_onInterrupt = self.__onInterrupt 367 self.OnInterrupt += self.__ref_onInterrupt 368 self.__stop_flag = False

369 - def __onInterrupt(self, task, cancel):

370 cancel() 371 self.__stop_flag = True

372 - def run(self):

373 QubX = qubx.global_namespace.QubX 374 qubx.task.Tasks.add_task(self) 375 self.status = 'Running...' 376 try: 377 build_styled_model(QubX.Modeling.Utils.optimizer, self.add_state, delta_ll=self.delta_ll, task=self, should_quit=lambda: self.__stop_flag) 378 finally: 379 try: 380 self.face.btnRun.set_sensitive(True) 381 except: 382 pass 383 gobject.idle_add(self.face.request_show) 384 gobject.idle_add(self.receiver) 385 qubx.task.Tasks.remove_task(self)

386

387 388 389 -def build_styled_model(optimizer, add_state, optimize=True, delta_ll=1.0, data_prep=None, task=None, should_quit=lambda: False):

390 """Builds up a model by, for all extant classes, repeatedly adding a state until LL improves by less than delta_ll. 391 Can be called from non-gobject Task threads if you provide the L{qubx.task.Task}. 392 393 @param optimizer: L{Optimizer} 394 @param add_state: func(class_index) called on gobject thread, operates on QubX.Models.file 395 @param optimize: False to evaluate model LL as-is 396 @param data_prep: optional prepared data to override onscreen 397 @param task: L{qubx.task.Task} if called from Task/robot thread 398 """ 399 QubX = qubx.global_namespace.QubX 400 model = QubX.Models.file 401 kept_tree = model.as_tree() 402 classes_present = sorted(list(set(state.Class for state in model.states))) 403 if task: 404 if optimize: 405 def get_score(): 406 return task.gui_call_recv(optimizer.optimize, data_prep=data_prep, wait=False, receiver=task)[0]

407 else: 408 def get_score(): 409 return task.gui_call_recv(optimizer.get_ll, data_prep=data_prep, wait=False, receiver=task)[0] 410 def gui_add_state(class_index): 411 task(add_state(class_index)) 412 def call_add_state(class_index): 413 return task.gui_call_recv(gui_add_state, class_index) 414 def gui_rewind(tree): 415 task(model.from_tree(tree)) 416 def rewind(tree): 417 task.gui_call_recv(gui_rewind, tree) 418 else: 419 if optimize: 420 def get_score(): 421 return optimizer.optimize(data_prep=data_prep) 422 else: 423 def get_score(): 424 return optimizer.get_ll(data_prep=data_prep) 425 def call_add_state(class_index): 426 return add_state(class_index) 427 rewind = model.from_tree 428 429 postLL = get_score() 430 preLL = postLL - 2 * delta_ll 431 for cls in classes_present: 432 while (postLL - preLL) > delta_ll: 433 if should_quit(): 434 return 435 kept_tree = model.as_tree() 436 call_add_state(cls) 437 preLL, postLL = postLL, get_score() 438 rewind(kept_tree) 439 preLL = postLL - 2 * delta_ll 440 get_score() 441

442 # star topology 443 444 -def AngleOfStates(m, s1, s2):

445 ss = m.states 446 dx = ss[s2, 'x'] - ss[s1, 'x'] 447 dy = ss[s2, 'y'] - ss[s1, 'y'] 448 mag = sqrt(dx*dx + dy*dy) 449 angle = asin(dy/mag) 450 if dx < 0: 451 if dy >= 0: 452 angle = pi - angle 453 else: 454 angle = -pi - angle 455 return angle

456

457 -def AddStarState(hub_index):

458 def add_state(cls): 459 QubX = qubx.global_namespace.QubX 460 model = QubX.Models.file 461 states = model.states 462 hub_ix = min(hub_index, len(states)-1) 463 464 angles = [ (AngleOfStates(model, hub_ix, i), i) for i in xrange(states.size) if i != hub_ix ] 465 angles.sort() 466 angles.append( (angles[0][0] + 2*pi, angles[0][1]) ) 467 # (ang. dist. from sorted_i to +1, state_ix, sorted state_ix) 468 anglediffs = [ (angles[i+1][0] - angles[i][0], angles[i][1], i) for i in range(len(angles) - 1) ] 469 anglediffs.sort() 470 insertafter = max(anglediffs) 471 insertaftersorted = insertafter[2] 472 abef, aaft = angles[ insertaftersorted ][0], angles[ insertaftersorted + 1 ][0] 473 insertangle = abef + (aaft - abef) / 2 474 if insertangle >= 2*pi: 475 insertangle = insertangle - 2*pi 476 x0, y0 = states[hub_ix, 'x'], states[hub_ix, 'y'] 477 dx, dy = cos(insertangle), sin(insertangle) 478 if dx != 0: 479 xmag = max( (100 - x0) / dx, -x0 / dx ) 480 else: 481 xmag = 1.0e10 482 if dy != 0: 483 ymag = max( (100 - y0) / dy, -y0 / dy ) 484 else: 485 ymag = 1.0e10 486 mag = 0.8 * min( xmag, ymag ) 487 488 states.append({'Class' : cls, 'x' : x0 + mag * dx, 'y' : y0 + mag * dy}) 489 model.rates.append({'From' : hub_ix, 'To' : len(states)-1}) 490 model.rates.append({'To' : hub_ix, 'From' : len(states)-1})

491 return add_state 492

493 # chain topology 494 495 -def AddChainState(cls):

496 QubX = qubx.global_namespace.QubX 497 model = QubX.Models.file 498 rates = model.rates 499 states = model.states 500 501 prev, statesOfCls, post = SubChain(model, cls) 502 503 s = len(states) 504 if len(statesOfCls) == 1 and prev is None: 505 states.append({'Class' : cls, 'x' : 5, 'y' : states[0, 'y']}) 506 rates.append({'From' : statesOfCls[0], 'To' : s}) 507 rates.append({'To' : statesOfCls[0], 'From' : s}) 508 elif len(statesOfCls) == 1 and post is None: 509 states.append({'Class' : cls, 'x' : 95, 'y' : states[0, 'y']}) 510 rates.append({'From' : statesOfCls[0], 'To' : s}) 511 rates.append({'To' : statesOfCls[0], 'From' : s}) 512 else: 513 if post: statesOfCls.append(post) 514 gaps = [ (states[statesOfCls[i+1], 'x'] - states[statesOfCls[i], 'x'], i) for i in range(len(statesOfCls)-1) ] 515 gaps.sort() 516 after = gaps[ len(gaps) - 1 ][1] 517 lh, rh = statesOfCls[after], statesOfCls[after+1] 518 states.append({'Class' : cls, 'x' : (states[lh,'x'] + states[rh,'x'])/2, 'y' : (states[lh,'y'] + states[rh,'y'])/2}) 519 520 rf = rb = None 521 for r in xrange(rates.size): 522 if (rates[r,'From'] == lh) and (rates[r,'To'] == rh): 523 rf = rates.get_row(r) 524 elif (rates[r, 'To'] == lh) and (rates[r,'From'] == rh): 525 rb = rates.get_row(r) 526 for r_del in reversed(sorted([x.Index for x in (rf, rb) if not (x is None)])): 527 rates.remove(r_del) 528 529 k0f = 10.0 if (rf is None) else (rf.k0 / 2) 530 k0b = 10.0 if (rb is None) else (rb.k0 / 2) 531 Lf = '' if (rf is None) else rf.Ligand 532 Lb = '' if (rb is None) else rb.Ligand 533 Vf = '' if (rf is None) else rf.Voltage 534 Vb = '' if (rb is None) else rb.Voltage 535 Pf = '' if (rf is None) else rf.Pressure 536 Pb = '' if (rb is None) else rb.Pressure 537 if rf is None: 538 k1f = .001 if Vf else 0.0 539 k2f = 1.0 if Pf else 0.0 540 else: 541 k1f = rf.k1 / 2 542 k2f = rf.k2 / 2 543 if rb is None: 544 k1b = .001 if Vb else 0.0 545 k2b = 1.0 if Pb else 0.0 546 else: 547 k1b = rb.k1 / 2 548 k2b = rb.k2 / 2 549 rates.append({'From' : lh, 'To' : s, 'k0' : k0f, 'k1' : k1f, 'k2' : k2f, 'Ligand' : Lf, 'Voltage' : Vf, 'Pressure' : Pf}) 550 rates.append({'From' : s, 'To' : lh, 'k0' : k0b, 'k1' : k1b, 'k2' : k2b, 'Ligand' : Lb, 'Voltage' : Vb, 'Pressure' : Pb}) 551 rates.append({'From' : rh, 'To' : s, 'k0' : k0b, 'k1' : k1b, 'k2' : k2b, 'Ligand' : Lb, 'Voltage' : Vb, 'Pressure' : Pb}) 552 rates.append({'From' : s, 'To' : rh, 'k0' : k0f, 'k1' : k1f, 'k2' : k2f, 'Ligand' : Lf, 'Voltage' : Vf, 'Pressure' : Pf})

553

554 -def Connections(model):

555 connections = collections.defaultdict(lambda: []) 556 for rate in model.rates: 557 connections[ rate.From ].append( rate.To ) 558 return connections

559

560 -def LeftmostState(model):

561 return min( [(model.states[i, 'x'], i) for i in range(model.states.size)] )[1]

562

563 -def ChainStateOrder(model):

564 conn = Connections(model) 565 at = LeftmostState(model) 566 order = [at] 567 while len(conn[at]): 568 next = conn[at][0] 569 del conn[at][0] 570 conn[next].remove(at) 571 at = next 572 order.append(at) 573 return order

574

575 -def SubChain(model, cls):

576 order = ChainStateOrder(model) 577 subchain = [(order[i], i) for i in range(len(order)) if model.states[order[i], 'Class'] == cls] # assumes all states of cls are contiguous in model 578 prev = None 579 post = None 580 581 ixInOrder = subchain[0][1] # index of first subchain item in order 582 if ixInOrder > 0: 583 prev = order[ ixInOrder - 1 ] 584 ixInOrder = subchain[ len(subchain) - 1 ][1] 585 if ixInOrder < (len(order) - 1): 586 post = order[ ixInOrder + 1 ] 587 588 subchain = [item[0] for item in subchain] 589 return prev, subchain, post

590 591 592 593 594 595 @Propertied(Property('trials_per_model', 3, 'Number of times to optimize each model'), 596 Property('no_loops', True, 'True to discard models with cycles'), 597 Property('balance_loops', True, 'True to add detailed balance constraints to all cycles'), 598 Property('max_connections', 10, 'Discards models with more connections between states'), 599 Property('max_exit_rates', 5, 'Discards models with excess connections to/from a single state'), 600 Property('initial_k0', 100.0, 'Starting pre-exponential rate constant'), 601 Property('random_rates', True, 'True to scale starting rates randomly from initial_k0'), 602 Property('random_from', 0.1, 'Minimum random scaling of starting rates'), 603 Property('random_to', 5.0, 'Maximum random scaling of starting rates'), 604 Property('output_name', 'Model Search', 'Name of output "Trials" table'), 605 Property('max_output', 0, 'If nonzero, sorts by AIC (or LL if AIC is 0), and outputs only the top entries.'), 606 Property('decoration_script', """# This script can adjust each model before optimizing. 607 # Work on search_model. search_number is the nauty class. For examples, turn on script recording (Admin:Scripts panel). 608 609 # e.g. adding pressure-dependence to opening rates: 610 #rates, states = search_model.rates, search_model.states 611 #for rate in rates: # rate is read-only; write to rates[index, field] 612 # if states[rate.From, 'Class'] < states[rate.To, 'Class']: 613 # rates[rate.Index, 'Pressure'] = "Pressure" 614 # rates[rate.Index, 'k2'] = 0.01 615 """, 'Script which is run for each trial, before optimizing'), 616 Property('post_script', """# This script runs after optimizing. Work on search_model. 617 """, 'Script which is run for each trial, after optimizing', in_tree=False))

618 -class ModelSearchFace(qubx.faces.Face):

619 __explore_featured = ['model_lists', 'txtScript', 'deferred_update', 'list_models', 'run']

620 - def __init__(self):

621 qubx.faces.Face.__init__(self, 'ModelSearch', 'QubX.Modeling.Utils.ModelSearch') 622 self.__ref = Reffer() 623 self.model_lists = qubx.tree.Open(os.path.join(qubx.global_namespace.app_path, 'msearch-models.qtr'), True) 624 self.model_lists.close() 625 self.propertied_connect_settings('Modeling.Utils.ModelSearch') 626 self.left_right = pack_item(gtk.HPaned(), self) 627 self.panOptions = gtk.HBox() 628 self.panOptions.show() 629 self.left_right.pack1(self.panOptions, False, True) 630 self.panOptions1 = pack_item(gtk.VBox(), self.panOptions, expand=True) 631 line = pack_item(gtk.HBox(True), self.panOptions1) 632 pack_label('Trials per model:', line) 633 self.txtTrialsPerModel = pack_item(qubx.GTK.NumEntry(self.trials_per_model, acceptIntGreaterThan(0), width_chars=6), line) 634 self.propertied_connect_NumEntry('trials_per_model', self.txtTrialsPerModel) 635 self.chkNoLoops = pack_check('No loops', self.panOptions1) 636 self.propertied_connect_check('no_loops', self.chkNoLoops) 637 self.chkBalanceLoops = pack_check('Balance loops', self.panOptions1) 638 self.propertied_connect_check('balance_loops', self.chkBalanceLoops) 639 line = pack_item(gtk.HBox(True), self.panOptions1) 640 pack_label('Max connections:', line) 641 self.txtMaxConnections = pack_item(qubx.GTK.NumEntry(self.max_connections, acceptIntGreaterThan(1), width_chars=6), line) 642 self.propertied_connect_NumEntry('max_connections', self.txtMaxConnections) 643 line = pack_item(gtk.HBox(), self.panOptions1) 644 pack_label('Max exit rates per state:', line) 645 line = pack_item(gtk.HBox(True), self.panOptions1) 646 pack_item(gtk.EventBox(), line) 647 self.txtMaxExitRates = pack_item(qubx.GTK.NumEntry(self.max_exit_rates, acceptIntGreaterThan(0), width_chars=6), line) 648 self.propertied_connect_NumEntry('max_exit_rates', self.txtMaxExitRates) 649 650 self.panOptions2 = pack_item(gtk.VBox(), self.panOptions, expand=True) 651 line = pack_item(gtk.HBox(True), self.panOptions2) 652 pack_label('Initial k0:', line) 653 self.txtInitialK0 = pack_item(qubx.GTK.NumEntry(self.initial_k0, acceptFloatGreaterThan(0.0), '%.5g', width_chars=6), line) 654 self.propertied_connect_NumEntry('initial_k0', self.txtInitialK0) 655 self.chkRandomRates = pack_check('Randomize starting rates', self.panOptions2) 656 self.propertied_connect_check('random_rates', self.chkRandomRates) 657 line = pack_item(gtk.HBox(True), self.panOptions2) 658 pack_label(' from:', line) 659 self.txtRandomFrom = pack_item(qubx.GTK.NumEntry(self.random_from, acceptFloatBetween(1e-5, 1.0), '%.5g', width_chars=6), line) 660 self.propertied_connect_NumEntry('random_from', self.txtRandomFrom) 661 line = pack_item(gtk.HBox(True), self.panOptions2) 662 pack_label(' to:', line) 663 self.txtRandomTo = pack_item(qubx.GTK.NumEntry(self.random_to, acceptFloatGreaterThan(1.0), '%.5g', width_chars=6), line) 664 self.propertied_connect_NumEntry('random_to', self.txtRandomTo) 665 pack_label('times initial k0', self.panOptions2) 666 line = pack_item(gtk.HBox(), self.panOptions2) 667 self.lblModelCount = pack_label('', line) 668 669 self.panScript = gtk.VBox() 670 self.panScript.show() 671 self.left_right.pack2(self.panScript, True, True) 672 lbl = pack_label('Per-trial model decoration script:', self.panScript) 673 ### scrolled script view; update self.decoration_script on exit textview 674 self.txtScript = gtk.TextView() 675 pack_scrolled(self.txtScript, self.panScript, size_request=(100, 40), expand=True) 676 self.txtScript.connect('focus_out_event', self.__onExitScript) 677 qubx.GTK.SetFixedWidth(self.txtScript) 678 self.txtScript.get_buffer().set_text(self.decoration_script) 679 680 line = pack_item(gtk.HBox(), self.panScript, at_end=True) 681 pack_label('Output Trials name:', line) 682 self.txtOutputName = pack_item(qubx.GTK.NumEntry(self.output_name), line, expand=True) 683 self.propertied_connect_NumEntry('output_name', self.txtOutputName) 684 self.mnuPresets = pack_item(qubx.settingsGTK.PresetsMenu('Modeling.Utils.ModelSearch', qubx.global_namespace.QubX.appname), line) 685 self.btnRun = pack_button('Run Search', line, self.__onClickRun) 686 687 QubX = qubx.global_namespace.QubX 688 self.deferred_update = qubx.pyenv.DeferredAction(delay_ms=5) 689 QubX.Models.OnSwitch += self.__ref(self.__onSwitchModel) 690 self.__model = None 691 self.__onSwitchModel(QubX.Models, QubX.Models.file)

692

693 - def onShow(self, showing):

694 qubx.faces.Face.onShow(self, showing) 695 if showing: 696 self.deferred_update(self, self.__update_model_count)

697 - def propertied_set(self, value, name):

698 super(ModelSearchFace, self).propertied_set(value, name) 699 if name in ('no_loops', 'max_connections', 'max_exit_rates'): 700 self.deferred_update(self, self.__update_model_count) 701 if name == 'decoration_script': 702 self.txtScript.get_buffer().set_text(value)

703 - def __onExitScript(self, txt, event):

704 buf = txt.get_buffer() 705 self.propertied_on_user_set('decoration_script', buf.get_text(buf.get_start_iter(), buf.get_end_iter()))

706 - def __update_model_count(self):

707 if self.__model and self.showing: 708 self.lblModelCount.set_text("(%i models)" % len(self.list_models(self.__model, self.no_loops, self.max_connections, self.max_exit_rates).nauty_codes))

709 - def __onSwitchModel(self, Models, file):

710 if self.__model == file: return 711 if self.__model: 712 self.__model.states.OnInsert -= self.__ref(self.__onChangeStates) 713 self.__model.states.OnRemoved -= self.__ref(self.__onChangeStates) 714 self.__model.states.OnSet -= self.__ref(self.__onSetState) 715 self.__model = file 716 if self.__model: 717 self.__model.states.OnInsert += self.__ref(self.__onChangeStates) 718 self.__model.states.OnRemoved += self.__ref(self.__onChangeStates) 719 self.__model.states.OnSet += self.__ref(self.__onSetState) 720 if self.showing: 721 self.deferred_update(self, self.__update_model_count)

722 - def __onChangeStates(self, *args):

723 if self.showing: 724 self.deferred_update(self, self.__update_model_count)

725 - def __onSetState(self, i, field, val, prev, undoing):

726 if self.showing and (field == 'Class'): 727 self.deferred_update(self, self.__update_model_count)

728 - def list_models(self, file=None, no_loops=None, max_connections=None, max_exit_rates=None):

729 mdl = self.__model if (file is None) else file 730 if not mdl: return Anon(nstate=0, nauty_codes=[]) 731 noloop = self.no_loops if (no_loops is None) else no_loops 732 maxconn = self.max_connections if (max_connections is None) else max_connections 733 maxexit = self.max_exit_rates if (max_exit_rates is None) else max_exit_rates 734 nauty_name, nauty_state_order = nauty_model_class(mdl) 735 list_node = self.model_lists[nauty_name] 736 if noloop: 737 list_node = list_node['NoLoops'] 738 nauty_codes = list_node.data[:] 739 nstate = len(mdl.states) 740 for i in reversed(xrange(len(nauty_codes))): 741 if not nauty_keep_model(nstate, nauty_codes[i], maxconn, maxexit): 742 del nauty_codes[i] 743 return Anon(nstate=nstate, nauty_codes=nauty_codes, state_order=nauty_state_order)

744 - def __onClickRun(self, btn):

745 qubx.pyenv.env.OnScriptable('%s.run()' % self.global_name) 746 self.run(wait=False)

747 - def run(self, receiver=lambda: None, wait=True):

748 models = self.list_models() 749 if not models.nauty_codes: 750 receiver() 751 return 752 decorator = qubx.pyenv.ScriptModule('<ModelSearch>', qubx.pyenv.env.globals) 753 postor = qubx.pyenv.ScriptModule('<ModelSearch>', qubx.pyenv.env.globals) 754 try: 755 decorator.preload(self.decoration_script) 756 except: 757 qubx.pyenvGTK.show_message(traceback.format_exc(), title='Model Search: Decoration script error', parent=self.parent_window) 758 receiver() 759 return 760 if decorator: 761 def decorate(model, number): 762 qubx.global_namespace.search_model = model 763 qubx.global_namespace.search_number = number 764 decorator(raise_exceptions=True)

765 else: 766 decorate = None 767 try: 768 postor.preload(self.post_script) 769 except: 770 qubx.pyenvGTK.show_message(traceback.format_exc(), title='Model Search: Post script error', parent=self.parent_window) 771 receiver() 772 return 773 if postor: 774 def post_script(model): 775 qubx.global_namespace.search_model = model 776 postor(raise_exceptions=True)

777 else: 778 post_script = None 779 #nauty_setup_model(self.__model, models.state_order, models.nauty_codes[int(round(random.uniform(-0.4, len(models.nauty_codes)-0.6)))], 780 # self.balance_loops, self.initial_k0, self.random_rates, self.random_from, self.random_to, decorate) ### 781 if wait: 782 return qubx.pyenv.call_async_wait(self.run, wait=False) 783 self.btnRun.set_sensitive(False) 784 task = ModelSearchTask(self.__model, models, self.trials_per_model, self.balance_loops, self.initial_k0, 785 self.random_rates, self.random_from, self.random_to, decorate, post_script, self.output_name, receiver) 786 task.OnException += task_exception_to_console 787 task.start() 788

789 -def nauty_model_class(model):

790 clazz = [model.states[i, 'Class'] for i in xrange(len(model.states))] 791 class_size = [0] * (max(clazz) + 1) 792 for c in clazz: 793 class_size[c] += 1 794 class_order = [] 795 for n in reversed(xrange(1, len(clazz))): 796 if n in class_size: 797 for i, nc in enumerate(class_size): 798 if nc == n: 799 class_order.append(i) 800 nauty_class_chars = ['N', str(len(clazz))] 801 state_order = [] 802 for i, sel_cls in enumerate(class_order): 803 for j, c in enumerate(clazz): 804 if c == sel_cls: 805 nauty_class_chars.append(str(i+1)) 806 state_order.append(j) 807 return ''.join(nauty_class_chars), state_order

808

809 -def nauty_keep_model(nstate, code, max_conn, max_exit):

810 counts = [0] * nstate 811 for i, conn in enumerate(nauty_connections(nstate, code)): 812 if i == max_conn: 813 return False 814 a, b = conn 815 counts[a] += 1 816 if counts[a] > max_exit: 817 return False 818 counts[b] += 1 819 if counts[b] > max_exit: 820 return False 821 return True

822

823 -def nauty_setup_model(model, state_order, code, balance_loops, initial_k0, random_rates, random_from, random_to, decorate):

824 model.rates.clear() 825 for a, b in nauty_connections(model.states.size, code): 826 model.rates.append({'From' : state_order[a], 'To' : state_order[b], 827 'k0' : random_rates and random_scaled(initial_k0, random_from, random_to) or initial_k0}) 828 model.rates.append({'To' : state_order[a], 'From' : state_order[b], 829 'k0' : random_rates and random_scaled(initial_k0, random_from, random_to) or initial_k0}) 830 if balance_loops: 831 model.balance_all_loops() 832 if decorate: 833 try: 834 decorate(model, code) 835 except: 836 qubx.pyenvGTK.show_message(traceback.format_exc(), title='Error in ModelSearch decoration script')

837 # reset channel count, amp/cond?

838 839 -def nauty_connections(nstate, code):

840 maxpath = nstate * (nstate - 1) / 2 841 pathsel = 1 << (maxpath - 1) 842 for a in xrange(nstate): 843 for b in xrange(a+1, nstate): 844 if code & pathsel: 845 yield(a, b) 846 pathsel >>= 1

847

848 -def random_scaled(center, lo_mult, hi_mult):

849 return center * exp(random.uniform(log(lo_mult), log(hi_mult)))

850

851 852 853 -class ModelSearchTask(qubx.task.Task):

854 - def __init__(self, model, models, trials_per_model, balance_loops, initial_k0, random_rates, random_from, random_to, decorate, post_script, output_name, receiver=lambda: None):

855 qubx.task.Task.__init__(self, 'ModelSearch') 856 self.__ref_onInterrupt = self.__onInterrupt 857 self.OnInterrupt += self.__ref_onInterrupt 858 self.__stop_flag = False 859 self.model = model 860 self.models = models 861 self.trials_per_model = trials_per_model 862 self.balance_loops = balance_loops 863 self.initial_k0 = initial_k0 864 self.random_rates = random_rates 865 self.random_from = random_from 866 self.random_to = random_to 867 self.decorate = decorate 868 self.post_script = post_script 869 self.output_name = output_name 870 self.receiver = receiver 871 QubX = qubx.global_namespace.QubX 872 QubX.Trials.get_trialset(output_name).clear() 873 QubX.Trials.trialset_by_name = output_name 874 self.optimizer = QubX.Modeling.Utils.optimizer

875 - def start_trial(self, code, label, receiver):

876 nauty_setup_model(self.model, self.models.state_order, code, self.balance_loops, self.initial_k0, 877 self.random_rates, self.random_from, self.random_to, self.decorate) 878 self.trial_label = label 879 self.optimizer.optimize(self.optimizer.get_model_prep(self.model), wait=False, receiver=receiver, on_trial=self.__onTrial)

880 - def end_trial(self, code, label, receiver):

881 if self.post_script: 882 try: 883 self.post_script(self.model) 884 except: 885 qubx.pyenvGTK.show_message(traceback.format_exc(), title='Error in ModelSearch post script') 886 receiver()

887 - def __onTrial(self, model_tree, row_dict):

888 row_dict['Label'] = self.trial_label 889 qubx.global_namespace.QubX.Trials.get_trialset(self.output_name).append_trial(model_tree, row_dict)

890 - def __onInterrupt(self, task, cancel):

891 cancel() 892 self.__stop_flag = True

893 - def run(self):

894 if self.trials_per_model > 1: 895 trial_label = lambda code, t: "%i-%i" % (code, t+1) 896 else: 897 trial_label = lambda code, t: code 898 qubx.task.Tasks.add_task(self) 899 self.status = 'Running...' 900 try: 901 for i, code in enumerate(self.models.nauty_codes): 902 for t in xrange(self.trials_per_model): 903 if self.__stop_flag: 904 raise KeyboardInterrupt 905 self.gui_call_recv(self.start_trial, code, trial_label(code, t), receiver=self) 906 self.gui_call_recv(self.end_trial, code, trial_label(code, t), receiver=self) 907 self.progress = (i*self.trials_per_model + t) * 100.0 / (len(self.models.nauty_codes * self.trials_per_model)) 908 finally: 909 QubX = qubx.global_namespace.QubX 910 QubX.Modeling.Utils.ModelSearch.btnRun.set_sensitive(True) 911 gobject.idle_add(QubX.Modeling.Utils.request_show) 912 gobject.idle_add(self.receiver) 913 gobject.idle_add(lambda: QubX.Tables.show_table(QubX.Tables.find_table('Trials'))) 914 qubx.task.Tasks.remove_task(self)

915 916 917 918 @Propertied(Property('optimize', False, 'False to evaluate LL without optimizing'), 919 Property('all_at_once', False, 'False to perturb and evaluate parameters one at a time, into the Trials table'), 920 Property('perturb_rates', True, ''), 921 Property('perturb_channelCount', False, ''), 922 Property('perturb_amp', False, ''), 923 Property('perturb_std', False, ''), 924 Property('perturb_cond', False, ''), 925 Property('perturb_condstd', False, ''), 926 Property('perturb_pct', 10.0, 'Perturbed parameter in gaussian(mean=original, std=original * (pct/100))'), 927 Property('trials', 1, 'Number of repeats'))

928 -class PerturbFace(qubx.faces.Face):

929 __explore_featured = ['run']

930 - def __init__(self):

931 qubx.faces.Face.__init__(self, 'Perturb', 'QubX.Modeling.Utils.Perturb') 932 self.propertied_connect_settings('Modeling.Utils.Perturb') 933 columns = pack_item(gtk.HBox(), self, expand=True) 934 col = pack_item(gtk.VBox(), columns) 935 pack_label('Adjust', col) 936 col = pack_item(gtk.VBox(), columns) 937 self.chkEach = pack_radio('each', col) 938 self.chkAll = pack_radio('all at once', col, group=self.chkEach) 939 self.propertied_connect_radios('all_at_once', [(False, self.chkEach), (True, self.chkAll)]) 940 col = pack_item(gtk.VBox(), columns) 941 self.chkPerturbRates = pack_check('rate constants', col) 942 self.propertied_connect_check('perturb_rates', self.chkPerturbRates) 943 self.chkPerturbChannelCount = pack_check('channel count', col) 944 self.propertied_connect_check('perturb_channelCount', self.chkPerturbChannelCount) 945 self.chkPerturbAmp = pack_check('Classes:Amp', col) 946 self.propertied_connect_check('perturb_amp', self.chkPerturbAmp) 947 self.chkPerturbStd = pack_check('Classes:Std', col) 948 self.propertied_connect_check('perturb_std', self.chkPerturbStd) 949 self.chkPerturbCond = pack_check('Classes:Cond', col) 950 self.propertied_connect_check('perturb_cond', self.chkPerturbCond) 951 self.chkPerturbCondStd = pack_check('Classes:CondStd', col) 952 self.propertied_connect_check('perturb_condstd', self.chkPerturbCondStd) 953 col = pack_item(gtk.VBox(), columns) 954 line = pack_item(gtk.HBox(), col) 955 pack_label(' randomly by ', line) 956 self.txtPct = pack_item(qubx.GTK.NumEntry(self.perturb_pct, acceptFloatGreaterThan(0.0), '%.0f', width_chars=4), line) 957 self.propertied_connect_NumEntry('perturb_pct', self.txtPct) 958 pack_label('%, then ', line) 959 col = pack_item(gtk.VBox(), columns) 960 chkEval = pack_radio('evaluate or', col) 961 chkOpt = pack_radio('optimize', col, group=chkEval) 962 self.propertied_connect_radios('optimize', [(False, chkEval), (True, chkOpt)]) 963 self.btnRun = pack_button('Run', col, self.__onClickRun, at_end=True) 964 line = pack_item(gtk.HBox(), col, at_end=True) 965 pack_label('Trials:', line) 966 self.txtTrials = pack_item(qubx.GTK.NumEntry(self.trials, acceptIntGreaterThan(0), width_chars=2), line, expand=True) 967 self.propertied_connect_NumEntry('trials', self.txtTrials)

968 - def __onClickRun(self, btn):

969 qubx.pyenv.env.OnScriptable('%s.run()' % self.global_name) 970 self.run(wait=False)

971 - def run(self, receiver=lambda: None, wait=True):

972 if wait: 973 return qubx.pyenv.call_async_wait(self.run, wait=False) 974 self.btnRun.set_sensitive(False) 975 task = PerturbTask(self.optimize, gen_perturbations(qubx.global_namespace.QubX.Models.file, 976 self.all_at_once, self.perturb_rates, self.perturb_channelCount, 977 self.perturb_amp, self.perturb_std, self.perturb_cond, self.perturb_condstd, 978 self.perturb_pct, self.trials), 979 receiver) 980 task.OnException += task_exception_to_console 981 task.start()

982

983 -def gen_perturbations(model, all_at_once, perturb_rates, perturb_channelCount, perturb_amp, perturb_std, perturb_cond, perturb_condstd, 984 perturb_pct, trials):

985 orig_tree = model.as_tree() 986 yield ('Initial', model) 987 def perturb(x): 988 dev_pct = random.gauss(0, perturb_pct) 989 return (x + (dev_pct / 100.0) * x, dev_pct)

990 991 rates = model.rates 992 classes = model.classes 993 states = model.states 994 classes_used = set(states[i, 'Class'] for i in xrange(len(states))) 995 for t in xrange(1, trials+1): 996 if all_at_once: 997 model.from_tree(orig_tree) 998 if perturb_rates: 999 for rate in rates: 1000 if not all_at_once: 1001 model.from_tree(orig_tree) 1002 rates[rate.Index, 'k0'], dev_pct = perturb(rate.k0) 1003 if not all_at_once: 1004 yield ("k0 %i->%i + %.1f%% (%i)" % (rate.From, rate.To, dev_pct, t), model) 1005 if rate.k1 and rate.Voltage: 1006 if not all_at_once: 1007 model.from_tree(orig_tree) 1008 rates[rate.Index, 'k1'], dev_pct = perturb(rate.k1) 1009 if not all_at_once: 1010 yield ("k1 %i->%i + %.1f%% (%i)" % (rate.From, rate.To, dev_pct, t), model) 1011 if rate.k2 and rate.Pressure: 1012 if not all_at_once: 1013 model.from_tree(orig_tree) 1014 rates[rate.Index, 'k2'], dev_pct = perturb(rate.k2) 1015 if not all_at_once: 1016 yield ("k2 %i->%i + %.1f%% (%i)" % (rate.From, rate.To, dev_pct, t), model) 1017 if perturb_channelCount: 1018 if not all_at_once: 1019 model.from_tree(orig_tree) 1020 cc, dev_pct = perturb(model.channelCount) 1021 model.channelCount = max(1, int(round(cc))) 1022 if not all_at_once: 1023 yield ("Channel count + %.1f%% (%i)" % (dev_pct, t), model) 1024 for flag, field in [(perturb_amp, 'Amp'), (perturb_std, 'Std'), (perturb_cond, 'Cond'), (perturb_condstd, 'CondStd')]: 1025 if flag: 1026 for c in sorted(classes_used): 1027 if not all_at_once: 1028 model.from_tree(orig_tree) 1029 classes[c, field], dev_pct = perturb(classes[c, field]) 1030 if not all_at_once: 1031 yield ("%s[%i] + %.1f%% (%i)" % (field, c, dev_pct, t), model) 1032 if all_at_once: 1033 yield ('%i'%t, model) 1034

1035 -class PerturbTask(qubx.task.Task):

1036 - def __init__(self, optimize, lbl_models, receiver=lambda: None):

1037 qubx.task.Task.__init__(self, 'Perturb') 1038 self.__ref_onInterrupt = self.__onInterrupt 1039 self.OnInterrupt += self.__ref_onInterrupt 1040 self.__stop_flag = False 1041 self.__first_trial = True # Initial model, always eval 1042 QubX = qubx.global_namespace.QubX 1043 QubX.Trials.get_trialset('Perturb').clear() 1044 QubX.Trials.trialset_by_name = 'Perturb' 1045 self.receiver = receiver 1046 self.optimizer = QubX.Modeling.Utils.optimizer 1047 self.optimize = optimize 1048 self.lbl_models = lbl_models # generator: (label, model) to be worked on gui thread

1049 - def __onInterrupt(self, task, cancel):

1050 cancel() 1051 self.__stop_flag = True

1052 - def run_next(self, **kw):

1053 def on_trial(model_tree, row): 1054 row['Label'] = label 1055 qubx.global_namespace.QubX.Trials.get_trialset('Perturb').append_trial(model_tree, row)

1056 try: 1057 label, model = self.lbl_models.next() 1058 if (not self.__first_trial) and self.optimize and not self.__first_trial: 1059 action = self.optimizer.optimize 1060 else: 1061 action = self.optimizer.get_ll 1062 self.__first_trial = False 1063 action(self.optimizer.get_model_prep(model), on_trial=on_trial, **kw) 1064 return True 1065 except StopIteration: 1066 return False 1067 except: 1068 traceback.print_exc() 1069 return False

1070 - def run(self):

1071 qubx.task.Tasks.add_task(self) 1072 self.status = 'Running...' 1073 try: 1074 while self.idle_wait(self.run_next): 1075 if self.__stop_flag: 1076 break 1077 finally: 1078 qubx.global_namespace.QubX.Modeling.Utils.Perturb.btnRun.set_sensitive(True) 1079 gobject.idle_add(qubx.global_namespace.QubX.Modeling.Utils.request_show) 1080 gobject.idle_add(self.receiver) 1081 qubx.task.Tasks.remove_task(self)

1082

Source Code for Module qubx.modeling_utils