Source Code for Module qubx.kalman

  1  """Kalman filter panel for Data. 
  2   
  3  Copyright 2014 Research Foundation State University of New York  
  4  This file is part of QUB Express.                                           
  5   
  6  QUB Express is free software; you can redistribute it and/or modify           
  7  it under the terms of the GNU General Public License as published by  
  8  the Free Software Foundation, either version 3 of the License, or     
  9  (at your option) any later version.                                   
 10   
 11  QUB Express is distributed in the hope that it will be useful,                
 12  but WITHOUT ANY WARRANTY; without even the implied warranty of        
 13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         
 14  GNU General Public License for more details.                          
 15   
 16  You should have received a copy of the GNU General Public License,    
 17  named LICENSE.txt, in the QUB Express program directory.  If not, see         
 18  <http://www.gnu.org/licenses/>.                                       
 19   
 20  """ 
 21   
 22  import gobject 
 23  import gtk 
 24  from gtk import gdk 
 25  import numpy 
 26  import os 
 27  import qubx.chop 
 28  import qubx.data_types 
 29  import qubx.faces 
 30  import qubx.fast.kalman 
 31  import qubx.GTK 
 32  import qubx.notebook 
 33  import qubx.pyenv 
 34  import traceback 
 35  from itertools import izip, count 
 36  from qubx.accept import * 
 37  from qubx.GTK import pack_item, pack_space, pack_hsep, pack_vsep, pack_label, pack_button, pack_check, pack_radio, pack_scrolled 
 38  from qubx.settings import Property, Propertied 
 39  from qubx.task import * 
 40  from qubx.util_types import WeakEvent, Reffer, bind, bind_with_args_before 
 41   
 42 -def task_exception_to_console(task, typ, val, tb): 
 43      traceback.print_exception(typ, val, tb) 
 44   
 45 -class KalmanTask(Task): 
 46 -    def __init__(self, receiver, **kw): 
 47          Task.__init__(self, 'Kalman') 
 48          self.__receiver = receiver 
 49          self.__kw = kw 
 50          QubX = qubx.global_namespace.QubX 
 51          self.dataview = QubX.Data.view 
 52          self.iSignal = QubX.DataSource.signal 
 53          self.datafits = self.dataview.file.fits[self.iSignal].idl 
 54          self.datafits.clear() 
 55          self.chunks = chunks = [] 
 56          segs = QubX.DataSource.get_segmentation() 
 57          for seg in segs: 
 58              first = True 
 59              for chunk in seg.chunks: 
 60                  if chunk.included: 
 61                      chunk_w_data = chunk.get_samples() 
 62                      chunk_w_data.is_first = first 
 63                      first = False 
 64                      chunk_w_data.src_seg = seg.index 
 65                      chunk_w_data.off_in_seg = chunk.f - seg.offset 
 66                      chunks.append(chunk_w_data) 
 67 -    def run(self): 
 68          try: 
 69              qubx.task.Tasks.add_task(self) 
 70              self.run_kalman(**self.__kw) 
 71          finally: 
 72              gobject.idle_add(self.__receiver) 
 73              qubx.task.Tasks.remove_task(self) 
 74 -    def run_kalman(self, process=1.0, process_std=1e-5, process_to_measurement=1.0, measurement_std=1e-3, process_init=UNSET_VALUE, show_stds=False): 
 75          # analyze each included chunk: 
 76          kf = qubx.fast.kalman.Filter1D(process, process_std**2, process_to_measurement, measurement_std**2) 
 77          for chunk in self.chunks: 
 78              try: 
 79                  if chunk.is_first: 
 80                      kf.reset() 
 81                      if process_init == UNSET_VALUE: 
 82                          kf.X[0] = chunk.samples[0] / process_to_measurement 
 83                      else: 
 84                          kf.X[0] = process_init 
 85                  arr = chunk.samples.astype('float64') 
 86                  del chunk.samples 
 87                  if show_stds: 
 88                      stds = kf.filter(arr, True) 
 89                      gobject.idle_add(self.__onResult, chunk, arr.astype('float32'), stds.astype('float32')) 
 90                  else: 
 91                      kf.filter(arr) 
 92                      gobject.idle_add(self.__onResult, chunk, arr.astype('float32')) 
 93              except: 
 94                  traceback.print_exc() 
 95          # print 'final std:', sqrt(kf.P[0,0]) 
 96          gobject.idle_add(self.__onDone) 
 97 -    def __onResult(self, chunk, fit, stds=None): 
 98          ff = numpy.arange(chunk.f, chunk.l+1, dtype='int32') 
 99          self.datafits.set_fit(ff, ff, fit, stds) 
100 -    def __onDone(self): 
101          self.dataview.file.OnChangeFits(self.dataview.file, self.iSignal) 
102   
103   
104  @Propertied(Property('process', 1.0, "process scalar"), 
105              Property('process_std', 1e-5, "process standard deviation"), 
106              Property('process_to_measurement', 1.0, "process to measurement scalar"), 
107              Property('measurement_std', 1e-3, "measurement standard deviation"), 
108              Property('process_init', UNSET_VALUE, "initial process value"), 
109              Property('show_stds', False, "Display +/- one standard deviation")) 
110 -class KalmanFace(qubx.faces.Face): 
111      __explore_featured = ['run'] 
112 -    def __init__(self, name="Kalman", global_name="QubX.Tools.Kalman"): 
113          qubx.faces.Face.__init__(self, name, global_name) 
114          self.set_size_request(300, 130) 
115          self.propertied_connect_settings('Kalman') 
116          self.__ref = Reffer() 
117           
118          v = gtk.VBox() 
119          self.scroll = pack_scrolled(v, self, with_vp=True, expand=True) 
120          vh = pack_item(gtk.HBox(True), v) 
121          pack_label('Process Std.:', vh) 
122          self.txtProcessStd = pack_item(qubx.GTK.NumEntry(self.process_std, acceptFloatGreaterThanOrEqualTo(0.0), '%.3g'), vh) 
123          self.propertied_connect_NumEntry('process_std', self.txtProcessStd) 
124          vh = pack_item(gtk.HBox(True), v) 
125          pack_label('Measurement Std.:', vh) 
126          self.txtMeasurementStd = pack_item(qubx.GTK.NumEntry(self.measurement_std, acceptFloatGreaterThanOrEqualTo(0.0), '%.3g'), vh) 
127          self.propertied_connect_NumEntry('measurement_std', self.txtMeasurementStd) 
128          vh = pack_item(gtk.HBox(True), v) 
129          pack_label('Start value:', vh) 
130          self.txtProcessInit = pack_item(qubx.GTK.NumEntry(self.process_init, acceptFloatOrUnset, formatFloatOrUnset('%.3g')), vh) 
131          self.txtProcessInit.set_tooltip_text('If this is blank, the start value comes from the first data point') 
132          self.propertied_connect_NumEntry('process_init', self.txtProcessInit) 
133          vh = pack_item(gtk.HBox(), v) 
134          self.chkShowStds = pack_check('Show +/- one std. (posterior distribution)', vh, self.show_stds) 
135          self.propertied_connect_check('show_stds', self.chkShowStds) 
136          vh = pack_item(gtk.HBox(True), v) 
137          self.__live = False 
138          self.__live_serial = 0 
139          self.chkLive = pack_check('live', vh, on_toggle=self.__onToggleLive) 
140          self.chkLive.set_tooltip_text('Re-run the filter automatically when settings change') 
141          vhh = pack_item(gtk.HBox(), vh) 
142          self.btnRun = pack_button('Run filter on Data Source', vhh, self.__ref(self.__onClickRun)) 
143          QubX = qubx.global_namespace.QubX 
144          QubX.DataSource.OnChangeSamples += self.__ref(self.try_live) 
145 -    def propertied_set(self, value, name): 
146          super(KalmanFace, self).propertied_set(value, name) 
147          self.try_live() 
148 -    def try_live(self): 
149          if self.__live: 
150              self.__live_serial += 1 
151              gobject.idle_add(self.run_live, self.__live_serial) 
152 -    def run_live(self, serial): 
153          if serial != self.__live_serial: 
154              return 
155          if self.btnRun.get_sensitive(): 
156              self.run(wait=False) 
157          else: 
158              gobject.timeout_add(300, self.try_live) 
159 -    def __onToggleLive(self, chk): 
160          self.__live = chk.get_active() 
161          self.try_live() 
162 -    def __onClickRun(self, btn): 
163          qubx.pyenv.env.OnScriptable('%s.run(process_std=%s, measurement_std=%s, process_init=%s, show_stds=%s)' % (self.global_name, self.process_std, self.measurement_std, formatFloatOrUnset(self.process_init), repr(self.show_stds))) 
164          self.run(wait=False) 
165 -    def run(self, process_std=None, measurement_std=None, process_init=None, show_stds=None, wait=True, receiver=lambda: None): 
166          if wait: 
167              return qubx.pyenv.call_async_wait(lambda rcv: self.run(process_std, measurement_std, process_init, show_stds, wait=False, receiver=rcv)) 
168          kw = {} 
169          kw['process_std'] = self.process_std if (process_std is None) else process_std 
170          kw['measurement_std'] = self.measurement_std if (measurement_std is None) else measurement_std 
171          kw['process_init'] = self.process_init if (process_init is None) else process_init 
172          kw['show_stds'] = self.show_stds if (show_stds is None) else show_stds 
173           
174          self.btnRun.set_sensitive(False) 
175          task = KalmanTask(lambda: self.__onDone(receiver), **kw) 
176          task.OnException += task_exception_to_console 
177          task.start() 
178 -    def __onDone(self, receiver): 
179          self.btnRun.set_sensitive(True) 
180          gobject.idle_add(receiver) 
181