1 """Provides a visual panel (L{qubx.faces.Face}) to idealize single-molecule data.
2
3 /* Copyright 2008-2014 Research Foundation State University of New York */
4
5 /* This file is part of QUB Express. */
6
7 /* QUB Express is free software; you can redistribute it and/or modify */
8 /* it under the terms of the GNU General Public License as published by */
9 /* the Free Software Foundation, either version 3 of the License, or */
10 /* (at your option) any later version. */
11
12 /* QUB Express is distributed in the hope that it will be useful, */
13 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
14 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
15 /* GNU General Public License for more details. */
16
17 /* You should have received a copy of the GNU General Public License, */
18 /* named LICENSE.txt, in the QUB Express program directory. If not, see */
19 /* <http://www.gnu.org/licenses/>. */
20
21 """
22
23 import itertools
24 import gc
25 import os
26 import sys
27 import time
28 import traceback
29 import gobject
30
31 from qubx.data_types import *
32 from qubx.util_types import *
33 from qubx.faces import *
34 from qubx.GTK import *
35 import qubx.notebook
36 import qubx.pyenv
37 import qubx.settings
38 import qubx.tree
39 from qubx.model import *
40 from qubx.settings import Property, Propertied
41 from qubx.task import *
42 import qubopt.amp
43 import qubopt.skm
44
45
46 @Propertied(Property('Method', 'Viterbi', 'Name of chosen idealizer'))
47 -class IdlFace(Face):
48 """Hosts idealization methods, many from qubx_single_molecule plugin. A method is provided by an Idealizer,
49 which is a gtk container widget with these properties:
50
51 - has .shared_controls, an empty container with space for "Method: [ whatever |v] [ Idealize ]"
52 - has .method (string) set by IdlFace upon selecting a method from the menu, to allow one Idealizer to provide multiple methods
53 - has .idealize(segments, model, on_finish), with
54 segments: list of L{qubx.data_types.SourceSeg}
55 model: L{qubx.model.QubModel}
56 on_finish: you must call gobject.idle_add(on_finish) upon completion (or failure)
57 You're expected to do work in another thread; see qubx.task.Task or qubx.task.Robot.
58
59 @ivar methods: dict(name : Idealizer)
60 @ivar idealizer: user's choice
61 """
62
63 __explore_featured = ['methods', 'idealizer', 'controls', 'register_method', 'unregister_method', 'idealize']
64
66 Face.__init__(self, 'Idealize', 'QubX.Modeling.Idealize')
67 self.__ref = Reffer()
68 self.propertied_connect_settings('Idealize')
69 self.__requested_method = self.Method
70 self.__running = False
71 self.set_size_request(200, 50)
72 self.QubX = qubx.pyenv.env.globals['QubX']
73
74 self.methods = {}
75 self.idealizer = None
76
77 qubx.notebook.Notebook.register_auto('Idealize.Messages', 'Messages, on Idealize', True)
78 qubx.notebook.Notebook.register_auto('Idealize.Model.Picture', 'Model picture, on Idealize', True)
79 qubx.notebook.Notebook.register_auto('Idealize.Model.PictureNO', 'Model picture, no overlays, on Idealize')
80 qubx.notebook.Notebook.register_auto('Idealize.Data.Picture', 'Data picture, on Idealize', True)
81 qubx.notebook.Notebook.register_auto('Idealize.Data.PictureNO', 'Data picture, no overlays, on Idealize')
82 qubx.notebook.Notebook.register_auto('Idealize.DurHist.Charts', 'DurHist charts, on Idealize')
83 qubx.notebook.Notebook.register_auto('Idealize.AmpHist.Charts', 'AmpHist charts, on Idealize')
84 qubx.notebook.Notebook.register_auto('Idealize.Model.Classes', 'Classes table, on Idealize')
85
86 row = self.controls = gtk.HBox()
87 self.mnuMethod = pack_item(StaticComboList([]), row, expand=True)
88 self.mnuMethod.OnChange += self.__ref(self.__onChangeMethod)
89 self.btnIdl = pack_button('Idealize', row, on_click=self.__onClickIdl)
90
92 """Registers an Idealizer plugin panel."""
93 if name in self.methods:
94 raise KeyError('Idealize method "%s" is already registered by %s' % (name, self.methods[name].__class__))
95 self.methods[name] = panel
96 self.mnuMethod.add_choice(name)
97 if (len(self.methods) == 1) or (self.__requested_method == name):
98 self.Method = name
100 """Removes a registered plugin panel."""
101 if (not (panel is None)) and (self.methods[name] != panel):
102 raise ValueError('Idealize method "%s" is not registered to this panel' % name)
103 self.mnuMethod.remove_choice(name)
104 del self.methods[name]
132 - def idealize(self, idealizer=None, segments=None, model=None, wait=True):
133 """Runs an idealizer. All params default to GUI values.
134
135 @param idealizer: one of the values in .methods
136 @param segments: list of L{qubx.data_types.SourceSeg}
137 @param model: L{qubx.model.QubModel}
138 @param wait: False to return immediately while idealization happens in background
139 """
140 if wait:
141 qubx.pyenv.call_async_wait(lambda recv, idl, segs, mdl: self.__start_idl(idl, segs, mdl, recv),
142 idealizer, segments, model)
143 else:
144 self.__start_idl(idealizer, segments, model)
145 - def __start_idl(self, idealizer, segs, model, receiver=lambda: None):
167 qubx.notebook.Notebook.send(self.QubX.Models.view.nbPicture, auto_id='Idealize.Model.Picture')
168 qubx.notebook.Notebook.send(self.QubX.Models.view.nbPictureNO, auto_id='Idealize.Model.PictureNO')
169 qubx.notebook.Notebook.send(self.QubX.Tables.find_view('Classes').nbTable, auto_id='Idealize.Model.Classes')
170 qubx.notebook.Notebook.send(self.QubX.Data.view.hires.nbPicture, auto_id='Idealize.Data.Picture')
171 qubx.notebook.Notebook.send(self.QubX.Data.view.hires.nbPictureNO, auto_id='Idealize.Data.PictureNO')
172
173
174 try:
175 self.QubX.Figures.index('DurHist')
176 self.QubX.Figures.DurHist.update(force=True)
177 self.QubX.Figures.DurHist.robot.do(self.__dhrobot_send_durhists)
178 except KeyError:
179 pass
180
181 self.QubX.Figures.AmpHist.update(force=True)
182 self.QubX.Figures.AmpHist.robot.do(self.__ahrobot_send_amphists)
183
184 if self.__out_list:
185 rcv = self.__receiver
186 rcv_meas = lambda tbl, off, ct: rcv()
187 if len(self.__segs) == len(self.__out_list):
188 self.QubX.Data.view.measure_list(self.__out_list, script_name='Idealized.py', wait=False, receiver=rcv_meas)
189 else:
190 self.QubX.Data.view.measure_segments(self.__segs, table=self.__out_list, row_offset=self.__out_list_ix,
191 script_name='Idealized.py', receiver=rcv_meas)
192 self.__receiver = None
193
194 self.btnIdl.set_sensitive(True)
195 if self.__receiver:
196 self.__receiver()
197 self.__receiver = self.__segs = self.__out_list = None
198
206 self.QubX.Figures.AmpHist.plots[0].controls.robot.sync()
207 qubx.notebook.Notebook.send(self.QubX.Figures.AmpHist.plots[0].nbChartFit, auto_id='Idealize.AmpHist.Charts')
208
209
210
211 @Propertied(Property('auto_init', 1, 'False to start from Classes:Amp'),
212 Property('auto_init_up', 1, 'False to open in the negative direction'),
213 Property('update_model', 1, 'True to update Classes:Amp and Std'))
215 __explore_featured = ['global_name', 'robot', 'txtStatus', 'outStatus', 'shared_controls', 'idealize']
217 gtk.HBox.__init__(self)
218 self.global_name = 'QubX.Modeling.Idealize.methods["Baum-Welch"]'
219 self.__ref = Reffer()
220 self.propertied_connect_settings('Idealize_BW')
221 self.robot = Robot('Baum-Welch', self.__ref(lambda: Tasks.add_task(self.robot)),
222 self.__ref(lambda: Tasks.remove_task(self.robot)))
223 self.robot.OnException += self.__ref(self.__onException)
224 self.robot.OnInterrupt += self.__ref(self.__onInterrupt)
225
226 self.set_size_request(200, 50)
227 self.txtStatus = gtk.TextView()
228 self.outStatus = TextViewAppender(self.txtStatus)
229 self.txtStatus.set_editable(False)
230 SetFixedWidth(self.txtStatus)
231 pack_scrolled(self.txtStatus, self, expand=True)
232 self.outStatus.write("""See http://en.wikipedia.org/wiki/Baum%E2%80%93Welch_algorithm
233
234 """)
235
236 column = pack_item(gtk.VBox(), self)
237 chk = pack_check('Auto-init', column)
238 self.propertied_connect_check('auto_init', chk)
239 chk = pack_check('Opening in + direction (up)', column)
240 self.propertied_connect_check('auto_init_up', chk)
241 chk = pack_check('Update model (Classes:Amp, Std)', column)
242 self.propertied_connect_check('update_model', chk)
243 self.shared_controls = pack_item(gtk.VBox(), column, at_end=True)
244
245 - def idealize(self, segments, model, on_finish):
262 - def robot_idealize(self, segments, config, model, update_model, on_finish):
270 - def __show_result(self, config, segments, model, update_model, result):
271 if update_model:
272 amp = result['amp'].data
273 sd = result['sd'].data
274 cc = update_model.classes
275 for i in xrange(min(len(amp), len(sd))):
276 cc[i, 'Amp'] = amp[i]
277 cc[i, 'Std'] = sd[i]
278
279 if segments:
280 idl_updater = segments[0].file.update_idl(segments[0].signal)
281 for seg, res in izip(segments, qubx.tree.children(result['DataSet'], 'Segment')):
282 ndwell = res.find('DwellCount').data
283 ndwell = ndwell[0] if len(ndwell) else 0
284 if ndwell:
285 idl_updater.set_dwells(ndwell, res['Firsts'].storage.data, res['Lasts'].storage.data, res['Classes'].storage.data,
286 res['amp'].data[:], res['sd'].data[:])
287 if segments:
288 idl_updater.done()
289
290 if self.__out_list:
291 lst = self.__out_list
292 off = self.__out_list_ix
293 for s, res in enumerate(qubx.tree.children(result['DataSet'], 'Segment')):
294 amp = res['amp'].data
295 sd = res['sd'].data
296 for i in xrange(len(amp)):
297 lst[s+off, 'BW Amp %i' % i] = amp[i]
298 lst[s, 'BW Std %i' % i] = sd[i]
299 if i:
300 lst[s+off, 'BW Delta Amp %i' % i] = (amp[i] - amp[0])
301 lst[s+off, 'BW Exc Std %i' % i] = sqrt(max(0.0, sd[i]**2 - sd[0]**2))
302
303 buf = cStringIO.StringIO()
304 buf.write("Baum-Welch (AMP) results [Class: mean +/- std]:\n\n")
305 for seg, res in izip(segments, qubx.tree.children(result["DataSet"], "Segment")):
306 if len(segments) > 1:
307 buf.write("Segment %i:\n" % (seg.index+1))
308 amp = res["amp"].data
309 sd = res["sd"].data
310 for i in xrange(min(len(amp), len(sd))):
311 buf.write(" %i: %.4g +/- %.4g\n" % (i, amp[i], sd[i]))
312 buf.write("\n")
313
314 amp = result["amp"].data
315 sd = result["sd"].data
316 if len(amp) > 1:
317 exc_amp = [a - amp[0] for a in amp]
318 exc_sd = [sqrt(max(0.0, s**2 - sd[0]**2)) for s in sd]
319 buf.write("Delta amp and excess std:\n\n")
320 for i in xrange(1, len(exc_amp)):
321 buf.write(" %i: %.4g +/- %.4g\n" % (i, exc_amp[i], exc_sd[i]))
322 buf.write("\n")
323 msgs = buf.getvalue()
324 self.outStatus.write(msgs)
325
326 qubx.notebook.Notebook.send(qubx.notebook.NbText("""Idealize: %s
327 model: %s
328 method: Baum-Welch
329 auto-init: %s
330
331 %s""" % (segments[0].file.path,
332 model.path,
333 bool(config['AutoInit'].data[0]),
334 msgs)),
335 auto_id='Idealize.Messages')
337 traceback.print_exception(typ, val, trace, file=self.outStatus)
341 gobject.idle_add(self.outStatus.write, "%s\n" % s)
342
343
344
345 @Propertied(Property('fixed_amp', [], 'List of classes to exempt from amp reestimation'),
346 Property('fixed_sd', [], 'List of classes to exempt from std dev. reestimation'),
347 Property('max_iterations', 10, 'maximum skm iterations'),
348 Property('fix_q', 0, 'True to hold rates (A matrix) constant'),
349 Property('conv_ll', 0.01, 'Convergence LL criterion'),
350 Property('drop_first_last', 0, 'True to drop first and last event, on the theory they are incomplete'),
351 Property('drop_if_class', 0, 'True to drop only first/last events of a certain class'),
352 Property('drop_class', 0, 'Class of first/last events to be dropped'))
354 __explore_featured = ['global_name', 'robot', 'txtStatus', 'outStatus', 'shared_controls', 'idealize', 'skm_controls']
356 gtk.HBox.__init__(self)
357 self.global_name = 'QubX.Modeling.Idealize.methods["SKM"]'
358 self.__ref = Reffer()
359 self.__method = 'SKM'
360 self.propertied_connect_settings('Idealize_SKM')
361 try:
362 len(self.fixed_amp)
363 except:
364 self.fixed_amp = [self.fixed_amp]
365 try:
366 len(self.fixed_sd)
367 except:
368 self.fixed_sd = [self.fixed_sd]
369 self.robot = Robot('SKM', self.__ref(lambda: Tasks.add_task(self.robot)),
370 self.__ref(lambda: Tasks.remove_task(self.robot)))
371 self.robot.OnException += self.__ref(self.__onException)
372 self.robot.OnInterrupt += self.__ref(self.__onInterrupt)
373
374 self.skm_controls = []
375 self.set_size_request(200, 50)
376 column = pack_item(gtk.VBox(), self, expand=True)
377 h = pack_item(gtk.HBox(), column)
378 self.skm_controls.append(h)
379 tip = 'To selectively prevent reestimation, enter a list of Index values (Classes table)'
380 pack_label('Fixed amps:', h)
381 txt = pack_item(qubx.GTK.NumEntry([], acceptIntList(acceptIntGreaterThanOrEqualTo(0)), formatList("%s"), width_chars=16), h)
382 txt.set_tooltip_text(tip)
383 self.propertied_connect_NumEntry('fixed_amp', txt)
384 pack_label('Fixed stds:', h)
385 txt = pack_item(qubx.GTK.NumEntry([], acceptIntList(acceptIntGreaterThanOrEqualTo(0)), formatList("%s"), width_chars=16), h)
386 txt.set_tooltip_text(tip)
387 self.propertied_connect_NumEntry('fixed_sd', txt)
388
389 self.txtStatus = gtk.TextView()
390 self.outStatus = TextViewAppender(self.txtStatus)
391 self.txtStatus.set_editable(False)
392 SetFixedWidth(self.txtStatus)
393 pack_scrolled(self.txtStatus, column, expand=True)
394 self.outStatus.write("""Based on "Restoration of single-channel currents using the segmental k-means method based on hidden Markov modeling," Qin F, Biophys J. 2004 Mar; 86(3):1488-501.
395
396 Before idealizing, make sure to "Grab" the amplitude distributions)
397 (right-click a state in the model).
398
399 """)
400
401 column = pack_item(gtk.VBox(), self)
402 chk = pack_check("Don't re-estimate rates", column)
403 self.skm_controls.append(chk)
404 self.propertied_connect_check('fix_q', chk)
405 h = pack_item(gtk.HBox(), column)
406 pack_label('Max iterations:', h)
407 txt = pack_item(qubx.GTK.NumEntry(10, acceptIntGreaterThan(0), width_chars=7), h, at_end=True)
408 self.propertied_connect_NumEntry('max_iterations', txt)
409 h = pack_item(gtk.HBox(), column)
410 pack_label('Conv LL:', h)
411 txt = pack_item(qubx.GTK.NumEntry(1.0, acceptFloatGreaterThan(0.0), '%.3g', width_chars=7), h, at_end=True)
412 self.propertied_connect_NumEntry('conv_ll', txt)
413 chk = pack_check('Drop first/last event', column)
414 self.propertied_connect_check('drop_first_last', chk)
415 h = pack_item(gtk.HBox(), column)
416 chk = pack_check(' only if Class ==', h)
417 self.propertied_connect_check('drop_if_class', chk)
418 txt = pack_item(qubx.GTK.NumEntry(10, acceptIntGreaterThanOrEqualTo(0), width_chars=7), h, at_end=True)
419 self.propertied_connect_NumEntry('drop_class', txt)
420
421 self.shared_controls = pack_item(gtk.VBox(), column, at_end=True)
422
424 self.__method = x
425 skm = self.__method == 'SKM'
426 for widget in self.skm_controls:
427 if skm:
428 widget.show()
429 else:
430 widget.hide()
431 method = property(lambda self: self.__method, lambda self, x: self.set_method(x))
432
433
434
435
436 - def idealize(self, segments, model, on_finish):
437 QubX = qubx.global_namespace.QubX
438 segs = self.__segs = QubX.DataSource.get_segmentation()
439 if ((QubX.DataSource.source == qubx.data_types.DATASOURCE_LIST) or
440 ( (QubX.DataSource.source == qubx.data_types.DATASOURCE_SCREEN) and
441 (len(segs) == 1) and
442 QubX.Data.file.list and (0 <= QubX.Data.file.list.find_selection(segs[0].f, segs[0].l)))):
443 self.__out_list = QubX.Data.file.list
444 self.__out_list_ix = QubX.Data.file.list.find_selection(segs[0].f, segs[0].l) if (QubX.DataSource.source == qubx.data_types.DATASOURCE_SCREEN) else 0
445 else:
446 self.__out_list = None
447
448 self.__idl_method = self.__method
449 if model.states and segments:
450 config = qubopt.skm.get_skm_config(model, self.fixed_amp, self.fixed_sd, self.method == 'SKM', self.max_iterations,
451 self.conv_ll, self.fix_q, self.drop_first_last, self.drop_if_class, self.drop_class)
452 self.__config = config
453 self.robot.do(self.robot_idealize, segments, config, model, on_finish)
454 else:
455 gobject.idle_add(on_finish)
465 if segments:
466 idl_updater = segments[0].file.update_idl(segments[0].signal)
467 if model.channelCount == 1:
468 class_order = sorted(list(set(st.Class for st in model.states)))
469 for seg, res in izip(segments, qubx.tree.children(result['DataSet'], 'Segment')):
470 ndwell = res.find('DwellCount').data
471 ndwell = ndwell[0] if len(ndwell) else 0
472 if ndwell:
473 ff = numpy.zeros(shape=(ndwell,), dtype='int32')
474 ll = numpy.zeros(shape=(ndwell,), dtype='int32')
475 dd = res['Durations'].storage.data
476 cc = res['Classes'].storage.data
477 amp = res['amp'].data[:]
478 sd = res['sd'].data[:]
479 if model.channelCount == 1:
480 for di in xrange(ndwell):
481 cc[di] = class_order[cc[di]]
482 nn = dd.reshape((ndwell,)) / (seg.file.sampling * 1e3)
483 nn += 0.5
484 nn = nn.astype('int32')
485 at = res.data[0]
486 for i in xrange(ndwell):
487 ff[i] = at
488 at += nn[i]
489 ll[i] = at - 1
490 if seg is segments[-1]:
491 Nclass = len(amp)
492 counts = numpy.zeros(shape=(Nclass,), dtype='float32')
493 for i in xrange(ndwell):
494 counts[cc[i]] += dd[i]
495 counts /= numpy.sum(counts)
496 ampplot = qubx.global_namespace.QubX.Figures.AmpHist.plots[0]
497 ampbot = ampplot.controls.robot
498 ampbot.set_curve(qubx.fit.Curves['GaussianUnit (%i)' % min(5, Nclass)])
499 for c in xrange(min(5, Nclass)):
500 ampbot.set_param('Amp%i'%(c+1), value=counts[c])
501 ampbot.set_param('Mean%i'%(c+1), value=amp[c])
502 ampbot.set_param('Std%i'%(c+1), value=sd[c])
503 idl_updater.set_dwells(ndwell, ff, ll, cc, amp, sd)
504
505 if segments:
506 idl_updater.done()
507
508 if self.__out_list and (self.__idl_method == 'SKM'):
509 lst = self.__out_list
510 off = self.__out_list_ix
511 for s, res in enumerate(qubx.tree.children(result['DataSet'], 'Segment')):
512 amp = res['amp'].data
513 sd = res['sd'].data
514 for i in xrange(len(amp)):
515 lst[s+off, 'SKM Amp %i' % i] = amp[i]
516 lst[s, 'SKM Std %i' % i] = sd[i]
517 if i:
518 lst[s+off, 'SKM Delta Amp %i' % i] = (amp[i] - amp[0])
519 lst[s+off, 'SKM Exc Std %i' % i] = sqrt(max(0.0, sd[i]**2 - sd[0]**2))
520
521 buf = cStringIO.StringIO()
522 if (self.__idl_method == 'SKM'):
523 buf.write("SKM results [Class: mean +/- std, observed A matrix]:\n\n")
524 for seg, res in izip(segments, qubx.tree.children(result["DataSet"], "Segment")):
525 if len(segments) > 1:
526 buf.write("Segment %i:\n" % (seg.index+1))
527 amp = res["amp"].data
528 sd = res["sd"].data
529 for i in xrange(min(len(amp), len(sd))):
530 buf.write(" %i: %.4g +/- %.4g\n" % (i, amp[i], sd[i]))
531 buf.write("A:%s\n" % res["A"].storage.data)
532 buf.write("\n")
533
534 qubx.global_namespace.skm_result = result
535
536 amp = result["amp"].data
537 sd = result["sd"].data
538 if len(amp) > 1:
539 exc_amp = [a - amp[0] for a in amp]
540 exc_sd = [sqrt(max(0.0, s**2 - sd[0]**2)) for s in sd]
541 buf.write("Delta amp and excess std:\n\n")
542 for i in xrange(1, len(exc_amp)):
543 buf.write(" %i: %.4g +/- %.4g\n" % (i, exc_amp[i], exc_sd[i]))
544 buf.write("\n")
545 msgs = buf.getvalue()
546 if msgs:
547 self.outStatus.write(msgs)
548
549 qubx.notebook.Notebook.send(qubx.notebook.NbText("""Idealize: %s
550 model: %s
551 method: %s
552 fixed amps: %s
553 fixed stds: %s
554 don't reestimate rates: %s
555 max iterations: %s
556 conv ll: %s
557 drop first/last: %s
558 drop only matching class: %s
559 drop class: %s
560
561 %s""" % (segments[0].file.path,
562 model.path,
563 self.__idl_method,
564 ", ".join(str(x) for x in config["FixAmp"].data[:]),
565 ", ".join(str(x) for x in config["FixSD"].data[:]),
566 repr(bool(config["FixQ"].data[0])),
567 config["MaxIterations"].data[0],
568 config["ConvLL"].data[0],
569 repr(bool(config["DropFirstLast"].data[0])),
570 repr(bool(config["DropIfClass"].data[0])),
571 config["DropClass"].data[0],
572 msgs)),
573 auto_id='Idealize.Messages')
575 traceback.print_exception(typ, val, trace, file=self.outStatus)
577 cancel()
578 self.__config["StopFlag"].data[0] = 1
580 gobject.idle_add(self.outStatus.write, "%s\n" % s)
581