Source Code for Module qubx_single.rates

  1  """Provides a visual panel (L{qubx.faces.Face}) to optimize single-molecule rate constants. 
  2   
  3  Uses math from U{qubx_single_molecule<http://www.qub.buffalo.edu/src/qub-express/api/qubx_single_molecule>}  
  4   
  5  /* Copyright 2008-2014 Research Foundation State University of New York   */ 
  6   
  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 itertools 
 26  import gc 
 27  import os 
 28  import sys 
 29  import time 
 30  import traceback 
 31   
 32  from qubx.data_types import * 
 33  from qubx.util_types import * 
 34  from qubx.faces import * 
 35  from qubx.GTK import * 
 36  import qubx.notebook 
 37  import qubx.optimize 
 38  import qubx.optimizeGTK 
 39  import qubx.pyenv 
 40  from qubx.settings import Propertied, Property 
 41  import qubx.tree 
 42  from qubx.trialset import * 
 43  from qubx.model import * 
 44  import qubx.fast.model 
 45  from qubx.fast.fast_utils import buffer_as_array_of_int, buffer_as_array_of_double, buffer_as_matrix_of_int, buffer_as_matrix_of_double 
 46  from qubx.task import * 
 47  from max_subi_ll import * 
 48  from max_inter_ll import inter_ll_or_0 
 49   
 50  import scipy.optimize 
 51  from numpy import * 
 52   
 53   
 54  MAXITER = 300 
 55  LBL_TDEAD = 'Tdead (ms)' 
 56   
 57  OPT_SIMPLEX = 0 
 58  OPT_DFP     = 1 
 59   
 60 -def cl_dimension(Ns): 
 61      dim = 4 
 62      while dim < Ns: 
 63          dim *= 2 
 64      return dim 
 65   
 66 -class ModelPrep(Anon): 
 67      """Represents a model ready for computation.""" 
 68 -    def __init__(self, model=None, modeltree=None): # one or the other, not both; modeltree is more efficient 
 69          Anon.__init__(self, source=model, modeltree=modeltree or model.as_tree(), model=QubModel()) 
 70          self.model.from_tree(self.modeltree) 
 71          self.K0, self.K1, self.K2, L, V, P = ModelToRateMatricesP(self.model) 
 72 -    def set_rates2(self, K0, K1, K2): 
 73          self.K0, self.K1, self.K2 = K0, K1, K2 
 74          K012ToModel(K0, K1, K2, self.model) 
 75   
 76 -class DataPrep(Anon): 
 77      """Represents data ready for computation.""" 
 78 -    def __init__(self, segs, ws): 
 79          Anon.__init__(self, segs=segs, ws=ws) 
 80           
 81 -class Workspace(Anon): 
 82      """Holds data, model, and results for one file.""" 
 83      pass 
 84   
 85 -def get_data_prep(segs, model_prep, wait=True, receiver=lambda dp: None): 
 86      if wait: 
 87          return qubx.pyenv.call_async_wait(lambda rcv: get_data_prep(segs, model_prep, wait=False, receiver=rcv)) 
 88      else: 
 89          QubX = qubx.global_namespace.QubX 
 90          view_of_file = dict([(view.file, view) for view in QubX.Data.views]) 
 91          view_segm = [] 
 92          vs_by_file = {} 
 93          for seg in segs: 
 94              if not (seg.file in vs_by_file): 
 95                  vs = (view_of_file[seg.file], []) 
 96                  view_segm.append(vs) 
 97                  vs_by_file[seg.file] = vs 
 98              vs_by_file[seg.file][1].append(seg) 
 99          if not view_segm: 
100              gobject.idle_add(receiver, DataPrep(segs, [])) 
101              return 
102          workspaces = [] 
103          data_prep = DataPrep(segs, workspaces) 
104          # stimulus idl could take a while, so its results come back on gobject.idle_add. 
105          # we'll use idle_add to alternate recursively between build_workspace and receive_stim, 
106          # until all workspaces are built. 
107          def build_workspace(): 
108              i = len(workspaces) 
109              view, segm = view_segm[i] 
110              ws = Workspace(view=view) 
111              workspaces.append(ws) 
112              if view: 
113                  for i in xrange(view.file.constants.size): 
114                      if LBL_TDEAD == view.file.constants.get(i, 'Name'): 
115                          ws.td_ms = view.file.constants.get(i, 'Value') 
116                          break 
117                  else: 
118                      ws.td_ms = view.file.sampling * 1e3 
119                      view.file.constants.append({'Name':LBL_TDEAD, 'Value':ws.td_ms}) 
120                      view.file.constants.select(view.file.constants.size-1, 'Value') 
121               
122                  ws.sampling_ms = view.file.sampling * 1e3 
123                  ws.sig_names = [view.file.signals.get(segm[0].signal, 'Name')] 
124                  ws.signals = [ [] ] 
125                  for seg in segm: 
126                      ff, ll, cc, dd = QubX.DataSource.get_idealization(seg, mark_excluded=True, get_fragments=True, get_durations=True) 
127                      aa = view.file.ideal[segm[0].signal].seg[seg.index].amp 
128                      ws.signals[0].append( (cc, dd, aa) ) 
129                  QubX.Modeling.Stimulus.call_with_stimulus_idl(receive_stim, [(seg.f, seg.l) for seg in segm], datafile=view.file, model=model_prep.source or model_prep.model) 
130              else: 
131                  receive_stim([], []) 
132          def receive_stim(stim_names, stimidl): 
133              i = len(workspaces) - 1 
134              ws = workspaces[i] 
135              ws.sig_names.extend(stim_names) 
136              for segs_flcda in stimidl: 
137                  ws.signals.append( [(cc, dd, aa) for ff,ll,cc,dd,aa in segs_flcda] ) 
138              try: 
139                  ws.segs, ws.mcls, ws.scls = ProcessSignals(ws.td_ms, ws.signals) 
140              except: 
141                  traceback.print_exc() 
142                  ws.segs, ws.mcls, ws.scls = [], [], [] 
143              del ws.signals 
144              if len(workspaces) == len(view_segm): 
145                  gobject.idle_add(receiver, data_prep) # exit from recursion, return workspaces via idle_add 
146              else: 
147                  gobject.idle_add(build_workspace) 
148          build_workspace() 
149   
150 -class MultiplexedSignals(Requestable): 
151 -    def __init__(self, datasrc, Data, stimulus): 
152          Requestable.__init__(self) 
153          self.datasrc = datasrc 
154          self.Data = Data 
155          self.stimulus = stimulus 
156          self.__ref = Reffer() 
157          self.__use_group = None 
158          self.__data = None 
159          self.OnChangeTdead = WeakEvent() 
160          self.datasrc.OnChangeIdealization += self.__ref(self.__onChangeIdl) 
161          self.Data.table.OnInsert += self.__ref(self.__onChangeFileList) 
162          self.Data.table.OnRemoved += self.__ref(self.__onChangeFileList) 
163          self.Data.table.OnSet += self.__ref(self.__onSetFileList) 
164          self.Data.OnSwitch += self.__ref(self.__onSwitchDataFile) 
165          self.__datafile = None 
166          self.__onSwitchDataFile(self.Data, self.Data.file) 
167          self.stimulus.OnInsert += self.__ref(self.__onChangeStimulus) 
168          self.stimulus.OnChange += self.__ref(self.__onChangeStimulus) 
169          self.stimulus.OnRemoving += self.__ref(self.__onChangeStimulus) 
170          self.invalidate() 
171      tdead = property(lambda self: self.get_tdead(), lambda self, x: self.set_tdead(x)) 
172 -    def get_tdead(self): 
173          if self.__datafile: 
174              for i in xrange(self.__datafile.constants.size): 
175                  if LBL_TDEAD == self.__datafile.constants.get(i, 'Name'): 
176                      return self.__datafile.constants.get(i, 'Value') 
177              td_ms = self.__datafile.sampling * 1e3 
178              self.__datafile.constants.append({'Name':LBL_TDEAD, 'Value':td_ms}) 
179              return td_ms 
180 -    def set_tdead(self, x, scriptable=False): 
181          if self.__datafile: 
182              qubx.pyenv.env.OnScriptable('QubX.Data.file.constants[QubX.Data.file.constants.index(%s), "Value"] = %s' % (repr(LBL_TDEAD), repr(x))) 
183              cc = self.__datafile.constants 
184              try: 
185                  cc[cc.index(LBL_TDEAD), 'Value'] = x 
186              except KeyError: 
187                  print 'surprise!' 
188                  cc.append({'Name':LBL_TDEAD, 'Value':x}) 
189 -    def set_use_group(self, x): 
190          if self.__use_group != x: 
191              self.__use_group = x 
192              self.invalidate() 
193      use_group = property(lambda self: self.__use_group, lambda self, x: self.set_use_group(x)) 
194 -    def __onChangeIdl(self, *args): 
195          self.invalidate() # TODO: on scroll, if use_group != None, no invalidate 
196 -    def __onChangeFileList(self, *args): 
197          if self.__use_group != None: 
198              self.invalidate() 
199 -    def __onSetFileList(self, i, field_name, val, prev, undoing): 
200          if (self.__use_group != None) and (field_name == 'Group'): 
201              self.invalidate() 
202 -    def __onChangeStimulus(self, *args): 
203          self.invalidate() 
204 -    def __onSwitchDataFile(self, Data, file): 
205          if file == self.__datafile: return 
206          if self.__datafile: 
207              self.__datafile.constants.OnSet -= self.__ref(self.__onSetConstant) 
208          self.__datafile = file 
209          if self.__datafile: 
210              self.__datafile.constants.OnSet += self.__ref(self.__onSetConstant) 
211          self.OnChangeTdead() 
212 -    def __onSetConstant(self, r, field, val, prev, undoing): 
213          if self.__datafile.constants[r,'Name'] == LBL_TDEAD: 
214              self.invalidate() 
215              self.OnChangeTdead() 
216 -    def get(self, receiver): 
217          if self.__use_group != None: 
218              view_segm = self.datasrc.get_segmentation_files(group=self.__use_group) 
219          else: 
220              view_segm = [(self.Data.view, self.datasrc.get_segmentation())] 
221          segs = [] 
222          for view, segm in view_segm: 
223              segs.extend(segm) 
224          get_data_prep(segs, self.model_prep, wait=False, receiver=receiver) 
225   
226   
227  QUBX_SINGLE_RATES_METHODS = ["MIL", "MSL"] 
228   
229  @Propertied(Property('use_Peq', True, "True to start each segment at equilibrium; False to use States:Pr"), 
230              Property('use_group', False, "True to fit a group of data files", in_tree=False), 
231              Property('use_group_index', 1, 'Desired Group in Data table'), 
232              Property('useGPU', True, "True to use OpenCL when available"), 
233              Property("method", QUBX_SINGLE_RATES_METHODS[0], "choice of algorithm in %s" % ', '.join(QUBX_SINGLE_RATES_METHODS)), 
234              Property("tdead_in_samples", True, "False to display Tdead in milliseconds") 
235              ) 
236 -class RatesFace(Face): 
237 -    def __init__(self): 
238          Face.__init__(self, 'MILRates', 'QubX.Modeling.MILRates') 
239          self.__ref = Reffer() 
240          self.propertied_connect_settings('Rates') 
241   
242          self.set_size_request(200, 50) 
243          self.QubX = qubx.pyenv.env.globals['QubX'] 
244          self.QubX.OnQuit += self.__ref(self.qubx_finalize) 
245          self.models = self.QubX.Models 
246          self.models.OnSwitch += self.__ref(self.__onSwitchModel), 'Rates.onSwitchModel' 
247          self.robot = Robot('MILRates') 
248          self.robot.OnException += self.__ref(self.__onException) 
249          self.robot.do(self.__robot_setup_cl_context) 
250          self.plexig = MultiplexedSignals(self.QubX.DataSource, self.QubX.Data, self.QubX.Modeling.Stimulus) 
251          self.plexig.use_group = self.use_group_index if self.use_group else None 
252          self.plexig.OnInvalidate += self.__ref(self.__onInvalidateData) 
253          self.plexig.OnChangeTdead += self.__ref(self.__onChangeTdead) 
254          self.model_prep = Product(self.__get_model_prep) 
255          self.model_prep.OnInvalidate += self.__ref(self.__onInvalidateModel) 
256          self.rates = Product(self.__get_rates) 
257          self.rates.OnInvalidate += self.__ref(self.__onInvalidateModel) 
258          self.LL = 0.0 
259          self.iterations = 0 
260          self.grads = [] 
261          self.__on_trial = None 
262   
263          qubx.optimize.setup_opt("qubx_single_rates_optimize", 
264                                  {"max step" : 0.1}) 
265           
266          self.OnLL = WeakEvent() 
267          self.OnIteration = WeakEvent() 
268          self.OnOptimized = WeakEvent() 
269   
270          qubx.notebook.Notebook.register_auto('Rates_eval.Messages', 'Messages, on MILRates: Get LL', False) 
271          qubx.notebook.Notebook.register_auto('Rates_eval.Model.Picture', 'Model picture, on MILRates: Get LL', False) 
272          qubx.notebook.Notebook.register_auto('Rates_eval.Model.PictureNO', 'Model picture, no overlays, on MILRates: Get LL', False) 
273          qubx.notebook.Notebook.register_auto('Rates_eval.Data.Picture', 'Data picture, on MILRates: Get LL', False) 
274          qubx.notebook.Notebook.register_auto('Rates_eval.Data.PictureNO', 'Data picture, no overlays, on MILRates: Get LL', False) 
275          qubx.notebook.Notebook.register_auto('Rates_eval.DurHist.Charts', 'DurHist charts, on MILRates: Get LL', False) 
276          qubx.notebook.Notebook.register_auto('Rates_eval.DurHist.TimeConst', 'DurHist time constants, on MILRates: Get LL', False) 
277          qubx.notebook.Notebook.register_auto('Rates_eval.Model.Rates', 'Rates table, on MILRates: Optimize', False) 
278   
279          qubx.notebook.Notebook.register_auto('Rates_opt.Messages', 'Messages, on MILRates: Optimize', True) 
280          qubx.notebook.Notebook.register_auto('Rates_opt.Model.Picture', 'Model picture, on MILRates: Optimize', True) 
281          qubx.notebook.Notebook.register_auto('Rates_opt.Model.PictureNO', 'Model picture, no overlays, on MILRates: Optimize', False) 
282          qubx.notebook.Notebook.register_auto('Rates_opt.Data.Picture', 'Data picture, on MILRates: Optimize', True) 
283          qubx.notebook.Notebook.register_auto('Rates_opt.Data.PictureNO', 'Data picture, no overlays, on MILRates: Optimize', False) 
284          qubx.notebook.Notebook.register_auto('Rates_opt.DurHist.Charts', 'DurHist charts, on MILRates: Optimize', True) 
285          qubx.notebook.Notebook.register_auto('Rates_opt.DurHist.TimeConst', 'DurHist time constants, on MILRates: Optimize', True) 
286          qubx.notebook.Notebook.register_auto('Rates_opt.Model.Rates', 'Rates table, on MILRates: Optimize', False) 
287          qubx.notebook.Notebook.register_auto('Rates_opt.Model.Constraints', 'Kinetic Constraints table, on MILRates: Optimize', True) 
288   
289          cols = pack_item(gtk.HBox(), self, expand=True) 
290          ctrls = pack_item(gtk.VBox(), cols, expand=True) 
291          row = pack_item(gtk.HBox(), ctrls) 
292          self.chkPeq = pack_check('start at equilibrium prob.   ', row) 
293          self.propertied_connect_check('use_Peq', self.chkPeq) 
294          pack_space(row, expand=True) 
295          self.chkUseSingle = pack_radio('Single file', row) 
296          self.chkUseGroup = pack_radio('Multi-file: Group', row, group=self.chkUseSingle) 
297          self.chkUseGroup.set_tooltip_text('Assign Groups in the Data table') 
298          self.propertied_connect_radios('use_group', [(False, self.chkUseSingle), 
299                                                       (True, self.chkUseGroup)]) 
300          self.txtUseGroup = pack_item(qubx.GTK.NumEntry(self.use_group_index, acceptIntGreaterThan(-1), width_chars=3), row) 
301          self.propertied_connect_NumEntry('use_group_index', self.txtUseGroup) 
302          pack_space(row, expand=True) 
303   
304          pack_label('', row).set_markup('t<sub>dead</sub>:') 
305          td = (1e-3*self.plexig.tdead/self.QubX.Data.file.sampling) if self.tdead_in_samples else self.plexig.tdead 
306          self.txtTdead = pack_item(qubx.GTK.NumEntry(td, acceptFloatGreaterThanOrEqualTo(0.0), '%.2g', width_chars=5), row) 
307          self.txtTdead.set_tooltip_text("""Note for multi-file: this is Tdead for the onscreen Data file only. 
308  To edit another file's Tdead, first select it in the Data panel.""") 
309          self.txtTdead.OnChange += self.__ref(self.__onEditTdead) 
310          self.chkTdeadMS = pack_radio('[ms]', row) 
311          self.chkTdeadSamples = pack_radio('[samples]', row, group=self.chkTdeadMS) 
312          self.propertied_connect_radios('tdead_in_samples', [(False, self.chkTdeadMS), (True, self.chkTdeadSamples)]) 
313          pack_space(row, expand=True) 
314   
315          self.txtStatus = gtk.TextView() 
316          self.outStatus = TextViewAppender(self.txtStatus) 
317          self.txtStatus.set_editable(False) 
318          SetFixedWidth(self.txtStatus) 
319          pack_scrolled(self.txtStatus, ctrls, expand=True) 
320          self.outStatus.write("""Based on "Estimating single-channel kinetic parameters from idealized patch-clamp data containing missed events," Qin F, Auerbach A, Sachs F, Biophys J. 1996 Jan; 70(1):264-80 
321  and "Maximum likelihood estimation of aggregated Markov processes," Qin F, Auerbach A, Sachs F, Proc Biol Sci. 1997 Mar 22; 264(1380):375-83. 
322   
323  """) 
324   
325          ctrls = pack_item(gtk.VBox(), cols) 
326          row = pack_item(gtk.HBox(), ctrls) 
327          pack_label('LL:', row) 
328          self.txtLL = pack_item(gtk.Entry(), row, expand=True) 
329          self.txtLL.set_editable(False) 
330          self.txtLL.set_width_chars(8) 
331          self.chkUseGPU = pack_check("use GPU", ctrls) 
332          self.chkUseGPU.set_tooltip_text("If an OpenCL-capable graphics processor is available") 
333          self.propertied_connect_check('useGPU', self.chkUseGPU) 
334          self.mnuMethod = pack_item(gtk.combo_box_new_text(), ctrls) 
335          for i, method_name in enumerate(QUBX_SINGLE_RATES_METHODS): 
336              self.mnuMethod.append_text(method_name) 
337              if method_name == self.method: 
338                  self.mnuMethod.set_active(i) 
339          self.mnuMethod.connect('changed', self.__onChangeMethod) 
340          self.btnGetLL = pack_button('Get LL', ctrls, on_click=self.__onClickGetLL) 
341          self.btnAddToTrials = pack_button('Add to Trials', ctrls, on_click=self.__onClickAddToTrials) 
342          self.btnOpt = pack_button('Optimize Rates', ctrls, on_click=self.__onClickOpt, at_end=True) 
343          self.optOptions = qubx.optimizeGTK.OptOptionsFace('MILRates', 'qubx_single_rates_optimize', 'QubX.About.MILRates') 
344           
345          self.__data = self.__model = None 
346          self.__optimizing = False 
347          self.__onSwitchModel(self.models, self.models.file) 
348          self.cl_context = None 
349           
350          self.QubX.Modeling.Utils.add_optimizer('MIL', ModelPrep, get_data_prep, self.get_ll_mil, self.optimize_mil, self.OnLL, self.OnIteration, self.OnOptimized) 
351          self.QubX.Modeling.Utils.add_optimizer('MSL', ModelPrep, get_data_prep, self.get_ll_msl, self.optimize_msl, self.OnLL, self.OnIteration, self.OnOptimized) 
352      optimizing = property(lambda self: self.__optimizing) 
353 -    def propertied_set(self, value, name): 
354          super(RatesFace, self).propertied_set(value, name) 
355          self.__gray_ll() 
356          if name in ('use_group', 'use_group_index'): 
357              self.plexig.use_group = self.use_group_index if self.use_group else None 
358          elif name == 'method': 
359              i = QUBX_SINGLE_RATES_METHODS.index(value) 
360              self.mnuMethod.set_active(i) 
361          elif name == 'tdead_in_samples': 
362              self.__onChangeTdead() 
363 -    def onShow(self, showing): 
364          if showing: 
365              self.QubX.About.append_face(self.optOptions) 
366              self.QubX.About.show_face('MILRates') 
367          else: 
368              self.QubX.About.remove_face(self.QubX.About.index('MILRates')) 
369 -    def qubx_finalize(self): 
370          self.robot.stop() 
371          self.data = None 
372          self.model = None 
373          #self.QubX.OnQuit -= self.__ref(self.qubx_finalize) 
374          #self.models.OnSwitch -= self.__ref(self.__onSwitchModel) 
375 -    def __onException(self, robot, typ, val, trace): 
376          traceback.print_exception(typ, val, trace, file=self.outStatus) 
377 -    def __onChangeTdead(self): 
378          td = self.plexig.tdead 
379          if self.tdead_in_samples: 
380              self.txtTdead.value = 1e-3*td / self.QubX.Data.file.sampling 
381          else: 
382              self.txtTdead.value = td 
383 -    def __onEditTdead(self, txt, val): 
384          if self.tdead_in_samples: 
385              self.plexig.set_tdead(val * 1e3 * self.QubX.Data.file.sampling, True) 
386          else: 
387              self.plexig.set_tdead(val, True) 
388 -    def __onSwitchModel(self, models, model): 
389          if self.__optimizing: 
390              self.robot.do(gobject.idle_add, self.set_model, model) 
391          else: 
392              self.model = model 
393 -    def set_model(self, model): 
394          if model == self.__model: return 
395          if self.__model: 
396              for event in (self.__model.states.OnSet, 
397                            self.__model.states.OnInsert, 
398                            self.__model.states.OnRemoved, 
399                            self.__model.rates.OnInsert, 
400                            self.__model.rates.OnRemoved, 
401                            self.__model.constraints_kin.OnEdit, 
402                            self.__model.OnChangeBalanceLoops, 
403                            self.__model.OnSetChannelCount): 
404                  event -= self.__ref(self.__onSetModel) 
405              self.__model.rates.OnSet -= self.__ref(self.__onSetRate) 
406          self.__model = model 
407          if self.__model: 
408              for event in (self.__model.states.OnSet, 
409                            self.__model.states.OnInsert, 
410                            self.__model.states.OnRemoved, 
411                            self.__model.rates.OnInsert, 
412                            self.__model.rates.OnRemoved, 
413                            self.__model.constraints_kin.OnEdit, 
414                            self.__model.OnChangeBalanceLoops, 
415                            self.__model.OnSetChannelCount): 
416                  event += self.__ref(self.__onSetModel) 
417              self.__model.rates.OnSet += self.__ref(self.__onSetRate) 
418              self.__onSetModel() 
419      model = property(lambda self: self.__model, lambda self, x: self.set_model(x)) 
420 -    def __onInvalidateData(self, product): 
421          self.model_prep.invalidate() 
422          self.rates.invalidate() 
423          self.robot.do(gobject.idle_add, self.__gray_ll) 
424 -    def __onInvalidateModel(self, product): 
425          if not self.__optimizing: 
426              self.robot.do(gobject.idle_add, self.__gray_ll) 
427 -    def __onSetModel(self, *args): 
428          self.model_prep.invalidate() 
429          self.rates.invalidate() 
430 -    def __onSetRate(self, r, field, val, prev, undoing): 
431          self.rates.invalidate() 
432 -    def __onChangeMethod(self, mnu): 
433          self.propertied_on_user_set('method', QUBX_SINGLE_RATES_METHODS[mnu.get_active()]) 
434 -    def __gray_ll(self, gray=True): 
435          color = gdk.color_parse(gray and "#999999" or "#000000") 
436          self.txtLL.modify_text(gtk.STATE_NORMAL, color) 
437 -    def __begin_task(self, receiver=None): 
438          self.__receiver = receiver 
439          self.txtStatus.get_buffer().set_text('') 
440          Tasks.add_task(self.robot) 
441          for ctrl in [self.btnGetLL, self.btnAddToTrials, self.btnOpt, self.chkPeq]: 
442              ctrl.set_sensitive(False) 
443 -    def __end_task(self, result): 
444          try: 
445              self.LL = result[0] 
446          except: 
447              self.LL = result 
448          Tasks.remove_task(self.robot) 
449          for ctrl in [self.btnGetLL, self.btnAddToTrials, self.btnOpt, self.chkPeq]: 
450              ctrl.set_sensitive(True) 
451          if self.method == 'MIL': 
452              self.request_show() 
453          self.__optimizing = False 
454          if self.__receiver: 
455              self.__receiver(result) 
456              del self.__receiver 
457 -    def __robot_setup_cl_context(self): 
458          self.cl_context = msl_accel_context(MSL_ACCEL_OPENCL) 
459 -    def __get_model_prep(self): 
460          if self.__model == None: 
461              return None 
462          return ModelPrep(self.__model) 
463 -    def __get_rates(self): 
464          if self.__model == None: 
465              return None, None 
466          signal_names = self.QubX.Modeling.Stimulus.active_names 
467          K0, K1, K2, L, V, P = ModelToRateMatricesP(self.__model, signal_names) 
468          return K0, K1, K2 
469 -    def __set_keys(self, **kw): 
470          for k, v in kw.iteritems(): 
471              try: 
472                  setattr(self, k, v) 
473              except: 
474                  print 'MIL/MSL optimizer has no option named %s' % k 
475 -    def __onClickGetLL(self, btn): 
476          qubx.pyenv.env.OnScriptable('QubX.Modeling.MILRates.get_ll()') 
477          self.get_ll(wait=False) 
478 -    def get_ll_mil(self, *args, **kw): 
479          self.method = 'MIL' 
480          return self.get_ll(*args, **kw) 
481 -    def get_ll_msl(self, *args, **kw): 
482          self.method = 'MSL' 
483          return self.get_ll(*args, **kw) 
484 -    def get_ll(self, model_prep=None, data_prep=None, wait=True, receiver=lambda LL: None, on_trial=None, **kw): 
485          if wait: 
486              return qubx.pyenv.call_async_wait(self.__start_get_ll, model_prep, data_prep, on_trial=on_trial, **kw)[0] 
487          else: 
488              self.__start_get_ll(receiver, model_prep, data_prep, on_trial=on_trial, **kw) 
489 -    def __start_get_ll(self, receiver=lambda LL: None, model_prep=None, data_prep=None, on_trial=None, **kw): 
490          self.__begin_task(receiver) 
491          self.__set_keys(**kw) 
492          mp = model_prep 
493          if mp is None: 
494              mp = self.model_prep.val 
495              mp.set_rates2(*self.rates.val) 
496          self.robot.do(self.robot_get_ll, self.use_Peq, mp, data_prep, on_trial) 
497 -    def __onGotLL(self, LL, on_trial): 
498          self.__end_task(LL) 
499          self.txtLL.set_text('%.6g' % LL) 
500          self.__gray_ll(False) 
501          qubx.notebook.Notebook.send(qubx.notebook.NbText('MILRates: Get LL (%s): %.6g' % (self.robot.method, LL)), auto_id='Rates_eval.Messages') 
502          qubx.notebook.Notebook.send(self.QubX.Models.view.nbPicture, auto_id='Rates_eval.Model.Picture') 
503          qubx.notebook.Notebook.send(self.QubX.Models.view.nbPictureNO, auto_id='Rates_eval.Model.PictureNO') 
504          qubx.notebook.Notebook.send(self.QubX.Data.view.hires.nbPicture, auto_id='Rates_eval.Data.Picture') 
505          qubx.notebook.Notebook.send(self.QubX.Data.view.hires.nbPictureNO, auto_id='Rates_eval.Data.PictureNO') 
506          qubx.notebook.Notebook.send(self.QubX.Tables.find_view('Rates').nbTable, auto_id='Rates_eval.Model.Rates') 
507          # publish durhist after any possible durhist.robot activity 
508          self.QubX.Figures.DurHist.update(force=True) 
509          self.QubX.Figures.DurHist.robot.do(self.__dhrobot_send_durhists, 'Rates_eval.DurHist.Charts', 'Rates_eval.DurHist.TimeConst') 
510          self.OnLL(LL) 
511          if on_trial: 
512              self.LL = LL 
513              on_trial(*self.make_trial()) 
514 -    def __dhrobot_send_durhists(self, auto_id, tc_id): 
515          gobject.idle_add(self.__send_durhists, auto_id, tc_id) 
516 -    def __send_durhists(self, auto_id, tc_id): 
517          qubx.notebook.Notebook.send(self.QubX.Figures.DurHist.nbCharts, auto_id=auto_id) 
518          qubx.notebook.Notebook.send(self.QubX.Figures.DurHist.nbTimeConst, auto_id=tc_id) 
519 -    def __onClickAddToTrials(self, btn): 
520          qubx.pyenv.env.OnScriptable('QubX.Modeling.MILRates.add_to_trials(*QubX.Modeling.MILRates.make_trial())') 
521          self.add_to_trials(*self.make_trial()) 
522 -    def add_to_trials(self, mtree, row): 
523          Trials = qubx.global_namespace.QubX.Trials 
524          ts = Trials.get_trialset(self.method, in_qsf=True) 
525          ts.append_trial(mtree, row) 
526          Trials.trialset = ts 
527 -    def make_trial(self, fields={}, initial_tree=None, model_prep=None, data_prep=None): 
528          if model_prep: 
529              model = model_prep.source or model_prep.model 
530          else: 
531              model = self.QubX.Models.file 
532          if data_prep and data_prep.segs: 
533              data = data_prep.segs[0].file 
534          else: 
535              data = self.QubX.Data.file 
536          mtree = model.as_tree() 
537          nParam = fields['nParam'] if ('nParam' in fields) else 0 
538          row = make_trial_row(model.path, data.path, self.LL, nParam, nData=self.robot.Nevent, 
539                               nEvent=self.robot.Nevent) 
540          if initial_tree: 
541              add_trial_rates(row, initial_tree, 'Initial ') 
542          add_trial_rates(row, mtree) 
543          for k,v in fields.iteritems(): 
544              row[k] = v 
545          return mtree, row 
546 -    def __onClickOpt(self, btn): 
547          qubx.pyenv.env.OnScriptable('QubX.Modeling.MILRates.optimize(on_trial=QubX.Modeling.MILRates.add_to_trials)') 
548          self.optimize(wait=False, on_trial=self.add_to_trials) 
549 -    def optimize_mil(self, *args, **kw): 
550          self.method = 'MIL' 
551          return self.optimize(*args, **kw) 
552 -    def optimize_msl(self, *args, **kw): 
553          self.method = 'MSL' 
554          return self.optimize(*args, **kw) 
555 -    def optimize(self, model_prep=None, data_prep=None, wait=True, receiver=lambda LL: None, on_trial=None, **kw): 
556          if wait: 
557              return qubx.pyenv.call_async_wait(self.__start_opt, model_prep, data_prep, on_trial=on_trial, **kw)[0] 
558          else: 
559              self.__start_opt(receiver, model_prep, data_prep, on_trial=on_trial, **kw) 
560 -    def __start_opt(self, receiver=lambda LL: None, model_prep=None, data_prep=None, ext_return=False, on_trial=None, **kw): 
561          self.iterations = 0 
562          self.iter_trials = [] 
563          self.grads = [] 
564          self.__on_trial = on_trial 
565          self.__optimizing = True 
566          self.__begin_task(receiver) 
567          self.__set_keys(**kw) 
568          mp = model_prep 
569          if mp is None: 
570              mp = self.model_prep.val 
571              if mp.source: 
572                  mp.source.undoStack.enable(False) 
573              mp.set_rates2(*self.rates.val) 
574          self.__mtree_pre_opt = mp.model.as_tree() 
575          Trials = qubx.global_namespace.QubX.Trials 
576          Trials.get_trialset("%s Iterations"%self.method, in_qsf=True).clear() 
577          self.robot.do(self.robot_opt, self.use_Peq, 
578                        qubx.optimize.get_opts("qubx_single_rates_optimize", 
579                                               on_start=lambda rec, opts: gobject.idle_add(self.optOptions.optimizers.onStartOpt, rec, opts), 
580                                               on_end=lambda rec, opts: gobject.idle_add(self.optOptions.optimizers.onEndOpt, rec, opts)), 
581                        mp, data_prep, ext_return) 
582          if self.method == 'MIL': 
583              self.QubX.Figures.DurHist.request_show() 
584          self.QubX.Models.view.optimizing = True 
585 -    def __onIteration(self, LL, K0, K1, K2, model_prep, grads): 
586          self.iterations += 1 
587          self.LL = LL 
588          if model_prep.source: 
589              K012ToModel(K0, K1, K2, model_prep.source) 
590          self.txtLL.set_text('%.6g' % LL) 
591          self.__gray_ll(False) 
592          self.outStatus.write('Iter: %i\tLL: %8.5g\tgrads: %s\n' % (self.iterations, LL, ', '.join("%8.5g"%g for g in grads))) 
593          self.OnIteration(self.iterations, LL, self.grads) 
594 -    def __onOptDone(self, LL, workspaces, model_prep, ext_return): 
595          self.QubX.Models.view.optimizing = False 
596          self.outStatus.write('Finished optimizing (%s).                                                 LL = %.6g\n' % (self.robot.method, LL)) 
597          try: 
598              if model_prep.source: 
599                  model_prep.source.undoStack.enable() 
600                  model_prep.source.undoStack.push_undo(bind(model_prep.source.from_tree, self.__mtree_pre_opt), 
601                                                        bind(model_prep.source.from_tree, model_prep.source.as_tree())) 
602              self.__set_model_confidence(model_prep, self.robot.K, self.robot.dK) 
603              if model_prep.source: 
604                  model_prep.source.undoStack.seal_undo('set rates') 
605              self.txtLL.set_text('%.6g' % LL) 
606              self.__gray_ll(False) 
607          finally: 
608              if ext_return: 
609                  self.__end_task((LL, self.iterations, self.grads, self.hessian)) 
610              else: 
611                  self.__end_task(LL) 
612          if self.__on_trial: 
613              mtree, trial_row = self.make_trial({'nParam' : len(self.grads), 'Iterations' : self.iterations}, 
614                                                 initial_tree=self.__mtree_pre_opt, model_prep=model_prep, data_prep=self.robot.data_prep) 
615              for i, g in enumerate(self.grads): 
616                  trial_row['Gradient %i' % i] = g 
617              self.__on_trial(mtree, trial_row) 
618              self.__on_trial = None 
619          qubx.notebook.Notebook.send(qubx.notebook.NbText("""MILRates: Finished optimizing %d events after %d iterations (%s). 
620  Data:      %s 
621  Model:     %s 
622  P0:        %s 
623  LL:        %.6g""" % (self.robot.Nevent, self.robot.iter, self.robot.method, 
624                        '\n           '.join(os.path.split(ws.view.file.path)[1] for ws in workspaces), 
625                        model_prep.source and os.path.split(model_prep.source.path)[1] or "(from script)", 
626                        self.robot.Peq and "equilibrium" or "fixed (States table)", 
627                        LL)), auto_id='Rates_opt.Messages') 
628          qubx.notebook.Notebook.send(self.QubX.Models.view.nbPicture, auto_id='Rates_opt.Model.Picture') 
629          qubx.notebook.Notebook.send(self.QubX.Models.view.nbPictureNO, auto_id='Rates_opt.Model.PictureNO') 
630          qubx.notebook.Notebook.send(self.QubX.Data.view.hires.nbPicture, auto_id='Rates_opt.Data.Picture') 
631          qubx.notebook.Notebook.send(self.QubX.Data.view.hires.nbPictureNO, auto_id='Rates_opt.Data.PictureNO') 
632          qubx.notebook.Notebook.send(self.QubX.Tables.find_view('Rates').nbTable, auto_id='Rates_opt.Model.Rates') 
633          qubx.notebook.Notebook.send(self.QubX.Tables.find_view('Kinetic Constraints').nbTable, auto_id='Rates_opt.Model.Constraints') 
634          # publish durhist after any possible durhist.robot activity 
635          self.QubX.Figures.DurHist.update(force=True) 
636          self.QubX.Figures.DurHist.robot.do(self.__dhrobot_send_durhists, 'Rates_opt.DurHist.Charts', 'Rates_opt.DurHist.TimeConst') 
637   
638          ts = qubx.global_namespace.QubX.Trials.get_trialset("%s Iterations"%self.method, in_qsf=True) 
639          for mtree, row in self.iter_trials: 
640              ts.append_trial(mtree, row) 
641   
642          self.OnOptimized(self.iterations, LL, self.grads, self.hessian) 
643 -    def __set_model_confidence(self, model_prep, K, dK): 
644          have_dK = any(dK) 
645          self.outStatus.write('%10s%7s' % ("Rate From", "To"), scroll=False) 
646          self.outStatus.write_bold('%14s' % "k0", scroll=False) 
647          if have_dK: 
648              self.outStatus.write('%11s' % "+/-", scroll=False) 
649          self.outStatus.write_bold('%14s' % 'k1', scroll=False) 
650          if have_dK: 
651              self.outStatus.write('%11s' % "+/-", scroll=False) 
652          self.outStatus.write_bold('%14s' % 'k2', scroll=False) 
653          if have_dK: 
654              self.outStatus.write('%11s' % "+/-", scroll=False) 
655          self.outStatus.write('\n') 
656          model = model_prep.source 
657          if not model: 
658              return 
659          k_ix = K_index(model.states.size) 
660          rr = model.rates 
661          for r in xrange(rr.size): 
662              a, b = [rr.get(r, f) for f in ('From', 'To')] 
663              i = k_ix.index((a,b)) 
664              self.outStatus.write('%10d%7d' % (a,b), scroll=False) 
665              self.outStatus.write_bold('%14.5g' % rr.get(r, 'k0'), scroll=False) 
666              k0_lo = k0_hi = k0_std = 0.0 
667              if dK[i]: 
668                  k0_lo = float(exp(K[i] - dK[i])) 
669                  k0_hi = float(exp(K[i] + dK[i])) 
670                  k0_std = float(exp(K[i]) * dK[i]) 
671                  rr.set(r, 'k0 lo', k0_lo) 
672                  rr.set(r, 'k0 hi', k0_hi) 
673                  rr.set(r, 'k0 std', k0_std) 
674              if have_dK: 
675                  self.outStatus.write('%11.3g' % k0_std, scroll=False) 
676              i += len(k_ix) 
677              self.outStatus.write_bold('%14.5g' % K[i], scroll=False) 
678              k1_lo = k1_hi = k1_std = 0.0 
679              if dK[i]: 
680                  k1_lo = float(K[i] - dK[i]) 
681                  k1_hi = float(K[i] + dK[i]) 
682                  k1_std = float(dK[i]) 
683                  rr.set(r, 'k1 lo', k1_lo) 
684                  rr.set(r, 'k1 hi', k1_hi) 
685                  rr.set(r, 'k1 std', k1_std) 
686              if have_dK: 
687                  self.outStatus.write('%11.3g' % k1_std, scroll=False) 
688              i += len(k_ix) 
689              self.outStatus.write_bold('%14.5g' % K[i], scroll=False) 
690              k2_lo = k2_hi = k2_std = 0.0 
691              if dK[i]: 
692                  k2_lo = float(K[i] - dK[i]) 
693                  k2_hi = float(K[i] + dK[i]) 
694                  k2_std = float(dK[i]) 
695                  rr.set(r, 'k2 lo', k2_lo) 
696                  rr.set(r, 'k2 hi', k2_hi) 
697                  rr.set(r, 'k2 std', k2_std) 
698              if have_dK: 
699                  self.outStatus.write('%11.3g' % k2_std, scroll=False) 
700              self.outStatus.write('\n') 
701 -    def robot_on_status(self, msg): 
702          gobject.idle_add(self.outStatus.write, msg+"\n") 
703 -    def robot_setup(self, Peq, model_prep, data_prep): 
704          self.robot.Peq = Peq 
705          self.robot.model_prep = model_prep 
706          self.robot.model = self.robot.model_prep.model 
707          self.robot.K = K012toK(self.robot.model_prep.K0, self.robot.model_prep.K1, self.robot.model_prep.K2) 
708          self.robot.clazz = ModelToClazz(self.robot.model) 
709          self.robot.method = self.method 
710          if self.robot.Peq: 
711              self.robot.P0 = None 
712          else: 
713              self.robot.P0 = ModelToP0(self.robot.model) 
714          if self.robot.clazz.shape[0] == 0: 
715              self.robot.data_prep = Anon(ws=[]) # no model -> do nothing 
716          elif data_prep: 
717              self.robot.data_prep = data_prep 
718          else: 
719              self.plexig.model_prep = model_prep 
720              self.robot.data_prep = self.robot.gui_call_recv(self.plexig.request, self.robot)[0] 
721          self.robot.workspaces = self.robot.data_prep.ws 
722          self.robot.Nevent = 0 
723          for ws in self.robot.workspaces: 
724              try: 
725                  ws.single_model.dispose() # work around callback-related memory leak 
726              except: 
727                  pass 
728              ws.single_model = ModelToFast(self.robot.model, ws.sig_names, self.robot_on_status) 
729              ws.multi_model = qubx.fast.model.MultiModel(ws.single_model) 
730              ws.multi_model.update() 
731              if self.robot.Peq: 
732                  ws.P0 = None 
733              else: 
734                  ws.P0 = buffer_as_array_of_double(ws.multi_model.multi_model[0].P0, ws.multi_model.multi_model[0].Nstate) 
735              if self.robot.method == 'MIL': 
736                  if len(ws.scls) > 1: 
737                      gobject.idle_add(self.outStatus.write, "MIL can't analyze data with a changing stimulus.  Please use MSL instead.\n") 
738                      return False 
739                  elif __builtins__['any'](mc[0] < 0 for mc in ws.mcls): 
740                      gobject.idle_add(self.outStatus.write, "MIL can't analyze data with blank or excluded regions.  Please extract, or use MSL instead.\n") 
741                      return False 
742              self.robot.Nevent += sum(len(durs) for clss, durs in ws.segs) 
743              if not self.cl_context: 
744                  self.__robot_setup_cl_context() 
745              if not ws.cl_data: 
746                  ws.cl_data = msl_accel_data(self.cl_context, ws.segs, ws.mcls, ws.scls) 
747              ws.cl_model = msl_accel_models(self.cl_context, cl_dimension(ws.multi_model.multi_model[0].Nstate), 1, len(ws.scls)) 
748          return True 
749 -    def robot_slice_data(self, lbl, segs, mcls, scls, td): 
750          self.robot.Nevent += sum(len(durs) for clss, durs in segs) 
751          gobject.idle_add(self.outStatus.write, '%s: %i events\n' % (lbl, self.robot.Nevent)) 
752          return [([seg], mcls, scls) for seg in segs if len(seg[0])] 
753 -    def robot_eval(self, K): 
754          if not self.robot.workspaces: 
755              gobject.idle_add(self.outStatus.write, "The model is empty.\n") 
756              return 0.0 
757          elif not self.robot.data_prep.ws: 
758              gobject.idle_add(self.outStatus.write, "There are no data files in this group (see the Data table).\n") 
759              return 0.0 
760          if self.robot.Nevent == 0: 
761              gobject.idle_add(self.outStatus.write, "To use this panel, first Idealize your single-molecule data.\nFor ensemble data, try MacRates.\n") 
762              return 0.0 
763          LL = 0.0 
764          K0, K1, K2 = KtoK012(K, self.robot.workspaces[0].single_model.K0) 
765          argses = [] 
766          iarr, imat, dmat = buffer_as_array_of_int, buffer_as_matrix_of_int, buffer_as_matrix_of_double 
767          for ws in self.robot.workspaces: 
768              ws.single_model.K0[:,:] = K0[:,:] 
769              ws.single_model.K1[:,:] = K1[:,:] 
770              ws.single_model.K2[:,:] = K2[:,:] 
771              ws.multi_model.update() 
772              mm = ws.multi_model.multi_model 
773              Ns = mm[0].Nstate 
774              method = self.cl_context.run_mil if (self.robot.method == 'MIL') else self.cl_context.run_msl 
775              ws.cl_model.reset() 
776              if ws.cl_model.setup_model(ws.cl_data, iarr(mm[0].clazz, Ns), ws.P0, dmat(mm[0].K0[0], Ns, Ns), dmat(mm[0].K1[0], Ns, Ns), dmat(mm[0].K2[0], Ns, Ns), 
777                                         imat(mm[0].L[0], Ns, Ns), imat(mm[0].V[0], Ns, Ns), imat(mm[0].P[0], Ns, Ns), 
778                                         ws.td_ms*1e-3, lambda msg, obj: gobject.idle_add(self.outStatus.write, msg+"\n")): 
779                  return -1e21 
780              if method(ws.cl_data, ws.cl_model, self.useGPU and MSL_ACCEL_OPENCL or MSL_ACCEL_OPENMP): 
781                  gobject.idle_add(self.disable_gpu) 
782                  if method(ws.cl_data, ws.cl_model, MSL_ACCEL_OPENMP): 
783                      return -1e20 
784              LL += ws.cl_model.ll[0] 
785          return LL 
786 -    def disable_gpu(self): 
787          self.useGPU = False 
788          self.chkGPU.set_sensitive(False) 
789          self.chkGPU.set_tooltip_text("OpenCL is not available; needs a better graphics card or new drivers.") 
790 -    def robot_get_ll(self, Peq, model_prep, data_prep, on_trial): 
791          if not self.robot_setup(Peq, model_prep, data_prep): 
792              gobject.idle_add(self.__onGotLL, 0.0, on_trial) 
793          else: 
794              try: 
795                  LL = self.robot_eval(self.robot.K) 
796                  gobject.idle_add(self.__onGotLL, LL, on_trial) 
797              except: 
798                  gobject.idle_add(self.__end_task, -1.0) 
799                  raise 
800          del self.robot.workspaces 
801 -    def robot_opt(self, Peq, optimize, model_prep, data_prep, ext_return): 
802          self.robot.optimize = optimize 
803          self.robot.progress = 0.0 
804          self.robot.iter = 0 
805          self.robot.LL = 0.0 
806          if not self.robot_setup(Peq, model_prep, data_prep): 
807              if ext_return: 
808                  gobject.idle_add(self.__end_task, (-1.0, 0, [], None)) 
809              else: 
810                  gobject.idle_add(self.__end_task, -1.0) 
811              return 
812          try: 
813              K0, K1, K2 = KtoK012(self.robot.K, self.robot.workspaces[0].single_model.K0) 
814              A, B = self.robot.model.constraints_kin.get_matrices_p(K0, K1, K2, 
815                                                                     self.robot.workspaces[0].single_model.L, self.robot.workspaces[0].single_model.V, 
816                                                                     self.robot.workspaces[0].single_model.P) 
817              #A, B = reduce_constraints(A, B) 
818              #self.robot.A, self.robot.B, self.robot.Ainv, pars = linear_constraints(A, B, self.robot.K) 
819              self.robot.Kr, AinR, BinR, self.robot.Kr_indices = qubx.fast.model.EliminateFixedParams(self.robot.K, A, B) 
820              A, B = qubx.fast.model.ReduceConstraints(AinR, BinR) 
821              self.robot.A, self.robot.B, self.robot.Ainv, pars = qubx.fast.model.SetupLinearConstraints(A, B, self.robot.Kr) 
822              self.robot.Asc, self.robot.Bsc, self.robot.Asci, pars = start_at_ones(pars) 
823              self.robot.pars = array([x for x in pars.A[:,0]]) 
824              self.robot.npar = len(self.robot.pars) 
825              if self.robot.npar > 0: 
826                  for ws in self.robot.workspaces: 
827                      ws.cl_models = msl_accel_models(self.cl_context, cl_dimension(ws.multi_model.multi_model[0].Nstate), 1+2*self.robot.npar, len(ws.scls)) 
828                  score, pars, iterations, grads, Hessian = self.robot.optimize(array(self.robot.pars), self.robot_func, self.robot_iter, self.robot_grad_func) 
829                  # make sure final rates are showing (if it overshot and gave up) 
830                  self.robot_iter(self.robot.pars, grads, iterations, - self.robot.LL) 
831              else: 
832                  score, pars, iterations, grads, Hessian = self.robot_eval(self.robot.K), [], 0, [], numpy.zeros(shape=(0,0)) 
833                  self.robot.LL = - score 
834              dKr = CalcStdOfK(self.robot.A, self.robot.Asc, Hessian) 
835              self.robot.dK = zeros(shape=self.robot.K.shape, dtype='float64') 
836              qubx.fast.model.Kr_to_K(dKr, self.robot.dK, self.robot.Kr_indices) 
837              self.iterations = iterations 
838              self.grads = grads 
839              self.hessian = Hessian 
840          except: 
841              traceback.print_exc() 
842              score = -1.0 
843              if ext_return: 
844                  gobject.idle_add(self.__end_task, -1.0, self.robot.iter, [], None) 
845              else: 
846                  gobject.idle_add(self.__end_task, -1.0) 
847              raise 
848          self.LL = score 
849          gobject.idle_add(self.__onOptDone, self.robot.LL, self.robot.workspaces, self.robot.model_prep, ext_return) 
850          del self.robot.workspaces 
851          del self.robot.model_prep 
852 -    def robot_func(self, pars): 
853          self.robot.Kr = ParsToK(parsToK(array(pars[:self.robot.npar]), self.robot.Asc, self.robot.Asci, self.robot.Bsc), self.robot.A, self.robot.Ainv, self.robot.B) 
854          qubx.fast.model.Kr_to_K(self.robot.Kr, self.robot.K, self.robot.Kr_indices) 
855          LL = self.robot_eval(self.robot.K) 
856          #print -LL,'of',self.robot.K 
857          return - LL 
858 -    def robot_grad_func(self, pars, grads, dx): 
859          iarr, imat, dmat = buffer_as_array_of_int, buffer_as_matrix_of_int, buffer_as_matrix_of_double 
860          LL = 0.0 
861          for ix in xrange(self.robot.npar): 
862              grads[ix] = 0.0 
863          valids = [] 
864          xpars = numpy.array(pars[:self.robot.npar], dtype='float64', copy=True) 
865          for ws in self.robot.workspaces: 
866              ws.cl_models.reset() 
867              for ix in xrange(0, self.robot.npar+1): 
868                  deltas = [0] if (ix == self.robot.npar) else [-dx, dx] 
869                  for d in deltas: 
870                      xpars[:self.robot.npar] = pars[:self.robot.npar] 
871                      if d: 
872                          xpars[ix] += d 
873   
874                      self.robot.Kr = ParsToK(parsToK(xpars, self.robot.Asc, self.robot.Asci, self.robot.Bsc), self.robot.A, self.robot.Ainv, self.robot.B) 
875                      qubx.fast.model.Kr_to_K(self.robot.Kr, self.robot.K, self.robot.Kr_indices) 
876                      #print 'q',self.robot.K 
877                      K0, K1, K2 = KtoK012(self.robot.K, self.robot.workspaces[0].single_model.K0) 
878                      ws.single_model.K0[:,:] = K0[:,:] 
879                      ws.single_model.K1[:,:] = K1[:,:] 
880                      ws.single_model.K2[:,:] = K2[:,:] 
881                      ws.multi_model.update() 
882                      mm = ws.multi_model.multi_model 
883                      Ns = mm[0].Nstate 
884                       
885                      valid = True 
886                      if ws.cl_models.setup_model(ws.cl_data, iarr(mm[0].clazz, Ns), ws.P0, dmat(mm[0].K0[0], Ns, Ns), dmat(mm[0].K1[0], Ns, Ns), dmat(mm[0].K2[0], Ns, Ns), 
887                                                  imat(mm[0].L[0], Ns, Ns), imat(mm[0].V[0], Ns, Ns), imat(mm[0].P[0], Ns, Ns), 
888                                                  ws.td_ms*1e-3, lambda msg, obj: gobject.idle_add(self.outStatus.write, msg+"\n")): 
889                          valid = False 
890                          if not d: # central params must succeed: 
891                              return 1e20 # usually positive eigenvalues 
892                      valids.append(valid) 
893              method = self.cl_context.run_mil if (self.robot.method == 'MIL') else self.cl_context.run_msl 
894              ws.cl_err = method(ws.cl_data, ws.cl_models, self.useGPU and MSL_ACCEL_OPENCL or MSL_ACCEL_OPENMP) 
895                  #print ws.cl_data.ll 
896                  #print ws.cl_models.ll 
897                  #print LL 
898              if ws.cl_err: 
899                  gobject.idle_add(self.disable_gpu) 
900                      # retry no opencl 
901                  ws.cl_err = method(ws.cl_data, ws.cl_models, MSL_ACCEL_OPENMP) 
902                  if ws.cl_err: 
903                      return 1e20 
904              LL += ws.cl_models.ll[ws.cl_models.Nmodel-1] 
905              for ix in xrange(self.robot.npar): 
906                  if valids[2*ix+1] and valids[2*ix]: 
907                      grads[ix] -= (ws.cl_models.ll[2*ix+1] - ws.cl_models.ll[2*ix]) / (2*dx) # ws.cl_models.Nmodel-1]) / dx #  
908                  elif valids[2*ix+1]: 
909                      grads[ix] -= (ws.cl_models.ll[2*ix+1] - ws.cl_models.ll[ws.cl_models.Nmodel-1]) / dx; 
910                  elif valids[2*ix]: 
911                      grads[ix] -= (ws.cl_models.ll[ws.cl_models.Nmodel-1] - ws.cl_models.ll[2*ix]) / dx; 
912          return -LL 
913 -    def robot_iter(self, pars, grads, iterations, score): 
914          self.robot.iter = iterations 
915          self.robot.LL = - score 
916          self.robot.pars[:] = pars[:self.robot.npar] 
917          self.robot.progress = (self.robot.iter * 100.0 / MAXITER) 
918          self.robot.Kr = ParsToK(parsToK(self.robot.pars, self.robot.Asc, self.robot.Asci, self.robot.Bsc), self.robot.A, self.robot.Ainv, self.robot.B) 
919          qubx.fast.model.Kr_to_K(self.robot.Kr, self.robot.K, self.robot.Kr_indices) 
920          K0, K1, K2 = KtoK012(self.robot.K, self.robot.workspaces[0].single_model.K0) 
921          gobject.idle_add(self.__onIteration, - score, K0, K1, K2, self.robot.model_prep, grads) 
922          mtree, row = self.make_trial({"Iteration": iterations}, model_prep=self.robot.model_prep, data_prep=self.robot.data_prep) 
923          self.iter_trials.append((mtree, row)) 
924   
925   
926 -def inter_ll_from_subi_args(args): 
927      use_lib, clazz, P0, K0, K1, K2, L, V, P, tdead_sec, segs, mcls, scls = args 
928      constants = scls[0] 
929      renumber = array([mc[0] for mc in mcls], dtype='int32') 
930      mil_segs = [(array(renumber[classes], copy=True, dtype='int32'), durations) for classes, durations in segs] 
931      mil_args = (clazz, P0, K0, K1, K2, L, V, P, constants, tdead_sec, mil_segs) 
932      return inter_ll_or_0(mil_args) 
933