qubx_single.dur

1 #/* Copyright 2008-2014 Research Foundation State University of New York */ 2 3 #/* This file is part of QUB Express. */ 4 5 #/* QUB Express is free software; you can redistribute it and/or modify */ 6 #/* it under the terms of the GNU General Public License as published by */ 7 #/* the Free Software Foundation, either version 3 of the License, or */ 8 #/* (at your option) any later version. */ 9 10 #/* QUB Express is distributed in the hope that it will be useful, */ 11 #/* but WITHOUT ANY WARRANTY; without even the implied warranty of */ 12 #/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ 13 #/* GNU General Public License for more details. */ 14 15 #/* You should have received a copy of the GNU General Public License, */ 16 #/* named LICENSE.txt, in the QUB Express program directory. If not, see */ 17 #/* <http://www.gnu.org/licenses/>. */ 18 19 from __future__ import with_statement 20 21 from cStringIO import StringIO 22 from qubx.hist import * 23 from qubx.model import * 24 from qubx_single.rates import * 25 import qubx.fit 26 import qubx.pyenv 27 import qubx.pyenvGTK 28 import qubx.table 29 import qubx.toolspace 30 import numpy 31 import qubx.fast.model 32 from qubx.fast.fast_utils import c_double, c_double_p, cdata, buffer_as_array_of_int 33 from qubx.settings import Property, Propertied 34 35 NBIN_DUR = 50 36 LBL_TDEAD = 'Tdead (ms)' 37 38 MAX_CHANNEL_COUNT = 8 39 MAX_DURHIST_CLASSES = 32 40 41 COLOR_TCRIT = ('qubx_single.dur_hist.tcrit', (1, .2, .2, .7)) 42 ColorInfo[COLOR_TCRIT[0]].label = 'Figures:DurHist Tcrit line' 43 44 45 @Propertied(Property('bins', 96, 'number of bins'), 46 Property('smooth', True, 'True to spread short events among adjacent bins, adjusting for sampling error'), 47 Property('sqrt', True, 'True to display sqrt(y axis)'))

48 -class DurHistFace(HistFace):

49 __explore_featured = ['lo', 'hi', 'auto', 'nbTimeConst', 'signals', 'controls', 'save_as', 50 'copy_to_clipboard', 'copy_image_to_clipboard', 'save_to', 'fit_all_exponentials']

51 - def __init__(self):

52 HistFace.__init__(self, 'DurHist', 'QubX.Figures.DurHist') 53 self.__ref = Reffer() 54 self.propertied_connect_settings('qubx_single.dur_hist') 55 self.__lo = 0.0 56 self.__hi = 1.0 57 self.__auto = True 58 self.__in_all_exp = False 59 self.__last_path = '' 60 self.file_i_bins = 0 61 self.file_i_bars = 1 62 self.file_right_edge = True 63 self.file_sqrt = True 64 self.file_norm = False 65 self.cancel_all_exp = False 66 self.OnChangePage += self.__ref(self.__onChangePage) 67 self.OnAddPlot += self.__ref(self.__onAddPlot) 68 self.OnRemovingPlot += self.__ref(self.__onAddPlot) 69 self.nbTimeConst = qubx.notebook.NbItems("All time constants", self.global_name and '%s.nbTimeConst'%self.global_name or '') 70 self.signals = SinglePlex(SingleSignals(self.QubX.DataSource, self.QubX.Data, self.QubX.Modeling.Stimulus)) 71 self.signals.OnInvalidate += self.__ref(self.__on_invalidata) 72 self.QubX.Modeling.MILRates.model_prep.OnInvalidate += self.__ref(self.__on_invalimodel) 73 self.QubX.Modeling.MILRates.rates.OnInvalidate += self.__ref(self.__on_invalirates) 74 75 self.controls = qubx.faces.Face("DurHist", global_name='QubX.About.DurHist') 76 self.controls.show() 77 line = gtk.HBox(True) 78 line.show() 79 lbl = gtk.Label('Bins:') 80 lbl.show() 81 line.pack_start(lbl) 82 self.txtBins = qubx.GTK.NumEntry(self.bins, acceptIntGreaterThan(1)) 83 self.propertied_connect_NumEntry('bins', self.txtBins) 84 self.txtBins.show() 85 line.pack_start(self.txtBins) 86 self.controls.pack_start(line, False, True) 87 line = gtk.HBox(True) 88 line.show() 89 line = gtk.HBox() 90 line.show() 91 self.chkAuto = gtk.CheckButton('auto') 92 self.chkAuto.set_active(self.__auto) 93 self.chkAuto.connect('toggled', self.__onToggleAuto) 94 self.chkAuto.show() 95 line.pack_start(self.chkAuto, False, True) 96 lbl = gtk.Label('lo:') 97 lbl.show() 98 line.pack_start(lbl, False, True) 99 self.txtLo = qubx.GTK.NumEntry(self.__lo, format='%.2g') 100 self.txtLo.set_width_chars(4) 101 self.txtLo.OnChange += self.__ref(self.__onChangeLo) 102 self.txtLo.show() 103 line.pack_start(self.txtLo, False, True) 104 lbl = gtk.Label('hi:') 105 lbl.show() 106 line.pack_start(lbl, False, True) 107 self.txtHi = qubx.GTK.NumEntry(self.__hi, format='%.2g') 108 self.txtHi.set_width_chars(4) 109 self.txtHi.OnChange += self.__ref(self.__onChangeHi) 110 self.txtHi.show() 111 line.pack_start(self.txtHi, False, True) 112 self.controls.pack_start(line, False, True) 113 self.chkSmooth = gtk.CheckButton('Smooth') 114 self.propertied_connect_check('smooth', self.chkSmooth) 115 self.chkSmooth.show() 116 self.controls.pack_start(self.chkSmooth, False, True) 117 self.chkSqrt = gtk.CheckButton('Sqrt ordinate') 118 self.propertied_connect_check('sqrt', self.chkSqrt) 119 self.chkSqrt.show() 120 self.controls.pack_start(self.chkSqrt, False, True) 121 line = gtk.HBox() 122 line.show() 123 self.btnSave = gtk.Button('Save...') 124 self.btnSave.connect('clicked', self.__onClickSave) 125 self.btnSave.show() 126 line.pack_start(self.btnSave, False, True) 127 self.btnCopy = gtk.Button('Copy') 128 self.btnCopy.connect('clicked', self.__onClickCopy) 129 self.btnCopy.show() 130 line.pack_start(self.btnCopy, False, True) 131 self.btnCopyImage = gtk.Button('Image...') 132 self.btnCopyImage.connect('clicked', self.__onClickCopyImage) 133 self.btnCopyImage.show() 134 line.pack_start(self.btnCopyImage, False, True) 135 self.controls.pack_start(line, False, True) 136 scr = gtk.ScrolledWindow() 137 scr.set_policy(gtk.POLICY_AUTOMATIC, gtk.POLICY_AUTOMATIC) 138 scr.show() 139 self.txtStim = gtk.TextView() 140 self.outStim = qubx.GTK.TextViewAppender(self.txtStim) 141 self.txtStim.show() 142 scr.add(self.txtStim) 143 self.controls.pack_end(scr, True, True) 144 145 self.__need_data = True 146 self.__need_model = True 147 self.__need_rates = True 148 self.robot.model = QubModel() 149 self.robot.OnInterrupt += self.__ref(self.__onInterrupt) 150 151 # generally short-running, so it takes more time to create and remove the progress bar: 152 self.robot.on_start_item = lambda: None 153 self.robot.on_finish_item = lambda: None 154 155 gobject.idle_add(self.update) 156 157 self.Ns = self.Nc = 0 158 self.path = None

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

160 super(DurHistFace, self).propertied_set(value, name) 161 if name in ('bins', 'smooth', 'sqrt'): 162 self.__on_invalidata()

163 - def set_auto(self, x):

164 redo = x != self.__auto 165 self.__auto = x 166 self.chkAuto.set_active(x) 167 if redo: 168 self.__need_data = True 169 self.update()

170 auto = property(lambda self: self.__auto, lambda self, x: self.set_auto(x))

171 - def set_lo(self, x, auto=False):

172 self.__lo = x 173 self.txtLo.setValue(x) 174 if not auto: 175 self.auto = False 176 self.__need_data = True 177 self.update()

178 lo = property(lambda self: self.__lo, lambda self, x: self.set_lo(x))

179 - def set_hi(self, x, auto=False):

180 self.__hi = x 181 self.txtHi.setValue(x) 182 if not auto: 183 self.auto = False 184 self.__need_data = True 185 self.update()

186 hi = property(lambda self: self.__hi, lambda self, x: self.set_hi(x))

187 - def __onToggleAuto(self, chk):

188 qubx.pyenv.env.OnScriptable('%s.auto = %s' % (self.global_name, repr(chk.get_active()))) 189 self.auto = chk.get_active()

190 - def __onChangeLo(self, txt, val):

191 qubx.pyenv.env.OnScriptable('%s.lo = %s' % (self.global_name, repr(val))) 192 self.lo = val

193 - def __onChangeHi(self, txt, val):

194 qubx.pyenv.env.OnScriptable('%s.hi = %s' % (self.global_name, repr(val))) 195 self.hi = val

196 - def __onClickSave(self, btn):

197 fname = qubx.GTK.SaveAs('Save as...', self.path or "", 'DurHist', self.save_as) 198 if fname: 199 self.path = os.path.split(fname)[0]

200 - def save_as(self, fname):

201 qubx.pyenv.env.OnScriptable('%s.save_as(%s)' % (self.global_name, repr(fname))) 202 self.save_to(open(fname, 'w'))

203 - def __onClickCopy(self, btn):

204 qubx.pyenv.env.OnScriptable('%s.copy_to_clipboard()' % self.global_name) 205 self.copy_to_clipboard()

206 - def copy_to_clipboard(self):

207 buf = cStringIO.StringIO() 208 self.save_to(buf) 209 gtk.clipboard_get('CLIPBOARD').set_text(buf.getvalue())

210 - def __onClickCopyImage(self, btn):

211 if gtk.RESPONSE_ACCEPT == TheDurHistCopyDialog.run(self.plots[0], len(self.plots)): 212 w = TheDurHistCopyDialog.copy_width 213 h = TheDurHistCopyDialog.copy_height 214 c = TheDurHistCopyDialog.copy_class 215 qubx.pyenv.env.OnScriptable('%s.copy_image_to_clipboard(w=%s, h=%s, c=%s)' % 216 (self.global_name, repr(w), repr(h), repr(c))) 217 self.copy_image_to_clipboard(w, h, c)

218 - def copy_image_to_clipboard(self, w, h, c):

219 pixmap = gdk.Pixmap(self.window, w, h, -1) 220 ctx = pixmap.cairo_create() 221 ctx.set_source_rgb(1,1,1) 222 ctx.paint() 223 self.plots[c].draw_to_context(ctx, w, h) 224 del ctx 225 pixbuf = gdk.Pixbuf(gdk.COLORSPACE_RGB, False, 8, w, h) 226 pixbuf.get_from_drawable(pixmap, gdk.colormap_get_system(), 0, 0, 0, 0, -1, -1) 227 clipboard = gtk.clipboard_get('CLIPBOARD') 228 clipboard.set_image(pixbuf) 229 del pixmap 230 del pixbuf 231 qubx.pyenv.env.gc_collect_on_idle() 232 self.plots[c].redraw_canvas() # re-calculate geometry for mouse

233 - def onShow(self, showing):

234 HistFace.onShow(self, showing) 235 if showing: 236 self.QubX.About.append_face(self.controls) 237 self.QubX.About.show_face('DurHist') 238 else: 239 self.QubX.About.remove_face( self.QubX.About.index('DurHist') )

240 - def qubx_finalize(self):

241 self.robot.stop()

242 - def __onInterrupt(self, robot, cancel):

243 if self.__in_all_exp: 244 self.cancel_all_exp = True 245 cancel()

246 - def __on_invalidata(self, *args):

247 self.__need_data = True 248 # self.robot.interrupt() 249 self.update()

250 - def __on_invalimodel(self, *args):

251 self.__need_model = True 252 # self.robot.interrupt() 253 self.update()

254 - def __on_invalirates(self, *args):

255 self.__need_rates = True 256 # self.robot.interrupt() 257 self.update()

258 - def save_to(self, f, scls=None):

259 sc = scls 260 if sc == None: 261 sc = self.page 262 for c,plot in enumerate(self.plots): 263 hist = self.histrec[sc][c] 264 f.write('Binmax\tlog10(Binmax)\t') 265 if self.sqrt: 266 f.write('sqrt(count%i/total)\t'%(c+1)) 267 else: 268 f.write('count%i/total\t'%(c+1)) 269 if hist.pdf != None: 270 f.write('PDF%i\t'%(c+1)) 271 f.write('\t'*len(hist.components)) 272 f.write('\n') 273 for b in xrange(self.bins): 274 for c,plot in enumerate(self.plots): 275 hist = self.histrec[sc][c] 276 f.write('%.8g\t' % hist.bins[b]) 277 f.write('%.8g\t' % log10(hist.bins[b])) 278 f.write('%.8g\t' % hist.bars[b]) 279 if hist.pdf != None: 280 f.write('%.8g\t' % hist.pdf[b]) 281 for m in hist.components: 282 f.write('%.8g\t' % m[b]) 283 f.write('\n')

284 - def __onChangePage(self, index):

285 if (0 <= index < len(self.histrec)) and len(self.histrec[index]): 286 self.txtStim.get_buffer().set_text('\n\n'.join([self.histrec[index][0].txt_stim] + [hh.txt_time_const for hh in self.histrec[index]]))

287 - def __onAddPlot(self, index):

288 plot = self.plots[index] 289 ref = plot.xref = Reffer() 290 plot.OnDraw += ref(bind_with_args_before(self.__onDrawPlot, index)) 291 if index == 0: 292 self.plots[0].subNotebook.items.append(self.nbTimeConst) 293 plot.nbTimeConst = qubx.notebook.NbDynText('Time constants', '%s.plots[%i].nbTimeConst' % (self.global_name, index), bind(self.__nb_time_const, index)) 294 plot.subNotebook.items.append(plot.nbTimeConst) 295 self.nbTimeConst.items.insert(index, plot.nbTimeConst) 296 plot.layAllExponentials = qubx.toolspace.Layer(x=-20, y=-2.5*qubx.fit_space.LINE_EMS, w=19, h=qubx.fit_space.LINE_EMS, cBG=plot.controls.cLayer) 297 plot.actAllExponentials = bind(self.__onClickAllExponentials, plot) 298 plot.layAllExponentials.add_sublayer(qubx.toolspace.SubLayer_Label('Fit all exponentials...', 0, 1, x=0, y=0, w=15, h=qubx.fit_space.LINE_EMS, 299 color=plot.controls.cButton, action=plot.actAllExponentials)) 300 plot.layAllExponentials.add_sublayer(qubx.toolspace.SubLayer_PNG(os.path.join(qubx.global_namespace.app_path, 'icons', 'folder.png'), x=16, y=0, w=qubx.fit_space.LINE_EMS, h=qubx.fit_space.LINE_EMS, tooltip="Open histogram from a file for fitting...", 301 action=ref(bind(self.__open_file, index)))) 302 plot.layAllExponentials.add_sublayer(qubx.toolspace.SubLayer_Label('?', 0, 1, x=18, y=0, w=1, h=qubx.fit_space.LINE_EMS, tooltip="About this algorithm...", 303 color=plot.controls.cButton, action=ref(self.__show_all_exponentials_info))) 304 plot.controls.add_layer(plot.layAllExponentials) 305 plot.layAllExpoLLR = qubx.toolspace.Layer(x=-16, y=-8*qubx.fit_space.LINE_EMS, w=15, h=qubx.fit_space.LINE_EMS, cBG=plot.controls.cLayer) 306 plot.layAllExpoLLR.add_sublayer(qubx.toolspace.SubLayer_Label('LLR: ', 1, 0, x=0, y=0, w=6, h=qubx.fit_space.LINE_EMS, color=plot.controls.cText)) 307 plot.subAllExpoLLR = qubx.toolspace.SubLayer_Label('', -1, 0, x=7, y=0, w=8, h=qubx.fit_space.LINE_EMS, color=plot.controls.cNumber) 308 plot.layAllExpoLLR.add_sublayer(plot.subAllExpoLLR) 309 plot.showingLLR = false 310 plot.layFileInfo = qubx.toolspace.Layer(x=-25, y=2, w=19, h=qubx.fit_space.LINE_EMS, cBG=plot.controls.cLayer) 311 plot.subFileName = qubx.toolspace.SubLayer_Label('', 0, 1, x=0, y=0, w=17, h=qubx.fit_space.LINE_EMS, 312 color=plot.controls.cText, 313 action=ref(bind(self.__open_file, index))) 314 plot.layFileInfo.add_sublayer(plot.subFileName) 315 plot.subFileClose = qubx.toolspace.SubLayer_Label('x', 0, 1, x=18, y=0, w=1, h=qubx.fit_space.LINE_EMS, 316 color=plot.controls.cButton, 317 action=ref(bind(self.__close_file, index))) 318 plot.layFileInfo.add_sublayer(plot.subFileClose) 319 plot.has_file = false

320 - def __onRemovingPlot(self, index):

321 del self.nbTimeConst.items[index]

322 - def __nb_time_const(self, index):

323 hist = self.histrec[self.page][index] 324 return hist.txt_time_const

325 - def __onDrawPlot(self, context, w, h, index):

326 hist = self.histrec[self.page][index] 327 if not hist.tcrit: 328 return 329 plot = self.plots[index] 330 if plot.layerset == plot.controls: 331 return 332 context.set_source_rgba(*plot.appearance.color(COLOR_TCRIT)) 333 context.set_line_width(2*plot.appearance.line_width) 334 for t in hist.tcrit: 335 x = plot.t2x(log10(1e3*t)) 336 context.move_to(x, 0) 337 context.line_to(x, h) 338 context.stroke()

339 - def robot_update(self):

340 #qubx.pyenv.env.stderr.write("%s\n"%('robot begin')) 341 self.robot.progress = 0 342 self.robot.new_data = self.__need_data 343 if self.__need_data: 344 segs, mcls, scls, names, sampling_ms, td_ms = self.robot.gui_call_recv(self.signals.request, self.robot) 345 shown_sc = max(0, min(int(self.page), len(scls)-1)) 346 self.robot_update_data(segs, mcls, scls, names, sampling_ms, td_ms, shown_sc) 347 self.robot_update_pdfs()

348 - def robot_update_data(self, segs, mcls, scls, names, sampling_ms, td_ms, shown_sc):

349 self.robot.names = names 350 sampling = sampling_ms * 1e-3 351 auto, lo, hi = self.auto, self.lo, self.hi 352 lo -= td_ms 353 hi -= td_ms 354 lo *= 1e-3 355 hi *= 1e-3 356 if not segs: return 357 self.robot.progress = 20 358 per_scls = events_per_scls(segs, mcls, scls) 359 if not per_scls: return 360 self.robot.progress = 30 361 if auto: 362 t_min = [len(du) and min(du) or 0.0 for cl, du in per_scls] 363 t_max = [len(du) and max(du) or 0.0 for cl, du in per_scls] 364 if len(t_min): 365 with self.robot.main_hold: 366 self.set_lo(t_min[shown_sc]*1e3+td_ms, auto=True) 367 self.set_hi(t_max[shown_sc]*1e3+td_ms, auto=True) 368 else: 369 t_min = len(scls) * [lo] 370 t_max = len(scls) * [hi] 371 self.robot.Nc = min(MAX_DURHIST_CLASSES, max(len(cl) and max(cl) or 0 for cl, du in per_scls) + 1) 372 stimulus = [ dict([(name, scval) for name, scval in izip(names, scvals)]) for scvals in scls ] 373 #Hist = self.smooth and DurHist_Smooth or DurHist 374 Hist = self.smooth and DurHist_Smooth_lib or DurHist_lib 375 cap_stim = [] 376 txt_stim = [] 377 for sc in xrange(len(per_scls)): 378 buf_stim = StringIO() 379 cap = "" 380 for i, nm in enumerate(names): 381 if i > 0: 382 buf_stim.write('%s: %.6g\n' % (nm, stimulus[sc][nm])) 383 if i > 1: 384 cap = cap + ', ' 385 cap = cap + '%.3g' % stimulus[sc][nm] 386 cap = cap and (' at '+cap) 387 txt_stim.append(buf_stim.getvalue()) 388 cap_stim.append(cap) 389 390 def make_rec(hist_out, sc, c): 391 bins, bars = hist_out 392 logbins = numpy.log10(bins) 393 if self.sqrt: 394 bars = numpy.sqrt(bars) 395 caption = 'DurHist of class %d%s'%(c, cap_stim[sc]) 396 return HistRec(caption, caption, logbins, bars, linbins=bins, 397 nb_headers=['log10_ms', 'count_per_n'], nb_xlabel='log10(ms)', nb_ylabel='count / n', 398 data_colorref=qubx.modelGTK.COLOR_CLASS(c), txt_stim=txt_stim[sc], txt_time_const="")

399 hist_sc_cl = [] 400 for sc, data in enumerate(per_scls): 401 bins = DurHist_Bins([data], td_ms, t_min[sc], t_max[sc], self.bins) 402 hist_of_sc = [] 403 for c in xrange(self.robot.Nc): 404 if (c < len(self.plots)) and self.plots[c].has_file: 405 hist_of_sc.append(self.histrec[sc][c]) 406 else: 407 hist_of_sc.append(make_rec(Hist(c, [data], td_ms, sampling, bins), sc, c)) 408 hist_sc_cl.append(hist_of_sc) 409 410 self.robot.progress = 45 411 with self.robot.main_hold: 412 self.stimulus = stimulus 413 self.sig_names = names 414 self.histrec = hist_sc_cl 415 self.Ns = len(scls) 416 self.td = td_ms 417 self.__need_data = False 418 # windows can't handle off-thread window creation 419 gobject.idle_add(self.set_page_count, self.Ns) 420 gobject.idle_add(self.update_plots)

421

422 - def robot_update_pdfs(self):

423 self.robot.progress = 50 424 if self.robot.new_data or self.__need_model or self.__need_rates: 425 if self.__need_model: 426 self.robot.model.undoStack.clear() 427 mt = self.QubX.Modeling.MILRates.model_prep.val.modeltree 428 if mt.find('States').find('State').isNull: 429 return 430 self.robot.model.from_tree( self.QubX.Modeling.MILRates.model_prep.val.modeltree ) 431 elif self.__need_rates: 432 self.__need_rates = False 433 rrates = self.robot.model.rates 434 K0, K1, K2 = self.QubX.Modeling.MILRates.rates.val 435 if not (K0 is None): 436 for i in xrange(rrates.size): 437 a, b = [rrates.get(i, field) for field in ['From', 'To']] 438 rrates.set(i, 'k0', K0[a,b]) 439 rrates.set(i, 'k1', K1[a,b]) 440 rrates.set(i, 'k2', K2[a,b]) 441 self.robot.progress = 70 442 self.robot.model.channelCount = min(MAX_CHANNEL_COUNT, self.robot.model.channelCount) 443 single_model = ModelToFast(self.robot.model, self.robot.names, printer) 444 multi_model = qubx.fast.model.MultiModel(single_model) 445 multi_model.update() 446 mm = multi_model.multi_model 447 clazz = buffer_as_array_of_int(mm[0].clazz, mm[0].Nstate) 448 449 classmap = sorted(list(set(st.Class for st in self.robot.model.states))) 450 451 stimclasses_ = numpy.zeros(shape=(self.Ns, len(self.robot.names)), dtype='float64') 452 for s in xrange(self.Ns): 453 for ni, nm in enumerate(self.robot.names): 454 stimclasses_[s,ni] = self.stimulus[s][nm] 455 #print stimclasses_ 456 stimclasses = cdata(stimclasses_, c_double_p) 457 ratesm = qubx.fast.model.StimRates(mm, self.Ns, len(self.robot.names), stimclasses, None, 0.0) # skip ratesm.AA generation 458 ratesm.setup_QQ() 459 for s in xrange(self.Ns): 460 Q = numpy.matrix(ratesm.QQ[s]) 461 #print Q 462 for c in xrange(self.plot_count): 463 if self.robot.model.channelCount == 1: 464 if c in classmap: 465 cl = classmap.index(c) 466 else: 467 continue 468 elif c in clazz: 469 cl = c 470 else: 471 continue 472 hist = self.histrec[s][c] 473 if (not (hist.bins is None)) and (len(hist.bins) > 1): 474 hist.pdf, hist.components, hist.time_const_tau, hist.time_const_amp, hist.tcrit = PDF(Q, clazz, cl, self.td, hist.linbins, not self.QubX.Modeling.MILRates.optimizing) 475 if self.sqrt: 476 hist.pdf = sqrt(hist.pdf) 477 hist.components = [sqrt(cc) for cc in hist.components] 478 hist.txt_time_const = "Class %i:\nTime Constants: (Tau [ms], Amp)\n" % c + \ 479 "\n".join(["%.4g, %.4g" % (1e3*tc_tau, tc_amp) for (tc_tau, tc_amp) in izip(hist.time_const_tau, hist.time_const_amp)]) 480 if hist.tcrit: 481 hist.txt_time_const += "\n\nTcrit:\n" + ", ".join("%.4g" % (1e3*crit) for crit in hist.tcrit) + "\n\n" 482 #qubx.pyenv.env.stderr.write("%s\n"%('robot model')) 483 single_model.dispose() # work around callback-related memory leak 484 self.__need_model = self.__need_rates = False 485 #qubx.pyenv.env.stderr.write("%s\n"%('robot done')) 486 gobject.idle_add(self.update_plots)

487 - def __onClickAllExponentials(self, plot):

488 if self.__in_all_exp: 489 return 490 qubx.pyenv.env.OnScriptable('QubX.Figures.DurHist.fit_all_exponentials(%i)' % self.plots.index(plot)) 491 self.fit_all_exponentials(self.plots.index(plot))

492 - def fit_all_exponentials(self, index):

493 if self.__in_all_exp: 494 return 495 self.__in_all_exp = True 496 self.cancel_all_exp = False 497 plot = self.plots[index] 498 plot.controls.remove_layer(plot.layAllExponentials) 499 ### todo: prompt, remember properties: 500 is_sqrt = plot.is_sqrt if plot.has_file else self.sqrt 501 self.robot.do(self.robot_all_exponentials, plot, is_sqrt, plot.controls.xx, plot.controls.yy, ALL_EXPONENTIALS_SIGNIFICANCE)

502 - def robot_all_exponentials(self, plot, is_sqrt, xx, yy, significance): # ignoring weight expression

503 pool = multiprocessing.Pool() 504 def map_func(func, args): 505 result = [x for x in pool.imap(func, args)] 506 return result 507 try: 508 expo = AllExponentials() 509 params = [1.0] * len(expo.params) 510 if is_sqrt: 511 yy = yy**2 512 for pct in expo.run_generator(xx, yy, params, map_func=map_func, significance=significance): 513 gobject.idle_add(self.__show_expo, plot, is_sqrt, params[:], pct > 99.5) 514 except: 515 traceback.print_exc() 516 pool.close() 517 gobject.idle_add(self.__done_expo, plot, expo.report_buf.getvalue(), expo.llr) ### plus e.g. N event required

518 - def __show_expo(self, plot, is_sqrt, params, can_fit):

519 N = int(round(params[0])) 520 robot = plot.controls.robot 521 robot.set_curve(AllExponentials_Display, "%i %sExponential%s, see Landowne et al, 2013"%(N, is_sqrt and "(sqrt) " or "", (N>1) and "s" or "")) 522 for i in xrange(N): 523 robot.set_param("tau%i"%(i+1), params[2*i+1], 1e-16, UNSET_VALUE, can_fit) 524 robot.set_param("area%i"%(i+1), params[2*i+2], 0.0, UNSET_VALUE, can_fit)

525 - def __done_expo(self, plot, report, llr):

526 if not plot.showingLLR: 527 plot.controls.add_layer(plot.layAllExpoLLR) 528 plot.subAllExpoLLR.label = '%.2g' % llr 529 # plot.controls.robot.fit() 530 plot.controls.add_layer(plot.layAllExponentials) 531 print report 532 self.__in_all_exp = False 533 pass ### show more info

534 - def __open_file(self, index):

535 qubx.GTK.Open("Open duration histogram...", self.__last_path, bind_with_args_before(self.__do_open, index), parent=self.parent_window)

536 - def __do_open(self, path, index):

537 plot = self.plots[index] 538 self.__last_path, fname = os.path.split(path) 539 headers, types, columns = qubx.table.read_table_text(open(path, 'r').read()) 540 i_bins, i_bars, right_edge, is_sqrt, is_norm = self.prompt_file_columns(headers, types, columns) 541 if i_bins is None: 542 return 543 self.__close_file(index, False) 544 plot.subFileName.label = fname if (len(fname) < 15) else (fname[:15]+'...') 545 if not plot.has_file: 546 plot.has_file = True 547 plot.controls.add_layer(plot.layFileInfo) 548 plot.is_sqrt = is_sqrt 549 bins = columns[i_bins] 550 bars = columns[i_bars] 551 if is_norm: 552 total = 0.0 553 for x in bars: 554 total += x 555 if total: 556 bars = [x/total for x in bars] 557 if not right_edge: 558 w = (bins[2] - bins[1]) 559 bins = [b+w for b in bins] 560 linbins = numpy.power(10, bins) 561 self.histrec[self.page][index] = HistRec(fname, fname, bins, bars, linbins=linbins) 562 self.update_plots()

563 - def __close_file(self, index, stayingShut=True):

564 plot = self.plots[index] 565 if plot.has_file and stayingShut: 566 plot.has_file = False 567 plot.controls.remove_layer(plot.layFileInfo) 568 self.__on_invalidata()

569 - def prompt_file_columns(self, headers, types, columns):

570 dlg = gtk.Dialog('Import duration histogram...', self.parent_window, gtk.DIALOG_MODAL) 571 h = pack_item(gtk.HBox(), dlg.vbox) 572 pack_label('Bin column: ', h) 573 print headers, types, len(columns) 574 colChoices = ['%s: %s, ...' % (name, ', '.join(['%.3g'%v for v in columns[i][:3]])) for i, name in enumerate(headers)] 575 print colChoices 576 mnuBin = pack_item(qubx.GTK.StaticComboList(colChoices), h) 577 mnuBin.active_i = self.file_i_bins 578 chkRight = pack_check('contains right edge of each bin (otherwise, the left edge)', dlg.vbox, self.file_right_edge) 579 h = pack_item(gtk.HBox(), dlg.vbox) 580 pack_label('Bar column: ', h) 581 mnuBar = pack_item(qubx.GTK.StaticComboList(colChoices), h) 582 mnuBar.active_i = self.file_i_bars 583 chkSqrt = pack_check('contains sqrt(count/total) (otherwise, count/total)', dlg.vbox, self.file_sqrt) 584 chkNorm = pack_check('normalize (divide by sum)', dlg.vbox, self.file_norm) 585 dlg.add_buttons('OK', gtk.RESPONSE_ACCEPT, 'Cancel', gtk.RESPONSE_REJECT) 586 response = dlg.run() 587 if response == gtk.RESPONSE_ACCEPT: 588 self.file_i_bins = mnuBin.active_i 589 self.file_i_bars = mnuBar.active_i 590 self.file_right_edge = chkRight.get_active() 591 self.file_sqrt = chkSqrt.get_active() 592 self.file_norm = chkNorm.get_active() 593 dlg.destroy() 594 if response == gtk.RESPONSE_ACCEPT: 595 return self.file_i_bins, self.file_i_bars, self.file_right_edge, self.file_sqrt, self.file_norm 596 return None, None, None, None, None

597 - def __show_all_exponentials_info(self, *args):

598 global ALL_EXPONENTIALS_MAX_N, ALL_EXPONENTIALS_AREA_CUTOFF, ALL_EXPONENTIALS_NEAR_THRESHOLD, ALL_EXPONENTIALS_SIGNIFICANCE, ALL_EXPONENTIALS_FINAL_CUTOFF 599 result = qubx.pyenvGTK.prompt_entries( 600 [ 601 ('Based on:', None, None, None), 602 ('Exponential Sum-Fitting of Dwell-Time Distributions without Specifying Starting Parameters', None, None, None), 603 ('David Landowne, Bin Yuan, and Karl L. Magleby', None, None, None), 604 ('Biophysical Journal Volume 104 June 2013 2383-2391', None, None, None), 605 ('Max components:', ALL_EXPONENTIALS_MAX_N, acceptIntGreaterThan(1), str), 606 ('Area cutoff:', ALL_EXPONENTIALS_AREA_CUTOFF, acceptFloatGreaterThan(0), '%.3g'), 607 ('Nearness threshold:', ALL_EXPONENTIALS_NEAR_THRESHOLD, acceptFloatGreaterThan(0), '%.3g'), 608 ('Significance:', ALL_EXPONENTIALS_SIGNIFICANCE, acceptFloatGreaterThan(0), '%.2f'), 609 ('Final area cutoff ("step 11"--enter 0 to skip):', ALL_EXPONENTIALS_FINAL_CUTOFF, acceptFloatGreaterThan(0), '%.3g') 610 ], 611 'About "Fit all exponentials"...' 612 ) 613 if not (result is None): 614 ALL_EXPONENTIALS_MAX_N = result[4] 615 ALL_EXPONENTIALS_AREA_CUTOFF = result[5] 616 ALL_EXPONENTIALS_NEAR_THRESHOLD = result[6] 617 ALL_EXPONENTIALS_SIGNIFICANCE = result[7] 618 ALL_EXPONENTIALS_FINAL_CUTOFF = result[8]

619

620 -def events_per_scls(segs, mcls, scls): # [(model_classes, durations)]

621 Ns = len(scls) 622 per_scls = [ ([], []) for x in scls ] 623 for cl, du in segs: 624 for c, d in izip(cl, du): 625 mc, sc = mcls[c] 626 if mc >= 0: 627 per_scls[sc][0].append(mc) 628 per_scls[sc][1].append(d) 629 return [(numpy.array(cc, dtype='int32'), numpy.array(dd, dtype='float32')) for cc, dd in per_scls] 630

631 632 -def DurHist_Bins(data, td_ms, t_min=None, t_max=None, nbin=NBIN_DUR):

633 tmax = t_max 634 if tmax == None: 635 tmax = max([ max([dd for cc,dd in itertools.izip(classes, durations) if cc==c]) for classes, durations in data ]) 636 tmax = tmax * 1e3 + td_ms 637 tmin = t_min 638 if tmin == None: 639 tmin = min([ min([dd for cc,dd in itertools.izip(classes, durations) if cc==c]) for classes, durations in data ]) 640 tmin = tmin * 1e3 + td_ms 641 N = nbin 642 bins = zeros(shape=(N,), dtype='float32') 643 if tmax == tmin: 644 return bins 645 dx = log10(tmax/tmin)/(N-1) 646 r = pow(10.0, dx) 647 bins[0] = tmin 648 for i in xrange(1,N): 649 bins[i] = bins[i-1] * r 650 for i in reversed(range(1,N)): 651 if bins[i] == bins[i-1]: 652 bins[i-1:N-1] = bins[i:N] 653 N -= 1 654 bins.resize(N) 655 return bins

656

657 -def DurHist(c, data, td_ms, sampling, bins):

658 N = len(bins) 659 bars = zeros(shape=(N,)) 660 tot_count = 0 661 for classes, durations in data: 662 for i in xrange(len(durations)): 663 if classes[i] != c: continue 664 dur = durations[i] * 1e3 + td_ms 665 last = 0.0 666 for b in xrange(N): 667 here = bins[b] 668 if last < dur <= here: 669 bars[b] += 1.0 670 last = here 671 tot_count += 1 672 if tot_count: 673 bars /= tot_count 674 return bins, bars

675

676 -def DurHist_lib(c, data, td_ms, sampling, bins):

677 cc, dd = data[0] 678 return bins, durhist_calhist(c, cc, dd*1e3, bins, td_ms)

679

680 -def DurHist_Smooth_lib(c, data, td_ms, sampling, bins):

681 cc, dd = data[0] 682 return bins, durhist_calhist_smooth(c, cc, dd*1e3, bins, td_ms, sampling)

683

684 -def PDF(Q, clazz, a, td, bins, wantTcrit=True):

685 together = zeros(shape=bins.shape) 686 components = [] 687 taus = [] 688 amps = [] 689 690 try: 691 Q = QtoQe(Q, clazz, td*1e-3) 692 A = clazz == a 693 Z = clazz != a 694 Qaa = Q[ix_(A, A)] 695 Qaz = Q[ix_(A, Z)] 696 Qza = Q[ix_(Z, A)] 697 Qzz = Q[ix_(Z, Z)] 698 W, V = linalg.eig( (Qaa.I*Qaz*Qzz.I*Qza).T ) 699 i_unit = argmin(abs(W-1)) 700 pe = V[:,i_unit].T 701 pe /= sum(pe) 702 703 W, V = linalg.eig(Qaa) 704 Vi = V.I 705 spectre = matrix(zeros(shape=Qaa.shape)) 706 for i in xrange(Qaa.shape[0]): 707 for j in xrange(Qaa.shape[0]): 708 for k in xrange(Qaa.shape[0]): 709 spectre[j,k] = V[j,i] * Vi[i,k] 710 coeff = sum(pe * spectre * Qaz) / W[i] 711 itau = W[i] 712 component = zeros(shape=bins.shape) 713 last = bins[0] * bins[0] / bins[1] 714 for b in xrange(len(bins)): 715 next = bins[b] 716 component[b] = coeff * (exp(itau*(next-td)*1e-3) - exp(itau*(last-td)*1e-3)) 717 together[b] += component[b] 718 last = next 719 components.append(component) 720 taus.append(-1.0/itau) 721 amps.append(-coeff) 722 except linalg.LinAlgError: 723 pass 724 sortix = [i for i in xrange(len(taus))] 725 sortix.sort(key=lambda i: taus[i]) 726 taus = [taus[i] for i in sortix] 727 amps = [amps[i] for i in sortix] 728 crits = SolveTCrits(taus, amps) if wantTcrit else [] 729 return together, components, taus, amps, crits

730

731 -def SolveTCrits(taus, amps):

732 #curve = qubx.fit.Curve("a1 * exp(-Tc / t1) - a2 * (1 - exp(-Tc / t2))") 733 #curve = qubx.fit.Curve("a1 * (exp(-Tc / t1) - 1) + a2 * (exp(-Tc / t2) - 1) - a1*t1") 734 curve = qubx.fit.Curve("a1 * exp(-Tc / t1) + a2 * exp(-Tc / t2) - a2") 735 iA1, iT1, iA2, iT2, iTc = [curve.params.index(x) for x in ["a1", "t1", "a2", "t2", "Tc"]] 736 for i in xrange(len(curve.params)): 737 if i != iTc: 738 curve.can_fit[i] = False 739 fitter = qubx.fit.SimplexFitter() 740 xx = yy = array([0.0], dtype='float32') 741 crits = [] 742 params = [0.0] * len(curve.params) 743 for i in xrange(len(taus)-1): 744 params[iA1] = amps[i] 745 params[iT1] = taus[i] 746 params[iA2] = amps[i+1] 747 params[iT2] = taus[i+1] 748 params[iTc] = (taus[i] + taus[i+1]) / 2 749 #curve.lo[iTc] = min(taus[i], taus[i+1]) 750 #curve.hi[iTc] = max(taus[i], taus[i+1]) 751 result, ssr, iterations, fit = fitter(curve, params, xx, yy) 752 crits.append(result[iTc]) 753 return crits

754 755 756 757 intersect = lambda a0, a1, b0, b1: (max(a0, b0), min(a1, b1)) 758 tri_rising_y = lambda peak, x: (x - peak + 1) 759 tri_falling_y = lambda peak, x: - ( x - peak - 1)

760 -def tri_rising_area(peak, fro, to):

761 x0, x1 = intersect(peak-1, peak, fro, to) 762 if x0 < x1: 763 y0 = tri_rising_y(peak, x0) 764 y1 = tri_rising_y(peak, x1) 765 return (0.5*(x1-x0)*(y0+y1)) 766 return 0.0

767 -def tri_falling_area(peak, fro, to):

768 x0, x1 = intersect(peak, peak+1, fro, to) 769 if x0 < x1: 770 y0 = tri_falling_y(peak, x0) 771 y1 = tri_falling_y(peak, x1) 772 return (0.5*(x1-x0)*(y0+y1)) 773 return 0.0

774 tri_area = lambda peak, fro, to: tri_rising_area(peak, fro, to) + tri_falling_area(peak, fro, to) 775 776 WAY_TOO_MANY_SAMPLES = 65536

777 778 -def spread_events(bins, bars, samples, count, sampling_ms):

779 if (count == 0) or (len(bins) == 0): return 780 a = 0.0 781 for i, b in enumerate(bins): 782 bars[i] += count * tri_area(samples, a/sampling_ms, b/sampling_ms) 783 a = b

784

785 -def DurHist_Smooth(c, data, td_ms, sampling, bins):

786 td = td_ms * 1e-3 787 sampling_ms = sampling * 1e3 788 N = len(bins) 789 bars = zeros(shape=(N,)) 790 # first bin by how many samples, 791 sampleCounts = [] 792 # but treat extra-long events separately so we don't need 200 million empty bins 793 tot_count = 0 794 for classes, durations in data: 795 for i in xrange(len(durations)): 796 if classes[i] != c: continue 797 tot_count += 1 798 samples = int(round((durations[i] + td) / sampling)) 799 if samples > WAY_TOO_MANY_SAMPLES: 800 spread_events(bins, bars, samples, 1, sampling_ms) 801 else: 802 while len(sampleCounts) <= samples: 803 sampleCounts.append(0) 804 sampleCounts[samples] += 1 805 # then distribute events among the log bins 806 # using empirical observation that events recorded as k samples long are actually between (k-1)dt and (k+1)dt, 807 # distributed with mean of k*dt and linear falloff to 0 at +/- dt 808 for samples, count in enumerate(sampleCounts): 809 spread_events(bins, bars, samples, count, sampling_ms) 810 if tot_count: 811 bars /= tot_count 812 return bins, bars

813

814 815 816 -class DurHistCopyDialog(qubx.GTK.CopyDialog):

817 - def __init__(self, parent=None):

818 qubx.GTK.CopyDialog.__init__(self, 'Copy histogram image', parent) 819 line = gtk.HBox(True) 820 line.show() 821 lbl = gtk.Label('Class:') 822 lbl.show() 823 line.pack_start(lbl) 824 self.txtClass = qubx.GTK.NumEntry(1, acceptInt) 825 self.txtClass.set_width_chars(11) 826 self.txtClass.show() 827 line.pack_start(self.txtClass) 828 self.vbox.pack_start(line, False, True)

829 - def run(self, view, Nclass):

830 self.txtClass.accept = acceptIntBetween(1, Nclass) 831 self.txtClass.value = min(Nclass, self.txtClass.value) 832 response = qubx.GTK.CopyDialog.run(self, view) 833 self.copy_class = self.txtClass.value - 1 834 return response

835 836 TheDurHistCopyDialog = DurHistCopyDialog()

837 838 839 -class SingleSignals(Requestable):

840 - def __init__(self, datasrc, Data, stimulus):

841 Requestable.__init__(self) 842 self.__ref = Reffer() 843 self.datasrc = datasrc 844 self.Data = Data 845 self.datasrc.OnChangeIdealization += self.__ref(self.__onChangeIdl) 846 self.__data = None 847 self.stimulus = stimulus 848 self.stimulus.OnInsert += self.__ref(self.__onChangeStimulus) 849 self.stimulus.OnChange += self.__ref(self.__onChangeStimulus) 850 self.stimulus.OnRemoving += self.__ref(self.__onChangeStimulus) 851 self.Data.OnSwitch += self.__ref(self.__onSwitchDataFile) 852 self.__datafile = None 853 self.__onSwitchDataFile(self.Data, self.Data.file) 854 self.__onChangeIdl()

855 - def qubx_finalize(self):

856 pass

857 - def __onChangeIdl(self, *args):

858 self.invalidate()

859 - def __onChangeStimulus(self, *args):

860 self.invalidate()

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

862 if file == self.__datafile: return 863 if self.__datafile: 864 self.__datafile.constants.OnSet -= self.__ref(self.__onSetConstant) 865 self.__datafile = file 866 if self.__datafile: 867 self.__datafile.constants.OnSet += self.__ref(self.__onSetConstant)

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

869 if self.__datafile.constants[r,'Name'] == LBL_TDEAD: 870 self.invalidate()

871 - def get(self, receiver):

872 data = qubx.pyenv.env.globals['QubX'].Data.view 873 if not data: 874 gobject.idle_add(receiver, None, None, None, None) # signals, names, sampling_ms, td_ms 875 return 876 877 for i in xrange(data.file.constants.size): 878 if LBL_TDEAD == data.file.constants.get(i, 'Name'): 879 td = data.file.constants.get(i, 'Value') 880 break 881 else: 882 td = data.file.sampling * 1e3 883 data.file.constants.append({'Name':LBL_TDEAD, 'Value':td}) 884 data.file.constants.select(data.file.constants.size-1, 'Value', sender=self) 885 886 segsrc = self.datasrc.get_segmentation() 887 sig_names = [data.file.signals.get(segsrc[0].signal, 'Name')] if segsrc else [] 888 signals = [ [] ] 889 for seg in segsrc: 890 ff, ll, cc, dd = self.datasrc.get_idealization(seg, mark_excluded=True, get_fragments=True, get_durations=True) 891 aa = data.file.ideal[segsrc[0].signal].seg[seg.index].amp 892 signals[0].append( (cc, dd, aa) ) 893 def receive_stim(stim_names, stimidl): 894 sig_names.extend(stim_names) 895 for segs_flcda in stimidl: 896 signals.append( [(cc, dd, aa) for ff,ll,cc,dd,aa in segs_flcda] ) 897 gobject.idle_add(receiver, signals, sig_names, qubx.pyenv.env.globals['QubX'].Data.file.sampling * 1e3, td)

898 self.stimulus.call_with_stimulus_idl(receive_stim, [(seg.f, seg.l) for seg in segsrc])

899

900 901 -class SinglePlex(Requestable):

902 - def __init__(self, signals):

903 Requestable.__init__(self) 904 self.__ref = Reffer() 905 self.signals = signals 906 self.signals.OnInvalidate += self.__ref(self.__onInvalidate)

907 - def __onInvalidate(self, signals):

908 self.invalidate()

909 - def get(self, receiver):

910 def receive_signals(signals, names, sampling_ms, td_ms): 911 segs, mcls, scls = ProcessSignals(td_ms, signals) 912 gobject.idle_add(receiver, segs, mcls, scls, names, sampling_ms, td_ms)

913 self.signals.request(receive_signals)

914 915 916 917 918 919 920 ALL_EXPONENTIALS_MAX_N = 20 921 ALL_EXPONENTIALS_AREA_CUTOFF = 1e-5 922 ALL_EXPONENTIALS_NEAR_THRESHOLD = 0.02 923 ALL_EXPONENTIALS_FIT_REPS = 5 924 ALL_EXPONENTIALS_FIT_ITERS = 30 925 ALL_EXPONENTIALS_SIGNIFICANCE = 0.05 926 ALL_EXPONENTIALS_FINAL_CUTOFF = 1e-5

927 928 -def llr_of_significance(sig):

929 return log((1.0-.0005)/sig)

930

931 -def ll_of_ssr(ssr, n):

932 s2 = ssr/n 933 return -(n/2)*(log(2*pi)+log(s2)) - ssr/(2*s2)

934

935 -class AllExponentials(qubx.fit.BaseCurve):

936 """This class implements a method to find all significant exponential components of a duration histogram, published in 937 938 Exponential Sum-Fitting of Dwell-Time Distributions without Specifying Starting Parameters 939 David Landowne,* Bin Yuan, and Karl L. Magleby* 940 Biophysical Journal Volume 104 June 2013 2383-2391 941 """

942 - def __init__(self, expr="0.0", locals={}):

943 qubx.fit.BaseCurve.__init__(self) 944 self.name = 'AllExponentials' 945 self.expr = "see Landowne et al, 2013" 946 self.params = ["N"] 947 for i in xrange(ALL_EXPONENTIALS_MAX_N): 948 self.params.extend(["t%i" % i, "a%i" % i]) 949 self.lo = [UNSET_VALUE]*len(self.params) 950 for i in xrange(ALL_EXPONENTIALS_MAX_N): 951 self.lo[2*i+1] = 1e-16 952 self.lo[2*i+2] = 0.0 953 self.hi = [UNSET_VALUE]*len(self.params) 954 self.can_fit = [False] * len(self.params) 955 self.param_defaults = [1.0] * len(self.params) 956 self.report_buf = StringIO() 957 self.llr = 0.0

958 - def report(self, *args):

959 self.report_buf.write('%s\n' % ' '.join(str(x) for x in args) )

960 - def eval(self, param_vals, xx, vvv=[]):

961 # area under amp * exp(t / -tau): 962 # amp * tau * (exp(-10**x{i-1} / -tau) - exp(-10**(r*x{i}) / tau)) 963 # where r = exp( log(10) * min(delta between two log10 x bins) ) 964 r = qubx.fit.calc_r(xx) 965 xp = numpy.power(10.0, xx) 966 xp /= r 967 yy_all = numpy.zeros(shape=xx.shape, dtype=xx.dtype) 968 for i in xrange(1, 1+2*max(0, min(ALL_EXPONENTIALS_MAX_N, int(round(param_vals[0])))), 2): 969 tau, area = param_vals[i:i+2] 970 if tau > 0.0: 971 amp = area / tau 972 xc = xp / (- tau) 973 yy = numpy.exp(xc) 974 xc *= r 975 yy -= numpy.exp(xc) 976 yy *= amp * tau 977 yy_all += yy 978 return yy_all

979 # 1. Log bin the data (done by DurHist panel) -- xx contains log10(left bound of each bin) 980 # 2. Generate the time constants for all possible exponentials

981 - def get_initial_params(self, xx, x0=None, x1=None):

982 """Adjusts can_fit, returns initial param values spaced evenly between min(xx) and max(xx) on log axis.""" 983 N = ALL_EXPONENTIALS_MAX_N 984 area = 1.0 / N 985 params = [area] * len(self.params) 986 params[0] = float(N) 987 if x0 is None: 988 X0, X1 = xx[0], xx[-1] 989 else: 990 X0, X1 = log10(x0), log10(x1) 991 t0 = X0 992 m = (X1 - X0) / (N - 1) 993 for i in xrange(N): 994 params[2*i+1] = 10**(t0 + i*m) 995 self.can_fit[2*i+1] = False 996 self.can_fit[2*i+2] = True 997 return params

998 # 3. Fit all areas with fixed time constants (done using a Fitter) 999 # 4. Delete the end exponentials with insignificant area, generate new evenly spaced taus, fit their areas

1000 - def get_params_minus_ends(self, params):

1001 """Adjusts can_fit, returns updated param values spaced evenly between smallest and largest legit expo in input params.""" 1002 N = ALL_EXPONENTIALS_MAX_N # int(round(params[0])) 1003 i0 = 0 1004 while (i0 < N) and (params[2*i0+2] < ALL_EXPONENTIALS_AREA_CUTOFF): 1005 self.report('4. Deleting initial insignificant component at',params[2*i0+1]) 1006 self.report(params) 1007 i0 += 1 1008 i1 = N-1 1009 while (i0 < i1) and (params[2*i1+2] < ALL_EXPONENTIALS_AREA_CUTOFF): 1010 self.report('4. Deleting final insignificant component at',params[2*i1+1]) 1011 self.report(params) 1012 i1 -= 1 1013 if (i0 != 0) or (i1 != (N-1)): 1014 return self.get_initial_params(None, params[2*i0+1], params[2*i1+1]) 1015 return params

1016 # 5. Delete all exponentials with negligible area

1017 - def get_params_minus_negligible(self, params, cutoff, stepno):

1018 """Adjusts can_fit, returns updated param values with the smallest expos removed.""" 1019 N = int(round(params[0])) 1020 params_out = [0] * len(params) 1021 i_out = 0 1022 for i_in in xrange(N): 1023 if params[2*i_in+2] >= cutoff: 1024 params_out[2*i_out+1], params_out[2*i_out+2] = params[2*i_in+1], params[2*i_in+2] 1025 i_out += 1 1026 else: 1027 self.report('%d. Deleting insignificant component at'%stepno,params[2*i_in+1]) 1028 self.report(params) 1029 params_out[0] = float(i_out) 1030 for i in xrange(i_out): 1031 self.can_fit[2*i+1] = True 1032 self.can_fit[2*i+2] = True 1033 for i in xrange(i_out, ALL_EXPONENTIALS_MAX_N): 1034 self.can_fit[2*i+1] = False 1035 self.can_fit[2*i+2] = False 1036 return params_out

1037 # 6. Fit remaining areas and taus (using a Fitter) 1038 # 7. Combine exponentials with very close time constants (in order of decreasing time const, re-fitting after each)

1039 - def get_params_combined(self, params, fitter, xx, yy, vvv, ww, ssr):

1040 """Adjusts can_fit, returns updated, fitted param values with near-identical time constants combined into one expo, and final ssr""" 1041 N = [int(round(params[0]))] 1042 params_out = params[:] 1043 # make sure they're in order of increasing tau 1044 expos = [(params[2*i+1], params[2*i+2]) for i in xrange(N[0])] 1045 expos.sort() 1046 for tau, area in expos: 1047 params_out[2*i+1], params_out[2*i+2] = tau, area 1048 1049 def combine(i): 1050 t1, a1, t2, a2 = params_out[2*i+1:2*i+5] 1051 params_out[2*i+1] = t1*a1/(a1+a2) + t2*a2/(a1+a2) 1052 params_out[2*i+2] = a1+a2 1053 for j in xrange(i+1,N[0]-1): 1054 params_out[2*j+1:2*j+3] = params_out[2*j+3:2*j+5] 1055 N[0] -= 1 1056 params_out[0] = float(N[0]) 1057 self.can_fit[2*N[0]+1] = self.can_fit[2*N[0]+2] = False 1058 params_out[2*N[0]+1] = params_out[2*N[0]+2] = 0.0

1059 for i in reversed(xrange(N[0]-1)): 1060 t0 = params_out[2*i+1] 1061 if abs((params_out[2*i+3] - t0) / t0) < ALL_EXPONENTIALS_NEAR_THRESHOLD: 1062 self.report('7. Combining close components at',t0,'and',params_out[2*i+3]) 1063 self.report(params) 1064 combine(i) 1065 params_out[:], ssr, iterations, ff = fitter(self, params_out, xx, yy, vvv, ww) 1066 return params_out, ssr

1067 # 8. Determine the best fit with one fewer exponential

1068 - def get_fit_with_fewer(self, params, fitter, xx, yy, vvv, ww, map_func):

1069 """One at a time, removes each exponential and re-fits; adjusts can_fit, returns reduced params, ssr of trial with max likelihood.""" 1070 N = int(round(params[0])) 1071 ssr_params = map_func(all_expo_fit_without, [(xx, yy, vvv, ww, params, self.can_fit, i) for i in xrange(N)]) 1072 ssr_without = [ssr for ssr, params_w in ssr_params] 1073 best_i = numpy.argmin(ssr_without) 1074 self.report('8. Removing component at',params[2*best_i+1]) 1075 self.report(params) 1076 params_out = ssr_params[best_i][1] 1077 # then we consolidate the output parameters 1078 N -= 1 1079 for j in xrange(best_i, N): 1080 params_out[2*j+1:2*j+3] = params_out[2*j+3:2*j+5] 1081 params_out[0] = float(N) 1082 params_out[2*N+1] = params_out[2*N+2] = 0.0 1083 self.can_fit[2*N+1] = self.can_fit[2*N+2] = False 1084 return params_out, ssr_without[best_i]

1085 # 9. Repeat step 8 until it makes the likelihood of the fit significantly worse

1086 - def try_fit_with_fewer(self, params, fitter, xx, yy, vvv, ww, ll0, map_func):

1087 params_fewer, ssr_fewer = self.get_fit_with_fewer(params, fitter, xx, yy, vvv, ww, map_func) 1088 ll = ll_of_ssr(ssr_fewer, len(xx)) 1089 llr = ll0 - ll 1090 return llr, params_fewer, ssr_fewer

1091 # the following method was incorporated into run_generator() for responsiveness (more yielding)

1092 - def get_minimal_fit(self, params, fitter, xx, yy, vvv, ww, ssr0, significance):

1093 """Adjusts can_fit, returns fully reduced params, ssr.""" 1094 target_llr = llr_of_significance(significance or ALL_EXPONENTIALS_SIGNIFICANCE) 1095 ll0 = ll_of_ssr(ssr0, len(xx)) 1096 if params[0] > 1: 1097 params_fewer, ssr = self.get_fit_with_fewer(params, fitter, xx, yy, vvv, ww) 1098 ll = ll_of_ssr(ssr, len(xx)) 1099 llr = ll0 - ll 1100 self.llr = llr 1101 if llr <= target_llr: 1102 return self.get_minimal_fit(params_fewer, fitter, xx, yy, vvv, ww, ssr0, significance) 1103 N = int(round(params[0])) 1104 self.can_fit[2*(N-1)+1] = self.can_fit[2*(N-1)+2] = True # re-enable params from failed get_fit_with_fewer 1105 return params, ssr0

1106 # 10. Exclude brief erroneous exponentials -- with tau less than the dead time, 1107 # and whose fitted area contributes < 1e-5 of the total fitted area 1108 # area * (exp(-10**xx[0] / -tau) - exp(-10**xx[-1] / tau))

1109 - def get_params_minus_brief_errs(self, params, xx):

1110 """Returns params, stripped of brief erroneous components.""" 1111 N = int(round(params[0])) 1112 params_out = params[:] 1113 left, right = 10**xx[0], 10**xx[-1] 1114 def fitted_area(i): 1115 tau, area = params[2*i+1], params[2*i+2] 1116 return area * (exp(-left / tau) - exp(-right / tau))

1117 i_out = 0 1118 for i in xrange(N): 1119 if (params[2*i+1] < left) and (fitted_area(i) < ALL_EXPONENTIALS_AREA_CUTOFF): 1120 self.report('10. Deleting brief erroneous component at',params[2*i+1]) 1121 self.report(params) 1122 N -= 1 1123 self.can_fit[2*N+1] = self.can_fit[2*N+2] = False 1124 params_out[2*N+1] = params_out[2*N+2] = 0.0 1125 else: 1126 params_out[2*i_out+1:2*i_out+3] = params[2*i+1:2*i+3] 1127 i_out += 1 1128 params_out[0] = float(i_out) 1129 return params_out

1130 - def run_generator(self, xx, yy, params, ww=None, fitter=None, map_func=map, significance=None):

1131 """Modifies params; yields percent complete after each step.""" 1132 self.report_buf = StringIO() 1133 self.llr = 0.0 1134 f = qubx.fit.LMFitter() if (fitter is None) else fitter 1135 def fit(*args): 1136 for i in xrange(ALL_EXPONENTIALS_FIT_REPS): 1137 params, ssr, iterations, ff = f(*args) 1138 if iterations < ALL_EXPONENTIALS_FIT_ITERS: 1139 break 1140 return params, ssr, iterations, ff

1141 vvv = [] 1142 params[:] = self.get_initial_params(xx) 1143 yield 20.0 1144 DurHist = qubx.global_namespace.QubX.Figures.DurHist 1145 if DurHist.cancel_all_exp: return 1146 params[:], ssr, iterations, ff = fit(self, params, xx, yy, vvv, ww) 1147 yield 30.0 1148 if DurHist.cancel_all_exp: return 1149 params[:] = self.get_params_minus_ends(params) 1150 params[:], ssr, iterations, ff = fit(self, params, xx, yy, vvv, ww) 1151 yield 40.0 1152 if DurHist.cancel_all_exp: return 1153 params[:] = self.get_params_minus_negligible(params, ALL_EXPONENTIALS_AREA_CUTOFF, 5) 1154 yield 50.0 1155 if DurHist.cancel_all_exp: return 1156 params[:], ssr, iterations, ff = fit(self, params, xx, yy, vvv, ww) 1157 yield 60.0 1158 if DurHist.cancel_all_exp: return 1159 params[:], ssr = self.get_params_combined(params, fit, xx, yy, vvv, ww, ssr) 1160 yield 70.0 1161 if DurHist.cancel_all_exp: return 1162 1163 # 9. Repeat step 8 until it makes the likelihood of the fit significantly worse 1164 target_llr = llr_of_significance(significance or ALL_EXPONENTIALS_SIGNIFICANCE) 1165 ll0 = ll_of_ssr(ssr, len(xx)) 1166 N = int(round(params[0])) 1167 while N > 1: 1168 llr, params_fewer, ssr_fewer = self.try_fit_with_fewer(params, fit, xx, yy, vvv, ww, ll0, map_func) 1169 if llr > target_llr: 1170 self.report('...failed, llr =',llr,', keeping component') 1171 break 1172 self.report('...success, llr =',llr) 1173 self.llr = llr 1174 params[:], ssr = params_fewer, ssr_fewer 1175 N = int(round(params[0])) 1176 yield 90 - N 1177 if DurHist.cancel_all_exp: return 1178 self.can_fit[2*(N-1)+1] = self.can_fit[2*(N-1)+2] = True # re-enable params from failed get_fit_with_fewer 1179 1180 yield 90.0 1181 if DurHist.cancel_all_exp: return 1182 params[:] = self.get_params_minus_brief_errs(params, xx) 1183 params[:], ssr, iterations, ff = fit(self, params, xx, yy, vvv, ww) 1184 yield 99.0 1185 if ALL_EXPONENTIALS_FINAL_CUTOFF: # "step 11" 1186 params[:] = self.get_params_minus_negligible(params, ALL_EXPONENTIALS_FINAL_CUTOFF, 11) 1187 params[:], ssr, iterations, ff = fit(self, params, xx, yy, vvv, ww) 1188 if DurHist.cancel_all_exp: return 1189 yield 100.0

1190 - def run(self, xx, yy, ww=None, fitter=None, map_func=map, significance=None):

1191 """Returns fitted list of (tau, area) of all significant exponentials.""" 1192 params = [1.0] * len(self.params) 1193 for pct in self.run_generator(xx, yy, params, ww, fitter, map_func, significance): 1194 pass 1195 result = [] 1196 for i in xrange(int(round(params[0]))): 1197 result.append(tuple(params[2*i+1:2*i+3])) 1198 return result

1199

1200 -class AllExponentials_Display(qubx.fit.BaseCurve):

1201 - def __init__(self, expr='1 (sqrt) Exponential, to use: click "Fit all exponentials" in the DurHist panel', locals={}):

1202 qubx.fit.BaseCurve.__init__(self) 1203 self.name = 'AllExponentials' 1204 self.params = [] 1205 self.lo = [] 1206 self.hi = [] 1207 self.can_fit = [] 1208 self.param_defaults = [] 1209 self.can_resample = False 1210 self.last_fit = [] 1211 1212 self.expr = expr # actually smuggling in N and sqrt 1213 self.N = int(re.search(r"([0-9]+) ", expr).group(1)) 1214 self.sqrt = 'sqrt' in expr 1215 Ninv = 1.0 / self.N 1216 for i in xrange(self.N): 1217 self.params.extend(["tau%i"%(i+1), "area%i"%(i+1)]) 1218 self.lo.extend([1e-16, 0.0]) 1219 self.hi.extend([UNSET_VALUE, UNSET_VALUE]) 1220 self.can_fit.extend([True, True]) 1221 self.param_defaults.extend([(i+1)*Ninv, Ninv])

1222 - def eval(self, param_vals, xx, vvv=[]):

1223 # area under amp * exp(t / -tau): 1224 # amp * tau * (exp(-10**x / -tau) - exp(-10**(r*x) / tau)) 1225 # where r = exp( log(10) * min(delta between two log10 x bins) ) 1226 r = qubx.fit.calc_r(xx) 1227 xp = numpy.power(10.0, xx) 1228 xp /= r 1229 yy_all = numpy.zeros(shape=xx.shape, dtype=xx.dtype) 1230 for i in xrange(self.N): 1231 tau, area = param_vals[2*i:2*i+2] 1232 if tau > 0.0: 1233 amp = area / tau 1234 xc = xp / (- tau) 1235 yy = numpy.exp(xc) 1236 xc *= r 1237 yy -= numpy.exp(xc) 1238 yy *= amp * tau 1239 yy_all += yy 1240 if self.sqrt: 1241 self.last_fit = numpy.sqrt(yy_all) 1242 else: 1243 self.last_fit = yy_all 1244 return self.last_fit

1245 qubx.fit.Curves['AllExponentials'] = AllExponentials_Display

1246 1247 -def all_expo_fit_without(args):

1248 try: 1249 # for more efficient test, we set area to 0.0 and can_fit to False 1250 xx, yy, vvv, ww, params, can_fit, i = args 1251 fitter = qubx.fit.LMFitter() # not generally possible to serialize a Fitter with multiprocessing, so ignoring caller's choice of fitter here 1252 p = params[:] 1253 p[2*i+2] = 0.0 1254 c = AllExponentials() 1255 c.can_fit[:] = can_fit[:] 1256 c.can_fit[2*i+1] = c.can_fit[2*i+2] = False 1257 fit_params, ssr, iterations, ff = fitter(c, p, xx, yy, vvv, ww) 1258 return ssr, fit_params 1259 except: 1260 traceback.print_exc() 1261 return 0.0

1262

Source Code for Module qubx_single.dur_hist