MRS MIDI II support added

akvo/gui/akvoGUI.py

@@ -120,6 +120,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         # Menu items #
         ##############
         self.ui.actionOpen_GMR.triggered.connect(self.openGMRRAWDataset)
+        self.ui.actionLoad_MIDI.triggered.connect(self.loadMIDI2Dataset)
         self.ui.actionSave_Preprocessed_Dataset.triggered.connect(self.SavePreprocess)
         self.ui.actionExport_Preprocessed_Dataset.triggered.connect(self.ExportPreprocess)
         self.ui.actionExport_Preprocessed_Dataset.setEnabled(False)
@@ -129,7 +130,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         ###########
         # Buttons #
         ###########
-        self.ui.loadDataPushButton.pressed.connect(self.loadRAW) 
+        #self.ui.loadDataPushButton.pressed.connect(self.loadRAW) 
         self.ui.sumDataGO.pressed.connect( self.sumDataChans )
         self.ui.bandPassGO.pressed.connect( self.bandPassFilter )
         self.ui.filterDesignPushButton.pressed.connect( self.designFilter )
@@ -884,7 +885,79 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         self.RAWDataProc.enableDSPTrigger.connect(self.enableDSP)
         self.RAWDataProc.doneTrigger.connect(self.doneStatus)
         self.RAWDataProc.updateProcTrigger.connect(self.updateProc)
+    
+    def loadMIDI2Dataset (self):
+        """ Opens a MIDI file and extracts header info
+        """
+        try:
+            with open('.midi2.last.path') as f: 
+                fpath = f.readline()  
+                pass
+        except IOError as e:
+            fpath = '.'
+         
+        self.headerstr = QtWidgets.QFileDialog.getExistingDirectory(self, 'Load MIDI 2 Directory', fpath)
+        self.ui.headerFileTextBrowser.clear()
+        self.ui.headerFileTextBrowser.append(self.headerstr)
+        
+        if len(self.headerstr) == 0:
+            return 
+
+        # TODO, should rename this class, use MIDI path  
+        self.connectGMRDataProcessor()
+        self.RAWDataProc.readMIDI2Header(str(self.headerstr))
 
+        # look in the directory for all the files 
+        self.ui.loadDataPushButton.pressed.connect(self.loadMIDI2) 
+
+        # If we got this far, enable all the widgets
+        self.ui.lcdNumberTauPulse1.setEnabled(True)
+        self.ui.lcdNumberNuTx.setEnabled(True)
+        self.ui.lcdNumberTuneuF.setEnabled(True)
+        self.ui.lcdNumberSampFreq.setEnabled(True)
+        self.ui.lcdNumberNQ.setEnabled(True)
+
+        self.ui.headerFileBox.setEnabled(True)
+        self.ui.headerFileBox.setChecked( True )
+        self.ui.headerBox2.setVisible(True) 
+        self.ui.inputRAWParametersBox.setEnabled(True)
+        self.ui.loadDataPushButton.setEnabled(True)
+         
+        # make plots as you import the dataset
+        self.ui.plotImportCheckBox.setEnabled(True)
+        self.ui.plotImportCheckBox.setChecked(True)
+         
+        # Update info from the header into the GUI
+        self.ui.pulseTypeTextBrowser.clear()
+        self.ui.pulseTypeTextBrowser.append(self.RAWDataProc.pulseType)
+        self.ui.lcdNumberNuTx.display(self.RAWDataProc.transFreq)
+        self.ui.lcdNumberTauPulse1.display(1e3*self.RAWDataProc.pulseLength[0])
+        self.ui.lcdNumberTuneuF.display(self.RAWDataProc.TuneCapacitance)
+        self.ui.lcdNumberSampFreq.display(self.RAWDataProc.samp)
+        self.ui.lcdNumberNQ.display(self.RAWDataProc.nPulseMoments)
+        self.ui.DeadTimeSpinBox.setValue(1e3*self.RAWDataProc.deadTime)
+        self.ui.CentralVSpinBox.setValue( self.RAWDataProc.transFreq )
+
+        # set the B0 field according to Tx as an initial guess
+        self.ui.intensitySpinBox.setValue( self.RAWDataProc.transFreq/GAMMAH )           
+ 
+        if self.RAWDataProc.pulseType != "FID":
+            self.ui.lcdNumberTauPulse2.setEnabled(1)
+            self.ui.lcdNumberTauPulse2.display(1e3*self.RAWDataProc.pulseLength[1])
+            self.ui.lcdNumberTauDelay.setEnabled(1)
+            self.ui.lcdNumberTauDelay.display(1e3*self.RAWDataProc.interpulseDelay)
+        
+        self.ui.FIDProcComboBox.clear() 
+        if self.RAWDataProc.pulseType == "4PhaseT1" or self.RAWDataProc.pulseType == "T1":
+            self.ui.FIDProcComboBox.insertItem(0, "Pulse 1") 
+            self.ui.FIDProcComboBox.insertItem(1, "Pulse 2") 
+            self.ui.FIDProcComboBox.insertItem(2, "Both")    
+            self.ui.FIDProcComboBox.setCurrentIndex (1)
+        elif self.RAWDataProc.pulseType == "FID":
+            self.ui.FIDProcComboBox.insertItem(0, "Pulse 1") 
+            self.ui.FIDProcComboBox.setCurrentIndex (0)
+    
+ 
     def openGMRRAWDataset(self):
         """ Opens a GMR header file
         """
@@ -913,6 +986,9 @@ class ApplicationWindow(QtWidgets.QMainWindow):
 
         self.connectGMRDataProcessor()
         self.RAWDataProc.readHeaderFile(str(self.headerstr))
+        
+        # make sure we will use GMR path 
+        self.ui.loadDataPushButton.pressed.connect(self.loadRAW) 
 
         # If we got this far, enable all the widgets
         self.ui.lcdNumberTauPulse1.setEnabled(True)
@@ -1082,6 +1158,12 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         INFO["nDAQVersion"] = self.RAWDataProc.nDAQVersion
         INFO["log"] = yaml.dump( self.YamlNode )  
 
+        # 1.6.4 and on 
+        INFO["Instrument"] = self.RAWDataProc.Instrument 
+        if self.RAWDataProc.Instrument == "MIDI 2":
+            INFO["MIDIGain"] = self.RAWDataProc.MIDIGain
+            INFO["datadir"] = self.RAWDataProc.datadir
+
         TXRX = []
         for ir in range(0, self.ui.txRxTable.rowCount() ):
             txrx = []
@@ -1096,12 +1178,14 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         if "Gate integrate" in self.YamlNode.Stacking.keys():
             INFO["GATED"] = self.RAWDataProc.GATED
 
-        print("META SAVE")
-        print("INFO log", INFO["log"])
+        #print("META SAVE")
+        #print("INFO log", INFO["log"])
 
         self.RAWDataProc.DATADICT["INFO"] = INFO 
 
         pickle.dump(self.RAWDataProc.DATADICT, save)
+        #pickle.dump(self.RAWDataProc, save) # doesn't work :-( 
+
         save.close()
 
     # Export XML file suitable for USGS ScienceBase Data Release
@@ -1296,7 +1380,109 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         self.RAWDataProc.doneTrigger.connect(self.doneStatus)
         
         self.enableAll()
+
+    def loadMIDI2(self):
+
+        #################################################
+        # Check to make sure we are ready to process
+
+        # Header
+        if self.RAWDataProc == None:
+            err_msg = "You need to load a header first."
+            reply = QtWidgets.QMessageBox.critical(self, 'Error', 
+                err_msg) #, QtWidgets.QMessageBox.Yes, QtWidgets.QMessageBox.No)
+            return
+        
+        # Stacks 
+        try:
+            self.procStacks = np.array(eval(str("np.r_["+self.ui.stacksLineEdit.text())+"]"))
+        except:
+            err_msg = "You need to set your stacks correctly.\n" + \
+                      "This should be a Python Numpy interpretable list\n" + \
+                      "of stack indices. For example 1:24 or 1:4,8:24"
+            QtWidgets.QMessageBox.critical(self, 'Error', err_msg) 
+            return
+
+        # Data Channels
+        #Chan = np.arange(0,9,1)
+        try:
+            self.dataChan = np.array(eval(str("np.r_["+self.ui.dataChanLineEdit.text())+"]"))
+        except:
+            #QMessageBox messageBox;
+            #messageBox.critical(0,"Error","An error has occured !");
+            #messageBox.setFixedSize(500,200);
+            #quit_msg = "Are you sure you want to exit the program?"
+            #reply = QtWidgets.QMessageBox.question(self, 'Message', 
+            #    quit_msg, QtWidgets.QMessageBox.Yes, QtWidgets.QMessageBox.No)
+            err_msg = "You need to set your data channels correctly.\n" + \
+                      "This should be a Python Numpy interpretable list\n" + \
+                      "of indices. For example 1 or 1:3 or 1:3 5\n\n" + \
+                      "valid GMR data channels fall between 1 and 8. Note that\n" +\
+                      "1:3 is not inclusive of 3 and is the same as 1,2 "
+            reply = QtWidgets.QMessageBox.critical(self, 'Error', 
+                err_msg) #, QtWidgets.QMessageBox.Yes, QtWidgets.QMessageBox.No)
+            return
+        #############################
+        # Reference Channels
+        # TODO make sure no overlap between data and ref channels
+        self.refChan = np.array( () )
+        if str(self.ui.refChanLineEdit.text()): # != "none":
+            try:
+                self.refChan = np.array(eval(str("np.r_["+self.ui.refChanLineEdit.text())+"]"))
+            except:
+                err_msg = "You need to set your reference channels correctly.\n" + \
+                      "This should be a Python Numpy interpretable list\n" + \
+                      "of indices. For example 1 or 1:3 or 1:3 5\n\n" + \
+                      "valid GMR data channels fall between 1 and 8. Note that\n" +\
+                      "1:3 is not inclusive of 3 and is the same as 1,2 "
+                QtWidgets.QMessageBox.critical(self, 'Error', err_msg) 
+                return
+
+        #####################################################
+        # Load data
+        
+        self.lock("loading MIDI 2 dataset")
+        self.procThread = thread.start_new_thread(self.RAWDataProc.loadMIDI2, \
+                (str(self.headerstr), self.procStacks, self.dataChan, self.refChan, \
+                 str(self.ui.FIDProcComboBox.currentText()), self.ui.mplwidget, \
+                1e-3 * self.ui.DeadTimeSpinBox.value( ), self.ui.plotImportCheckBox.isChecked() )) #, self)) 
+
+        self.YamlNode.Import["MIDI Header"] = self.headerstr
+        self.YamlNode.Import["opened"] = datetime.datetime.now().isoformat() 
+        self.YamlNode.Import["pulse Type"] = str(self.RAWDataProc.pulseType) 
+        self.YamlNode.Import["stacks"] = self.procStacks.tolist() 
+        self.YamlNode.Import["data channels"] = self.dataChan.tolist()  
+        self.YamlNode.Import["reference channels"] = self.refChan.tolist() 
+        self.YamlNode.Import["pulse records"] = str(self.ui.FIDProcComboBox.currentText())  
+        self.YamlNode.Import["instrument dead time"] = (1e-3 * self.ui.DeadTimeSpinBox.value( ))    
+
+        self.Log (  )     
+        
+        # enable META tab 
+        self.ui.METATab.setEnabled(1)
+        self.ui.siteBox.setEnabled(1)
+
+        # should be already done
+#        QtCore.QObject.connect(self.RAWDataProc, QtCore.SIGNAL("updateProgress(int)"), self.updateProgressBar)
+#        QtCore.QObject.connect(self.RAWDataProc, QtCore.SIGNAL("enableDSP()"), self.enableDSP)
+#        QtCore.QObject.connect(self.RAWDataProc, QtCore.SIGNAL("doneStatus()"), self.doneStatus)
+
+        #self.ui.ProcessedBox.setEnabled(True)
+        self.ui.lcdNumberFID1Length.setEnabled(1)
+        self.ui.lcdNumberFID2Length.setEnabled(1)
+        self.ui.lcdNumberResampFreq.setEnabled(1)
+        self.ui.lcdTotalDeadTime.setEnabled(1)
+
+        self.ui.lcdTotalDeadTime.display( self.ui.DeadTimeSpinBox.value( ) )
+        #self.ui.lcdTotalDeadTime.display( round(1e3*(self.RAWDataProc.DATADICT["Pulse 1"]["TIMES"][0]-self.RAWDataProc.DATADICT["Pulse 1"]["PULSE_TIMES"][-1]), 3) )
+        
+        #self.ui.lcdNumberFID1Length.display(0)
+        #self.ui.lcdNumberFID2Length.display(0)
+        #self.ui.lcdNumberResampFreq.display( self.RAWDataProc.samp )
  
+        self.mpl_toolbar = NavigationToolbar2QT(self.ui.mplwidget, self.ui.mplwidget)
+        self.ui.mplwidget.draw()
+
     def loadRAW(self):
 
         #################################################
@@ -1305,8 +1491,8 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         # Header
         if self.RAWDataProc == None:
             err_msg = "You need to load a header first."
-            reply = QtGui.QMessageBox.critical(self, 'Error', 
-                err_msg) #, QtGui.QMessageBox.Yes, QtGui.QMessageBox.No)
+            reply = QtWidgets.QMessageBox.critical(self, 'Error', 
+                err_msg) #, QtWidgets.QMessageBox.Yes, QtWidgets.QMessageBox.No)
             return
         
         # Stacks 
@@ -1316,7 +1502,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
             err_msg = "You need to set your stacks correctly.\n" + \
                       "This should be a Python Numpy interpretable list\n" + \
                       "of stack indices. For example 1:24 or 1:4,8:24"
-            QtGui.QMessageBox.critical(self, 'Error', err_msg) 
+            QtWidgets.QMessageBox.critical(self, 'Error', err_msg) 
             return
 
         # Data Channels
@@ -1328,15 +1514,15 @@ class ApplicationWindow(QtWidgets.QMainWindow):
             #messageBox.critical(0,"Error","An error has occured !");
             #messageBox.setFixedSize(500,200);
             #quit_msg = "Are you sure you want to exit the program?"
-            #reply = QtGui.QMessageBox.question(self, 'Message', 
-            #    quit_msg, QtGui.QMessageBox.Yes, QtGui.QMessageBox.No)
+            #reply = QtWidgets.QMessageBox.question(self, 'Message', 
+            #    quit_msg, QtWidgets.QMessageBox.Yes, QtWidgets.QMessageBox.No)
             err_msg = "You need to set your data channels correctly.\n" + \
                       "This should be a Python Numpy interpretable list\n" + \
                       "of indices. For example 1 or 1:3 or 1:3 5\n\n" + \
                       "valid GMR data channels fall between 1 and 8. Note that\n" +\
                       "1:3 is not inclusive of 3 and is the same as 1,2 "
-            reply = QtGui.QMessageBox.critical(self, 'Error', 
-                err_msg) #, QtGui.QMessageBox.Yes, QtGui.QMessageBox.No)
+            reply = QtWidgets.QMessageBox.critical(self, 'Error', 
+                err_msg) #, QtWidgets.QMessageBox.Yes, QtWidgets.QMessageBox.No)
             return
         #############################
         # Reference Channels
@@ -1351,13 +1537,13 @@ class ApplicationWindow(QtWidgets.QMainWindow):
                       "of indices. For example 1 or 1:3 or 1:3 5\n\n" + \
                       "valid GMR data channels fall between 1 and 8. Note that\n" +\
                       "1:3 is not inclusive of 3 and is the same as 1,2 "
-                QtGui.QMessageBox.critical(self, 'Error', err_msg) 
+                QtWidgets.QMessageBox.critical(self, 'Error', err_msg) 
                 return
 
         #####################################################
         # Load data
-
         self.lock("loading RAW GMR dataset")
+        
 
         if self.RAWDataProc.pulseType == "FID":
             self.procThread = thread.start_new_thread(self.RAWDataProc.loadFIDData, \
@@ -1466,6 +1652,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         # Adaptive filtering 
         self.ui.adaptBox.setEnabled(True)
         self.ui.adaptBox.setChecked(True)
+        self.ui.plotRLS.setEnabled(True)
         
         # FD Adaptive filtering 
         self.ui.adaptFDBox.setEnabled(True)
@@ -1648,6 +1835,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
                 self.ui.adaptTruncateSpinBox.value(), \
                 self.ui.adaptMuSpinBox.value(), \
                 str(self.ui.PCAComboBox.currentText()), \
+                self.ui.plotRLS.isChecked(), \
                 self.ui.mplwidget))
 
     def sumDataChans(self): 
@@ -1661,7 +1849,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
 
         self.dataChan = [self.dataChan[0]]
         self.ui.sumDataBox.setEnabled(False)    
-        thread.start_new_thread( self.RAWDataProc.sumData, ( self.ui.mplwidget, 7 ) )
+        thread.start_new_thread( self.RAWDataProc.sumData, ( self.ui.mplwidget, self.ui.sumType.currentText(), self.ui.sumAll.isChecked() ) )
  
     def adaptFilterFD(self):
         self.lock("FD noise cancellation filter")

akvo/gui/main.ui

@@ -1576,7 +1576,7 @@ background: dark grey;
                <bool>false</bool>
               </property>
               <property name="sizePolicy">
-               <sizepolicy hsizetype="Fixed" vsizetype="Expanding">
+               <sizepolicy hsizetype="Preferred" vsizetype="Expanding">
                 <horstretch>0</horstretch>
                 <verstretch>0</verstretch>
                </sizepolicy>
@@ -1608,7 +1608,7 @@ background: dark grey;
                 </widget>
                </item>
                <item row="0" column="1">
-                <widget class="QComboBox" name="comboBox">
+                <widget class="QComboBox" name="sumType">
                  <item>
                   <property name="text">
                    <string>sum</string>
@@ -1646,6 +1646,22 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
+               <item row="1" column="0">
+                <widget class="QCheckBox" name="sumAll">
+                 <property name="sizePolicy">
+                  <sizepolicy hsizetype="Preferred" vsizetype="Fixed">
+                   <horstretch>0</horstretch>
+                   <verstretch>0</verstretch>
+                  </sizepolicy>
+                 </property>
+                 <property name="toolTip">
+                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;In case of more than two data channels, would you like to sum all of the channels together?&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+                 </property>
+                 <property name="text">
+                  <string>Sum All</string>
+                 </property>
+                </widget>
+               </item>
               </layout>
              </widget>
             </item>
@@ -1947,13 +1963,6 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
-               <item row="15" column="0">
-                <widget class="Line" name="line_7">
-                 <property name="orientation">
-                  <enum>Qt::Horizontal</enum>
-                 </property>
-                </widget>
-               </item>
                <item row="11" column="0">
                 <widget class="QLabel" name="label_27">
                  <property name="text">
@@ -2156,6 +2165,13 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
+               <item row="15" column="0">
+                <widget class="Line" name="line_7">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
               </layout>
              </widget>
             </item>
@@ -2195,6 +2211,69 @@ background: dark grey;
                <bool>true</bool>
               </property>
               <layout class="QGridLayout" name="gridLayout_12">
+               <item row="2" column="1">
+                <widget class="QDoubleSpinBox" name="adaptTruncateSpinBox">
+                 <property name="toolTip">
+                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;This filter is a time-domain filter that takes some time to get going. Time-domain filters do a better job compared to frequency-domain filters in the presence of non-stationary noise. &lt;/p&gt;&lt;p&gt;The filter is run backwards, so often the late times will not be cancelled as well. You may trim records off the back using this input. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+                 </property>
+                 <property name="whatsThis">
+                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;This filter is a time-domain filter that takes some time to get going. Time-domain filters do a better job compared to frequency-domain filters in the presence of non-stationary noise.  &lt;/p&gt;&lt;p&gt;The filter is run backwards, so often the late times will not be cancelled as well. You may trim records off the back using this input. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+                 </property>
+                 <property name="maximum">
+                  <double>1000.000000000000000</double>
+                 </property>
+                 <property name="value">
+                  <double>800.000000000000000</double>
+                 </property>
+                </widget>
+               </item>
+               <item row="1" column="1">
+                <widget class="QDoubleSpinBox" name="adaptLambdaSpinBox">
+                 <property name="toolTip">
+                  <string>Forgetting factor, how quickly does the filter adapt.</string>
+                 </property>
+                 <property name="decimals">
+                  <number>3</number>
+                 </property>
+                 <property name="minimum">
+                  <double>0.950000000000000</double>
+                 </property>
+                 <property name="maximum">
+                  <double>1.000000000000000</double>
+                 </property>
+                 <property name="value">
+                  <double>0.990000000000000</double>
+                 </property>
+                </widget>
+               </item>
+               <item row="2" column="0">
+                <widget class="QLabel" name="label_46">
+                 <property name="text">
+                  <string>Truncate [ms]</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="0" column="2">
+                <widget class="QLabel" name="label_51">
+                 <property name="text">
+                  <string>Mu</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="1" column="0">
+                <widget class="QLabel" name="label_44">
+                 <property name="text">
+                  <string>Forgetting factor (λ)</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="1" column="2">
+                <widget class="QLabel" name="label_52">
+                 <property name="text">
+                  <string>PCA on ref</string>
+                 </property>
+                </widget>
+               </item>
                <item row="0" column="0">
                 <widget class="QLabel" name="label_43">
                  <property name="text">
@@ -2215,10 +2294,19 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
-               <item row="0" column="2">
-                <widget class="QLabel" name="label_51">
+               <item row="4" column="3">
+                <widget class="QPushButton" name="adaptGO">
+                 <property name="styleSheet">
+                  <string notr="true">#adaptGO {
+    background: green;
+}
+
+#adaptGO:disabled{
+    background: black;
+}</string>
+                 </property>
                  <property name="text">
-                  <string>Mu</string>
+                  <string>GO</string>
                  </property>
                 </widget>
                </item>
@@ -2241,39 +2329,6 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
-               <item row="1" column="0">
-                <widget class="QLabel" name="label_44">
-                 <property name="text">
-                  <string>Forgetting factor (λ)</string>
-                 </property>
-                </widget>
-               </item>
-               <item row="1" column="1">
-                <widget class="QDoubleSpinBox" name="adaptLambdaSpinBox">
-                 <property name="toolTip">
-                  <string>Forgetting factor, how quickly does the filter adapt.</string>
-                 </property>
-                 <property name="decimals">
-                  <number>3</number>
-                 </property>
-                 <property name="minimum">
-                  <double>0.950000000000000</double>
-                 </property>
-                 <property name="maximum">
-                  <double>1.000000000000000</double>
-                 </property>
-                 <property name="value">
-                  <double>0.990000000000000</double>
-                 </property>
-                </widget>
-               </item>
-               <item row="1" column="2">
-                <widget class="QLabel" name="label_52">
-                 <property name="text">
-                  <string>PCA on ref</string>
-                 </property>
-                </widget>
-               </item>
                <item row="1" column="3">
                 <widget class="QComboBox" name="PCAComboBox">
                  <property name="toolTip">
@@ -2294,42 +2349,26 @@ background: dark grey;
                  </item>
                 </widget>
                </item>
-               <item row="2" column="0">
-                <widget class="QLabel" name="label_46">
-                 <property name="text">
-                  <string>Truncate [ms]</string>
+               <item row="3" column="0" colspan="4">
+                <widget class="Line" name="line_15">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
                  </property>
                 </widget>
                </item>
-               <item row="2" column="1">
-                <widget class="QDoubleSpinBox" name="adaptTruncateSpinBox">
-                 <property name="toolTip">
-                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;This filter is a time-domain filter that takes some time to get going. Time-domain filters do a better job compared to frequency-domain filters in the presence of non-stationary noise. &lt;/p&gt;&lt;p&gt;The filter is run backwards, so often the late times will not be cancelled as well. You may trim records off the back using this input. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
-                 </property>
-                 <property name="whatsThis">
-                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;This filter is a time-domain filter that takes some time to get going. Time-domain filters do a better job compared to frequency-domain filters in the presence of non-stationary noise.  &lt;/p&gt;&lt;p&gt;The filter is run backwards, so often the late times will not be cancelled as well. You may trim records off the back using this input. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
-                 </property>
-                 <property name="maximum">
-                  <double>1000.000000000000000</double>
+               <item row="4" column="1">
+                <widget class="QCheckBox" name="plotRLS">
+                 <property name="enabled">
+                  <bool>false</bool>
                  </property>
-                 <property name="value">
-                  <double>800.000000000000000</double>
+                 <property name="text">
+                  <string>Plot</string>
                  </property>
-                </widget>
-               </item>
-               <item row="2" column="3">
-                <widget class="QPushButton" name="adaptGO">
-                 <property name="styleSheet">
-                  <string notr="true">#adaptGO {
-    background: green;
-}
-
-#adaptGO:disabled{
-    background: black;
-}</string>
+                 <property name="checkable">
+                  <bool>true</bool>
                  </property>
-                 <property name="text">
-                  <string>GO</string>
+                 <property name="checked">
+                  <bool>true</bool>
                  </property>
                 </widget>
                </item>
@@ -4143,8 +4182,8 @@ background: dark grey;
               <rect>
                <x>0</x>
                <y>0</y>
-               <width>100</width>
-               <height>30</height>
+               <width>96</width>
+               <height>26</height>
               </rect>
              </property>
              <attribute name="label">
@@ -4525,6 +4564,7 @@ p, li { white-space: pre-wrap; }
      <string>File</string>
     </property>
     <addaction name="actionOpen_GMR"/>
+    <addaction name="actionLoad_MIDI"/>
     <addaction name="separator"/>
     <addaction name="actionOpen_Preprocessed_Dataset"/>
     <addaction name="separator"/>
@@ -4614,6 +4654,11 @@ p, li { white-space: pre-wrap; }
     <string>Inversion</string>
    </property>
   </action>
+  <action name="actionLoad_MIDI">
+   <property name="text">
+    <string>Load MIDI 2 Data </string>
+   </property>
+  </action>
  </widget>
  <customwidgets>
   <customwidget>

akvo/tressel/calcAkvoKernel--save.py

@@ -0,0 +1,138 @@
+
+import os, sys
+import numpy as np
+from ruamel import yaml
+
+import pyLemma.LemmaCore as lc 
+import pyLemma.Merlin as mrln 
+import pyLemma.FDEM1D as em1d 
+
+import numpy as np
+
+#import matplotlib.pyplot as plt 
+#import seaborn as sns
+#sns.set(style="ticks")
+#import cmocean 
+#from SEGPlot import *
+#from matplotlib.ticker import FormatStrFormatter
+#import matplotlib.ticker as plticker
+
+
+# Converts Lemma/Merlin/Akvo serialized Eigen arrays into numpy ones for use by Python 
+class VectorXr(yaml.YAMLObject):
+    """
+    Converts Lemma/Merlin/Akvo serialized Eigen arrays into numpy ones for use by Python 
+    """
+    yaml_tag = u'VectorXr'
+    def __init__(self, array):
+        self.size = np.shape(array)[0]
+        self.data = array.tolist()
+    def __repr__(self):
+        # Converts to numpy array on import 
+        return "np.array(%r)" % (self.data)
+
+class AkvoData(yaml.YAMLObject):
+    """
+    Reads an Akvo serialized dataset into a standard python dictionary 
+    """
+    yaml_tag = u'AkvoData'
+    def __init__(self, array):
+        pass
+        #self.size = np.shape(array)[0]
+        #self.Imp = array.tolist()
+    def __repr__(self):
+        # Converts to a dictionary with Eigen vectors represented as Numpy arrays 
+        return self
+
+def loadAkvoData(fnamein):
+    """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be 
+        corrected in future kernel calculations. The current was reported but not the pulse length. 
+    """
+    fname = (os.path.splitext(fnamein)[0])
+    with open(fnamein, 'r') as stream:
+        try:
+            AKVO = (yaml.load(stream, Loader=yaml.Loader))
+        except yaml.YAMLError as exc:
+            print(exc)
+    return AKVO 
+
+
+def main():
+
+    if len(sys.argv) < 3:
+        print ("usage  python calcAkvoKernel.py   AkvoDataset.yaml  Coil1.yaml  " )
+        exit()
+
+    AKVO = loadAkvoData(sys.argv[1])
+
+    B_inc = AKVO.META["B_0"]["inc"]  
+    B_dec = AKVO.META["B_0"]["dec"]  
+    B0 = AKVO.META["B_0"]["intensity"]  
+
+    gamma = 2.67518e8
+    fT = AKVO.transFreq
+    #B0 = (fL*2.*np.pi) /gamma * 1e9
+ 
+    Coil1 = em1d.PolygonalWireAntenna.DeSerialize( sys.argv[2] )
+    Coil1.SetNumberOfFrequencies(1)
+    Coil1.SetFrequency(0, fT) 
+    Coil1.SetCurrent(1.)
+
+    lmod = em1d.LayeredEarthEM() 
+    lmod.SetNumberOfLayers(4)
+    lmod.SetLayerThickness([15.49, 28.18])
+    lmod.SetLayerConductivity([0.0, 1./16.91, 1./24.06, 1./33.23])
+
+    lmod.SetMagneticFieldIncDecMag( B_inc, B_dec, B0, lc.NANOTESLA )
+    
+    exit()
+
+    Kern = mrln.KernelV0()
+    Kern.PushCoil( "Coil 1", Coil1 )
+    Kern.SetLayeredEarthEM( lmod );
+    Kern.SetIntegrationSize( (200,200,200) )
+    Kern.SetIntegrationOrigin( (0,0,0) )
+    Kern.SetTolerance( 1e-9 )
+    Kern.SetMinLevel( 3 )
+    Kern.SetHankelTransformType( lc.FHTKEY201 )
+    Kern.AlignWithAkvoDataset( sys.argv[1] )
+
+    thick = np.geomspace(.5, 10,num=40)
+    iface = np.cumsum(thick)
+    Kern.SetDepthLayerInterfaces(iface)
+    #Kern.SetDepthLayerInterfaces(np.geomspace(1, 110, num=40))
+    #Kern.SetDepthLayerInterfaces(np.linspace(1, 110, num=50))
+    #Kern.SetDepthLayerInterfaces(np.geomspace(1, 110, num=40))
+ 
+    # autAkvoDataNode = YAML::LoadFile(argv[4]);
+    # Kern->AlignWithAkvoDataset( AkvoDataNode );
+
+    #Kern.SetPulseDuration(0.040)
+    #Kern.SetPulseCurrent(  [1.6108818092452406, 1.7549935078885168, 1.7666319459646016, 1.9270787752430283,
+    #    1.9455431806179229, 2.111931346726564, 2.1466747256211747, 2.3218217392379588,
+    #    2.358359967649008, 2.5495654202189058, 2.5957289164577992, 2.8168532605800802,
+    #    2.85505242699187, 3.1599429539069774, 3.2263673040205068, 3.6334182368296544,
+    #    3.827985200119751, 4.265671313014058, 4.582237014873297, 5.116839616183394,
+    #    5.515173073160611, 6.143620383280934, 6.647972282096122, 7.392577402979211,
+    #    8.020737177449933, 8.904435233295793, 9.701975105606063, 10.74508217792577,
+    #    11.743887525923592, 12.995985956061467, 14.23723766879807, 15.733870137824457,
+    #    17.290155933625808, 19.07016662950366, 21.013341340455703, 23.134181634845618,
+    #    25.570925414182238, 28.100862178905476, 31.13848909847073, 34.16791099558486,
+    #    37.95775984680512, 41.589619321873165, 46.327607251605286, 50.667786337299205,
+    #    56.60102493062895, 61.81174065797068, 69.23049946198458, 75.47409803238031,
+    #    84.71658869065816, 92.1855007134236, 103.77129947551164, 112.84577430578537,
+    #    127.55127257092909, 138.70199812969176, 157.7443764728878, 171.39653462998626]
+    #)
+
+    Kern.CalculateK0( ["Coil 1"], ["Coil 1"], False )
+
+    yml = open('akvoK3-' + str(Kern.GetTolerance()) + '.yaml', 'w')
+    print(Kern, file=yml)
+
+    K0 = Kern.GetKernel()
+    
+    #plt.matshow(np.abs(K0))
+    #plt.show()
+
+if __name__ == "__main__":
+    main()

akvo/tressel/decay-old.py

@@ -0,0 +1,729 @@
+import numpy, array #,rpy2
+from matplotlib import pyplot as plt
+import numpy as np
+from scipy.optimize import least_squares
+from rpy2.robjects.packages import importr
+
+import rpy2.robjects as robjects
+import rpy2.robjects.numpy2ri
+
+#import notch
+from numpy.fft import fft, fftfreq
+
+# We know/can calculate frequency peak, use this to guess where picks will be.
+# maybe have a sliding window that reports peak values.
+def peakPicker(data, omega, dt):
+
+    # compute window based on omega and dt
+    # make sure you are not aliased, grab every other peak
+    window = (2*numpy.pi) / (omega*dt)
+
+    data = numpy.array(data) 
+    peaks = []
+    troughs = []
+    times = []
+    times2 = []
+    indices = []
+    ws = 0
+    we = window
+    ii = 0
+    for i in range((int)(len(data)/window)):
+        
+        # initially was just returning this I think avg is better
+        #times.append( (ws + numpy.abs(data[ws:we]).argmax()) * dt )
+    
+        peaks.append(numpy.max(data[ws:we]))
+        times.append( (ws + data[ws:we].argmax()) * dt )
+        indices.append( ii + data[ws:we].argmax() )        
+
+        troughs.append(numpy.min(data[ws:we]))
+        times2.append( (ws + (data[ws:we]).argmin()) * dt )
+        indices.append( ii + data[ws:we].argmin() )        
+
+        ws += window
+        we += window
+        ii += (int)(we-ws)
+    
+    #return numpy.array(peaks), numpy.array(times)
+    
+    # Averaging peaks does a good job of removing bias in noise
+    return (numpy.array(peaks)-numpy.array(troughs))/2., \
+        (numpy.array(times)+numpy.array(times2))/2., \
+        indices           
+
+
+#################################################
+# Regress for T2 using rpy2 interface
+def regressCurve(peaks,times,sigma2=1,intercept=True):
+
+    # TODO, if regression fails, it might be because there is no exponential
+    # term, maybe do a second regression then on a linear model. 
+    b1  = 0                  # Bias
+    b2  = 0                  # Linear 
+    rT2 = 0.3                # T2 regressed
+    r   = robjects.r         
+
+    # Variable shared between R and Python
+    robjects.globalenv['b1'] = b1
+    robjects.globalenv['b2'] = b2
+    robjects.globalenv['rT2'] = rT2
+    robjects.globalenv['sigma2'] = sigma2
+    value = robjects.FloatVector(peaks)
+    times = robjects.FloatVector(numpy.array(times))
+    
+#    my_weights = robjects.RVector(value/sigma2)
+#    robjects.globalenv['my_weigts'] = my_weights
+
+#    if sigma2 != 0:
+#        print "weighting"
+#        tw = numpy.array(peaks)/sigma2 
+#        my_weights = robjects.RVector( tw/numpy.max(tw) )
+#    else:
+#        my_weights = robjects.RVector(numpy.ones(len(peaks))) 
+
+#    robjects.globalenv['my_weights'] = my_weights
+    
+    if (intercept):
+        my_list = robjects.r('list(b1=50, b2=1e2, rT2=0.03)')
+        my_lower = robjects.r('list(b1=0, b2=0, rT2=.005)')
+        my_upper = robjects.r('list(b1=20000, b2=2000, rT2=.700)')
+    else:
+        my_list = robjects.r('list(b2=1e2, rT2=0.3)')
+        my_lower = robjects.r('list(b2=0, rT2=.005)')
+        my_upper = robjects.r('list(b2=2000, rT2=.700)')
+
+    my_cont = robjects.r('nls.control(maxiter=1000, warnOnly=TRUE, printEval=FALSE)')
+
+    
+    if (intercept):
+        #fmla = robjects.RFormula('value ~ b1 + exp(-times/rT2)')
+        fmla = robjects.Formula('value ~ b1 + b2*exp(-times/rT2)')
+        #fmla = robjects.RFormula('value ~ b1 + b2*times + exp(-times/rT2)')
+    else:
+        fmla = robjects.Formula('value ~ b2*exp(-times/rT2)')
+
+    env = fmla.getenvironment()
+    env['value'] = value
+    env['times'] = times
+    
+    # ugly, but I get errors with everything else I've tried
+    my_weights = robjects.r('rep(1,length(value))')
+    for ii in range(len(my_weights)):
+        my_weights[ii] *= peaks[ii]/sigma2
+    Error = False
+    #fit = robjects.r.nls(fmla,start=my_list,control=my_cont,weights=my_weights)
+    if (sigma2 != 1):
+        print("SIGMA 2")
+        #fit = robjects.r.tryCatch(robjects.r.suppressWarnings(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port", \
+        #                     weights=my_weights)), 'silent=TRUE')
+        fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont))#, \
+                            # weights=my_weights))
+    else:
+        try:
+            fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port"))#,lower=my_lower,upper=my_upper))
+        except:
+            print("regression issue pass")
+            Error = True
+    # If failure fall back on zero regression values   
+    if not Error:
+        #Error = fit[3][0]
+        report =  r.summary(fit)
+    b1 = 0
+    b2 = 0 
+    rT2 = 1
+    if (intercept):
+        if not Error:
+            b1  =  r['$'](report,'par')[0]
+            b2  =  r['$'](report,'par')[1]
+            rT2 =  r['$'](report,'par')[2]
+            #print  report
+            #print  r['$'](report,'convergence')
+            #print  r['convergence'] #(report,'convergence')
+            #print  r['$'](report,'par')[13]
+            #print  r['$'](report,'par')[14]
+        else:
+            print("ERROR DETECTED, regressed values set to default")
+            b1 = 1e1
+            b2 = 1e-2
+            rT2 = 1e-2
+            #print r['$'](report,'par')[0]
+            #print r['$'](report,'par')[1]
+            #print r['$'](report,'par')[2]
+        return [b1,b2,rT2] 
+    else:
+        if not Error:
+            rT2 =  r['$'](report,'par')[1]
+            b2  =  r['$'](report,'par')[0]
+        else:
+            print("ERROR DETECTED, regressed values set to default")
+        return [b2, rT2] 
+
+#################################################
+# Regress for T2 using rpy2 interface
+def regressCurve2(peaks,times,sigma2=[None],intercept=True):
+
+    if sigma2[0] != None:
+        my_weights = robjects.FloatVector( sigma2 )
+
+    # TODO, if regression fails, it might be because there is no exponential
+    # term, maybe do a second regression then on a linear model. 
+    b1  = 0                  # Bias
+    b2  = 0                  # Linear 
+    bb2  = 0                 # Linear 
+    rT2 = 0.3                # T2 regressed
+    rrT2 = 1.3               # T2 regressed
+    r   = robjects.r         
+
+    # Variable shared between R and Python
+    robjects.globalenv['b1'] = b1
+    robjects.globalenv['b2'] = b2
+    robjects.globalenv['rT2'] = rT2
+    
+    robjects.globalenv['bb2'] = b2
+    robjects.globalenv['rrT2'] = rT2
+    
+    #robjects.globalenv['sigma2'] = sigma2
+    value = robjects.FloatVector(peaks)
+    times = robjects.FloatVector(numpy.array(times))
+    
+    
+    if (intercept):
+        my_list = robjects.r('list(b1=.50, b2=1e2, rT2=0.03, bb2=1e1, rrT2=1.3)')
+        my_lower = robjects.r('list(b1=0, b2=0, rT2=.005, bb2=0, rrT2=.005 )')
+        my_upper = robjects.r('list(b1=2000, b2=2000, rT2=.700, bb2=2000, rrT2=1.3 )')
+    else:
+        my_list  = robjects.r('list(b2=.5, rT2=0.3,  bb2=.5, rrT2=1.3)')
+        my_lower = robjects.r('list(b2=0,  rT2=.005, bb2=0,  rrT2=.005)')
+        my_upper = robjects.r('list(b2=1,  rT2=2.6,    bb2=1,  rrT2=2.6)')
+
+    my_cont = robjects.r('nls.control(maxiter=1000, warnOnly=TRUE, printEval=FALSE)')
+
+    
+    if (intercept):
+        #fmla = robjects.RFormula('value ~ b1 + exp(-times/rT2)')
+        fmla = robjects.Formula('value ~ b1 + b2*exp(-times/rT2) + bb2*exp(-times/rrT2)')
+        #fmla = robjects.RFormula('value ~ b1 + b2*times + exp(-times/rT2)')
+    else:
+        fmla = robjects.Formula('value ~ b2*exp(-times/rT2) + bb2*exp(-times/rrT2)')
+
+    env = fmla.getenvironment()
+    env['value'] = value
+    env['times'] = times
+    
+    # ugly, but I get errors with everything else I've tried
+    Error = False
+    #fit = robjects.r.nls(fmla,start=my_list,control=my_cont,weights=my_weights)
+    if (sigma2[0] != None):
+        #print("SIGMA 2")
+        #fit = robjects.r.tryCatch(robjects.r.suppressWarnings(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port", \
+        #                     weights=my_weights)), 'silent=TRUE')
+        fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm='port',weights=my_weights,lower=my_lower,upper=my_upper))#, \
+                            # weights=my_weights))
+    else:
+        try:
+            fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port"))#,lower=my_lower,upper=my_upper))
+        except:
+            print("regression issue pass")
+            Error = True
+    # If failure fall back on zero regression values   
+    if not Error:
+        #Error = fit[3][0]
+        report =  r.summary(fit)
+    b1 = 0
+    b2 = 0 
+    rT2 = 1
+    if (intercept):
+        if not Error:
+            b1  =  r['$'](report,'par')[0]
+            b2  =  r['$'](report,'par')[1]
+            rT2 =  r['$'](report,'par')[2]
+            #print  report
+            #print  r['$'](report,'convergence')
+            #print  r['convergence'] #(report,'convergence')
+            #print  r['$'](report,'par')[13]
+            #print  r['$'](report,'par')[14]
+        else:
+            print("ERROR DETECTED, regressed values set to default")
+            b1 = 1e1
+            b2 = 1e-2
+            rT2 = 1e-2
+            #print r['$'](report,'par')[0]
+            #print r['$'](report,'par')[1]
+            #print r['$'](report,'par')[2]
+        return [b1,b2,rT2, bb2, rrT2] 
+    else:
+        if not Error:
+            rT2 =  r['$'](report,'par')[1]
+            b2  =  r['$'](report,'par')[0]
+            rrT2 =  r['$'](report,'par')[3]
+            bb2  =  r['$'](report,'par')[2]
+        else:
+            print("ERROR DETECTED, regressed values set to default")
+        return [b2, rT2, bb2, rrT2] 
+
+def fun(x, t, y):
+    """ Cost function for regression, single exponential, no DC term 
+        x[0] = A0
+        x[1] = zeta 
+        x[2] = df
+        x[3] = T2
+    """
+    # concatenated real and imaginary parts  
+    pre =  np.concatenate((-x[0]*np.sin(2.*np.pi*x[2]*t + x[1])*np.exp(-t/x[3]), \
+                           +x[0]*np.cos(2.*np.pi*x[2]*t + x[1])*np.exp(-t/x[3])))  
+    return y-pre
+
+def fun2(x, t, y):
+    """ Cost function for regression, single exponential, no DC term 
+        x[0] = A0
+        x[1] = zeta 
+        x[2] = T2
+    """
+    # concatenated real and imaginary parts  
+    pre =  np.concatenate((x[0]*np.cos(x[1])*np.exp(-t/x[2]), \
+                       -1.*x[0]*np.sin(x[1])*np.exp(-t/x[2])))  
+    return y-pre
+
+
+def quadratureDetect2(X, Y, tt, x0="None"): 
+    """ Pure python quadrature detection using Scipy.  
+        X = real part of NMR signal 
+        Y = imaginary component of NMR signal 
+        tt = time 
+    """
+    print("Pure Python Quad Det", "TODO look at loss functions and method")
+    # Loss functions, linear, soft_l1, huber, cauchy, arctan 
+    # df
+    loss = 'cauchy'  #  'soft_l1'
+    method = 'trf'   # trf, dogbox, lm 
+    if x0=="None":
+        x0 = np.array( [1., 0., 0., .2] ) # A0, zeta, df, T2 
+        res_lsq = least_squares(fun, x0, args=(tt, np.concatenate((X, Y))), loss=loss, f_scale=1.0,\
+            bounds=( [1., -np.pi, -5, .005] , [1000., np.pi, 5, .800] ),
+            method=method 
+            )
+        x = res_lsq.x 
+        print ("df", x[0], x[1], x[2], x[3])
+    else:
+        res_lsq = least_squares(fun, x0, args=(tt, np.concatenate((X, Y))), loss=loss, f_scale=1.0,\
+            bounds=( [1., -np.pi, -5, .005] , [1000., np.pi, 5, .800] ),
+            method=method 
+            )
+
+        #bounds=( [0., 0, -20, .0] , [1., np.pi, 20, .6] ))
+
+    x = res_lsq.x 
+    return res_lsq.success, x[0], x[2], x[1], x[3]
+    
+    # no df
+    #x = np.array( [1., 0., 0.2] )
+    #res_lsq = least_squares(fun2, x, args=(tt, np.concatenate((X, Y))), loss='soft_l1', f_scale=0.1)
+    #x = res_lsq.x 
+    #return conv, E0,df,phi,T2
+    #return res_lsq.success, x[0], 0, x[1], x[2]
+
+def quadratureDetect(X, Y, tt, CorrectFreq=False, BiExp=False, CorrectDC=False):
+ 
+    r   = robjects.r        
+
+    if CorrectDC:
+        robjects.r(''' 
+             Xc1 <- function(E01, df, tt, phi, T2_1, DC) {
+	                DC + E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1)
+            }
+    
+            Yc1 <- function(E01, df, tt, phi, T2_1, DC) {
+	                DC - E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1)
+            } 
+            ''')
+    else:   
+        robjects.r(''' 
+             Xc1 <- function(E01, df, tt, phi, T2_1) {
+	                E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1)
+            }
+    
+            Yc1 <- function(E01, df, tt, phi, T2_1) {
+	                -E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1)
+            } 
+            ''')
+
+    # bi-exponential 
+    if CorrectDC:
+        robjects.r(''' 
+             Xc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2, DC) {
+	               DC + E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1) + 
+	                DC + E02*cos(2*pi*df*tt + phi) * exp(-tt/T2_2)
+            }
+
+            Yc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2, DC) {
+	                DC - E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1) + 
+	                DC - E02*sin(2*pi*df*tt + phi) * exp(-tt/T2_2)
+            } 
+            ''')
+    else:   
+        robjects.r(''' 
+             Xc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2) {
+	               E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1) + 
+	               E02*cos(2*pi*df*tt + phi) * exp(-tt/T2_2)
+            }
+
+            Yc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2) {
+	                -E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1) + 
+	                -E02*sin(2*pi*df*tt + phi) * exp(-tt/T2_2)
+            } 
+            ''')
+
+    # Make 0 vector 
+    Zero = robjects.FloatVector(numpy.zeros(len(X)))
+    
+    # Fitted Parameters
+    E01 = 0.
+    E02 = 0.
+    df = 0.
+    phi = 0.
+    T2_1 = 0.
+    T2_2 = 0.
+    DC = 0.
+    robjects.globalenv['DC'] = DC
+    robjects.globalenv['E01'] = E01
+    robjects.globalenv['E02'] = E02
+    robjects.globalenv['df'] = df
+    robjects.globalenv['phi'] = phi
+    robjects.globalenv['T2_1'] = T2_1
+    robjects.globalenv['T2_2'] = T2_2
+    XY = robjects.FloatVector(numpy.concatenate((X,Y)))
+    
+    # Arrays
+    tt = robjects.FloatVector(numpy.array(tt))
+    X = robjects.FloatVector(numpy.array(X))
+    Y = robjects.FloatVector(numpy.array(Y))
+    Zero = robjects.FloatVector(numpy.array(Zero))
+
+    
+
+    if BiExp:
+        if CorrectDC:
+            fmla = robjects.Formula('XY ~ c(Xc2( E01, E02, df, tt, phi, T2_1, T2_2, DC ), Yc2( E01, E02, df, tt, phi, T2_1, T2_2, DC ))')
+            if CorrectFreq:    
+                start = robjects.r('list(E01=.100, E02=.01,   df=0,    phi=0.    ,  T2_1=.100, T2_2=.01, DC=0.0)')
+                lower = robjects.r('list(E01=1e-6, E02=1e-6,  df=-50,  phi=-3.14 ,  T2_1=.001, T2_2=.001, DC=0.0)')
+                upper = robjects.r('list(E01=1.00, E02=1.0,   df=50,   phi=3.14  ,  T2_1=.800, T2_2=.8, DC=0.5)')
+            else:
+                start = robjects.r('list(E01=.100, E02=.01,   phi=0.9   ,  T2_1=.100, T2_2=.01,  DC=0.0)')
+                lower = robjects.r('list(E01=1e-6, E02=1e-6,  phi=-3.14 ,  T2_1=.001, T2_2=.001, DC=0.0)')
+                upper = robjects.r('list(E01=1.00, E02=1.0,   phi=3.14  ,  T2_1=.800, T2_2=.8,   DC=0.5)')
+        else:
+            fmla = robjects.Formula('XY ~ c(Xc2( E01, E02, df, tt, phi, T2_1, T2_2 ), Yc2( E01, E02, df, tt, phi, T2_1, T2_2))')
+            if CorrectFreq:    
+                start = robjects.r('list(E01=.100, E02=.01,   df=0,    phi=0.    ,  T2_1=.100, T2_2=.01)')
+                lower = robjects.r('list(E01=1e-6, E02=1e-6,  df=-50,  phi=-3.14 ,  T2_1=.001, T2_2=.001)')
+                upper = robjects.r('list(E01=1.00, E02=1.0,   df=50,   phi=3.14  ,  T2_1=.800, T2_2=.8)')
+            else:
+                start = robjects.r('list(E01=.100, E02=.01,   phi=0.9   ,  T2_1=.100, T2_2=.01)')
+                lower = robjects.r('list(E01=1e-6, E02=1e-6,  phi=-3.14 ,  T2_1=.001, T2_2=.001)')
+                upper = robjects.r('list(E01=1.00, E02=1.0,   phi=3.14  ,  T2_1=.800, T2_2=.8)')
+    else: 
+        if CorrectDC:
+            fmla = robjects.Formula('XY ~ c(Xc1( E01, df, tt, phi, T2_1, DC), Yc1( E01, df, tt, phi, T2_1,DC))')
+            if CorrectFreq:    
+                start = robjects.r('list(E01=.100, df=0   , phi=0.   , T2_1=.100, DC=0.0)')
+                lower = robjects.r('list(E01=1e-6, df=-50., phi=-3.14, T2_1=.001, DC=0.0)')
+                upper = robjects.r('list(E01=1.00, df=50. , phi=3.14 , T2_1=.800, DC=0.5)')
+            else:
+                start = robjects.r('list(E01=.100, phi= 0.  , T2_1=.100, DC=0.0)')
+                lower = robjects.r('list(E01=1e-6, phi=-3.13, T2_1=.001, DC=0.0)')
+                upper = robjects.r('list(E01=1.00, phi= 3.13, T2_1=.800, DC=0.5)')
+        else:
+            fmla = robjects.Formula('XY ~ c(Xc1( E01, df, tt, phi, T2_1), Yc1( E01, df, tt, phi, T2_1))')
+            if CorrectFreq:    
+                start = robjects.r('list(E01=.100, df=0     , phi=0.   ,  T2_1=.100)')
+                lower = robjects.r('list(E01=1e-6, df=-50. , phi=-3.14 ,  T2_1=.001)')
+                upper = robjects.r('list(E01=1.00, df=50.  , phi=3.14  ,  T2_1=.800)')
+            else:
+                start = robjects.r('list(E01=.100, phi= 0.  , T2_1=.100)')
+                lower = robjects.r('list(E01=1e-6, phi=-3.13, T2_1=.001)')
+                upper = robjects.r('list(E01=1.00, phi= 3.13, T2_1=.800)')
+
+    env = fmla.getenvironment()
+    env['Zero'] = Zero
+    env['X'] = X
+    env['Y'] = Y
+    env['XY'] = XY 
+    env['tt'] = tt
+
+    cont = robjects.r('nls.control(maxiter=10000, warnOnly=TRUE, printEval=FALSE)')
+    
+    fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=start, control=cont, lower=lower, upper=upper, algorithm='port')) #, \
+    #fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=start, control=cont)) #, \
+    report =  r.summary(fit)
+
+    conv = r['$'](fit,'convergence')[0]
+    #if conv:
+    #    print (report)
+    #    print ("conv", conv)
+    print ("Conv",  r['$'](fit,'convergence'))  # T2
+    print (report)
+    
+    if BiExp:
+        if CorrectFreq:    
+            E0   =  r['$'](report,'par')[0]   # E01
+            E0  +=  r['$'](report,'par')[1]   # E02
+            df  =  r['$'](report,'par')[2]   # offset
+            phi =  r['$'](report,'par')[3]   # phase 
+            T2  =  r['$'](report,'par')[4]   # T2
+        else:
+            E0   =  r['$'](report,'par')[0]   # E01
+            E0  +=  r['$'](report,'par')[1]   # E02
+            phi =  r['$'](report,'par')[2]   # phase 
+            T2  =  r['$'](report,'par')[3]   # T2
+    else:
+        if CorrectFreq:    
+            E0   =  r['$'](report,'par')[0]   # E01
+            df  =  r['$'](report,'par')[1]   # offset
+            phi =  r['$'](report,'par')[2]   # phase 
+            T2  =  r['$'](report,'par')[3]   # T2
+        else:
+            E0   =  r['$'](report,'par')[0]   # E01
+            phi =  r['$'](report,'par')[1]   # phase 
+            T2  =  r['$'](report,'par')[2]   # T2
+    #phi = 0.907655876627
+    #phi = 0
+    #print ("df", df)# = 0
+    return conv, E0,df,phi,T2
+    
+
+#################################################
+# Regress for T2 using rpy2 interface
+def regressSpec(w, wL, X): #,sigma2=1,intercept=True):
+
+    # compute s
+    s = -1j*w
+
+    # TODO, if regression fails, it might be because there is no exponential
+    # term, maybe do a second regression then on a linear model. 
+    a   = 0                  # Linear 
+    rT2 = 0.1                # T2 regressed
+    r   = robjects.r         
+
+    # Variable shared between R and Python
+    robjects.globalenv['a'] = a
+    robjects.globalenv['rT2'] = rT2
+    robjects.globalenv['wL'] = wL
+    robjects.globalenv['nb'] = 0
+
+    s = robjects.ComplexVector(numpy.array(s))
+    XX = robjects.ComplexVector(X)
+    Xr = robjects.FloatVector(numpy.real(X))
+    Xi = robjects.FloatVector(numpy.imag(X))
+    Xa = robjects.FloatVector(numpy.abs(X))
+    Xri = robjects.FloatVector(numpy.concatenate((Xr,Xi)))
+    
+    #my_lower = robjects.r('list(a=.001, rT2=.001, nb=.0001)')
+    my_lower = robjects.r('list(a=.001, rT2=.001)')
+    #my_upper = robjects.r('list(a=1.5, rT2=.300, nb =100.)')
+    my_upper = robjects.r('list(a=1.5, rT2=.300)')
+     
+    #my_list = robjects.r('list(a=.2, rT2=0.03, nb=.1)')
+    my_list = robjects.r('list(a=.2, rT2=0.03)')
+    my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
+    
+    #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
+    ##fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
+    #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
+    #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
+    fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))') # complex
+ 
+    env = fmla.getenvironment()
+    env['s'] = s
+    env['Xr'] = Xr
+    env['Xa'] = Xa
+    env['Xi'] = Xi
+    env['Xri'] = Xri
+    env['XX'] = XX
+     
+    #fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
+    fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
+    report =  r.summary(fit)
+    #print report 
+    #print  r.warnings()
+ 
+    a  =  r['$'](report,'par')[0]
+    rT2 =  r['$'](report,'par')[1]
+    nb =  r['$'](report,'par')[2]
+    
+    return a, rT2, nb
+
+#################################################
+# Regress for T2 using rpy2 interface
+def regressSpecComplex(w, wL, X): #,sigma2=1,intercept=True):
+
+    # compute s
+    s = -1j*w
+
+    # TODO, if regression fails, it might be because there is no exponential
+    # term, maybe do a second regression then on a linear model. 
+    a   = 1                  # Linear 
+    rT2 = 0.1                # T2 regressed
+    r   = robjects.r         
+    phi2 = 0                 # phase
+    wL2 = wL
+
+    # Variable shared between R and Python
+    robjects.globalenv['a'] = a
+    robjects.globalenv['rT2'] = rT2
+    robjects.globalenv['wL'] = wL
+    robjects.globalenv['wL2'] = 0
+    robjects.globalenv['nb'] = 0
+    robjects.globalenv['phi2'] = phi2
+
+    s = robjects.ComplexVector(numpy.array(s))
+    XX = robjects.ComplexVector(X)
+    Xr = robjects.FloatVector(numpy.real(X))
+    Xi = robjects.FloatVector(numpy.imag(X))
+    Xa = robjects.FloatVector(numpy.abs(X))
+    Xri = robjects.FloatVector(numpy.concatenate((X.real,X.imag)))
+
+    robjects.r(''' 
+        source('kernel.r')
+    ''')   
+    #Kw = robjects.globalenv['Kwri']
+     
+    #print (numpy.shape(X))
+    
+    #my_lower = robjects.r('list(a=.001, rT2=.001, nb=.0001)')
+    #my_lower = robjects.r('list(a=.001, rT2=.001)') # Working
+    my_lower = robjects.r('list(a=.001, rT2=.001, phi2=-3.14, wL2=wL-5)')
+    #my_upper = robjects.r('list(a=1.5, rT2=.300, nb =100.)')
+    my_upper = robjects.r('list(a=3.5, rT2=.300, phi2=3.14, wL2=wL+5)')
+     
+    #my_list = robjects.r('list(a=.2, rT2=0.03, nb=.1)')
+    my_list = robjects.r('list(a=.2, rT2=0.03, phi2=0, wL2=wL)')
+    my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
+    
+    #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
+    #fmla = robjects.Formula('Xi   ~   Im(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 ))') # envelope
+    #fmla = robjects.Formula('Xri ~ c(Re(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )), Im(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )))') # envelope
+    
+    #fmlar = robjects.Formula('Xr ~ (Kwr(a, phi2, s, rT2, wL)) ') # envelope
+    #fmlai = robjects.Formula('Xi ~ (Kwi(a, phi2, s, rT2, wL)) ') # envelope
+    fmla = robjects.Formula('Xri ~ c(Kwr(a, phi2, s, rT2, wL2), Kwi(a, phi2, s, rT2, wL2) ) ') # envelope
+    #fmla = robjects.Formula('Xri ~ (Kwri(a, phi2, s, rT2, wL)) ') # envelope
+    
+    #fmla = robjects.Formula('Xa ~ (abs(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )))') # envelope
+    #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
+    #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
+    
+    #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
+    
+    #        self.Gw[iw, iT2] = ((np.sin(phi2) *  (alpha + 1j*self.w[iw]) + self.wL*np.cos(phi2)) / \
+    #                               (self.wL**2 + (alpha+1.j*self.w[iw])**2 ))
+    #        self.Gw[iw, iT2] = ds * self.sc*((np.sin(phi2)*( alpha + 1j*self.w[iw]) + self.wL*np.cos(phi2)) / \
+    #                               (self.wL**2 + (alpha+1.j*self.w[iw])**2 ))
+    
+    # Works Amplitude Only!
+    #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))') # complex
+ 
+    env = fmla.getenvironment()
+    env['s'] = s
+    env['Xr'] = Xr
+    env['Xa'] = Xa
+    env['Xi'] = Xi
+    env['Xri'] = Xri
+    env['XX'] = XX
+     
+    fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
+    #fitr = robjects.r.tryCatch(robjects.r.nls(fmlar, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
+    
+    #env = fmlai.getenvironment()
+    #fiti = robjects.r.tryCatch(robjects.r.nls(fmlai, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
+    
+    #reportr =  r.summary(fitr)
+    #reporti =  r.summary(fiti)
+    report =  r.summary(fit)
+    #print( report )
+    #exit()
+    #print( reportr )
+    #print( reporti  )
+    #exit()
+    #print  r.warnings()
+ 
+    #a   =  (r['$'](reportr,'par')[0] + r['$'](reporti,'par')[0]) / 2.
+    #rT2 =  (r['$'](reportr,'par')[1] + r['$'](reporti,'par')[1]) / 2.
+    #nb  =  (r['$'](reportr,'par')[2] + r['$'](reporti,'par')[2]) / 2.
+    a   =  r['$'](report,'par')[0] 
+    rT2 =  r['$'](report,'par')[1] 
+    nb  =  r['$'](report,'par')[2] #phi2 
+
+    print ("Python wL2", r['$'](report,'par')[3] )   
+    print ("Python zeta", r['$'](report,'par')[2] )   
+ 
+    return a, rT2, nb
+
+
+
+###################################################################
+###################################################################
+###################################################################
+if __name__ == "__main__":
+
+    dt    = .0001
+    T2    = .1
+    omega = 2000.*2*numpy.pi
+    phi   = .0
+    T     = 8.*T2
+    
+    t = numpy.arange(0, T, dt)
+
+    # Synthetic data, simple single decaying sinusoid 
+    # with a single decay parameter and gaussian noise added 
+    data = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) + numpy.random.normal(0,.05,len(t)) \
+                         + numpy.random.randint(-1,2,len(t))*numpy.random.exponential(.2,len(t)) 
+    cdata = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) #+ numpy.random.normal(0,.25,len(t))
+    #data = numpy.random.normal(0,.25,len(t))
+
+    sigma2 = numpy.std(data[::-len(data)/4])
+    #sigma2 = numpy.var(data[::-len(data)/4])
+    print("sigma2", sigma2)    
+    
+    [peaks,times,indices] = peakPicker(data, omega, dt)
+    
+    [b1,b2,rT2] = regressCurve(peaks,times)
+    print("rT2 nonweighted", rT2)
+    
+    [b1,b2,rT2] = regressCurve(peaks,times,sigma2)
+    print("rT2 weighted", rT2)
+
+    envelope   =  numpy.exp(-t/T2)
+    renvelope  =  numpy.exp(-t/rT2)
+
+    #outf = file('regress.txt','w')
+    #for i in range(len(times)):
+    #    outf.write(str(times[i]) + "   " +  str(peaks[i]) + "\n")  
+    #outf.close()
+
+    plt.plot(t,data, 'b')
+    plt.plot(t,cdata, 'g', linewidth=1)
+    plt.plot(t,envelope, color='violet', linewidth=4)
+    plt.plot(t,renvelope, 'r', linewidth=4)
+    plt.plot(times, numpy.array(peaks), 'bo', markersize=8, alpha=.25)
+    plt.legend(['noisy data','clean data','real envelope','regressed env','picks'])
+    plt.savefig("regression.pdf")
+
+
+    # FFT check
+    fourier = fft(data)
+    plt.figure()
+    freq = fftfreq(len(data), d=dt)
+    plt.plot(freq, (fourier.real))
+    
+    plt.show()
+
+    # TODO do a bunch in batch mode to see if T2 estimate is better with or without 
+    # weighting and which model is best.
+
+    # TODO try with real data
+
+    # TODO test filters (median, FFT, notch)
+
+    # It looks like weighting is good for relatively low sigma, but for noisy data
+    # it hurts us. Check

akvo/tressel/gateIntegrate4.py

@@ -0,0 +1,392 @@
+from __future__ import division 
+
+import matplotlib as mpl
+mpl.use('pdf')
+
+#from rasterize import rasterize_and_save
+
+import matplotlib.patches as mpatches
+from pwctime import pwcTime
+from logbarrier import * 
+from perlin import perlin
+from scipy import stats
+import cmocean
+import sys 
+import numpy as np 
+import seaborn as sns
+
+def bootstrapWindows(N, nboot, isum, adapt=False):
+    """ Bootstraps noise as a function of gate width
+        N = input noise signal 
+        nboot = number of boostrap windows to perform 
+        isum = length of windows (L_i)
+        adapt = reduce nboot as window size increases
+    """
+    nc = np.shape(N)[0]
+    Means = {}
+
+    if adapt:
+        Means = -9999*np.ones((len(isum), nboot//isum[0])) # dummy value
+        for ii, nwin in enumerate(isum):  
+            for iboot in range(nboot//isum[ii]):
+                cs = np.random.randint(0,nc-nwin)
+                Means[ii,iboot] = np.mean( N[cs:cs+nwin] )
+        Means = np.ma.masked_less(Means, -9995)
+
+    else:
+        Means = np.zeros((len(isum), nboot))
+        for ii, nwin in enumerate(isum):  
+            for iboot in range(nboot):
+                cs = np.random.randint(0,nc-nwin)
+                Means[ii,iboot] = np.mean( N[cs:cs+nwin] )
+
+    return Means, np.array(isum)
+
+def gateIntegrate(T2D, T2T, gpd, sigma, stackEfficiency=2.):
+    """ Gate integrate the signal to gpd, gates per decade
+        T2D = the time series to gate integrate, complex 
+        T2T = the abscissa values 
+        gpd = gates per decade 
+        sigma = estimate of standard deviation for theoretical gate noise 
+        stackEfficiency = exponential in theoretical gate noise, 2 represents ideal stacking
+    """
+    
+    # use artificial time gates so that early times are fully captured
+    T2T0 = T2T[0]
+    T2TD = T2T[0] - (T2T[1]-T2T[0])
+    T2T -= T2TD
+    
+    #####################################
+    # calculate total number of decades #
+    # windows edges are approximate until binning but will be adjusted to reflect data timing, this 
+    # primarily impacts bins with a few samples  
+    nd = np.log10(T2T[-1]/T2T[0])               
+    tdd = np.logspace( np.log10(T2T[0]), np.log10(T2T[-1]), (int)(gpd*nd)+1, base=10, endpoint=True) 
+    tdl = tdd[0:-1]                 # approximate window left edges
+    tdr = tdd[1::]                  # approximate window right edges
+    td = (tdl+tdr) / 2.             # approximate window centres
+
+
+    Vars = np.zeros( len(td) ) 
+    htd = np.zeros( len(td), dtype=complex )
+    isum = np.zeros( len(td), dtype=int )  
+
+    ii = 0
+    for itd in range(len(T2T)):
+        if ( T2T[itd] > tdr[ii] ):
+            ii += 1
+            # correct window edges to centre about data 
+            tdr[ii-1] = (T2T[itd-1]+T2T[itd])*.5 
+            tdl[ii  ] = (T2T[itd-1]+T2T[itd])*.5
+        isum[ii] += 1
+        htd[ii] += T2D[ itd ]
+        Vars[ii] += sigma**2
+        
+    td = (tdl+tdr) / 2.             # actual window centres
+    sigma2 = np.sqrt( Vars * ((1/(isum))**stackEfficiency) ) 
+
+    # Reset abscissa where isum == 1 
+    # when there is no windowing going on 
+    td[isum==1] = T2T[0:len(td)][isum==1]
+
+    tdd = np.append(tdl, tdr[-1])
+
+    htd /= isum # average
+    T2T += T2TD # not used  
+    return td+T2TD, htd, tdd+T2TD, sigma2, isum  # centre abscissa, data, window edges, error 
+
+PhiD = []
+def invert(Time, t, v, sig, lambdastar):
+    """ helper function that simply calls logBarrier, here to allow for drop in repacement  
+    """
+    #model = logBarrier(Time.Genv, 1e-2*v, Time.T2Bins, MAXITER=5000, sigma=1e-2*sig, alpha=1e6, smooth="Both") 
+    model = logBarrier(Time.Genv, 1e-2*v, Time.T2Bins, lambdastar, MAXITER=750, sigma=1e-2*sig, alpha=1e6, smooth="Smallest") 
+    PhiD.append(model[2])
+    return model
+
+def gateTest(vc, vgc, pperlin, boot, lamdastar):
+    """ Performs gate integration and adds random noise 
+        vc = clean data (dense)
+        vgc = clean data at gates 
+        boot = if "boot" then bootstrap the gate noise 
+        lambdastar = l-curve or discrepency principle
+        pperlin = percent perlin noise, noise floor is maintained at 2.00 PU 
+    """
+
+    t = np.arange(2e-4, .3601, 2e-4)
+    zeta = np.pi / 3.
+    
+    v = np.copy(vc) # important!  
+
+    # Scaling factors to keep noise floor constant with increasing levels of 
+    # Perlin noise. These were determined using populations of 5,000 and hold to 
+    # two significant digits (i.e, 2.00 PU) 
+    PF = {0.0:0,\
+          2.5:.450,\
+          5.0:.6125,\
+          7.5:.765,\
+          10.0:.87375,\
+          12.5:.9725,\
+          15.0:1.05,\
+          17.5:1.1275,\
+          20.0:1.20,\
+          22.5:1.265,\
+          25.0:1.325}
+ 
+    # random noise
+    np.random.seed() # necessary for thread pool, otherwise all threads can get same numbers 
+    sigma = 2.*(1.-1e-2*pperlin)
+    eps = np.random.normal(0, sigma, len(t)) + \
+       1j*np.random.normal(0, sigma, len(t)) 
+    eps += PF[pperlin] * perlin(len(t), L=.3601, sigma_f=.005, sigma_r=0.72) # 1 PU std
+    v += eps
+
+    # Noise based on residual  
+    sigmahat = np.std( v.imag )
+    
+    gt, gd, we, err, isum = gateIntegrate(v.real, t, 20, sigmahat)
+    ge = np.copy(err)
+    if boot=="boot":
+        Means, isum2 = bootstrapWindows(v.imag, 20000, isum[isum!=1], adapt=True)
+        # STD 
+        #err[isum!=1] = np.ma.std(Means, axis=1, ddof=1)[isum!=1]
+        # MAD, only for windows > 1 
+        c = stats.norm.ppf(3./4.)
+        err[isum!=1] = np.ma.median(np.ma.abs(Means), axis=1) / c 
+    if boot=="uniform":
+        err = sigmahat
+ 
+    return gt, gd, gd-vgc, err, v.real, t, isum
+
+if __name__ == "__main__":
+
+    import multiprocessing 
+    import itertools 
+
+    from GJIPlot import *
+    import numpy as np 
+    import matplotlib.pyplot as plt 
+    from matplotlib.ticker import FormatStrFormatter
+    from matplotlib import ticker
+    from collections import OrderedDict
+
+    if len(sys.argv)<4:
+        print ( "Python script for generating plots used in GJI publication")
+        print ( "useage:")
+        print ( "python gateIntegrate4.py   NoiseType   Sigma_i  Lambda*   " )
+        exit()
+
+    if sys.argv[1] not in ['0.0','2.5','5.0','7.5','10.0','12.5','15.0','17.5','20.0','22.5','25.0']: 
+        print ( "PercentPerlin: [0.0,2.5,5.0...25.0] ", "got", sys.argv[1])
+        exit(1)
+
+    if sys.argv[2] != "gauss" and sys.argv[2] != "boot" and sys.argv[2] != "uniform":
+        print ( "Sigma_i: gauss | boot | uniform")
+        exit(1)
+    
+    if sys.argv[3] != "lcurve" and sys.argv[3] != "discrepency": 
+        print ( "Lambda*: lcurve | discrepency ")
+        exit(1)
+
+    #offwhite = (.98,.98,.98)
+    offwhite = (1.,1.,1.) 
+    mDarkBrown = '#eb811b' # alert colour 
+    mDarkTeal = '#23373b'
+    mLightBrown= "#EB811B"
+    mLightGreen = "#14B03D"
+
+    # Time series plot
+    fig = plt.figure(figsize=(pc2in(18),pc2in(18)), facecolor=offwhite)
+    ax2  = fig.add_axes([.195, .175, .750, .75], facecolor=offwhite)  # time
+    
+    # Main plot 
+    fig2 = plt.figure(figsize=(pc2in(20),pc2in(2*.5*20)))
+    ax1  = fig2.add_axes([.175, .410*1.5, .6,   .225*1.5])
+    ax1c = fig2.add_axes([.800, .410*1.5, .025, .225*1.5])
+    ax3  = fig2.add_axes([.175, .100*1.5, .495, .225*1.5], facecolor='None')
+    ax3r  = fig2.add_axes([.175, .100*1.5, .495, .225*1.5], facecolor='None', rasterized=True, sharex=ax3, sharey=ax3)
+    ax3b = fig2.add_axes([.825, .100*1.5, .1,   .225*1.5])
+
+    SIG = []
+    ER = []
+    GD = []
+    GT = []
+    V = []
+    MOD = []
+    CONV = []
+    PHID = []
+    PHIM = []
+    LSTAR = []
+    ns = 10000 #10000  #10000 # number of realizations for PDF 
+    ni = 5000  #5000  #1000  # number of inversions to plot 
+    t = np.arange(2e-4, .3601, 2e-4) # CMR sampling
+ 
+    #CMAP = cmocean.cm.solar
+    CMAP = cmocean.cm.gray_r
+    #CMAP = cmocean.cm.haline
+    #CMAP = cmocean.cm.tempo
+
+    ##############################################
+    # set up model 
+    lowT2 = .001
+    hiT2 = 1.0 
+    nT2 = 30
+    spacing = "Log_10"
+    Time = pwcTime()
+    Time.setT2(lowT2, hiT2, nT2, spacing) 
+    Time.setSampling( np.arange(2e-4, .3601, 2e-4) )
+    Time.generateGenv()
+    tmod = np.zeros(nT2)
+    tmod [8] = .15  # distribution centres...to be smoothed
+    tmod [20] = .1
+    for i in range(2):
+        tmod = np.convolve(tmod, np.array([.0625,.125,.1875,.25,.1875,.125,.0625]), 'same') 
+
+    vc = 100. * np.dot(Time.Genv, tmod) + 0j # in PU
+    gt, gd, we, err, isum = gateIntegrate(vc, t, 20, 3)
+
+    ##############################################
+    # Set up inversion 
+    Time = pwcTime()
+    Time.setT2(lowT2, hiT2, nT2, spacing) 
+    Time.setSampling( gt ) 
+    Time.generateGenv()
+    vgc = 100.*np.dot(Time.Genv, tmod) + 0j # in PU
+
+    # make the Pool of workers
+    print("pool gate integrate")
+    with multiprocessing.Pool() as pool: 
+        results = pool.starmap(gateTest, zip(np.tile(vc, (ns, 1)), np.tile(vgc, (ns,1)), itertools.repeat(eval(sys.argv[1])), \
+                                                                                         itertools.repeat(sys.argv[2]), \
+                                                                                         itertools.repeat(sys.argv[3])))
+    print("done pool gate integrate")
+   
+    # parse out results 
+    for i in range(ns):
+        gt,gd,ge,err,v,vt,isum = results[i] 
+        V.append(v.real)
+        GT.append(gt.real)
+        GD.append(gd.real)
+        ER.append( ge.real / err.real )
+        SIG.append( err.real )
+
+    print("pool inversions")
+    with multiprocessing.Pool() as pool: 
+        invresults = pool.starmap(invert, zip(itertools.repeat(Time), GT[0:ni], GD[0:ni], SIG[0:ni], itertools.repeat(sys.argv[3]) )) 
+    #print("done pool inversions",results[:][0])
+
+    # Parse results 
+    for i in range(ns):
+        #print("Sym %", round(100.*i/(float)(1)/(float)(ns)))
+        # invert
+        if i < ni:
+            #mod, conv, phid = invert(Time, gt, gd.real, err)
+            if sys.argv[3] == "discrepency":
+                mod, conv, phid_final = invresults[i] 
+            else:
+                mod, conv, phid_final, phim, phid, lstar = invresults[i] 
+            MOD.append(mod)
+            CONV.append(conv)
+            PHID.append(phid)
+            PHIM.append(phim)
+            LSTAR.append(lstar)
+
+    PHIM = np.array(PHIM)
+    PHID = np.array(PHID)
+    ER = np.array(ER)
+    MOD = np.array(MOD)
+    GD = np.array(GD)
+
+    ####################
+    # Time series plot #       
+    ax2.plot( 1e3*vt, V[0], color=mDarkTeal, label="$V_N$", linewidth=1, zorder=-32) #, rasterized=True) 
+    ax2.errorbar( 1e3*gt, GD[0], yerr=SIG[0], fmt='.', markersize=6, color=mLightBrown, label="$V_G$")
+    ax2.set_ylim([-10,30])
+    leg1 = ax2.legend( labelspacing=0.2, scatterpoints=1, numpoints=1, frameon=True )   
+    fixLeg(leg1)
+    ax2.set_xscale("log", nonposx='clip')  
+    ax2.set_ylabel(r"$V_N$ (PU)")
+    ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+    ax2.set_xlabel("time (ms)")
+    deSpine(ax2)
+    fig.savefig( sys.argv[1] + "-" + sys.argv[2] + "-" + sys.argv[3] + "-ts.pdf", dpi=400, facecolor=offwhite,edgecolor=offwhite)
+
+    # histogram of error statistic
+    bins = np.linspace( -3, 3, 40, endpoint=True )
+    HIST = []
+    for i in range(0,np.shape(ER)[1]):
+        hist, edges = np.histogram(ER[:,i], bins=bins, density=False)        
+        HIST.append(hist)
+    HIST =  np.array(HIST)/(float)(ns) # normalize 
+        
+    im = ax1.pcolor(1e3*we, edges, HIST.T, cmap=CMAP, vmin=0, vmax=.1, rasterized=True)
+    im.set_edgecolor('face')
+    cb = plt.colorbar(im, ax1c, label=r"probability density", format=FormatStrFormatter('%1.2f'))
+    cb.solids.set_rasterized(True) 
+    tick_locator = ticker.MaxNLocator(nbins=4)
+    cb.locator = tick_locator
+    cb.update_ticks()
+
+    ax1.set_xscale("log", nonposx='clip')
+    ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+    ax1.set_xlabel("time (ms)")
+    ax1.set_ylabel(r"gate error $\left( \left( {V_G - V_T} \right) / {\tilde{\sigma_i}} \right)$")
+
+    LMT2 = [] 
+    THETA = []
+    MODERR = []
+       
+    # plot a random sample of ns instead?  
+    for i in range(ni):
+        # plot log mean and amplitude
+        model = MOD[i] 
+        theta = np.sum( model )
+        LogMeanT2 = np.exp(np.sum( model * np.log( Time.T2Bins ) ) / theta ) 
+        LMT2.append(LogMeanT2)
+        THETA.append( np.sum(model)  )
+        MODERR.append( np.linalg.norm(model-tmod) )
+                
+    CONV = np.array(CONV)
+    THETA = np.array(THETA)
+    MOD = np.array(MOD)
+    MODERR = np.array(MODERR)
+        
+    #############################
+    # plot all models, 1 colour #
+    ires = ax3r.plot( 1e3*np.tile(Time.T2Bins, (np.sum(np.array(CONV)) ,1)).T , 1e2*MOD.T, color=mDarkTeal, alpha=.01, lw=.5, label="$\mathbf{f}_I$", zorder=0, rasterized=True) 
+    lns2, = ax3r.plot(1e3*Time.T2Bins, 1e2*tmod, color=mLightBrown, linewidth=2, label="$\mathbf{f}_T$")
+
+    handles, labels = ax3r.get_legend_handles_labels()
+    by_label = OrderedDict(zip(labels, handles))
+    leg3 = ax3r.legend(by_label.values(), by_label.keys(), labelspacing=0.2, scatterpoints=1, numpoints=1, frameon=True , loc="upper right")            
+    for line in leg3.get_lines():
+        line.set_linewidth(1)
+    for lh in leg3.legendHandles: 
+        lh.set_alpha(1)
+    fixLeg(leg3)
+
+    ###########################
+    # Error histogram on side #
+    ax3b.hist( 1e2*MODERR, bins='auto', orientation="horizontal", color=mDarkTeal, stacked=True, density=True, range=(0,20)) 
+    ax3b.axhline(1e2*np.mean(MODERR), linewidth=1.25, color=mLightBrown) #, color=CMAP(0.7), zorder=1)
+    deSpine(ax3b)
+    ax3b.set_xscale("log", nonposx='clip')  
+    ax3b.set_ylabel(r"$\Vert \mathbf{f}_I -\mathbf{f}_T \Vert$") # %(m$^3$/m$^3$)") #, color="C0")   
+    ax3b.set_xlabel("log probability\ndensity") #, color="C0")   
+    ax3.set_xlim( (1e3*Time.T2Bins[0], 1e3*Time.T2Bins[-1]) ) 
+    ax3.set_ylim( (0,5) )  
+    ax3.set_xlabel("$T_2$ (ms)")   
+    ax3.set_ylabel("partial water content (PU)") #, color="C0")   
+    ax3.set_xscale("log", nonposx='clip')  
+    ax3.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+    plt.setp(ax3r.get_xticklabels(), visible=False)
+    plt.setp(ax3r.get_yticklabels(), visible=False)
+    deSpine(ax3)
+    deSpine(ax3r)
+    
+    np.save("pperlin" + str(round(1e1*eval(sys.argv[1]))) + "-" + sys.argv[2] + "-" + sys.argv[3] + "-err", MODERR) 
+
+    plt.savefig("pperlin" + str(round(1e1*eval(sys.argv[1]))) + "-" + sys.argv[2] + "-" + sys.argv[3] + ".pdf", dpi=600, facecolor=offwhite, edgecolor=offwhite)
+    
+    plt.show()
+

akvo/tressel/invertTA.py

@@ -447,8 +447,8 @@ def main():
             ax1.set_ylim( ifaces[-1], ifaces[0] )
 
             ax2 = ax1.twiny()
-            ax2.plot( np.sum(DINV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='red' )
-            ax2.set_xlabel(u"total water (m$^3$/m$^3$)")
+            ax2.plot( np.sum(DINV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='mediumblue' )
+            ax2.set_xlabel(u"total water (m$^3$/m$^3$)", color='mediumblue')
             ax2.set_ylim( ifaces[-1], ifaces[0] )
             ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )   
             ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
@@ -506,14 +506,17 @@ def main():
 
 
     ax2 = ax1.twiny()
-    ax2.plot( np.sum(INV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='red' )
-    ax2.set_xlabel(u"total water (m$^3$/m$^3$)")
+    ax2.plot( np.sum(INV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='mediumblue' )
+    ax2.set_xlabel(u"NMR total water (m$^3$/m$^3$)", color='mediumblue')
     ax2.set_ylim( ifaces[-1], ifaces[0] )
     ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )   
     ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
     #ax2.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed'  )
     if CalcDOI:
         ax2.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed'  )
+        
+    ax2.tick_params(axis='x', colors='mediumblue')
+    plt.setp(ax2.get_xticklabels(), color="mediumblue")
 
     plt.savefig("akvoInversion.pdf")
 
@@ -596,8 +599,8 @@ def main():
         #ax1.xaxis.set_label_position('top') 
 
         ax2 = ax1.twiny()
-        ax2.plot( np.sum(INVc, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='red' )
-        ax2.set_xlabel(u"total water (m$^3$/m$^3$)")
+        ax2.plot( np.sum(INVc, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='mediumblue' )
+        ax2.set_xlabel(u"NMR total water (m$^3$/m$^3$)", color='mediumblue')
         ax2.set_ylim( ifaces[-1], ifaces[0] )
         ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )   
         ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
@@ -606,6 +609,9 @@ def main():
             ax2.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed'  )
         #ax2.xaxis.set_label_position('bottom') 
         #fig.suptitle("Non linear inversion")
+        ax2.tick_params(axis='x', colors='mediumblue')
+        plt.setp(ax2.get_xticklabels(), color="mediumblue")
+
         plt.savefig("akvoInversionNL.pdf")
 
 

akvo/tressel/logbarrier-nnls.py

@@ -0,0 +1,325 @@
+from __future__ import division
+import numpy as np
+from scipy.sparse.linalg import iterative  as iter
+from scipy.sparse.linalg import spilu  as ilu
+import scipy.sparse.linalg as spla
+from scipy.sparse import eye as seye
+import pylab 
+import pprint 
+from scipy.optimize import nnls 
+
+import matplotlib.pyplot as plt
+
+from akvo.tressel.SlidesPlot import * 
+
+def PhiB(mux, minVal, x):
+    phib = mux * np.abs( np.sum(np.log( x-minVal)) )
+    return phib
+
+def curvaturefd(x, y, t):
+    x1 = np.gradient(x,t) 
+    x2 = np.gradient(x1,t) 
+    y1 = np.gradient(y,t) 
+    y2 = np.gradient(y1,t) 
+    return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)
+
+def curvatureg(x, y):
+    from scipy.ndimage import gaussian_filter1d
+    #first and second derivative
+    x1 = gaussian_filter1d(x, sigma=1, order=1)#, mode='constant', cval=x[-1])
+    x2 = gaussian_filter1d(x1, sigma=1, order=1)#, mode='constant', cval=y[-1])
+    y1 = gaussian_filter1d(y, sigma=1, order=1)#, mode='constant', cval=x1[-1])
+    y2 = gaussian_filter1d(y1, sigma=1, order=1)#, mode='constant', cval=y1[-1])
+    return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)
+
+def logBarrier(A, b, T2Bins, lambdastar, x_0=0, xr=0, alpha=10, mu1=10, mu2=10, smooth=False, MAXITER=70, fignum=1000, sigma=1, callback=None):
+    """Impliments a log barrier Tikhonov solution to a linear system of equations 
+        Ax = b  s.t.  x_min < x < x_max. A log-barrier term is used for the constraint
+    """
+    # TODO input
+    minVal = 0.0
+    #maxVal = 1e8
+
+    Wd =    (np.eye(len(b)) / (sigma))        # Wd = eye( sigma )
+    WdTWd = (np.eye(len(b)) / (sigma**2))     # Wd = eye( sigma )
+
+    ATWdTWdA = np.dot(A.conj().transpose(), np.dot( WdTWd, A ))     # TODO, implicit calculation instead?
+    N = np.shape(A)[1]                        # number of model
+    M = np.shape(A)[0]                        # number of data
+    SIGMA   = .25 # .25 # lower is more aggresive relaxation of log barrier  
+    EPSILON = 1e-25 #1e-35 
+    #SIGMA   = .05 # .25 # lower is more aggresive relaxation of log barrier  
+    #EPSILON = 1e-10 #1e-35 
+
+    # reference model
+    if np.size(xr) == 1:
+        xr =  np.zeros(N)     
+    
+    # initial guess
+    if np.size(x_0) == 1:
+        x = 1e-10 + np.zeros(N)     
+    else:
+        x = 1e-10 + x_0
+        
+    # Construct model constraint base   
+    Phim_base = np.zeros( [N , N] ) 
+    a1 = .05     # smallest too
+    
+    # calculate largest term            
+    D1 = 1./abs(T2Bins[1]-T2Bins[0])
+    D2 = 1./abs(T2Bins[2]-T2Bins[1])
+    #a2 = 1. #(1./(2.*D1+D2))    # smooth
+    
+    if smooth == "Both":
+        #print ("Both small and smooth model")
+        for ip in range(N):
+            D1 = 0.
+            D2 = 0.
+            if ip > 0:
+                #D1 = np.sqrt(1./abs(T2Bins[ip]-T2Bins[ip-1]))**.5
+                D1 = (1./abs(T2Bins[ip]-T2Bins[ip-1])) #**2
+            if ip < N-1:
+                #D2 = np.sqrt(1./abs(T2Bins[ip+1]-T2Bins[ip]))**.5
+                D2 = (1./abs(T2Bins[ip+1]-T2Bins[ip])) #**2
+            if ip > 0:
+                Phim_base[ip,ip-1] =   -(D1)      
+            if ip == 0:
+                Phim_base[ip,ip  ] = 2.*(D1+D2)  
+            elif ip == N-1:
+                Phim_base[ip,ip  ] = 2.*(D1+D2) 
+            else:
+                Phim_base[ip,ip  ] = 2.*(D1+D2)
+            if ip < N-1:
+                Phim_base[ip,ip+1] =   -(D2)  
+        Phim_base /= np.max(Phim_base)            # normalize 
+        Phim_base += a1*np.eye(N)
+
+    elif smooth == "Smooth":
+        #print ("Smooth model")
+        for ip in range(N):
+            if ip > 0:
+                Phim_base[ip,ip-1] = -1    # smooth in log space
+            if ip == 0:
+                Phim_base[ip,ip  ] = 2.05   # Encourage a little low model
+            elif ip == N-1:
+                Phim_base[ip,ip  ] = 2.5   # Penalize long decays
+            else:
+                Phim_base[ip,ip  ] = 2.1   # Smooth and small
+            if ip < N-1:
+                Phim_base[ip,ip+1] = -1    # smooth in log space
+
+    elif smooth == "Smallest":
+        for ip in range(N):
+            Phim_base[ip,ip  ] = 1.
+    else: 
+        print("non valid model constraint:", smooth)
+        exit()
+    
+    Phi_m =  alpha*Phim_base
+    WmTWm = Phim_base # np.dot(Phim_base, Phim_base.T)            
+    b_pre = np.dot(A, x)
+    phid = np.linalg.norm( np.dot(Wd, (b-b_pre)) )**2
+    phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
+
+    mu2 = phim
+    phib = PhiB(mu1, 0, x) 
+    mu1 = 1e-4* ((phid + alpha*phim) / phib)
+
+    PHIM = []
+    PHID = []
+    MOD = []
+
+    ALPHA = []
+    ALPHA.append(alpha)
+    #ALPHA = np.linspace( alpha, 1, MAXITER  )
+    for i in range(MAXITER):
+        #alpha = ALPHA[i]
+
+        Phi_m =  alpha*Phim_base
+        
+        # reset mu1 at each iteration 
+        # Calvetti -> No ; Li -> Yes   
+        # without this, non monotonic convergence occurs...which is OK if you really trust your noise 
+        mu1 = 1e-4* ((phid + alpha*phim) / phib) 
+
+        WmTWm = Phim_base # np.dot(Phim_base, Phim_base.T)            
+        phid_old = phid
+        inner = 0
+
+        First = True # guarantee entry 
+
+        xp = np.copy(x) # prior step x 
+            
+        b2a = np.dot(A.conj().transpose(), np.dot(WdTWd, b-b_pre) ) - alpha*np.dot(WmTWm,(x-xr))
+        xg = nnls(ATWdTWdA + Phi_m, b2a)
+        x = xg[0]
+        
+        #b_pre = np.dot(A, x)
+        #phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+        #phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
+
+        while ( (phib / (phid+alpha*phim)) > EPSILON  or First==True ):
+        #while ( False ):
+
+            First = False
+            # Log barrier, keep each element above minVal
+            X1 = np.eye(N) * (x-minVal)**-1           
+            X2 = np.eye(N) * (x-minVal)**-2         
+            
+            # Log barrier, keep sum below maxVal TODO normalize by component. Don't want to push all down  
+            #Y1 = np.eye(N) * (maxVal - np.sum(x))**-1           
+            #Y2 = np.eye(N) * (maxVal - np.sum(x))**-2         
+            
+            AA = ATWdTWdA + mu1*X2 + Phi_m 
+            M = np.eye( N ) * (1./np.diag(ATWdTWdA + mu1*X2 + Phi_m))
+            ##M = seye( N ).dot(1./np.diag(ATWdTWdA + mu1*X2 + Phi_m))
+
+            #print("Incomplete LU",flush=True)
+            #M2 = ilu(AA) #, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None, diag_pivot_thresh=None, relax=None, panel_size=None, options=None)     
+            #ilu(AA, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None, 
+                #diag_pivot_thresh=None, relax=None, panel_size=None, options=None)[source]      
+            #M = spla.LinearOperator(np.shape(AA), M2.solve) 
+            #print("dun ILU", flush=True)
+ 
+            # Solve system (newton step) (Li)
+            b2 = np.dot(A.conj().transpose(), np.dot(WdTWd, b-b_pre) ) + 2.*mu1*np.diag(X1) - alpha*np.dot(WmTWm,(x-xr))
+            ztilde = iter.cg(AA, b2, M=M) 
+            h = (ztilde[0].real) 
+
+            # TODO move nnls outside of this loop...use as starting point maybe?
+            #print("nnls", flush=True)
+            #xg = nnls(AA, b2)
+            #print("nnls done", flush=True)
+            #print("calling cg", flush=True)
+            #print("out of cg", flush=True)
+            #print(xg[0], ztilde[0] , flush=True)
+            #print(np.linalg.norm(xg[0] - ztilde[0]), flush=True )
+
+            
+            # Solve system (direct solution) (Calvetti) 
+            #b2 = np.dot(A.conj().transpose(), np.dot(WdTWd, b)) + 2.*mu1*np.diag(X1) - alpha*np.dot(WmTWm,(x-xr))
+            #ztilde = iter.cg(AA, b2, M=M, x0=xg[0]) 
+            #h = (ztilde[0].real - x) 
+
+            # step size
+            d = np.min( (1, 0.95 * np.min(x/np.abs(h+1e-120))) )
+            
+            ##########################################################
+            # Update and fix any over/under stepping
+            x += d*h
+        
+            # Determine mu steps to take
+            s1 = mu1 * (np.dot(X2, ztilde[0].real) - 2.*np.diag(X1))
+            #s2 = mu2 * (np.dot(Y2, ztilde[0].real) - 2.*np.diag(Y1))
+
+            # determine mu for next step
+            mu1 = SIGMA/N * np.abs(np.dot(s1, x))
+            #mu2 = SIGMA/N * np.abs(np.dot(s2, x))
+            
+            b_pre = np.dot(A, x)
+            phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+            phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
+            phib = PhiB(mu1, minVal, x)
+            inner += 1
+ 
+        PHIM.append(phim)      
+        PHID.append(phid)      
+        MOD.append(np.copy(x))  
+
+        # determine alpha
+        scale = 1.5*(len(b)/phid)
+        #alpha *= np.sqrt(scale)
+        alpha *= min(scale, .95) # was .85...
+        #print("alpha", min(scale, 0.99))
+        #alpha *= .99 # was .85...
+        ALPHA.append(alpha)
+        #alpha = ALPHA[i+1]
+            
+        print("inversion progress", i, alpha, np.sqrt(phid/len(b)), phim, flush=True)       
+        
+
+#         if np.sqrt(phid/len(b)) < 0.97: 
+#             ibreak = -1
+#             print ("------------overshot--------------------", alpha, np.sqrt(phid/len(b)), ibreak)
+#             alpha *= 2. #0
+#             x -= d*h
+#             b_pre = np.dot(A, x)
+#             phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+#             phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )#**2
+#             mu1 = ((phid + alpha*phim) / phib)
+        if lambdastar == "discrepency": 
+            if np.sqrt(phid/len(b)) < 1.00 or alpha < 1e-5: 
+                ibreak = 1
+                print ("optimal solution found", alpha, np.sqrt(phid/len(b)), ibreak)
+                break
+        # slow convergence, bail and use L-curve 
+        # TI- only use L-curve. Otherwise results for perlin noise are too spurious for paper.  
+        if lambdastar == "lcurve":
+            if i > 4: 
+                kappa = curvaturefd(np.log(np.array(PHIM)), np.log(np.array(PHID)), ALPHA[0:i+1])#ALPHA[0:-1])
+                #kappa = curvatureg(np.log(np.array(PHIM)), np.log(np.array(PHID)))
+                print("max kappa", np.argmax(kappa), "distance from", i-np.argmax(kappa)) 
+            if i > 4 and (i-np.argmax(kappa)) > 4: # ((np.sqrt(phid_old/len(b))-np.sqrt(phid/len(b))) < 1e-4) : 
+            #if np.sqrt(phid/len(b)) < 3.0 and ((np.sqrt(phid_old/len(b))-np.sqrt(phid/len(b))) < 1e-3): 
+                ibreak = 1
+                MOD = np.array(MOD)
+                print ("###########################") #slow convergence", alpha, "phid_old", np.sqrt(phid_old/len(b)), "phid", np.sqrt(phid/len(b)), ibreak)
+                print ("Using L-curve criteria") 
+                #kappa = curvaturefd(np.log(np.array(PHIM)), np.log(np.array(PHID)), ALPHA[0:-1])
+                #kappa2 = curvatureg(np.log(np.array(PHIM)), np.log(np.array(PHID)))
+                #kappa = curvature( np.array(PHIM), np.array(PHID))
+                x = MOD[ np.argmax(kappa) ]
+                b_pre = np.dot(A, x)
+                phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+                phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
+                mu1 = ((phid + alpha*phim) / phib) 
+                print ("L-curve selected", alpha, "phid_old", np.sqrt(phid_old/len(b)), "phid", np.sqrt(phid/len(b)), ibreak)
+                print ("###########################")
+                if np.sqrt(phid/len(b)) <= 1:
+                    ibreak=0
+
+                fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
+                ax1 = fig.add_axes( [.2,.15,.6,.7] )
+                #plt.plot( (np.array(PHIM)),  np.log(np.array(PHID)/len(b)), '.-')
+                #plt.plot(  ((np.array(PHIM))[np.argmax(kappa)]) , np.log( (np.array(PHID)/len(b))[np.argmax(kappa)] ), '.', markersize=12)
+                #plt.axhline()
+                lns1 = plt.plot( np.log(np.array(PHIM)),  np.log(np.sqrt(np.array(PHID)/len(b))), '.-', label="L curve")
+                lns2 = plt.plot( np.log(np.array(PHIM))[np.argmax(kappa)], np.log(np.sqrt(np.array(PHID)/len(b))[np.argmax(kappa)]), '.', markersize=12, label="$\lambda^*$")
+                ax2 = plt.twinx()
+                lns3 = ax2.plot( np.log(np.array(PHIM)), kappa, color='orange', label="curvature" )
+
+                # Single legend 
+                lns = lns1+lns3
+                labs = [l.get_label() for l in lns]
+                ax2.legend(lns, labs, loc=0)
+
+                ax1.set_xlabel("$\phi_m$")
+                ax1.set_ylabel("$\phi_d$")
+                
+                ax2.set_ylabel("curvature")
+
+                plt.savefig('lcurve.pdf')
+                break
+
+    PHIM = np.array(PHIM)
+    PHID = np.array(PHID)
+
+    if (i == MAXITER-1 ):
+        ibreak = 2
+        print("Reached max iterations!!", alpha, np.sqrt(phid/len(b)), ibreak)
+        kappa = curvaturefd(np.log(np.array(PHIM)), np.log(np.array(PHID)), ALPHA[0:-1])
+        x = MOD[ np.argmax(kappa) ]
+        b_pre = np.dot(A, x)
+        phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+        phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
+        mu1 = ((phid + alpha*phim) / phib) 
+
+    if lambdastar == "lcurve":
+        return x, ibreak, np.sqrt(phid/len(b)), PHIM, PHID/len(b), np.argmax(kappa)
+    else:
+        return x, ibreak, np.sqrt(phid/len(b))
+
+
+
+if __name__ == "__main__":
+    print("Test")

akvo/tressel/mrsurvey.py

@@ -6,8 +6,9 @@ import sys
 import scipy
 from scipy import stats
 import copy
-import struct
+import struct, glob
 from scipy.io.matlab import mio
+import pandas as pd
 from numpy import pi
 from math import floor
 import matplotlib as mpl
@@ -243,6 +244,7 @@ class GMRDataProcessor(SNMRDataProcessor):
         self.nTransVersion   = HEADER[8]           # Transmitter version
         self.nDAQVersion     = HEADER[9]           # DAQ software version 
         self.nInterleaves    = HEADER[10]          # num interleaves
+        self.Instrument      = "GMR" 
 
         self.gain()
         
@@ -315,7 +317,9 @@ class GMRDataProcessor(SNMRDataProcessor):
                 istack = 0
                 for sstack in self.DATADICT["stacks"]:
                     if self.pulseType == "FID" or pulse == "Pulse 2":
-                        if floor(self.nDAQVersion) < 2:
+                        if self.Instrument == "MIDI 2":
+                            mod = 1
+                        elif self.Instrument == "GMR" and floor(self.nDAQVersion) < 2:
                             mod = 1
                         else:
                             mod = (-1.)**(ipm%2) * (-1.)**(sstack%2)
@@ -950,27 +954,36 @@ class GMRDataProcessor(SNMRDataProcessor):
         self.doneTrigger.emit() 
 
 
-    def sumData(self, canvas, fred):
+    def sumData(self, canvas, plusminus, sumAll):
+        print("In sumData", plusminus, sumAll)
+        if plusminus == "sum":
+            splus = "+"
+        else:
+            splus = "-"
+
         chans = copy.deepcopy(self.DATADICT[self.DATADICT["PULSES"][0]]["chan"]) #= np.array( ( self.DATADICT[pulse]["chan"][0], ) )
         nchan = len(chans)
         # Sum permutations of two channel combos
         for ich in range(nchan-1):
             for ch in chans[ich+1:]:
-                chsum = chans[ich] + "+" + ch
+                chsum = chans[ich] + splus + ch
                 for pulse in self.DATADICT["PULSES"]:
                     self.DATADICT[pulse][chsum] = {} 
                     for ipm in range(self.DATADICT["nPulseMoments"]):
                         self.DATADICT[pulse][chsum][ipm] = {} 
                         for istack in self.DATADICT["stacks"]:
-                            self.DATADICT[pulse][chsum][ipm][istack] = self.DATADICT[pulse][chans[ich]][ipm][istack] - self.DATADICT[pulse][ch][ipm][istack] 
-                    if chsum == "1+2":
+                            if plusminus == "sum":
+                                self.DATADICT[pulse][chsum][ipm][istack] = self.DATADICT[pulse][chans[ich]][ipm][istack] + self.DATADICT[pulse][ch][ipm][istack] 
+                            else:
+                                self.DATADICT[pulse][chsum][ipm][istack] = self.DATADICT[pulse][chans[ich]][ipm][istack] - self.DATADICT[pulse][ch][ipm][istack] 
+                    #if chsum == "1+2":
                         #self.DATADICT[pulse]["rchan"].pop()
                         #self.DATADICT[pulse]["rchan"].pop()
-                        self.DATADICT[pulse]["chan"].append(chsum)
+                    self.DATADICT[pulse]["chan"].append(chsum)
 
         # Sum all channels 
-        sumall = False
-        if sumall:
+        #sumall = False
+        if sumAll:
             chsum = ""
             for ch in chans:
                 chsum += ch + "+" 
@@ -1735,7 +1748,9 @@ class GMRDataProcessor(SNMRDataProcessor):
     def enableDSP(self):
         self.enableDSPTrigger.emit() 
 
-    def adaptiveFilter(self, M, flambda, truncate, mu, PCA, canvas):
+    def adaptiveFilter(self, M, flambda, truncate, mu, PCA, plot, canvas):
+
+        #plot = False
 
         canvas.reAx2(shx=False, shy=False)
         # ax1 is top plot of filter taps 
@@ -1762,13 +1777,14 @@ class GMRDataProcessor(SNMRDataProcessor):
         for istack in self.DATADICT["stacks"]:
             for ipm in range(self.DATADICT["nPulseMoments"]):
                 for pulse in self.DATADICT["PULSES"]:
-                    canvas.softClear()
-                    mmax = 0
-                    for ichan in self.DATADICT[pulse]["chan"]:
-                        canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9* self.DATADICT[pulse][ichan][ipm][istack], alpha=.5) 
-                        mmax = max(mmax, np.max(1e9*self.DATADICT[pulse][ichan][ipm][istack])) 
-                    canvas.ax2.set_ylim(-mmax, mmax) 
-                    canvas.ax2.set_prop_cycle(None)
+                    if plot:
+                        canvas.softClear()
+                        mmax = 0
+                        for ichan in self.DATADICT[pulse]["chan"]:
+                            canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9* self.DATADICT[pulse][ichan][ipm][istack], alpha=.5) 
+                            mmax = max(mmax, np.max(1e9*self.DATADICT[pulse][ichan][ipm][istack])) 
+                        canvas.ax2.set_ylim(-mmax, mmax) 
+                        canvas.ax2.set_prop_cycle(None)
                     for ichan in self.DATADICT[pulse]["chan"]:
                         #H = np.zeros(M)
                         RX = []
@@ -1810,42 +1826,46 @@ class GMRDataProcessor(SNMRDataProcessor):
                                                         RX,\
                                                         M, mu, PCA, flambda, H[pulse][ichan])
                                 iloop += 1
-                            print("Recalled ", iloop, "times with norm=", np.linalg.norm(hm1-H[pulse][ichan]))
+                            #print("Recalled ", iloop, "times with norm=", np.linalg.norm(hm1-H[pulse][ichan]))
                         else:
                             [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
                                                         RX,\
                                                         M, mu, PCA, flambda, H[pulse][ichan])
                         # replace
                         if truncate:
-                            canvas.ax2.plot( self.DATADICT[pulse]["TIMES"][0:itrunc], 1e9* e[::-1][0:itrunc],\
-                                label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan="  + str(ichan))
+                            if plot:
+                                canvas.ax2.plot( self.DATADICT[pulse]["TIMES"][0:itrunc], 1e9* e[::-1][0:itrunc],\
+                                    label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan="  + str(ichan))
                             self.DATADICT[pulse][ichan][ipm][istack] = e[::-1][0:itrunc]
                         else:
-                            canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9* e[::-1],\
-                                label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan="  + str(ichan))
+                            if plot:
+                                canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9* e[::-1],\
+                                    label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan="  + str(ichan))
                             self.DATADICT[pulse][ichan][ipm][istack] = e[::-1]
-                           
-     
-                        #canvas.ax1.plot( H[pulse][ichan].reshape(-1, len(RX)) ) # , label="taps") 
-                        canvas.ax1.plot( H[pulse][ichan][::-1].reshape(M, len(RX), order='F' ) ) #.reshape(-1, len(RX)) ) # , label="taps") 
 
-                        canvas.ax2.legend(prop={'size':10}, loc='upper right')
-                        #canvas.ax2.legend(prop={'size':6}, loc='upper right')
+                                               
+                        if plot:
+                            #canvas.ax1.plot( H[pulse][ichan].reshape(-1, len(RX)) ) # , label="taps") 
+                            canvas.ax1.plot( H[pulse][ichan][::-1].reshape(M, len(RX), order='F' ) ) #.reshape(-1, len(RX)) ) # , label="taps") 
 
-                        mh = np.max(np.abs( H[pulse][ichan] ))
-                        canvas.ax1.set_ylim( -mh, mh )
+                            canvas.ax2.legend(prop={'size':10}, loc='upper right')
+                            #canvas.ax2.legend(prop={'size':6}, loc='upper right')
 
-                        canvas.ax2.set_xlabel(r"time (s)", fontsize=10)
-                        canvas.ax2.set_ylabel(r"signal (nV)", fontsize=10)
+                            mh = np.max(np.abs( H[pulse][ichan] ))
+                            canvas.ax1.set_ylim( -mh, mh )
 
-                        canvas.ax1.set_xlabel(r"filter tap index", fontsize=10)
-                        canvas.ax1.set_ylabel(r"tap amplitude", fontsize=10)
+                            canvas.ax2.set_xlabel(r"time (s)", fontsize=10)
+                            canvas.ax2.set_ylabel(r"signal (nV)", fontsize=10)
 
-                    canvas.fig.tight_layout()
-                    deSpine(canvas.ax1)
-                    deSpine(canvas.ax2)
+                            canvas.ax1.set_xlabel(r"filter tap index", fontsize=10)
+                            canvas.ax1.set_ylabel(r"tap amplitude", fontsize=10)
 
-                    canvas.draw()
+                    if plot:
+                        canvas.fig.tight_layout()
+                        deSpine(canvas.ax1)
+                        deSpine(canvas.ax2)
+
+                        canvas.draw()
 
                     # truncate the reference channels too, in case you still need them for something. 
                     # Otherwise they are no longer aligned with the data 
@@ -2484,6 +2504,187 @@ class GMRDataProcessor(SNMRDataProcessor):
                 self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] = DATA[:,2][nps:nps+npul] * invCGain
             nds += ndr
             nps += ndr 
+    
+    def readMIDI2Header(self, directory):
+        """ Reads the header of the FID_Q1_D0_R1.dat which should be in very MIDI directory
+        """
+        print("Searching ", directory, "for Q files")
+        self.QFiles = []
+        for file in glob.glob(directory+'/FID_Q*R1.dat'):
+            self.QFiles.append(file)
+
+        fidname = self.QFiles[0] # "/FID_Q1_D0_R1.dat"
+
+        HEADER = {}
+        with open(fidname, 'rb') as FID:
+            #print(FID.name)
+            headerLine = 0 
+            for i in range(21):
+                tags  = FID.readline().split(b']')
+                tags[0] = tags[0].strip(b'[')
+                tags[1] = tags[1].decode().strip(  )
+                HEADER[ ''.join(map(chr, tags[0])) ] = tags[1]
+        #print(HEADER)
+
+        pulseTypeDict = {
+            0 : lambda: "FID",
+            2 : lambda: "T1",
+            3 : lambda: "SPINECHO",
+            4 : lambda: "4PhaseT1"
+        }
+        
+        pulseLengthDict = {
+            1 : lambda x: np.ones(1) * x,
+            2 : lambda x: np.ones(2) * x,
+            3 : lambda x: np.array([x, 2.*x]),
+            4 : lambda x: np.ones(2) * x
+        }
+
+        if HEADER["P2"] == "0":
+            self.pulseType       = "FID"
+        else:
+            self.pulseType       = "T1"
+
+        self.transFreq       = float(HEADER[ 'Excitation frequency /Hz' ])
+        self.maxBusV         = '24 V'
+        self.pulseLength     = [ float(HEADER['P1'])/self.transFreq, float(HEADER['P2'])/self.transFreq] 
+        self.interpulseDelay = 1e-3*float(HEADER["Delay /ms"])       # for T2, Spin Echo
+        self.repetitionDelay = float(HEADER['Pause /s'])             # delay between q pulses 
+        self.nPulseMoments   = len(self.QFiles)                           # Number of pulse moments per stack
+        self.TuneCapacitance = 0                                     # tuning capacitance in uF
+        self.nTransVersion   = "MIDI 2"                              # Transmitter version
+        self.nDAQVersion     = HEADER["Software Revision"]           # DAQ software version 
+        self.nInterleaves    = 0                                     # num interleaves
+        self.Instrument = "MIDI 2" 
+        self.datadir = directory 
+        self.MIDIGain = HEADER["Gains"].split()
+        # default 
+        self.samp                   = float(HEADER["Data Rate /Hz"])        # sampling frequency
+        self.dt                     = 1./self.samp                          # sampling rate 
+
+
+
+    def loadMIDI2(self, directory, procStacks, chanin, rchanin, FIDProc, canvas, deadTime, plot):
+        """Reads a MRS MIDI2 experiment. 
+        """
+        canvas.reAx3(True,False)
+
+        chan = []
+        for ch in chanin:
+            chan.append(str(ch)) 
+        
+        rchan = []
+        for ch in rchanin:
+            rchan.append(str(ch)) 
+
+        #print(self.QFiles)
+        # Set up the same structure as GMR 
+        PULSES = [FIDProc]
+        PULSES = ["Pulse 1"]
+
+        self.DATADICT = {}
+        self.DATADICT["nPulseMoments"] = self.nPulseMoments
+        self.DATADICT["stacks"] = procStacks
+        self.DATADICT["PULSES"] = PULSES
+        for pulse in PULSES: 
+            self.DATADICT[pulse] = {}
+            self.DATADICT[pulse]["chan"] = chan        # TODO these should not be a subet of pulse! for GMR all 
+            self.DATADICT[pulse]["rchan"] = rchan      #      data are consistent 
+            self.DATADICT[pulse]["CURRENT"] = {} 
+            for ichan in np.append(chan,rchan):
+                self.DATADICT[pulse][ichan] = {}
+                for ipm in range(self.nPulseMoments):
+                    self.DATADICT[pulse][ichan][ipm] = {} 
+                    self.DATADICT[pulse]["CURRENT"][ipm] = {} 
+                    for istack in procStacks:
+                        pass 
+                        #print("pulse", pulse, "ichan",type(ichan), ichan, "ipm", type(ipm), ipm, "istack",type(istack), istack)
+                        #
+                        #self.DATADICT[pulse][ichan][ipm][istack] = np.zeros(3)
+                        #self.DATADICT[pulse]["CURRENT"][ipm][istack] = np.zeros(3) 
+
+        iistack = 0
+        idead = int( (self.pulseLength[0]+deadTime) / self.dt)
+
+        for iq in range(self.nPulseMoments):
+
+            fidbase = self.datadir #+ self.QFiles[iq][0:-5]
+            
+            for istack in procStacks:
+                #fidname = fidbase + str(iq).zfill(2) + ".dat"
+                fidname = fidbase + "/FID_Q" + str(iq) + "_D0_R" + str(istack) + ".dat"
+ 
+                with open(fidname, 'rb') as FID:
+                    #print(FID.name)
+                    headerLine = 0 
+                    for i in range(100):
+                        line = FID.readline().strip()
+                        headerLine += 1
+                        if line == b'[Begin Data]':
+                            break 
+
+                    # use numpy for now, consider pandas for faster read? 
+                    DATA = np.genfromtxt(FID, skip_header=0, skip_footer=1 )
+                    #DATA = pd.read_csv(fidname, skiprows=headerLine, skipfooter=1, sep='\t', encoding='ascii')
+           
+                    for ichan in np.append(chan,rchan):
+                 
+                        if int(ichan) <= 3: 
+                            self.DATADICT["Pulse 1"][ichan][iq][istack] = DATA[idead:,int(ichan)] / float(self.MIDIGain[int(ichan)-1])
+
+                        elif int(ichan) > 3:     
+                            self.DATADICT["Pulse 1"][ichan][iq][istack] = DATA[idead:,int(ichan)+1] / float(self.MIDIGain[int(ichan)-1])
+                   
+                        # truncate after dead time  
+
+                        self.DATADICT["Pulse 1"]["TIMES"] = DATA[idead:,0]
+
+                        # truncate until dead time
+                        ipulse = int(self.pulseLength[0] / self.dt)
+                        self.DATADICT["Pulse 1"]["PULSE_TIMES"] = DATA[0:ipulse,0]
+                        self.DATADICT["Pulse 1"]["CURRENT"][iq][istack] = DATA[0:ipulse,4]
+
+                        if plot: 
+                            canvas.softClear()                           
+
+                            for ichan in chan:
+                                canvas.ax1.plot(self.DATADICT["Pulse 1"]["PULSE_TIMES"], self.DATADICT["Pulse 1"]["CURRENT"][iq][istack] , color='black')
+                                canvas.ax3.plot(self.DATADICT["Pulse 1"]["TIMES"],       self.DATADICT["Pulse 1"][ichan][iq][istack], label="Pulse 1 FID data ch. "+str(ichan)) #, color='blue')
+
+                            for ichan in rchan:
+                                canvas.ax2.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][iq][istack], label="Pulse 1 FID ref ch. "+str(ichan)) #, color='blue')
+
+                            # reference axis
+                            canvas.ax2.tick_params(axis='both', which='major', labelsize=10)
+                            canvas.ax2.tick_params(axis='both', which='minor', labelsize=10)
+                            #canvas.ax2.xaxis.set_ticklabels([])
+                            plt.setp(canvas.ax2.get_xticklabels(), visible=False)
+                            canvas.ax2.legend(prop={'size':10}, loc='upper right')
+                            canvas.ax2.set_title("stack "+str(istack)+" pulse index " + str(iq), fontsize=10)
+                            canvas.ax2.set_ylabel("RAW signal [V]", fontsize=10)
+
+                            canvas.ax1.set_ylabel("Current (A)", fontsize=10) 
+                            canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
+                            canvas.ax1.set_xlabel("time (s)", fontsize=10)
+
+                            canvas.ax3.legend(prop={'size':10}, loc='upper right')
+                            canvas.ax3.set_ylabel("RAW signal [V]", fontsize=10)
+                    
+                            canvas.fig.tight_layout()
+                            canvas.draw()
+
+                percent = (int) (1e2*((float)(iistack) / (len(procStacks)*self.nPulseMoments)))
+                self.progressTrigger.emit(percent) 
+                iistack += 1
+
+#                percent = (int) (1e2*((float)((iistack*self.nPulseMoments+ipm+1))  / (len(procStacks)*self.nPulseMoments)))
+#                self.progressTrigger.emit(percent) 
+#                iistack += 1
+
+        self.enableDSP()    
+        self.doneTrigger.emit()
+
+
 
     def loadGMRASCIIT1( self, rawfname, istack ):
         """Based on the geoMRI instrument manufactured by VistaClara. Imports 

setup.py

@@ -21,7 +21,7 @@ with open("README.md", "r") as fh:
     long_description = fh.read()
 
 setup(name='Akvo',     
-      version='1.6.3', 
+      version='1.7.0', 
       python_requires='>3.7.0', # due to pyLemma 
       description='Surface nuclear magnetic resonance workbench',
       long_description=long_description,