AkvoFluo/GUI/main.py

#!/usr/bin/env python
from __future__ import division

import sys, os
sys.path.append( '../.' )

import matplotlib
matplotlib.use('Qt4Agg')
matplotlib.rcParams['backend.qt4']='PySide'
#matplotlib.rcParams['mathtext.fontset']='stixsans' # sans-serif in plots
from MRProc  import MRProc
from pwctimeWhite import pwcTime
from logbarrier import *

from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
from pylab import meshgrid
from matplotlib.ticker import MaxNLocator

#from pyqtgraph.widgets.MatplotlibWidget import MatplotlibWidget
from MatplotlibWidget import MatplotlibWidget

#from PyQt5 import QtCore, QtGui 

from PySide import QtCore, QtGui 
from borehole2 import Ui_MainWindow
import brewer2mpl

import numpy as np
from matplotlib.figure import Figure
import matplotlib.pyplot as plt

from multiprocessing import Lock, Process, Queue, Pool, current_process 

#from pylasdev import * 

# thread is good for GUI/process control, but use multiprocess if you want to exploit many cores 
# in traditional HPC type application
try:    
    import thread 
except ImportError:
    import _thread as thread #Py3K changed it.

class MyForm(QtGui.QMainWindow): #, threading.Thread):

    def __init__(self, parent=None):

        #QtGui.QWidget.__init__(self, parent)
        super(MyForm, self).__init__(parent)
        self.ui = Ui_MainWindow()
        self.ui.setupUi(self)

        # Add progressbar to statusbar
        self.ui.barProgress =  QtGui.QProgressBar()
        self.ui.statusbar.addPermanentWidget(self.ui.barProgress, 0);
        self.ui.barProgress.setMaximumSize(100, 16777215);
        self.ui.barProgress.hide();        

        # signals and slots
        QtCore.QObject.connect(self.ui.actionLoadRecord, QtCore.SIGNAL("triggered()"), self.openSingleRecord )
        QtCore.QObject.connect(self.ui.actionLoadRecords, QtCore.SIGNAL("triggered()"), self.openMultipleRecords )
        QtCore.QObject.connect(self.ui.actionLoadLogSeries, QtCore.SIGNAL("triggered()"), self.openVCLogSeries )
        QtCore.QObject.connect(self.ui.actionOpen_CMR_Log, QtCore.SIGNAL("triggered()"), self.openCMRLog )
        
        # preprocess
        QtCore.QObject.connect(self.ui.WindowStackGO, QtCore.SIGNAL("clicked()"), self.windowAndStack )
        QtCore.QObject.connect(self.ui.envelopeGO, QtCore.SIGNAL("clicked()"), self.envelope )
        QtCore.QObject.connect(self.ui.phaseGO, QtCore.SIGNAL("clicked()"), self.phase )
        QtCore.QObject.connect(self.ui.gateGO, QtCore.SIGNAL("clicked()"), self.gate )
        QtCore.QObject.connect(self.ui.batchLoadDataPushButton, QtCore.SIGNAL("clicked()"), self.batch )
        
        # inversion
        QtCore.QObject.connect(self.ui.monoGO, QtCore.SIGNAL("clicked()"), self.mono )
        QtCore.QObject.connect(self.ui.multiGO, QtCore.SIGNAL("clicked()"), self.multi )
        
        # Diffusion
        QtCore.QObject.connect(self.ui.JgGO, QtCore.SIGNAL("clicked()"), self.jgFit )
        QtCore.QObject.connect(self.ui.diffusionGO, QtCore.SIGNAL("clicked()"), self.diffusion )
        QtCore.QObject.connect( self.ui.jgComboBox, QtCore.SIGNAL("currentIndexChanged(int)"), self.jgc ) #this, SLOT(comboBoxIndexChanged()) )       
        
        # forward modelling
        QtCore.QObject.connect(self.ui.modelGO, QtCore.SIGNAL("clicked()"), self.model )
        
        # \kappa estimates
        QtCore.QObject.connect(self.ui.sdrGO, QtCore.SIGNAL("clicked()"), self.sdr )
        QtCore.QObject.connect(self.ui.tcGO, QtCore.SIGNAL("clicked()"), self.tc )
 

        #QtCore.QObject.connect(self.ui.windowFilterGO, QtCore.SIGNAL("clicked()"), self.windowFilter )
        #QtCore.QObject.connect(self.ui.despikeGO, QtCore.SIGNAL("clicked()"), self.despikeFilter )
        #QtCore.QObject.connect(self.ui.adaptGO, QtCore.SIGNAL("clicked()"), self.adaptFilter )
        #QtCore.QObject.connect(self.ui.adaptFDGO, QtCore.SIGNAL("clicked()"), self.adaptFilterFD )

        ##########################################################################
        # modelling Table 
        self.ui.modelTableWidget.setRowCount(10)       
        self.ui.modelTableWidget.setColumnCount(5)      
        #Labels = QtCore.QStringList( ) # ["A", "B", "C"] )
        self.ui.modelTableWidget.setHorizontalHeaderLabels( ["Top [m]","Bottom [m]","T_2 [ms]","Porosity [%]", "noise std"] )
        self.ui.modelTableWidget.cellChanged.connect(self.cellChanged) 

        self.ui.modelTableWidget.setDragDropOverwriteMode(False)
        self.ui.modelTableWidget.setDragEnabled(True)
        self.ui.modelTableWidget.setDragDropMode(QtGui.QAbstractItemView.InternalMove)
        #item->setFlags(item->flags() & ~(Qt::ItemIsDropEnabled))

        # Only for CMR data
        self.burst = False

        # Dictionary of procesors
        self.ORS = {}
        
        # Do we have a whole hole dataset?
        self.DepthProfile = False

    def jgc(self):
        print("JCD Docc")
        if ( self.ui.jgComboBox.currentText() == "Constant" ):
            self.ui.jgcSpinBox.setEnabled(True)
        else:
            self.ui.jgcSpinBox.setEnabled(False)

    def ORSConnect(self, ff): 
        self.ORS[ff] = MRProc(  )  # Fast
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("updateProgress(int)"), self.updateProgressBar)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("enableDSP()"), self.enableDSP)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("enableINV()"), self.enableINV)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotRAW()"), self.plotRAW)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotWIN()"), self.plotWIN)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotENV()"), self.plotENV)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotLOG10ENV()"), self.plotLOG10ENV)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotMONO()"), self.plotMONO)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotBI()"), self.plotBI)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("plotDIST(int)"), self.plotDIST)
        QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("doneStatus()"), self.doneStatus)
        QtCore.QObject.connect(self,         QtCore.SIGNAL("doneStatus2()"), self.doneStatus)

    def doneStatus(self): #, enable): # unlocks GUI
        self.ui.statusbar.clearMessage ( )
        self.ui.barProgress.hide()
        #self.updateProc()       
        #self.enableAll()

    def lock(self, string):
        self.ui.statusbar.showMessage ( string )
        self.ui.barProgress.show()
        self.ui.barProgress.setValue(0)
        self.disable()

    def unlock(self):
        self.ui.statusbar.clearMessage ( )
        self.ui.barProgress.hide()
        self.enableAll()

    def done(self):
        self.ui.statusbar.showMessage ( "" )


    def updateProgressBar(self, percent):
        self.ui.barProgress.setValue(percent)


    def disable(self):
        self.ui.WindowStackBox.setEnabled(False)
        self.ui.envelopeGroupBox.setEnabled(False)
        self.ui.phaseGroupBox.setEnabled(False)
        self.ui.gateGroupBox.setEnabled(False)
        #self.ui.adaptFDBox.setEnabled(False)
        #self.ui.qCalcGroupBox.setEnabled(False)    
        #self.ui.FDSmartStackGroupBox.setEnabled(False)    


    def enableAll(self):
        self.enableDSP()
        #self.enableQC()

    def enableDSP(self):

        # Bandpass filter
        self.ui.WindowStackBox.setEnabled(True)
        self.ui.WindowStackBox.setChecked(True)
        #self.ui.WindowStackGO.setEnabled(False) # default on need to design first

        # downsample 
        self.ui.envelopeGroupBox.setEnabled(True)
        self.ui.envelopeGroupBox.setChecked(True)
 
        # Phase
        self.ui.phaseGroupBox.setEnabled(True)
        self.ui.phaseGroupBox.setChecked(True)

        # FD Adaptive filtering 
        #self.ui.adaptFDBox.setEnabled(True)
        #self.ui.adaptFDBox.setChecked(True)

        #self.enableQC() 
        #QtCore.QObject.connect(self.ORS, QtCore.SIGNAL("updateProc()"), self.updateProc)

    def enableINV(self):
        
        # time gating 
        self.ui.gateGroupBox.setEnabled(True)
        self.ui.gateGroupBox.setChecked(True)

        # mono
        self.ui.monoGroupBox.setEnabled(True)
        self.ui.monoGroupBox.setChecked(True)
        
        # bi
        #self.ui.biGroupBox.setEnabled(True)
        #self.ui.biGroupBox.setChecked(True)
        
        # dist 
        self.ui.multiGroupBox.setEnabled(True)
        self.ui.multiGroupBox.setChecked(True)

    def windowAndStack(self, ff=None):

        if ff == None:
            ff = self.headerstr[0]

        #self.ORS.WindowAndStack()                  # Fast
        self.lock("window and stack records")
        #self.ORS = MRProc(  )  # Fast
        #self.ORS.loadORSFile( self.headerstr )

        rawThread = thread.start_new_thread(self.ORS[ff].WindowAndStack, \
                (str(self.ui.fTypeComboBox.currentText()), self.ui.winTrimLeftSpinBox.value(), self.ui.winTrimRightSpinBox.value() )) 
        self.ui.logTextBrowser.append( "Windowed and stack: " + str(self.ui.fTypeComboBox.currentText() ))
    
    def envelope(self, ff=None):
        
        if ff == None:
            ff = self.headerstr[0]
        
        self.lock("window and stack records")
        rawThread = thread.start_new_thread(self.ORS[ff].T2EnvelopeDetect, (self.ui.offsetSpinBox.value(),) ) 
        self.ui.logTextBrowser.append( "Envelope Detect: " + str(self.ui.fTypeComboBox.currentText() ))
        # Now we can do inversions or gate integraion
        #self.enableINV( )

    def phase(self, ff=None):

        if ff == None:
            ff = self.headerstr[0]

        self.lock("window and stack records")
        rawThread = thread.start_new_thread(self.ORS[ff].CorrectedAmplitude, (0,0) ) 
        self.ui.logTextBrowser.append( "Phase correct") 
    
    def mono(self):

        self.lock("mono-exponential fit")
        # TODO thread this
        if self.ui.expComboBox.currentText() == "mono":
            for ff in self.headerstr:
                #rawThread = thread.start_new_thread(self.ORS[ str(ff) ].MonoFit, (self.ui.maskNSpinBox_2.value(), str(self.ui.interceptComboBox.currentText())) ) 
                self.ORS[ff].MonoFit(self.ui.maskNSpinBox_2.value(), str(self.ui.interceptComboBox.currentText())) 
                self.ui.logTextBrowser.append( "mono-exponential fit") 
            self.plotMONO()
        elif self.ui.expComboBox.currentText() == "bi":
            for ff in self.headerstr:
                #rawThread = thread.start_new_thread(self.ORS[ str(ff) ].MonoFit, (self.ui.maskNSpinBox_2.value(), str(self.ui.interceptComboBox.currentText())) ) 
                self.ORS[ff].BiFit(self.ui.maskNSpinBox_2.value(), str(self.ui.interceptComboBox.currentText())) 
                self.ui.logTextBrowser.append( "bi-exponential fit") 
            self.plotBI()


    def multiMP(self, q):
        # TODO consider use of pool here instead! Easiest way to limit number of processes
        #      Pool does not like to use class member functions. There are workarounds but I can't be bothered. 
        processes = []
#        Smooth = False
#        if str(self.ui.DistConstraint.currentText()) == "Smooth":
#            Smooth = True
#        elif str(self.ui.DistConstraint.currentText()) == "Smallest":
#            Smooth = False
#        else:
#            print ("Smooth operator my ass")
#            exit(3)
        for ff in self.headerstr:
            p = Process( target=self.ORS[ff].DistFitMP, args= (q, ff, 0, self.ui.maskNSpinBox_4.value(), self.ui.nT2SpinBox.value(), \
                                 1e-3*self.ui.lowT2SpinBox.value(), 1e-3*self.ui.hiT2SpinBox.value(), self.ui.distT2Box.currentText(), \
                                 self.ui.DistConstraint.currentText(), self.ui.betaScale.value()
                         ) )
            processes.append(p)
        return processes
 
    def multi(self):
        import pickle
        import multiprocessing

        self.lock("multi-exponential fit")
            
        q = Queue()
        procs = self.multiMP(q) 

        zombies = []
        while len(procs) > 0:  
            if  len(multiprocessing.active_children()) < multiprocessing.cpu_count():
                p = procs.pop()
                p.start()
                zombies.append(p)

        # kill the zombies...well manage them again
        # joining like this is sort of problematic as it locks main thread
        # better to use some kind of signal.
        for p in zombies:
            p.join()
        
        # We are finished        
        self.doneStatus()
    
        self.lock("parsing results")
        for ff in self.headerstr:
            tag = q.get() #['tag']) #q.get()
            Dict = pickle.load( open( str(tag)+".p", "rb" ) )
            #print (Dict)
            self.ORS[Dict['tag']].DistRes = Dict
        self.plotDISTPar() 

    def plotDISTPar(self):

        self.plotDIST( self.ui.maskNSpinBox_4.value() )
        self.ui.qplot_2.draw()  
        self.ui.qplot_3.draw()  
         
        if (self.DepthProfile == True):
            dp = []
            mod = []
            for ff in self.headerstr:
                dp.append( self.ORS[ff].depth )
                mod.append( self.ORS[ff].DistRes['mod'] )
            #self.msp0.plot(self.ORS[ff].DistRes['Time'].T2Bins, self.ORS[ff].DistRes['mod'], linewidth=2, label="recovered")
            X,Y = meshgrid(  self.ORS[ff].DistRes['Time'].T2Bins, dp )   
            #self.mp.matshow(Data, aspect = 'auto')
            self.mp.pcolor(X,Y, np.array(mod), cmap = 'hot_r')
            #self.mp.set_ylim( self.mp.get_ylim()[::-1] )
            self.ui.qplot_4.draw()
        
    def cellChanged(self):
        # TODO consider building the model whenever this is called. Would be nice to be able to 
        # do that. Would require instead dist of T2 I guess. 
        jj = self.ui.modelTableWidget.currentColumn()
        ii = self.ui.modelTableWidget.currentRow()
            
        #if self.ui.modelTableWidget.item(ii, jj) == None:
        #    return

        try:
            eval (str( self.ui.modelTableWidget.item(ii, jj).text() ))
        except:
            #print  ("cell changed ERROR",   str( self.ui.modelTableWidget.item(ii, jj).text() ) )
            #Error = QtGui.QErrorMessage()
            Error = QtGui.QMessageBox()
            #Error.showMessage("NOT A NUMERIC VALUE")
            Error.setWindowTitle("Error!")
            Error.setText("Non-numeric value encountered")
            Error.setDetailedText("Modelling parameters must be able to be cast into numeric values.")
            Error.exec_()

    def tc(self):
        """ Computes the permeabiliy based on Timur-Coates equation 
            k = c phi^m (ffv/bfv)^n  
            k = permeability
            phi = nmr porosity
            m = porosity exponent, 2-4 is standard. 
            ffv = free fluid volume (T2 above cutoff)
            bfv = bound fluid volume (T2 below cutoff)
            n = fluid exponent 2 is standard 
            c = prefactor
            cutoff = bfv : ffv cutoff, 33 ms is standard 
        """
        if not self.ui.saveDatBox_3.isChecked():
            self.ui.qplot_5.getFigure().clf()
            
            #self.dd = self.ui.qplot_5.getFigure().add_subplot(141)
            #self.kp = self.ui.qplot_5.getFigure().add_subplot(142, sharey = self.dd)
            #self.kt = self.ui.qplot_5.getFigure().add_subplot(143, sharey = self.dd)
            #self.kk = self.ui.qplot_5.getFigure().add_subplot(144, sharey = self.dd)
            #self.ui.qplot_5.getFigure().tight_layout() # TODO tweak spacing tight is almost good. 
            #plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
            
            # add_axes()
            self.dd = self.ui.qplot_5.getFigure().add_axes( [  .075, .15, .2, .8] , facecolor='black' )
            self.kp = self.ui.qplot_5.getFigure().add_axes( [  .300, .15, .2, .8] , sharey=self.dd )
            self.kt = self.ui.qplot_5.getFigure().add_axes( [  .525, .15, .2, .8] , sharey=self.dd )
            self.kk = self.ui.qplot_5.getFigure().add_axes( [  .750, .15, .2, .8] , sharey=self.dd )
            
            self.dc = self.ui.qplot_5.getFigure().add_axes( [  .075, .05, .2, .02]                  )

            #self.kkp = self.kp.twiny()
            self.kk.set_xlabel(r"$\kappa$ [mD]") 
            self.kp.set_xlabel(r"$\phi$ [m$^3$ / m$^3$]") 
            self.kt.set_xlabel(r"$T_2$ [s]") 
            self.dd.set_xlabel(r"time [s]") 
            self.dd.set_ylabel(r"depth [m]")
            #self.kt.set_ylabel(r"depth [m]")
            #self.kk.set_ylabel(r"depth [m]")
            
            # permeability and decay time in log10
            self.kk.set_xscale('log', basex=10) 
            self.kt.set_xscale('log', basex=10) 
            
            self.dd.set_title(r"CPMG time series")
            self.kp.set_title(r"NMR water (c,b,f)")
            self.kt.set_title(r"NMR $T_2$")
            self.kk.set_title(r"$\kappa_{NMR}$")


            #self.kt.xaxis.set_major_locator(MaxNLocator(5)) 
            #self.kk.xaxis.set_major_locator(MaxNLocator(5)) 

       
        c = self.ui.TC_c.value() 
        m = self.ui.TC_m.value()
        n =  self.ui.TC_n.value()
        cutoff = self.ui.TC_cutoff.value() * 1e-3
        cutoff2 = .003
            
        dp = []
        phi = []
        FFV = []
        BFV = []
        CBFV = []
        LogMeanT2 = []       
        ii = 0

        dat = open("NMR-log.dat","w")        
        dat.write( "#depth,phi,T_2ML,FFV,BFV,CBFV\n" )
        
        results = file("inversion.csv",'w') 

        Data = np.zeros( (len(self.ORS), len(self.ORS[0].T2T)) )
        
        times = self.ORS[0].T2T
        depTol = np.abs(self.ORS[self.headerstr[0]].depth - self.ORS[self.headerstr[-1]].depth) 
        t2scale = 100*depTol / len(self.ORS)  #TODO adjust based on total depth, was 20

        for ff in self.headerstr:
            nT2 = len(self.ORS[ff].DistRes['Time'].T2Bins)
            dp.append( self.ORS[ff].depth )
            phi.append( np.sum( self.ORS[ff].DistRes['mod'] ) )
            LogMeanT2.append( np.exp(np.sum( self.ORS[ff].DistRes['mod'] * np.log( self.ORS[ff].DistRes['Time'].T2Bins) ) / phi[ii] ) )     
            FFV.append( np.sum( self.ORS[ff].DistRes['mod'][ self.ORS[ff].DistRes['Time'].T2Bins>=cutoff])  )       
            BFV.append( np.sum( self.ORS[ff].DistRes['mod'][  np.logical_and( self.ORS[ff].DistRes['Time'].T2Bins<cutoff, self.ORS[ff].DistRes['Time'].T2Bins>=cutoff2) ] + 1e-6 )  ) # add small epsilon value       
            CBFV.append( np.sum( self.ORS[ff].DistRes['mod'][ self.ORS[ff].DistRes['Time'].T2Bins<cutoff2]) + 1e-6 )  # epsilon  
            self.kt.plot( self.ORS[ff].DistRes['Time'].T2Bins, dp[ii] - t2scale*self.ORS[ff].DistRes['mod'] , color='blue' ) 
            self.kt.fill_between( self.ORS[ff].DistRes['Time'].T2Bins, dp[ii] - t2scale*self.ORS[ff].DistRes['mod'], dp[ii]*np.ones( nT2 ), alpha=1, color='LimeGreen')
            Data[ff,:] = self.ORS[ff].T2D.imag

            for iid, val in enumerate(self.ORS[ff].DistRes['mod']):
                if iid == 0:
                    for it2, t2b in enumerate( self.ORS[ff].DistRes['Time'].T2Bins ):
                        if it2 == 0:
                            results.write( "%1.4f" %(t2b)  )
                        else:
                            results.write( ",%1.4f" %(t2b)  )
                    results.write("\n") 
                    results.write( "%1.4f" %(val)  )
                else:
                    results.write( ",%1.4f" %(val)  )
            results.write("\n") 

            dat.write( "%1.4f,%1.4f,%1.4f,%1.4f,%1.4f,%1.4f\n" %(dp[-1],phi[-1],LogMeanT2[-1],FFV[-1],BFV[-1],CBFV[-1]))
            #if np.abs( phi[-1] - (FFV[-1] + BFV[-1] + CBFV[-1]) ) > 1e-6:
            #    print("mismatch", phi[-1] - (FFV[-1] + BFV[-1] + CBFV[-1]), CBFV[-1] ) 
            #else:
            #    print("nope", CBFV[-1])
            ii += 1
        results.close()
        dat.close()
        #k =  ((c*np.array(phi))**m) * ( (np.array(FFV)/np.array(BFV))**n)
        k = c * 1e4 *  (np.array(phi)**m) * ((np.array(FFV)/(np.array(BFV)+np.array(CBFV)))**n)

        ############################################################
        # porosity logs        
        
        # TI, TODO these fill_between cause errors with Javelin data? Why?
        #self.kp.plot( phi, dp , label = r"$\phi_{NMR}$", color='blue' )

        #try:
        #    self.kp.fill_betweenx(np.array(dp), np.array(phi), alpha=1., color='blue' ) #, label = r"$BFV_{TC}$" )
        #except:
        #    print ("shapes of dp and phi", np.shape(np.array(dp)), np.shape(np.array(phi)))
        FFV = np.array(FFV)
        BFV = np.array(BFV)
        CBFV = np.array(CBFV)

        # Plot cumulitive of all pore sizes below, that way the visible portion corresponds with the PWC
        self.kp.plot( FFV+BFV+CBFV , dp , label = r"$FFV_{TC}$", color='blue' ) # don't plot implicit
        self.kp.fill_betweenx(dp, FFV+BFV+CBFV, alpha=1, color='blue' ) #, label = r"$BFV_{TC}$" )
       
        self.kp.plot( BFV+CBFV, dp , label = r"$BFV_{TC}$" , color='green')
        self.kp.fill_betweenx(dp, BFV+CBFV, alpha=1, color='green' )  #, label = r"$BFV_{TC}$" )
        
        self.kp.plot( CBFV, dp , label = r"$CBFV_{TC}$" , color='brown')
        self.kp.fill_betweenx(dp, CBFV, alpha=1, color='brown' )  #, label = r"$BFV_{TC}$" )

        #ax.fill_between(xs, ys, where=ys>=d, interpolate=True, color='blue')
        self.kp.set_ylim( [ max(self.kp.get_ylim()), min(self.kp.get_ylim()) ]) 
        self.kp.set_xlim( [ 0, self.kp.get_xlim()[1] ] )
        
        ####################################################################################
        # Data Plot       
        #self.dd.matshow()
        #X,Y = meshgrid( np.log10(times), dp )  
        # pcolor cuts one spot, need to extrapolate
        ma = self.ui.maskNSpinBox_4.value()
 
        dp2 = np.append( dp, dp[-1]+(dp[-1]-dp[-2]) )
        times2 = np.append( times, times[-1]+(times[-1]-times[-2]) )
        times2 = times2[ma:]
        #times2 = np.append( times[ma:], times[ma:-1]+(times[-1]-times[-2]) )
        dp2 = dp2 - (dp[1]-dp[0]) *.5 

        X,Y = meshgrid( times2, dp2 )   
        #self.mp.matshow(Data, aspect = 'auto')
        #pc = self.dd.pcolor(X,Y, Data, cmap = 'coolwarm', vmin = 0)
        pc = self.dd.pcolor(X,Y, Data[:,ma:], cmap = 'viridis') #, vmin = -np.max(Data[:,ma:]))
        #pc = self.dd.pcolor(X,Y, Data, cmap = 'seismic_r') #, vmin=-.5, vmax=.5) #, vmin = 0)
        #self.dd.set_ylim( self.dd.get_ylim()[::-1] )

        np.savez("cmrdata",times2=times2,dp2=dp2,X=X, Y=Y, data=Data[:,ma:])       
 
        self.dd.set_ylim( np.max(dp2), np.min(dp2) )
        self.dd.set_xlim( times2[0], times2[-1] )
        #self.dd.set_xlim( np.log10(times[0]), np.log10(times[-1]) )

        self.dd.set_xscale("log", nonposx='clip')           
 
        ############################################################
        # decay logs
        #self.kt.plot( LogMeanT2, dp , label = r"$T_{2ML}$" )
        #self.kt.set_ylim( [ max(self.kp.get_ylim()), min(self.kp.get_ylim()) ]) # use same range as others 
        
        # Don't plot huge K's 
        self.kk.set_xlim( [ min(self.kk.get_xlim()),  1e6 ]) # use same range as others 
        
        # permeability logs
        self.kk.plot( k, dp , label = r"$\kappa_{TC}$", color='red' )
        #self.kk.set_ylim( [ max(self.kk.get_ylim()), min(self.kk.get_ylim()) ]) 
        self.kk.set_xlim( [ 1e-1, self.kk.get_xlim()[1] ] )
        self.kk.grid() 

        # don't plot y axis labels, except on data plot (far left) 
        #self.kk.axes.get_yaxis().set_ticks([]) 
        #self.kp.axes.get_yaxis().set_ticks([]) 
        #self.kt.axes.get_yaxis().set_ticks([]) 
        # above kills dd labels too
        plt.setp( self.kk.axes.get_yticklabels(), visible=False)     
        plt.setp( self.kp.axes.get_yticklabels(), visible=False)     
        plt.setp( self.kt.axes.get_yticklabels(), visible=False)     

        #h1, l1 = self.kp.get_legend_handles_labels()
        #h2, l2 = self.kkp.get_legend_handles_labels()
        #self.kp.legend(h1+h2, l1+l2, numpoints=1)
       
        # try to avoid legends, they detract from logs 
        #self.kp.legend() 
        #self.kk.legend() 
        #self.kt.legend() 
            
        self.kp.xaxis.set_major_locator(MaxNLocator(3)) 
        #self.dd.xaxis.set_major_locator(MaxNLocator(3)) 

        cb = self.ui.qplot_5.getFigure().colorbar(pc, self.dc, orientation='horizontal')
        cb.solids.set_edgecolor("face")
        #tick_locator = MaxNLocator(nbins=5)
        cb.locator = MaxNLocator(5) #tick_locator
        cb.update_ticks() 

        self.ui.qplot_5.draw()
    
    def sdr(self):
        """ Estimates permeability based on Schlumberger-Doll Research (SDR) equation
                \kappa = c \phi^m T_{2ML}^n
        """
        if not self.ui.saveDatBox_3.isChecked():
            self.ui.qplot_5.getFigure().clf()
            self.kp = self.ui.qplot_5.getFigure().add_subplot(131)
            self.kt = self.ui.qplot_5.getFigure().add_subplot(132)
            self.kk = self.ui.qplot_5.getFigure().add_subplot(133)
 
            self.kk.set_xlabel(r"$\kappa$ [mD]") 
            self.kp.set_xlabel(r"$\phi$ [m$^3$ / m$^3$]") 
            self.kt.set_xlabel(r"$T_2$ [s]") 
            self.kp.set_ylabel(r"depth [m]")
            #self.kt.set_ylabel(r"depth [m]")
            #self.kk.set_ylabel(r"depth [m]")
            
            # permeability and decay time in log10
            self.kk.set_xscale('log', basex=10) 
            self.kt.set_xscale('log', basex=10) 
            
            self.kp.set_title(r"NMR water (b,t)")
            self.kt.set_title(r"NMR $T_2$")
            self.kk.set_title(r"$\kappa_{NMR}$")

            self.kp.xaxis.set_major_locator(MaxNLocator(4)) 
            #self.dd.xaxis.set_major_locator(MaxNLocator(4)) 
            #self.kt.xaxis.set_major_locator(MaxNLocator(5)) 
            #self.kk.xaxis.set_major_locator(MaxNLocator(5)) 

            #self.kp.set_title("permeability")
        
        c = self.ui.SDR_c.value() 
        m = self.ui.SDR_m.value()
        n =  self.ui.SDR_n.value()
        
        dp = []
        phi = []
        LogMeanT2 = []
        ii = 0
        for ff in self.headerstr:
            dp.append( self.ORS[ff].depth )
            phi.append( np.sum( self.ORS[ff].DistRes['mod'] ) )
            #theta = np.sum( self.ORS[ff].DistRes['mod'] )
            LogMeanT2.append( np.exp(np.sum( self.ORS[ff].DistRes['mod'] * np.log( self.ORS[ff].DistRes['Time'].T2Bins) ) / phi[ii] ) )        
            ii += 1        
        
        #print ("LogMeanT2", LogMeanT2) 
        # 1e3 converts to ms where c is designed for 
        k = c*np.array(phi)**m * ((1e3 * np.array(LogMeanT2))**n)

        self.kp.plot( phi, dp , label = r"$\phi_{NMR}$" )
        self.kt.plot( LogMeanT2, dp , linewidth=2, color='blue', label = r"$T_{2ML}$" )
        self.kk.plot( k, dp , label = r"$\kappa_{SDR}$", color='blue' )
        
        self.kp.set_ylim( [ max(self.kp.get_ylim()), min(self.kp.get_ylim()) ])
        self.kk.set_ylim( [ max(self.kp.get_ylim()), min(self.kp.get_ylim()) ])
        self.kt.set_ylim( [ max(self.kp.get_ylim()), min(self.kp.get_ylim()) ])

        #h1, l1 = self.kp.get_legend_handles_labels()
        #h2, l2 = self.kkp.get_legend_handles_labels()
        #self.kp.legend(h1+h2, l1+l2, numpoints=1)
            
        #self.kkp.grid() 
        self.ui.qplot_5.draw()

    def model(self):
        from pylab import meshgrid

        # TODO add noise std column which will be interesting for resolution analysis 
        #self.ui.modelTableWidget.setColumnCount(4)      
        #for ii in range( self.ui.modelTableWidget.rowCount() ):
        #    for jj in range(  self.ui.modelTableWidget.columnCount() ):
        #        if self.ui.modelTableWidget.item(ii, jj) != None:
        #            print( str(ii) +  "\t" + str(jj) + "\t" +  str( self.ui.modelTableWidget.item(ii, jj).text() ) )

        self.DepthProfile = True

        Bottom =  self.ui.doubleSpinBoxBottom.value() 
        Top = self.ui.doubleSpinBoxTop.value() 
        Dz = self.ui.doubleSpinBoxDz.value() 
        T = self.ui.doubleSpinBoxTTrain.value() 
        TauE = 1e-3 * self.ui.doubleSpinBoxTauE.value()  
        Mean = self.ui.doubleSpinBoxMean.value() 
        Std = self.ui.doubleSpinBoxStd.value() 

        times = np.arange(TauE, T, TauE)
        dp = np.arange(Bottom, Top, -Dz)
        Data = np.zeros( (len(dp), len(times)) )
        Noise = np.zeros( (len(dp), len(times)) )
        Mod = np.zeros( (len(dp) )) + 1e-6
        Por = np.zeros( (len(dp) )) 
        Por = np.zeros( (len(dp) )) 
        for ii in range( self.ui.modelTableWidget.rowCount() ):
            a = type(self.ui.modelTableWidget.item(ii, 0))
            if str(a) == "<class 'NoneType'>":  # ugly hack needed by PySide for some strange reason.
                pass #print ("NONE")
            else: #if type(self.ui.modelTableWidget.item(ii, 0)) != None:
                #print ("type", type(self.ui.modelTableWidget.item(ii, 0)))
                for iid in range(len(dp)):
                    #print self.ui.modelTableWidget.item(ii, 0).text(),  dp[iid],  self.ui.modelTableWidget.item(ii, 1).text()
                    try:
                        if dp[iid] > eval(str(self.ui.modelTableWidget.item(ii, 0).text()) ) and dp[iid] <= eval(str(self.ui.modelTableWidget.item(ii, 1).text()) ):
                            Mod[iid] = 1e-3 * eval( str(self.ui.modelTableWidget.item(ii, 2).text() ))
                            Por[iid] = 1e-2 * eval(str(self.ui.modelTableWidget.item(ii, 3).text() ))
                            # Append the data for multi-exponential  
                            Data[iid, :] += Por[iid] * np.exp( -times/ Mod[iid] ) + np.random.normal( 0, eval(str( self.ui.modelTableWidget.item(ii,4).text())), len(times) ) 
                            Noise[iid,:] += np.random.normal( 0, eval(str( self.ui.modelTableWidget.item(ii,4).text())), len(times) ) 
                    except:
                        pass
        self.ORS.clear() # = {} 
        self.headerstr = [] #.clear() # = []
        for iid in range(len(dp)):
            self.headerstr.append( iid )
            self.ORS[iid] = MRProc(  )  # Fast
            self.ORS[iid].NE = len(times)
            self.ORS[iid].TAUE = TauE 
            self.ORS[iid].NS = 1  # TODO consider what to do here  
            self.ORS[iid].T2T = times 
            self.ORS[iid].T2D = Noise[iid,:] +  1j*Data[iid,:]
            self.ORS[iid].T2N = Noise[iid,:] 
            self.ORS[iid].sigma = np.std( Noise[iid,:] ) 
            self.ORS[iid].depth = dp[iid] 
            
            # Doom, why
            #QtCore.QObject.connect(self.ORS[iid], QtCore.SIGNAL("enableDSP()"), self.enableDSP)
            #QtCore.QObject.connect(self.ORS[iid], QtCore.SIGNAL("enableINV()"), self.enableINV)
            #QtCore.QObject.connect(self.ORS[iid], QtCore.SIGNAL("doneStatus()"), self.doneStatus)
       

        # enable gate integrate
        self.enableDSP() 
        self.ui.bulkProcessGroupBox.setEnabled(True)
        self.ui.batchLoadDataPushButton.setEnabled(True)
        self.ui.gateGroupBox.setEnabled(True)

        # disable buttons, do flow 
        
        self.ui.WindowStackBox.setChecked(False)   
        self.ui.envelopeGroupBox.setChecked(False)
        self.ui.phaseGroupBox.setChecked(False)
        self.ui.gateGroupBox.setChecked(True)
        
        #self.ui.WindowStackGO.checked(False)
        #self.ui.envelopeGO.checked(False)
        #self.ui.phaseGO.checked(False)
        self.ui.gateGO.setEnabled(False)

        # Enable Inversions
        self.enableINV()
 
        #for iz in range(len(dp)):
            #Data[iz,:]  = 1e-2*dp[iz] + np.exp(-times/.324)
        #    Data[iz,:]  = Por[iz] * np.exp( -times/ Mod[iz] )

        # OK, so now I need to get these into data container class
 
        # for now just make data and plot it
        if not self.ui.saveDatBox_2.isChecked():
            self.ui.qplot_4.getFigure().clf()
            self.mp = self.ui.qplot_4.getFigure().add_subplot(111) 
            self.mp.set_title("synthetic data")

        X,Y = meshgrid( times, dp )   
        #self.mp.matshow(Data, aspect = 'auto')
        self.mp.pcolor(X,Y, Data, cmap = 'hot_r')
        self.mp.set_ylim( self.mp.get_ylim()[::-1] )
        self.ui.qplot_4.draw()

    def jgFit(self):
        T2Bw = 1e-3 * self.ui.T2bwSpinBox.value()
        #T2Bw = 2.
        Times = {}
        Times2 = np.empty(0)
        TauE = []
        Data = np.empty(0) #{}
        Datap = np.empty(0) #{}
        Sigma = np.empty(0) #{}
        mask = self.ui.maskNSpinBox_5.value()
        for ff in self.headerstr:
            Times[self.ORS[ff].TAUE] =  self.ORS[ff].T2T[mask::]   
            Times2 = np.append( Times2, self.ORS[ff].T2T[mask::] )  
            TauE.append(self.ORS[ff].TAUE)
            Data = np.append( Data,  np.imag(self.ORS[ff].T2D)[mask::] )
            Datap = np.append( Datap,  np.imag(self.ORS[ff].T2D)[mask::] / np.exp(-self.ORS[ff].T2T[mask::]/(self.ui.A0SpinBox.value()*T2Bw) ) )
            Sigma = np.append( Sigma,  self.ORS[ff].sigma[mask::] )
 
        # Plot up data, TODO need a separate thread for this, unresponsive GUI 
        #self.plotDIFF(Times2, Data, Sigma)
        #self.plotDIFF2(Times2, Data, Datap, Sigma)

        #############################################
        # Generate Kernel for inversion to J(G)
        Times = {}
        Times2 = np.empty(0)
        TauE = []
        Data = np.empty(0) #{}
        Sigma = np.empty(0) #{}
        mask = self.ui.maskNSpinBox_5.value()
        for ff in self.headerstr:
            Times[self.ORS[ff].TAUE] = self.ORS[ff].T2T[mask::]   
            Times2 = np.append( Times2, self.ORS[ff].T2T[mask::] )  
            TauE.append(self.ORS[ff].TAUE)
            Data = np.append( Data,  np.imag(self.ORS[ff].T2D)[mask::] )
            Sigma = np.append( Sigma,  self.ORS[ff].sigma[mask::] )
        self.GBins = np.linspace(self.ui.gminSpinBox.value(), self.ui.gmaxSpinBox.value(), num=self.ui.ngSpinBox.value(), endpoint=True) 
        #self.GBins = np.logspace(np.log10(self.ui.gminSpinBox.value()), np.log10(self.ui.gmaxSpinBox.value()), num=self.ui.ngSpinBox.value(), endpoint=True) 
        INV = pwcTime()        
 
        INV.generateJGKernel(TauE, Times, 1e-5*self.ui.dwSpinBox.value() , self.GBins, self.ui.DTau.value())
        self.jgModel  = logBarrier(INV.JG, Datap, MAXITER=500, sigma=Sigma, alpha=1e20) #, smooth=True) # 

        pre = np.dot(INV.JG, self.jgModel)
        
        #self.plotDIFF(Times2, Data, Sigma, Datap) #pre)
        self.plotDIFF2(Times2, Data, Datap, Sigma, pre)

        # Plot model
        self.ui.qplot_7.getFigure().clf()
        self.msp0 = self.ui.qplot_7.getFigure().add_axes([.15,.2,.75,.65]) 

        self.msp0.plot(self.GBins, self.jgModel, color='black', linewidth=3)
        self.msp0.grid()
        self.msp0.set_xlabel("$J(G)$ [G/cm]")
        self.msp0.set_ylabel("intensity")
        self.ui.qplot_7.draw()  
        #plt.show()

    def diffusion(self):
        from pylab import meshgrid
        #print ("Diffusion Inversion")
        Times = {}
        Times2 = np.empty(0)
        TauE = []
        Data = np.empty(0) #{}
        Sigma = np.empty(0) #{}
        mask = self.ui.maskNSpinBox_5.value()
        for ff in self.headerstr:
            Times[self.ORS[ff].TAUE] = self.ORS[ff].T2T[mask::]   
            Times2 = np.append( Times2, self.ORS[ff].T2T[mask::] )  
            TauE.append(self.ORS[ff].TAUE)
            Data = np.append( Data,  np.imag(self.ORS[ff].T2D)[mask::] )
            Sigma = np.append( Sigma,  self.ORS[ff].sigma[mask::] )


 
        self.INV = pwcTime()        
        #self.D = np.linspace(1e-9*self.ui.Dmin.value(), 1e-5*self.ui.Dmax.value(), num=self.ui.nDBins.value(), endpoint=True)  
        self.D = np.logspace(np.log10(1e-9*self.ui.Dmin.value()), np.log10(1e-5*self.ui.Dmax.value()), num=self.ui.nDBins.value(), endpoint=True)  
        
        self.INV.setT2( 1e-3*self.ui.lowT2SpinBox.value(), 1e-3*self.ui.hiT2SpinBox.value(), self.ui.nT2SpinBox.value(), self.ui.distT2Box.currentText() )       
        
        if self.ui.jgComboBox.currentText() == "Distribution":
            print("Using J(G) Distribution")
            self.INV.generateGDenv( TauE, Times, self.D, self.GBins, self.jgModel, self.ui.DTau.value() )
        else:
            print("Using constant G")
            self.INV.generateGDenv( TauE, Times, self.D, np.array([self.ui.jgcSpinBox.value()]), np.array([1.]) )
            

        # Plot up data, TODO need a separate thread for this, unresponsive GUI 
        self.plotDIFF(Times2, Data, Sigma)

        # OK now invert TODO, need to speed this up. Sparse matrices and/or implicit ATWdTWdA? 
        model  = logBarrier(self.INV.GD, Data, MAXITER=500, sigma=Sigma, alpha=1e10, callback=None) #self.plotDmod) #, smooth=True) # 
        pre = np.dot(self.INV.GD, model)
        
        self.plotDIFF(Times2, Data, Sigma, pre)
        self.plotDmod(model)

    def plotDmod(self, model):

        # plot models and fit
        self.DMod = np.reshape(model, (len(self.D), len(self.INV.T2Bins))) 

        # Save to file for analysis
        print("saving D mod to dmod.bin")
        res = open('dmod.bin', 'wb')
        np.save(res, self.D)
        np.save(res, self.INV.T2Bins)
        np.save(res, model)
    
        #if not self.ui.saveModBox.isChecked():
        self.ui.qplot_7.getFigure().clf()
        self.msp0 = self.ui.qplot_7.getFigure().add_axes([.15,.2,.65,.65]) 
        self.csp0 = self.ui.qplot_7.getFigure().add_axes([.85,.2,.025,.65]) 
        #self.msp0.set_xlabel(r"$T_2$ [s]", color='black')
        #self.msp0.set_ylabel(r"$A_0$ [rku]", color='black') 

        #for ff in self.headerstr: 
        #X,Y = meshgrid( np.concatenate( (np.array([.999*self.INV.T2Bins[0]]), self.INV.T2Bins) ), \
        #                np.concatenate( (np.array([.999*self.D[0]]), self.D)) ) # pcolor needs an extra bin  
        X,Y = meshgrid( self.INV.T2Bins, self.D ) # pcolor needs an extra bin  
        #cmap = plt.get_cmap('hot_r') #'
        imsa = self.msp0.pcolor(X, Y, self.DMod, cmap='hot_r', edgecolors='w') #, rasterized=True)
        self.msp0.set_xlabel("$T_2$ [s]") 
        self.msp0.set_ylabel("D  [cm$^2$/s]")
        self.msp0.set_yscale('log', basey=10) 
        self.msp0.set_xscale('log', basex=10)

        plt.colorbar(imsa, self.csp0) #self.msp0)

        #cb1 = mpl.colorbar.ColorbarBase(self.cbax, cmap=wmap,
        #                               norm=nnorm,
        #                               orientation='horizontal', 
        #                               extend='both',
        #                               format = r'%i')
 
        #self.msp0.imshow( self.DMod, interpolation='bicubic' ) #self.ORS[ff].MonoRes['rt20'], self.ORS[ff].MonoRes['b0'], 'o', markersize=6, label="mono")   
        #    self.msp0.legend(numpoints=1)
        #plt.ylabel("D [cm$^2$/s]")
        #plt.xlabel("$T_2$ [s]")
        self.ui.qplot_7.draw()   

        #print ("FINISHED DIFFUSION")

    def gate(self, ff=None):

        #print ("GATING")

        if ff == None:
            ff = self.headerstr[0]

        self.ui.logTextBrowser.append( "Gate integrate correct" + "GPD") 
        self.lock("time gating")
        rawThread = thread.start_new_thread(self.ORS[ff].gateIntegrate, (self.ui.gpdSpinBox.value(),self.ui.stackEfficiency.value()) ) 
        #rawThread = thread.start_new_thread(self.ORS[ff].CorrectedAmplitude, (0,0) ) 
        #self.ui.logTextBrowser.append( "Envelope Detect: " + str(self.ui.fTypeComboBox.currentText() ))

    def openCMRLog(self):
        print("CMR log process")
        import scipy.io as sio
        
        try:
            with open('.akvo.cmr.last.path') as f: 
                fpath = f.readline()  
                pass
        except IOError as e:
            fpath = '.'

        self.DepthProfile == True
        self.burst = False # True
        self.headerstrRAW = [QtGui.QFileDialog.getOpenFileName(self, 'Open Schlumberger CMR Log Series', fpath,  r"CMR MATLAB (*.mat)" )]
        self.headerstr = [os.path.normpath( str(self.headerstrRAW[0][0]) )]
    
        self.ui.headerFileTextBrowser.clear()
        self.ui.headerFileTextBrowser.append( self.headerstr[0] )
        
        # clear the processing log? TODO do we want to do this? Not sure what about multiple inputs?
        self.ui.logTextBrowser.clear()

        path,filen=os.path.split(str(self.headerstr[0] ))
        fileName, fileExtension = os.path.splitext( str(self.headerstr[0]) )

        f = open('.akvo.cmr.last.path', 'w')
        if str(self.headerstr[0]):
            f.write( str(self.headerstr[0]) ) # prompt last file      

        self.ui.headerFileBox.setEnabled(True)

        #########################################
        # LOAD the data and do some basic work
        self.ORSConnect( str(self.headerstr[0]) )       
        CMR_PHI_CONV = 0.00068 # this is temperature dependent, should be in files, but wasn't!  

        self.lock("loading RAW record")
        matfile = self.headerstr[0] 
        if fileExtension == ".mat": 
            
            Data =  sio.loadmat(self.headerstr[0])

            # the Schlumberger data has a burst mode at the end of the record which is designed only to capture 
            # fast decays. I don't see a record of the number in the files, is often 30          
            ne = np.shape(Data["echo_amp_x"].T[:,0])[0] - 30
            nb = 30  # number of burst echoes 

            T2T = 2e-4 * np.arange(1, ne+1, 1) # .2 ms seems to be CMR standard, not in data files
            T2Tb = 2e-4 * np.arange(1, nb+1, 1) # burst times 
#            T2T = np.concatenate( (T2T,T2Tb) )

            TAUE =  T2T[3] - T2T[2]
            NS = 1
            self.ORS.clear()     # = {} 
            self.headerstr = []  #.clear() # = []
            
            # TODO query for depth step and data quality factor 
            iiid = 0
            
            #################################
            dmin = 7672  # Morrow           #
            dmax = 7707  # Morrow Bottom    #
            #################################
            #dmin = 7672 - 150  # Morrow           #
            #dmax = 7707 - 120  # Morrow Bottom    #
            #dmin = 7672+10.25  # Morrow           #
            #dmax = 7707-10.25  # Morrow Bottom    #
            #dmin = 0 
            #dmax = 7510
            #dmax = 10000
            for iid in range(len(Data['depth'])):
                dep =  float(Data['depth'][iid])
                #if True:
                if dep > dmin and dep <= dmax:
                    self.headerstr.append( iiid )
                    self.ORS[iiid] = MRProc(  )  # Fast
                    self.ORS[iiid].NE = ne #len(Data['time'][:,0])
                    self.ORS[iiid].TAUE = TAUE 
                    self.ORS[iiid].NS = 1  # TODO consider what to do here  
                    self.ORS[iiid].T2T  = T2T
                    self.ORS[iiid].T2Tb = T2Tb
                    self.ORS[iiid].burst = False #True
                
                    # long record
                    self.ORS[iiid].T2D =  CMR_PHI_CONV*Data["echo_amp_r"].T[0:ne,iid] - CMR_PHI_CONV*1j*Data["echo_amp_x"].T[0:ne,iid] 
                    
                    # burst data
                    self.ORS[iiid].T2Db =  CMR_PHI_CONV*Data["echo_amp_r"].T[ne:,iid] - CMR_PHI_CONV*1j*Data["echo_amp_x"].T[ne:,iid] 
                    
                    # both
                    #self.ORS[iiid].T2D =  CMR_PHI_CONV*Data["echo_amp_r"].T[:,iid] - CMR_PHI_CONV*1j*Data["echo_amp_x"].T[:,iid] 
                    
                    self.ORS[iiid].depth = float(Data['depth'][iid]) * 0.3048 # Schlumberger logs in feet
                    iiid += 1
           
            self.enableDSP() 
            
            # disable buttons, do flow 
            self.ui.WindowStackGO.setEnabled(False)
            self.ui.envelopeGO.setEnabled(False)
            self.ui.phaseGO.setEnabled(False)
        
            self.ui.WindowStackBox.setChecked(False)   
            self.ui.envelopeGroupBox.setChecked(False)
            self.ui.phaseGroupBox.setChecked(True)
            self.ui.gateGroupBox.setChecked(True)
            self.ui.gateGroupBox.setEnabled(True)
        
            self.ui.gateGO.setEnabled(False)
            self.ui.bulkProcessGroupBox.setEnabled(True)
            self.ui.batchLoadDataPushButton.setEnabled(True)
        
        self.ui.logTextBrowser.append( "opened header file " + matfile )

    def openVCLogSeries(self):
        import scipy.io as sio

        try:
            with open('.akvo.last.vcl.path') as f: 
                fpath = f.readline()  
                pass
        except IOError as e:
            fpath = '.'
        
        self.DepthProfile == True

        self.headerstrRAW = [QtGui.QFileDialog.getOpenFileName(self, 'Open VistaClara Log Series', fpath,  r"CMR (*.mat)" )]
        self.headerstr = [os.path.normpath( str(self.headerstrRAW[0][0]) )]
    
        self.ui.headerFileTextBrowser.clear()
        self.ui.headerFileTextBrowser.append( self.headerstr[0] )
        
        # clear the processing log? TODO do we want to do this? Not sure what about multiple inputs?
        self.ui.logTextBrowser.clear()

        path,filen=os.path.split(str(self.headerstr[0] ))
        fileName, fileExtension = os.path.splitext( str(self.headerstr[0]) )

        f = open('.akvo.last.vcl.path', 'w')
        if str(self.headerstr[0]):
            f.write( str(self.headerstr[0]) ) # prompt last file      

        self.ui.headerFileBox.setEnabled(True)

        #########################################
        # LOAD the data and do some basic work
        self.ORSConnect( str(self.headerstr[0]) )       

        self.lock("loading RAW record")
        matfile = self.headerstr[0] #.clear() # = []
        if fileExtension == ".mat": 
            
            Data =  sio.loadmat(self.headerstr[0])
            T2T = Data['time'][:,0] 
            TAUE =  T2T[3] - T2T[2]
            NS = 1
            self.ORS.clear() # = {} 
            self.headerstr = [] #.clear() # = []
           
            # TODO query for depth step and data quality factor 
 
            for iid in range(len(Data['depth'])):
                self.headerstr.append( iid )
                self.ORS[iid] = MRProc(  )  # Fast
                self.ORS[iid].NE = len(Data['time'][:,0])
                self.ORS[iid].TAUE = TAUE 
                self.ORS[iid].NS = 1  # TODO consider what to do here  
                self.ORS[iid].T2T = T2T
                # What about moving window std.dev?  
                self.ORS[iid].sigma =  np.std(1e-2*Data['se_vector_wc'][:,iid][-100::]) # slowly varying at end?  
                self.ORS[iid].T2N = np.random.normal( 0, self.ORS[iid].sigma, len(T2T) )    
                self.ORS[iid].T2D = self.ORS[iid].T2N +  1j*1e-2*Data['se_vector_wc'][:,iid]

                #self.ORS[iid].depth = Data['depth'][iid][0] 
                # TODO pass into GUI for weird Javelin data that is not recorded correctly
                self.ORS[iid].depth = 2.25- .25 * Data['depth'][iid][0] 
               
           
            self.enableDSP() 
            
            # disable buttons, do flow 
            self.ui.WindowStackGO.setEnabled(False)
            self.ui.envelopeGO.setEnabled(False)
            self.ui.phaseGO.setEnabled(False)
        
            self.ui.WindowStackBox.setChecked(False)   
            self.ui.envelopeGroupBox.setChecked(False)
            self.ui.phaseGroupBox.setChecked(False)
            self.ui.gateGroupBox.setChecked(True)
            self.ui.gateGroupBox.setEnabled(True)
        
            self.ui.gateGO.setEnabled(False)
            self.ui.bulkProcessGroupBox.setEnabled(True)
            self.ui.batchLoadDataPushButton.setEnabled(True)
        
        self.ui.logTextBrowser.append( "opened header file " + matfile )

    def openMultipleRecords(self):

        try:
            with open('.akvo.last.path') as f: 
                fpath = f.readline()  
                pass
        except IOError as e:
            fpath = '.'
        
        self.headerstrRAW = QtGui.QFileDialog.getOpenFileNames(self, 'Open multiple record', fpath,   r"ORS RAW files (*.ors *.ors2);;Vista Clara Mid Field (*.vcm)")
        # ;;Vista Clara Javelin Processed (*.mat)" )
        self.headerstr = [] #[QtGui.QFileDialog.getOpenFileName(self, 'Open single record', fpath,  r"ORS RAW files (*.ors *.ors2 *.vcm);;VC Mid field (*.vcm)" )]

        # Windoze and \ in paths
        #for iff in range(len(self.headerstr)) :
        #    self.headerstr[iff] = os.path.normpath( str(self.headerstr[iff]) )

        # fix Windows stupid conventions
        for iff in range(len(self.headerstrRAW[0])) :
            #self.headerstr[iff] = os.path.normpath( str(self.headerstr[iff][0]) )
            self.headerstr.append( os.path.normpath( str(self.headerstrRAW[0][iff]) ))
        
        f = open('.akvo.last.path', 'w')
        if str(self.headerstr[0]):
            f.write( str(self.headerstr[0]) ) # prompt last file      

        # Enable all the DSP stuff, but disable individual GO buttons
        self.enableDSP() 
        self.ui.bulkProcessGroupBox.setEnabled(True)
        self.ui.batchLoadDataPushButton.setEnabled(True)
        self.ui.gateGroupBox.setEnabled(True)

        # disable buttons, do flow 
        self.ui.WindowStackGO.setEnabled(False)
        self.ui.envelopeGO.setEnabled(False)
        self.ui.phaseGO.setEnabled(False)
        self.ui.gateGO.setEnabled(False)

        for ff in self.headerstr:
            self.ORS[ff] = MRProc(  )  # Fast
            QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("enableDSP()"), self.enableDSP)
            QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("enableINV()"), self.enableINV)
            QtCore.QObject.connect(self.ORS[ff], QtCore.SIGNAL("doneStatus()"), self.doneStatus)
            #self.ORSConnect( ff )        
        
        self.ui.diffusionGroupBox.setEnabled(True)
        self.ui.diffusionGroupBox.setChecked(True)

        
        path,filen=os.path.split(str(self.headerstr[0] ))
        fileName, fileExtension = os.path.splitext( str(self.headerstr[0]) )

        if fileExtension == ".vcm": 

            self.ui.phaseGroupBox.setChecked(False)
            self.ui.envelopeGroupBox.setChecked(False)
            self.ui.WindowStackBox.setChecked(False)


    def hangon(self): 
        self.plotLOG10ENV()

        # reopen processing flow
        # Enable all the DSP stuff, but disable individual GO buttons
        self.enableDSP() 
        self.ui.bulkProcessGroupBox.setEnabled(True)
        self.ui.batchLoadDataPushButton.setEnabled(True)
        self.ui.gateGroupBox.setEnabled(True)

        # disable buttons, do flow 
        self.ui.WindowStackGO.setEnabled(False)
        self.ui.envelopeGO.setEnabled(False)
        self.ui.phaseGO.setEnabled(False)
        self.ui.gateGO.setEnabled(False)

        # Enable Inversions
        self.enableINV( )

    def batch(self):

        import threading
        # Plot all data
        #self.ui.qplot.getFigure().clf()
        #sp0 = self.ui.qplot.getFigure().add_subplot(211) 
        #sp1 = self.ui.qplot.getFigure().add_subplot(212) 
        
        #########################################
        # LOAD the data and do some basic work
        #self.unlock()
        self.lock("batch processing records")
        
        pars = {  'ftype' : str(self.ui.fTypeComboBox.currentText()),   \
                  'winTrimLeft' : self.ui.winTrimLeftSpinBox.value(),   \
                  'winTrimRight' : self.ui.winTrimRightSpinBox.value(), \
                  'offset' :  self.ui.offsetSpinBox.value(), \
                  'gpd' : self.ui.gpdSpinBox.value(),
                  'stackEfficiency' : self.ui.stackEfficiency.value()     }
        
        total = len(self.headerstr)
        iif = 1
        for ff in self.headerstr:
            
            rawThread = threading.Thread( target=self.ORS[ff].batchThread, args=(pars, ff, os.path.isfile(str(ff)), self.ui.WindowStackBox.isChecked(), \
                            self.ui.envelopeGroupBox.isChecked(),  self.ui.phaseGroupBox.isChecked(), self.ui.gateGroupBox.isChecked()) )
            self.updateProgressBar(1e2*iif/total)
            rawThread.start()
            rawThread.join()
            iif += 1
            #rawThread = thread.start_new_thread(self.ORS[ff].batchThread, (pars, ff, os.path.isfile(ff), self.ui.WindowStackBox.isChecked(), \
            #                self.ui.envelopeGroupBox.isChecked(),  self.ui.phaseGroupBox.isChecked(), self.ui.gateGroupBox.isChecked()) )

            #self.ORS[ff].batchThread(pars, ff, os.path.isfile(ff), self.ui.WindowStackBox.isChecked(), \
            #                self.ui.envelopeGroupBox.isChecked(),  self.ui.phaseGroupBox.isChecked(), self.ui.gateGroupBox.isChecked()) 
        
        self.unlock()
        self.hangon()
        #rawThread = thread.start_new_thread(self.ORS.loadORSFile, \
        #        (str(ff),))
        #rawThread = thread.start_new_thread(self.batchThread, (None,))
        #self.batchThread(ff)


    def openSingleRecord(self):

        try:
            with open('.akvo.last.path') as f: 
                fpath = f.readline()  
                pass
        except IOError as e:
            fpath = '.'
            

        self.headerstrRAW = [QtGui.QFileDialog.getOpenFileName(self, 'Open single record', fpath,  r"ORS RAW files (*.ors *.ors2 *.vcm);;VC Mid field (*.vcm)")]
        #;;VC Javelin (*.mat)" )]
        self.headerstr = [os.path.normpath( str(self.headerstrRAW[0][0]) )]
    
        #print(self.headerstr)
        #exit()            
    
        self.ui.headerFileTextBrowser.clear()
        self.ui.headerFileTextBrowser.append( self.headerstr[0] )
        
        # clear the processing log? TODO do we want to do this? Not sure what about multiple inputs?
        self.ui.logTextBrowser.clear()

        path,filen=os.path.split(str(self.headerstr[0] ))
        fileName, fileExtension = os.path.splitext( str(self.headerstr[0]) )

        f = open('.akvo.last.path', 'w')
        if str(self.headerstr[0]):
            f.write( str(self.headerstr[0]) ) # prompt last file      

        self.ui.headerFileBox.setEnabled(True)

        #########################################
        # LOAD the data and do some basic work
        self.ORSConnect( str(self.headerstr[0]) )       
       

        self.lock("loading RAW record")
        if fileExtension == ".ors": 
            rawThread = thread.start_new_thread(self.ORS[self.headerstr[0]].loadORSFile, \
                (str(self.headerstr[0]),)) 
        elif fileExtension == ".vcm":
            rawThread = thread.start_new_thread(self.ORS[self.headerstr[0]].loadVCMFile, \
                (str(self.headerstr[0]),)) 
            self.ui.gateGroupBox.setEnabled(True)

            self.ui.phaseGroupBox.setChecked(False)
            self.ui.envelopeGroupBox.setChecked(False)
            self.ui.WindowStackBox.setChecked(False)
            
            # disable buttons, do flow 
            self.ui.WindowStackGO.setEnabled(False)
            self.ui.envelopeGO.setEnabled(False)
            self.ui.phaseGO.setEnabled(False)
            self.ui.gateGO.setEnabled(True)
        
        self.ui.logTextBrowser.append( "opened header file " + self.headerstr[0] )
        self.enableDSP( )

    #################################################################################################
    #################  PLOTTING ROUTINES  ###########################################################
    #################################################################################################

    def plotRAW(self, ff=None):

        if ff == None:
            ff = self.headerstr[0]     

        self.ui.statusbar.showMessage ( "plotting RAW" )

        #########################################
        # plot
        self.ui.qplot.getFigure().clf()
        self.ui.qplot.draw()
      
        sp0 = self.ui.qplot.getFigure().add_subplot(211) 
        sp1 = self.ui.qplot.getFigure().add_subplot(212) 
        #DATASTACK = np.average(ORS.DATA, axis=0)
       
        # Can we instead concatinate some of these in order to only call plot 1 time? 
        nt = len(self.ORS[ff].TIMES)
        for iss in range(0, self.ORS[ff].NS):
            times = np.empty(nt*self.ORS[ff].NE)
            dati = np.empty(nt*self.ORS[ff].NE)
            datr = np.empty(nt*self.ORS[ff].NE)
            #COLOUR = np.empty((nt*self.ORS[ff].NE, 3))
            for ie in range(0, self.ORS[ff].NE):
                #COLOUR[nt*ie:nt*(ie+1),:] = np.array([ie/self.ORS[ff].NE, 0., 1.-ie/self.ORS[ff].NE]) # for plots
                times[nt*ie:nt*(ie+1)] = self.ORS[ff].TAUE*(float)(ie) + self.ORS[ff].TIMES
                dati[nt*ie:nt*(ie+1)] =  np.imag(self.ORS[ff].DATA[iss,ie]) 
                datr[nt*ie:nt*(ie+1)] =  np.real(self.ORS[ff].DATA[iss,ie]) 
                #sp0.plot( self.ORS[ff].TAUE*(float)(ie) + self.ORS[ff].TIMES, np.imag(self.ORS[ff].DATA[iss,ie]),  color=COLOUR)
                #sp1.plot( self.ORS[ff].TAUE*(float)(ie) + self.ORS[ff].TIMES, np.real(self.ORS[ff].DATA[iss,ie]),  color=COLOUR)
            sp0.plot(times, dati, '.', color='blue', markersize=1) #COLOUR) doesn't work
            sp1.plot(times, datr, '.', color='blue', markersize=1)
        sp0.set_title("imaginary")
        sp1.set_title("real")
        self.ui.qplot.getFigure().suptitle("RAW quadrature data")
        sp1.set_xlabel("time [s]")
        plt.setp(sp0.get_xticklabels(), visible=False)
        self.ui.qplot.draw()

    def plotWIN(self, ff=None):

        if ff == None:
            ff = self.headerstr[0]        
        
        #########################################
        # plot
        self.ui.qplot.getFigure().clf()
        self.ui.qplot.draw()
      
        sp0 = self.ui.qplot.getFigure().add_subplot(211) 
        sp1 = self.ui.qplot.getFigure().add_subplot(212) 
        
        nt = len(self.ORS[ff].TIMES)
        times = np.empty(nt*self.ORS[ff].NE)
        dati = np.empty(nt*self.ORS[ff].NE)
        datr = np.empty(nt*self.ORS[ff].NE)
        for ie in range(0, self.ORS[ff].NE):
            times[nt*ie:nt*(ie+1)] = self.ORS[ff].TAUE*(float)(ie) + self.ORS[ff].TIMES
            dati[nt*ie:nt*(ie+1)] =  np.imag(self.ORS[ff].DATASTACK[ie]) 
            datr[nt*ie:nt*(ie+1)] =  np.real(self.ORS[ff].DATASTACK[ie]) 
            #COLOUR = np.array([ie/self.ORS[ff].NE, 0., 1.-ie/self.ORS[ff].NE]) # for plots
            #sp0.plot( self.ORS[ff].TAUE*(float)(ie) + self.ORS[ff].TIMES, np.imag(self.ORS[ff].DATASTACK[ie]),  color=COLOUR)
            #sp1.plot( self.ORS[ff].TAUE*(float)(ie) + self.ORS[ff].TIMES, np.real(self.ORS[ff].DATASTACK[ie]),  color=COLOUR)
        sp0.plot(times, dati, '.', color='blue', markersize=1) #COLOUR) doesn't work
        sp1.plot(times, datr, '.', color='blue', markersize=1)
        
        sp0.set_title("imaginary")
        sp1.set_title("real")
        self.ui.qplot.getFigure().suptitle("RAW quadrature data")
        sp1.set_xlabel("time [s]")
        plt.setp(sp0.get_xticklabels(), visible=False)
        self.ui.qplot.draw()

    def plotENV(self, ff=None):
        
        if ff == None:
            ff = self.headerstr[0]     


        #########################################
        # plot
        self.ui.qplot.getFigure().clf()
        self.ui.qplot.draw()
        sp0 = self.ui.qplot.getFigure().add_subplot(211) 
        sp1 = self.ui.qplot.getFigure().add_subplot(212) 
        
        sp0.errorbar(self.ORS[ff].T2T, np.imag(self.ORS[ff].T2D), fmt='.', yerr=self.ORS[ff].sigma, label="data imag")
        sp1.errorbar(self.ORS[ff].T2T, np.real(self.ORS[ff].T2D), fmt='.', yerr=np.std(self.ORS[ff].T2D.real), label="data real")
        if self.ORS[ff].NS > 2:
            sp1.errorbar(self.ORS[ff].T2T, np.imag(self.ORS[ff].T2N), fmt='.', yerr=self.ORS[ff].sigma, label="noise imag")
       
        sp0.axhline(linewidth=1, color='#d62728') 
        sp1.axhline(linewidth=1, color='#d62728') 

        sp0.set_title("imaginary")
        sp1.set_title("real")
        self.ui.qplot.getFigure().suptitle("$T_2$ envelope dude")
        sp1.set_xlabel("time [s]")
        plt.setp(sp0.get_xticklabels(), visible=False)
        self.ui.qplot.draw()
    
    def plotLOG10ENV(self):
        #Valid names are: ['Blues', 'BuGn', 'BuPu', 'GnBu', 'Greens', 'Greys', 'OrRd', 'Oranges', 'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']
        nc = max(3, len(self.headerstr))
        if len(self.headerstr) > 9:
            nc = 9
        colours = brewer2mpl.get_map("Greens", "Sequential", nc).mpl_colors# [::-1]
        
        self.ui.qplot.getFigure().clf()

        if self.burst:
            sp0  = self.ui.qplot.getFigure().add_subplot(221, facecolor=[.4,.4,.4]) 
            sp0b = self.ui.qplot.getFigure().add_subplot(222, facecolor=[.4,.4,.4]) 
            sp1  = self.ui.qplot.getFigure().add_subplot(223, facecolor=[.4,.4,.4], sharex=sp0, sharey=sp0) 
            sp1b = self.ui.qplot.getFigure().add_subplot(224, facecolor=[.4,.4,.4], sharex=sp0b, sharey=sp0b)
            #plt.tight_layout(self.ui.qplot.getFigure()) 
        else:
            sp0 = self.ui.qplot.getFigure().add_subplot(211, facecolor=[.4,.4,.4]) 
            sp1 = self.ui.qplot.getFigure().add_subplot(212, facecolor=[.4,.4,.4], sharex=sp0, sharey=sp0) 
        
        if len(self.headerstr) <= 2:
            iif = 1
        else:
            iif = 0

        # why to -1?
        #for ff in self.headerstr[::-1]:
        for ff in self.headerstr:
            self.plotSingleLOG10ENV( self.ORS[ff], sp0, sp1, colours[iif%9], ff )
            if self.burst:
                self.plotSingleLOG10ENVBurst( self.ORS[ff], sp0b, sp1b, colours[iif%9], ff )
            iif += 1

        NLOG = True #True # True
        if NLOG:
            # write out noise log for GJI article 
            nlog = open("noise.dat", 'w')
            nlog2 = open("noise2.dat", 'w')
            first = False
            for ff in self.headerstr[::-1]:
                if not first:
                    times = self.ORS[ff].T2T
                    ifirst = True
                    for tt in times: 
                        if ifirst:
                            nlog.write( str(tt)  )
                            nlog2.write( str(tt)  )
                            ifirst = False
                        else:
                            nlog.write( r"," + str(tt)  )
                            nlog2.write( r"," + str(tt)  )
                    nlog.write( "\n" )
                    nlog2.write( "\n" )
                    first = True
                dat = self.ORS[ff].T2D.real
                ifirst = True
                for d in dat: 
                    if ifirst:
                        nlog.write( str(d)  )
                        ifirst = False
                    else:
                        nlog.write( r"," + str(d)  )
                nlog.write( "\n" )
                dat = self.ORS[ff].T2D.imag
                ifirst = True
                for d in dat: 
                    if ifirst:
                        nlog2.write( str(d)  )
                        ifirst = False
                    else:
                        nlog2.write( r"," + str(d)  )
                nlog2.write( "\n" )
            nlog.close()
            nlog2.close()

        sp0.set_title("imaginary")
        sp1.set_title("real")
        sp1.set_xlabel("time (s)")
        sp1.set_ylabel("signal (PU)")
        sp0.set_ylabel("signal (PU)")
        plt.setp(sp0.get_xticklabels(), visible=False)
        sp0.set_xscale('log', basex=10)
        sp1.set_xscale('log', basex=10)
        sp0.xaxis.grid(b=True, which='major', color='black', linestyle='-')
        sp1.xaxis.grid(b=True, which='major', color='black', linestyle='-')

        if self.burst:
            sp0b.set_title("burst im")
            sp1b.set_title("$\Re$burst")
            sp1b.set_xlabel("time [s]")
            plt.setp(sp0b.get_xticklabels(), visible=False)
            sp0b.set_xscale('log', basex=10)
            sp1b.set_xscale('log', basex=10)
            sp0b.xaxis.grid(b=True, which='major', color='black', linestyle='-')
            sp1b.xaxis.grid(b=True, which='major', color='black', linestyle='-')

        # making the tickmarks white, longer, and wider 
        #sp0.tick_params(axis='both', which='both', color='w')#,length=3,width=2)
        #sp1.tick_params(axis='both', which='both', color='w')#,length=3,width=2)

        self.ui.qplot.draw()
        self.doneStatus()

    def plotSingleLOG10ENV(self, ORS, sp0, sp1, colour, ff):

        #########################################
        # plot
        #self.ui.qplot.draw()
        #sp0 = self.ui.qplot.getFigure().add_subplot(211) 
        #sp1 = self.ui.qplot.getFigure().add_subplot(212)

        #dataf = open("javelin_" + str(ff) + ".dat", "w")
        #for ii in range(len(ORS.T2D)):
        #    dataf.write( "%f\t%f\t%f\n" %(ORS.T2T[ii],np.imag(ORS.T2D[ii]),ORS.sigma[ii])  )
       
        #print ("sigma", ORS.sigma) 
        sp0.errorbar(ORS.T2T, np.imag(ORS.T2D), fmt='.', color=colour, yerr=ORS.sigma, label="data imag")
        sp1.errorbar(ORS.T2T, np.real(ORS.T2D), fmt='.', color=colour, yerr=ORS.sigma, label="data real")
        if ORS.NS > 2:
            sp1.errorbar(ORS.T2T, np.imag(ORS.T2N), fmt='.', yerr=ORS.sigma, label="noise imag")
    
    def plotSingleLOG10ENVBurst(self, ORS, sp0, sp1, colour, ff):
        sp0.errorbar(ORS.T2Tb, np.imag(ORS.T2Db), fmt='.', color=colour, yerr=ORS.sigmaBurst, label="data imag")
        sp1.errorbar(ORS.T2Tb, np.real(ORS.T2Db), fmt='.', color=colour, yerr=ORS.sigmaBurst, label="data real")
        #if ORS.NS > 2:
        #    sp1.errorbar(ORS.T2Tb, np.imag(ORS.T2Nb), fmt='.', yerr=ORS.sigma, label="noise imag")
        

    def plotMONO(self):

        # TODO for masked data, show differently in fit

        #########################################
        # plot
        if not self.ui.saveDatBox.isChecked():
            self.ui.qplot_2.getFigure().clf()
            self.dsp0 = self.ui.qplot_2.getFigure().add_subplot(211) 
            self.dsp1 = self.ui.qplot_2.getFigure().add_subplot(212) 
        
            self.dsp0.set_title("data")
            self.dsp1.set_title("residuals with 1 $\sigma$ errorbars")
            #self.ui.qplot_2.getFigure().suptitle("multi-exponential fit")
            self.dsp1.set_xlabel("time [s]")
            plt.setp(self.dsp0.get_xticklabels(), visible=False)

            self.dsp0.set_xscale('log', basex=10)
            self.dsp1.set_xscale('log', basex=10)
     
        for ff in self.headerstr:   
            self.dsp0.errorbar(self.ORS[ff].T2T, np.imag(self.ORS[ff].T2D), fmt='o', yerr=self.ORS[ff].sigma, label="data imag")
            self.dsp0.plot(self.ORS[ff].MonoRes['t'], self.ORS[ff].MonoRes['env'], 'o', label=self.ORS[ff].MonoRes['rfit'])
            self.dsp1.errorbar(self.ORS[ff].MonoRes['t'], self.ORS[ff].MonoRes['env'] - self.ORS[ff].T2D.imag , fmt='o', yerr=self.ORS[ff].sigma, label="mono")
        
        #sp0.legend(numpoints=1)
        #sp1.legend(numpoints=1)
        self.ui.qplot_2.draw()
        
        #########################################
        # plot model

        if not self.ui.saveModBox.isChecked():
            self.ui.qplot_3.getFigure().clf()
            self.msp0 = self.ui.qplot_3.getFigure().add_subplot(111) 
            self.msp0.set_xlabel(r"$T_2$ [s]", color='black')
            self.msp0.set_ylabel(r"$A_0$ [rku]", color='black') 

        for ff in self.headerstr:   
            self.msp0.plot( self.ORS[ff].MonoRes['rt20'], self.ORS[ff].MonoRes['b0'], 'o', markersize=6, label="mono")   
            self.msp0.legend(numpoints=1)
            self.ui.qplot_3.draw()  

    def plotBI(self):

        # TODO for masked data, show differently in fit

        #########################################
        # plot
        if not self.ui.saveDatBox.isChecked():
            self.ui.qplot_2.getFigure().clf()
            self.dsp0 = self.ui.qplot_2.getFigure().add_subplot(211) 
            self.dsp1 = self.ui.qplot_2.getFigure().add_subplot(212) 
        
            self.dsp0.set_title("data")
            self.dsp1.set_title("residuals with 1 $\sigma$ errorbars")
            #self.ui.qplot_2.getFigure().suptitle("multi-exponential fit")
            self.dsp1.set_xlabel("time [s]")
            plt.setp(self.dsp0.get_xticklabels(), visible=False)

            self.dsp0.set_xscale('log', basex=10)
            self.dsp1.set_xscale('log', basex=10)
     
        for ff in self.headerstr:   
            self.dsp0.errorbar(self.ORS[ff].T2T, np.imag(self.ORS[ff].T2D), fmt='o', yerr=self.ORS[ff].sigma, label="data imag")
            self.dsp0.plot(self.ORS[ff].BiRes['t'], self.ORS[ff].BiRes['env'], 'o', label=self.ORS[ff].BiRes['rfit'])
            self.dsp1.errorbar(self.ORS[ff].BiRes['t'], self.ORS[ff].BiRes['env'] - self.ORS[ff].T2D.imag , fmt='o', yerr=self.ORS[ff].sigma, label="mono")
        
        #sp0.legend(numpoints=1)
        #sp1.legend(numpoints=1)
        self.ui.qplot_2.draw()
        
        #########################################
        # plot model

        if not self.ui.saveModBox.isChecked():
            self.ui.qplot_3.getFigure().clf()
            self.msp0 = self.ui.qplot_3.getFigure().add_subplot(111) 
            self.msp0.set_xlabel(r"$T_2$ [s]", color='black')
            self.msp0.set_ylabel(r"$A_0$ [rku]", color='black') 

        for ff in self.headerstr:   
            self.msp0.plot( self.ORS[ff].BiRes['rt20'], self.ORS[ff].BiRes['b0'], 'o', markersize=6, label="bi1")   
            self.msp0.plot( self.ORS[ff].BiRes['rrt20'], self.ORS[ff].BiRes['bb0'], 'o', markersize=6, label="bi2")   
            self.msp0.legend(numpoints=1)
            self.ui.qplot_3.draw()  

    def plotDIFF(self, Times2, Data, Sigma, Fit=None):
        
        #########################################
        # Plot up Data and Errorbars
        if not self.ui.saveDatBox.isChecked():
            self.ui.qplot_6.getFigure().clf()
            self.dsp0 = self.ui.qplot_6.getFigure().add_subplot(211) 
            self.dsp1 = self.ui.qplot_6.getFigure().add_subplot(212) 
        
            self.dsp0.set_title("data")
            self.dsp1.set_title("residuals with 1 $\sigma$ errorbars")
            #self.ui.qplot_2.getFigure().suptitle("multi-exponential fit")
            self.dsp1.set_xlabel("time [s]")
            plt.setp(self.dsp0.get_xticklabels(), visible=False)

            self.dsp0.set_xscale('log', basex=10)
            self.dsp1.set_xscale('log', basex=10)
     
        #for ff in self.headerstr:   
        self.dsp0.errorbar(Times2, Data, fmt='o', yerr=Sigma, label="data imag")
        if Fit != None:
            self.dsp0.plot(Times2, Fit, 'o', label="fit")
            #self.dsp0.plot(self.ORS[ff].T2T, self.ORS[ff].MonoRes['env'], 'o', label=self.ORS[ff].MonoRes['rfit'])
            self.dsp1.errorbar(Times2, Data - Fit , fmt='o', yerr=Sigma, label="mono")
        
        #sp0.legend(numpoints=1)
        #sp1.legend(numpoints=1)
        self.ui.qplot_6.draw() 

    def plotDIFF2(self, Times2, Data, Datap, Sigma, Fit=None):
        
        #########################################
        # Plot up Data and Errorbars
        if not self.ui.saveDatBox.isChecked():
            self.ui.qplot_6.getFigure().clf()
            self.dsp0 = self.ui.qplot_6.getFigure().add_subplot(211) 
            self.dsp1 = self.ui.qplot_6.getFigure().add_subplot(212) 
        
            self.dsp0.set_title("data")
            self.dsp1.set_title("residuals with 1 $\sigma$ errorbars")
            #self.ui.qplot_2.getFigure().suptitle("multi-exponential fit")
            self.dsp1.set_xlabel("time [s]")
            plt.setp(self.dsp0.get_xticklabels(), visible=False)

            self.dsp0.set_xscale('log', basex=10)
            self.dsp1.set_xscale('log', basex=10)
     
        #for ff in self.headerstr:   
        self.dsp0.errorbar(Times2, Data,  fmt='o', yerr=Sigma, label="data imag")
        self.dsp0.errorbar(Times2, Datap, fmt='o', yerr=Sigma, label="data imag")
        if Fit != None:
            self.dsp0.plot(Times2, Fit, 'o', label="fit")
            #self.dsp0.plot(self.ORS[ff].T2T, self.ORS[ff].MonoRes['env'], 'o', label=self.ORS[ff].MonoRes['rfit'])
            self.dsp1.errorbar(Times2, Datap - Fit , fmt='o', yerr=Sigma, label="mono")
        
        #sp0.legend(numpoints=1)
        #sp1.legend(numpoints=1)
        self.ui.qplot_6.draw() 

    def plotDIST(self, mask):

        #########################################
        # plot data
        if not self.ui.saveDatBox.isChecked():
            self.ui.qplot_2.getFigure().clf()
            self.dsp0 = self.ui.qplot_2.getFigure().add_subplot(211) 
            self.dsp1 = self.ui.qplot_2.getFigure().add_subplot(212) 
        
            self.dsp0.set_title("data")
            self.dsp1.set_title("residuals with 1 $\sigma$ errorbars")
            #self.ui.qplot_2.getFigure().suptitle("multi-exponential fit")
            self.dsp1.set_xlabel("time [s]")
            plt.setp(self.dsp0.get_xticklabels(), visible=False)

            self.dsp0.set_xscale('log', basex=10)
            self.dsp1.set_xscale('log', basex=10)
       
        for ff in self.headerstr: 
            self.dsp0.errorbar(self.ORS[ff].T2T, np.imag(self.ORS[ff].T2D), fmt='o', color='blue', yerr=self.ORS[ff].sigma, label="data")
            self.dsp0.plot(self.ORS[ff].T2T[mask:], self.ORS[ff].DistRes['env'], '.', color='red', label='dist')#, label=self.ORS[ff].MonoRes['rfit'])
            self.dsp1.errorbar(self.ORS[ff].T2T[mask:], self.ORS[ff].DistRes['env'] - self.ORS[ff].T2D.imag[mask:] , fmt='o', yerr=self.ORS[ff].sigma[mask:], label="dist")
        
        #self.dsp0.legend(numpoints=1)
        #self.dsp1.legend(numpoints=1)
        
        #self.ui.qplot_2.draw()  
        
        nc = max(3, len(self.headerstr))
        if len(self.headerstr) > 9:
            nc = 9
        colours = brewer2mpl.get_map("Greens", "Sequential", nc).mpl_colors# [::-1]
 
        #########################################
        # plot model

        if not self.ui.saveModBox.isChecked():
            self.ui.qplot_3.getFigure().clf()
            #self.msp0 = self.ui.qplot_3.getFigure().add_subplot(111, facecolor=[.4,.4,.4]) 
            self.msp0 = self.ui.qplot_3.getFigure().add_axes([.155,.2,.7,.65], facecolor=[.4,.4,.4]) 
            self.msp0.set_xlabel(r"$T_2$ [s]", color='black')
            self.msp0.set_ylabel(r"PWC (m$^3$/m$^3$)", color='black') 
            self.msp1 = self.msp0.twinx()     
            self.msp1.set_ylabel(r"total water (%)", color='black') 
        
        self.msp0.set_xscale('log', basex=10)
        self.msp0.xaxis.grid(b=True, which='major', color='black', linestyle='-')
        
        if len(self.headerstr) <= 2:
            iif = 1
        else:
            iif = 0
            
        print("%s,%s,%s"%("theta","LogMeanT2","filename")) 
        for ff in self.headerstr[::-1]: 
            #############################
            # Plot dist. model
            # consider normalizing by bin width? If you shift nBins it appears that water content varies
            self.msp0.plot(self.ORS[ff].DistRes['Time'].T2Bins, self.ORS[ff].DistRes['mod'], '-', color=colours[iif%nc], linewidth=2, label="recovered")
                    
            ################################
            # Log mean and total water
            theta = np.sum( self.ORS[ff].DistRes['mod'] )
            LogMeanT2 = np.exp(np.sum( self.ORS[ff].DistRes['mod'] * np.log( self.ORS[ff].DistRes['Time'].T2Bins) ) / theta )  

#            print("%1.5e,%1.5e,%s"%(theta,LogMeanT2, os.path.basename(ff))) 

            self.msp1.plot(LogMeanT2, 1e2*theta, 'o', markersize=6, color=colours[iif%nc], label="$T_{2ML}$")

            iif += 1
            

        #self.msp0.xaxis.grid(b=True, which='minor', color=[.2,.2,.2], linestyle='-')

        #self.dsp1.set_xscale('log', basex=10)
        #self.msp0.legend(numpoints=1)
        self.ui.qplot_3.draw()  

if __name__ == "__main__":
    
    import time
    import matplotlib.image as mpimg 

    app = QtGui.QApplication(sys.argv)
    
    # splash a little advertising
    #pixmap = QtGui.QPixmap("XRI_LOGO2.jpg")
    #splash = QtGui.QSplashScreen(pixmap, QtCore.Qt.WindowStaysOnTopHint)
    #splash.show()
    #splash.showMessage("Loaded modules"); 
    
    # Build main app
    myapp = MyForm()

    # dummy plot
    img=mpimg.imread('./lemma.png')
    for ax in [ myapp.ui.qplot, myapp.ui.qplot_2, myapp.ui.qplot_3 ]:
        subplot = ax.getFigure().add_subplot(111)
        ax.getFigure().patch.set_facecolor( None )
        ax.getFigure().patch.set_alpha( .0 )
        subplot.imshow(img) 
        subplot.xaxis.set_major_locator(plt.NullLocator()) 
        subplot.yaxis.set_major_locator(plt.NullLocator()) 
        ax.draw()
        #subplot = myapp.ui.qplot.getFigure().add_subplot(111)
        #myapp.ui.qplot.getFigure().patch.set_facecolor( None )
        #myapp.ui.qplot.getFigure().patch.set_alpha( .0 )
        #subplot.imshow(img) 
        #subplot.xaxis.set_major_locator(plt.NullLocator()) 
        #subplot.yaxis.set_major_locator(plt.NullLocator()) 
        #myapp.ui.qplot.draw()
    
    myapp.show()
    #splash.finish(myapp);

    sys.exit(app.exec_())