Addition of DOI metric based on incomplete resolution matrix construction

akvo/gui/akvoGUI.py

         invDict["T2Bins"]["low"] = T2lo
         invDict["T2Bins"]["high"] = T2hi
         invDict["T2Bins"]["number"] = NT2       
+        invDict["NonLinearRefinement"] = self.ui.NLButton.isChecked()
+        invDict["CalcDOI"] = self.ui.DOIButton.isChecked()
  
         node = yaml.YAML()
         kpo = open( "invert.yml", 'w' )

akvo/gui/main.ui

               <rect>
                <x>0</x>
                <y>0</y>
-               <width>96</width>
-               <height>26</height>
+               <width>100</width>
+               <height>30</height>
               </rect>
              </property>
              <attribute name="label">
               <rect>
                <x>0</x>
                <y>0</y>
-               <width>96</width>
-               <height>26</height>
+               <width>411</width>
+               <height>67</height>
               </rect>
              </property>
              <attribute name="label">
                <string>QT Inversion</string>
               </property>
               <layout class="QGridLayout" name="gridLayout_19">
-               <item row="2" column="1">
-                <widget class="QComboBox" name="invChan">
-                 <item>
-                  <property name="text">
-                   <string>Chan. 1</string>
-                  </property>
-                 </item>
-                 <item>
-                  <property name="text">
-                   <string>Chan. 2</string>
-                  </property>
-                 </item>
-                 <item>
-                  <property name="text">
-                   <string>Chan. 3</string>
-                  </property>
-                 </item>
-                 <item>
-                  <property name="text">
-                   <string>Chan. 4</string>
-                  </property>
-                 </item>
+               <item row="5" column="0">
+                <widget class="QLabel" name="label_81">
+                 <property name="text">
+                  <string>initial alpha </string>
+                 </property>
                 </widget>
                </item>
                <item row="4" column="1">
                  </item>
                 </widget>
                </item>
-               <item row="3" column="0">
-                <widget class="QLabel" name="label_80">
-                 <property name="text">
-                  <string>T2* model objective</string>
-                 </property>
-                </widget>
-               </item>
-               <item row="4" column="0">
-                <widget class="QLabel" name="label_82">
-                 <property name="text">
-                  <string>Depth model objective</string>
-                 </property>
-                </widget>
-               </item>
                <item row="0" column="1">
                 <widget class="QTextBrowser" name="dataText">
                  <property name="maximumSize">
                  </property>
                 </widget>
                </item>
-               <item row="6" column="0" colspan="2">
+               <item row="1" column="0">
+                <widget class="QPushButton" name="invKernelButton">
+                 <property name="text">
+                  <string>Select Kernel</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="3" column="0">
+                <widget class="QLabel" name="label_80">
+                 <property name="text">
+                  <string>T2* model objective</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="4" column="0">
+                <widget class="QLabel" name="label_82">
+                 <property name="text">
+                  <string>Depth model objective</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="11" column="0" colspan="2">
                 <widget class="QPushButton" name="invertButton">
                  <property name="styleSheet">
                   <string notr="true">#invertButton {
                  </property>
                 </widget>
                </item>
-               <item row="2" column="0">
-                <widget class="QLabel" name="label_79">
-                 <property name="text">
-                  <string>Data channel</string>
-                 </property>
+               <item row="2" column="1">
+                <widget class="QComboBox" name="invChan">
+                 <item>
+                  <property name="text">
+                   <string>Chan. 1</string>
+                  </property>
+                 </item>
+                 <item>
+                  <property name="text">
+                   <string>Chan. 2</string>
+                  </property>
+                 </item>
+                 <item>
+                  <property name="text">
+                   <string>Chan. 3</string>
+                  </property>
+                 </item>
+                 <item>
+                  <property name="text">
+                   <string>Chan. 4</string>
+                  </property>
+                 </item>
                 </widget>
                </item>
-               <item row="0" column="0">
-                <widget class="QPushButton" name="invDataButton">
+               <item row="6" column="0">
+                <widget class="QRadioButton" name="NLButton">
+                 <property name="toolTip">
+                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;The dephased (aka quadrature detected, aka corrected amplitude) inversion introduces non-linearity to the otherwise linear forward operator in that out of phase signals cannot destructively interefere with a real forward operator. Normally this impact is pretty minimal and can be ignored. However, in conductive media this is not always the case. Akvo can perform a non-linear inversion after the linear one to correct for this. The result of the linear inversion is used as an initial model in the non-linear inversion and the same level of regularisation is used. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+                 </property>
                  <property name="text">
-                  <string>Select Data</string>
+                  <string>Non Linear refinement</string>
+                 </property>
+                 <property name="checked">
+                  <bool>true</bool>
+                 </property>
+                 <property name="autoExclusive">
+                  <bool>false</bool>
                  </property>
                 </widget>
                </item>
-               <item row="5" column="0">
-                <widget class="QLabel" name="label_81">
+               <item row="2" column="0">
+                <widget class="QLabel" name="label_79">
                  <property name="text">
-                  <string>initial alpha </string>
+                  <string>Data channel</string>
                  </property>
                 </widget>
                </item>
                  </property>
                 </widget>
                </item>
-               <item row="1" column="0">
-                <widget class="QPushButton" name="invKernelButton">
-                 <property name="text">
-                  <string>Select Kernel</string>
-                 </property>
-                </widget>
-               </item>
                <item row="1" column="1">
                 <widget class="QTextEdit" name="kernelText">
                  <property name="maximumSize">
                  </property>
                 </widget>
                </item>
+               <item row="0" column="0">
+                <widget class="QPushButton" name="invDataButton">
+                 <property name="text">
+                  <string>Select Data</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="6" column="1">
+                <widget class="QRadioButton" name="DOIButton">
+                 <property name="toolTip">
+                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Geophysical inverse problems do not have a true, unambigous depth of investigation. &lt;/p&gt;&lt;p&gt;&lt;br/&gt;&lt;/p&gt;&lt;p&gt;Akvo makes a calculation of the DOI based on a resolution analysis approach. A pseudo 1D total water point spread function is calculated by inverting kronecker delta models at the midpoint of the &lt;span style=&quot; font-style:italic;&quot;&gt;T&lt;/span&gt;&lt;span style=&quot; font-style:italic; vertical-align:sub;&quot;&gt;2&lt;/span&gt;&lt;span style=&quot; font-style:italic; vertical-align:super;&quot;&gt;*&lt;/span&gt; distribution at he same regularisation level as the inversion. The recovered models are summed across the &lt;span style=&quot; font-style:italic;&quot;&gt;T&lt;/span&gt;&lt;span style=&quot; font-style:italic; vertical-align:sub;&quot;&gt;2&lt;/span&gt;&lt;span style=&quot; font-style:italic; vertical-align:super;&quot;&gt;*&lt;/span&gt; space in order to reduce variables in the analysis. hese models are then used as rows in the construction of the point spread function (Model Resolution Matrix). The point spread function is used to identify the point at which the peak depth in the recovered models deviates from the input model by at least 10 %. This is an arbitrary threshold, but is the same one used in[1] &lt;/p&gt;&lt;p&gt;&lt;br/&gt;&lt;/p&gt;&lt;p&gt;[1] Müller-Petke, M., &amp;amp; Yaramanci, U. (2008). Resolution studies for magnetic resonance sounding (MRS) using the singular value decomposition. &lt;span style=&quot; font-style:italic;&quot;&gt;Journal of Applied Geophysics&lt;/span&gt;, &lt;span style=&quot; font-style:italic;&quot;&gt;66&lt;/span&gt;(3-4), 165-175.   &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+                 </property>
+                 <property name="text">
+                  <string>Calculate DOI</string>
+                 </property>
+                 <property name="checked">
+                  <bool>true</bool>
+                 </property>
+                 <property name="autoExclusive">
+                  <bool>false</bool>
+                 </property>
+                </widget>
+               </item>
               </layout>
              </widget>
             </item>

akvo/tressel/invertTA.py

 from akvo.tressel.logbarrier import *
 import yaml,os
 
+import multiprocessing 
+import itertools
+ 
+from scipy.linalg import svd 
+
+from matplotlib.backends.backend_pdf import PdfPages
 from matplotlib.colors import LogNorm
 from matplotlib.colors import LightSource
 from matplotlib.ticker import ScalarFormatter
         K0[0] = np.concatenate( (K0[0].T, K0[ik].T) ).T
     K0 = K0[0]
 
+    #np.save("ifaces", ifaces)
+    #exit()
+
     #plt.matshow(np.real(K0))
     #plt.show()
     #exit()
 
-    ###################    
-    # VERY Simple DOI #
+    ##########################################################    
+    # VERY Simple Sensitivity based calc. of noise per layer #
     maxq = np.argmax(np.abs(K0), axis=1)
     maxK = .1 *  np.abs(K0)[ np.arange(0,len(ifaces)-1), maxq ] # 10% water is arbitrary  
     SNR = maxK / (VS[0][0])
 
-    #SNR[SNR>1] = 1
-    SNRidx = len(ifaces)-3 
-    while SNR[SNRidx] > SNR[SNRidx+1] and SNRidx > 2:
-        SNRidx -= 1
-
-    #print("IDX", SNRidx)
-    #plt.plot(ifaces[0:-1], SNR)
-    #plt.plot(ifaces[0:-1][SNRidx], SNR[SNRidx], '.',markersize=12)
-    #plt.gca().axhline(y=VS[0][0], xmin=0, xmax=ifaces[-1], color='r')
-    #plt.gca().axhline(y=1, xmin=0, xmax=ifaces[-1], color='r')
-    #K0T = np.dot(K0, K0.T) 
-    #K0T = np.dot(K0, np.dot( VS[0][0]* np.eye(len(ifaces)-1,1), K0.T) )
-    #K0T = np.dot(K0, np.dot( 1/(VS[0][0]**2) * np.eye(np.shape(K0)[1]), K0.T) )
-    #plt.matshow(0.05 * np.sqrt(K0T))
-    #plt.colorbar()
-    #plt.plot(ifaces[0:-1], np.diag( 0.01* np.sqrt(K0T))) 
-    #print(np.shape(K0T), len(ifaces)-1)
-    #plt.show()
-    #exit()
-
     ###############################################
     # Build full kernel
     ############################################### 
-    T2Bins = np.logspace( np.log10(cont["T2Bins"]["low"]), np.log10(cont["T2Bins"]["high"]), cont["T2Bins"]["number"], endpoint=True, base=10)  
-    KQT = np.real(buildKQT(np.abs(K0),tg,T2Bins))
 
-    # model resolution matrix 
-    np.linalg.svd(KQT)
+    T2Bins = np.logspace( np.log10(cont["T2Bins"]["low"]), np.log10(cont["T2Bins"]["high"]), cont["T2Bins"]["number"], endpoint=True, base=10) 
+    T2Bins2 = np.append( T2Bins, T2Bins[-1] + (T2Bins[-1]-T2Bins[-2]) )
+    NT2 = len(T2Bins)
 
-    exit()
+    KQT = np.real(buildKQT(np.abs(K0),tg,T2Bins))
 
     ###############################################
     # Linear Inversion 
     print("Calling inversion", flush=True)
     inv, ibreak, errn, phim, phid, mkappa, Wd, Wm, alphastar = logBarrier(KQT, np.ravel(V), T2Bins, "lcurve", MAXITER=150, sigma=np.ravel(VS), alpha=1e6, smooth="Smallest" ) 
 
-    ###############################################
-    # Non-linear refinement! 
-    ###############################################   
- 
-    KQTc = buildKQT(K0, tg, T2Bins)
-    prec = np.abs(np.dot(KQTc, inv))
-    phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
-    #PREc = np.reshape( prec, np.shape(V)  )
-    print("PHID linear=", errn, "PHID complex=", phidc/len(np.ravel(V)))
-    
-    res = nl.nonlinearinversion(inv, Wd, KQTc, np.ravel(V), Wm, alphastar )   
-    if res.success == True:    
-        INVc = np.reshape(res.x, (len(ifaces)-1,cont["T2Bins"]["number"]) )
-        prec = np.abs(np.dot(KQTc, res.x))
-        phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
-        #PREc = np.reshape( prec, np.shape(V)  )
-        print("PHID linear=", errn, "PHID nonlinear=", phidc/len(np.ravel(V)))
-   
-    # Perform second fit around results of first  
-    res = nl.nonlinearinversion(res.x, Wd, KQTc, np.ravel(V), Wm, alphastar )   
-    if res.success == True:    
-        INVc = np.reshape(res.x, (len(ifaces)-1,cont["T2Bins"]["number"]) )
-        prec = np.abs(np.dot(KQTc, res.x))
-        phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
-        #PREc = np.reshape( prec, np.shape(V)  )
-        print("PHID linear=", errn, "PHID nonlinear=", phidc/len(np.ravel(V)))
-
-    #plt.matshow(INVc)
-    #KQTc = buildKQT(K0,tg,T2Bins)
-
-    #plt.matshow(PREc, cmap='Blues')
-    #plt.gca().set_title("complex predicted")
-    #plt.colorbar()
-
-    ###############################################
-    # Appraise
-    ###############################################
-
 
-
-
- 
+    ####################
+    # Summary plots    #
+    ####################
+    
     pre = np.dot(KQT,inv) 
     PRE = np.reshape( pre, np.shape(V)  )
     plt.matshow(PRE, cmap='Blues')
-    plt.gca().set_title("predicted")
+    plt.gca().set_title("linear predicted")
     plt.colorbar()
 
     DIFF = (PRE-V) / VS
     md = np.max(np.abs(DIFF))
     plt.matshow(DIFF, cmap=cmocean.cm.balance, vmin=-md, vmax=md)
-    plt.gca().set_title("misfit / $\widehat{\sigma}$")
+    plt.gca().set_title("linear misfit / $\widehat{\sigma}$")
     plt.colorbar()
     
     plt.matshow(V, cmap='Blues')
     plt.gca().set_title("observed")
     plt.colorbar()
 
+    ###############################################
+    # Non-linear refinement! 
+    ###############################################   
+
+    nonLinearRefinement = cont['NonLinearRefinement']
+    if nonLinearRefinement: 
+
+        KQTc = buildKQT(K0, tg, T2Bins)
+        prec = np.abs(np.dot(KQTc, inv))
+        phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
+        print("PHID forward linear=", errn, "PHID forward nonlinear=", phidc/len(np.ravel(V)))
+    
+        res = nl.nonlinearinversion(inv, Wd, KQTc, np.ravel(V), Wm, alphastar )   
+        if res.success == True:    
+            INVc = np.reshape(res.x, (len(ifaces)-1,cont["T2Bins"]["number"]) )
+            prec = np.abs(np.dot(KQTc, res.x))
+            phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
+            PREc = np.reshape( prec, np.shape(V)  )
+            print("PHID linear=", errn, "PHID nonlinear=", phidc/len(np.ravel(V)))
+
+        while phidc/len(np.ravel(V)) > errn:  
+            phidc_old = phidc/len(np.ravel(V))
+            #alphastar *= .9 
+            res = nl.nonlinearinversion(res.x, Wd, KQTc, np.ravel(V), Wm, alphastar )   
+            if res.success == True:    
+                INVc = np.reshape(res.x, (len(ifaces)-1,cont["T2Bins"]["number"]) )
+                prec = np.abs(np.dot(KQTc, res.x))
+                phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
+                PREc = np.reshape( prec, np.shape(V)  )
+                print("PHID linear=", errn, "PHID nonlinear=", phidc/len(np.ravel(V)))
+            else:
+                break
+
+            if phidc_old - phidc/len(np.ravel(V)) < 0.005:
+                print("Not making progress reducing misfit in nonlinear refinement")
+                break
+    
+        plt.matshow(PREc, cmap='Blues')
+        plt.gca().set_title("nonlinear predicted")
+        plt.colorbar()
+
+        DIFFc = (PREc-V) / VS
+        md = np.max(np.abs(DIFF))
+        plt.matshow(DIFFc, cmap=cmocean.cm.balance, vmin=-md, vmax=md)
+        plt.gca().set_title("nonlinear misfit / $\widehat{\sigma}$")
+        plt.colorbar()
+ 
+    ###############################################
+    # Appraise DOI using simplified MRM 
+    ###############################################
+
+    CalcDOI = cont['CalcDOI']
+
+    if CalcDOI:
+    
+        pdf = PdfPages('resolution_analysis' + '.pdf' )
+        MRM = np.zeros((len(ifaces)-1, len(ifaces)-1))
+
+    #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]) )) 
+ 
+        # This could be parallelized 
+        for ilay in range(len(ifaces)-1):
+            iDeltaT2 = len(T2Bins)//2
+            deltaMod = np.zeros( (len(ifaces)-1, len(T2Bins)) )
+            deltaMod[ilay][iDeltaT2] = 0.3
+            dV = np.dot(KQT, np.ravel(deltaMod)) 
+            # invert 
+            dinv, dibreak, derrn = logBarrier(KQT, dV, T2Bins, "single", MAXITER=1, sigma=np.ravel(VS), alpha=alphastar, smooth="Smallest" ) 
+            print("Sum dinv from", str(ifaces[ilay]), "to", str(ifaces[ilay+1]), "=", np.sum(dinv))
+    
+            DINV = np.reshape(dinv, (len(ifaces)-1,cont["T2Bins"]["number"]) )
+            MRM[ilay,:] = np.sum(DINV, axis=1)
+
+            Y,X = meshgrid( ifaces, T2Bins2 )
+            fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
+            ax1 = fig.add_axes( [.2,.15,.6,.7] )
+            im = ax1.pcolor(X, Y, DINV.T, cmap=cmocean.cm.tempo, shading='auto')
+            ax1.plot( T2Bins[iDeltaT2], (ifaces[ilay]+ifaces[ilay+1])/2, 's', markersize=6, markeredgecolor='black') #, markerfacecolor=None )  
+            im.set_edgecolor('face')
+            ax1.set_xlabel(u"$T_2^*$ (ms)")
+            ax1.set_ylabel(u"depth (m)")
+            ax1.set_xlim( T2Bins2[0], T2Bins2[-1] )
+            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.set_ylim( ifaces[-1], ifaces[0] )
+            ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )   
+            ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
+
+            pdf.savefig(facecolor=[0,0,0,0])
+            plt.close(fig)
+
+        np.save("MRM", MRM)
+        centres = (ifaces[0:-1]+ifaces[1:])/2
+        X,Y = np.meshgrid(ifaces,ifaces)
+
+        fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
+        ax1 = fig.add_axes( [.2,.15,.6,.7] )
+        ax1.pcolor(X,Y,MRM, cmap = cmocean.cm.ice)
+        ax1.set_ylim(ifaces[-1], ifaces[0])
+        maxDepth = np.argmax(MRM, axis=0)
+
+        plt.plot(centres[maxDepth], centres, color='white')
+
+        # Determine DOI 
+        DOIMetric =  centres[maxDepth]/centres #> 0.9 
+        DOI = ifaces[ np.where(DOIMetric < 0.9 ) ][0]
+        plt.axhline(y=DOI, color='white', linestyle='-.')
+
+        ax1.set_ylim( ifaces[-1], ifaces[0] )
+        ax1.set_xlim( ifaces[0], ifaces[-1] )
+        ax1.set_xlabel(u"depth (m)")
+        ax1.set_ylabel(u"depth (m)")
+
+        pdf.close()
+
+
 
-    T2Bins = np.append( T2Bins, T2Bins[-1] + (T2Bins[-1]-T2Bins[-2]) )
     INV = np.reshape(inv, (len(ifaces)-1,cont["T2Bins"]["number"]) )
 
     #alphas = np.tile(SNR, (len(T2Bins)-1,1))
 
     ##############  LINEAR RESULT   ##########################
 
-    Y,X = meshgrid( ifaces, T2Bins )
+    Y,X = meshgrid( ifaces, T2Bins2 )
     fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
     ax1 = fig.add_axes( [.2,.15,.6,.7] )
     im = ax1.pcolor(X, Y, INV.T, cmap=cmocean.cm.tempo) #cmap='viridis')
     #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
     #ax1.axhline( y=ifaces[SNRidx], xmin=T2Bins[0], xmax=T2Bins[-1], color='black'  )
     im.set_edgecolor('face')
-    ax1.set_xlim( T2Bins[0], T2Bins[-1] )
+    ax1.set_xlim( T2Bins[0], T2Bins2[-1] )
     ax1.set_ylim( ifaces[-1], ifaces[0] )
     cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
     cb.locator = MaxNLocator( nbins = 4)
     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'  )
+    #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.xaxis.set_label_position('bottom') 
 
     plt.savefig("akvoInversion.pdf")
 
     ##############  NONLINEAR RESULT   ##########################
 
-    Y,X = meshgrid( ifaces, T2Bins )
-    fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
-    ax1 = fig.add_axes( [.2,.15,.6,.7] )
-    im = ax1.pcolor(X, Y, INVc.T, cmap=cmocean.cm.tempo) #cmap='viridis')
-    #im = ax1.pcolor(X[0:SNRidx,:], Y[0:SNRidx,:], INV.T[0:SNRidx,:], cmap=cmocean.cm.tempo) #cmap='viridis')
-    #im = ax1.pcolor(X[SNRidx::,:], Y[SNRidx::,:], INV.T[SNRidx::,:], cmap=cmocean.cm.tempo, alpha=.5) #cmap='viridis')
-    #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
-    #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
-    #ax1.axhline( y=ifaces[SNRidx], xmin=T2Bins[0], xmax=T2Bins[-1], color='black'  )
-    im.set_edgecolor('face')
-    ax1.set_xlim( T2Bins[0], T2Bins[-1] )
-    ax1.set_ylim( ifaces[-1], ifaces[0] )
-    cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
-    cb.locator = MaxNLocator( nbins = 4)
-    cb.ax.yaxis.set_offset_position('left')                         
-    cb.update_ticks()
+    if nonLinearRefinement: 
+        Y,X = meshgrid( ifaces, T2Bins )
+        fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
+        ax1 = fig.add_axes( [.2,.15,.6,.7] )
+        im = ax1.pcolor(X, Y, INVc.T, cmap=cmocean.cm.tempo) #cmap='viridis')
+        #im = ax1.pcolor(X[0:SNRidx,:], Y[0:SNRidx,:], INV.T[0:SNRidx,:], cmap=cmocean.cm.tempo) #cmap='viridis')
+        #im = ax1.pcolor(X[SNRidx::,:], Y[SNRidx::,:], INV.T[SNRidx::,:], cmap=cmocean.cm.tempo, alpha=.5) #cmap='viridis')
+        #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
+        #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
+        #ax1.axhline( y=ifaces[SNRidx], xmin=T2Bins[0], xmax=T2Bins[-1], color='black'  )
+        im.set_edgecolor('face')
+        ax1.set_xlim( T2Bins[0], T2Bins[-1] )
+        ax1.set_ylim( ifaces[-1], ifaces[0] )
+        cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
+        cb.locator = MaxNLocator( nbins = 4)
+        cb.ax.yaxis.set_offset_position('left')                         
+        cb.update_ticks()
  
-    ax1.set_xlabel(u"$T_2^*$ (ms)")
-    ax1.set_ylabel(u"depth (m)")
+        ax1.set_xlabel(u"$T_2^*$ (ms)")
+        ax1.set_ylabel(u"depth (m)")
     
-    ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
-    ax1.get_yaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
-    ax1.xaxis.set_major_locator( MaxNLocator(nbins = 4) )   
-
-    #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.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'  )
-    #ax2.xaxis.set_label_position('bottom') 
-    fig.suptitle("Non linear inversion")
-    plt.savefig("akvoInversionNL.pdf")
-
-
+        ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+        ax1.get_yaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+        ax1.xaxis.set_major_locator( MaxNLocator(nbins = 4) )   
+
+        #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.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.xaxis.set_label_position('bottom') 
+        fig.suptitle("Non linear inversion")
+        plt.savefig("akvoInversionNL.pdf")
 
     #############
     # water plot#
     ax = fig2.add_axes( [.2,.15,.6,.7] )
     
     # Bound water cutoff 
-    Bidx = T2Bins[0:-1]<33.0
+    Bidx = T2Bins<33.0
     twater = np.sum(INV, axis=1)
     bwater = np.sum(INV[:,Bidx], axis=1)
     
     ax.set_xlim( 0, ax.get_xlim()[1] )
     
     #ax.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed'  )
+    if CalcDOI:
+        ax.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed'  )
     
     plt.savefig("akvoInversionWC.pdf")
     plt.legend()  
  
 
     # Report results into a text file 
-    fr = pd.DataFrame( INV, columns=T2Bins[0:-1] )
+    fr = pd.DataFrame( INV, columns=T2Bins ) #[0:-1] )
     fr.insert(0, "layer top", ifaces[0:-1] )
     fr.insert(1, "layer bottom", ifaces[1::] )
     fr.insert(2, "NMR total water", np.sum(INV, axis=1) )
     fr.insert(3, "NMR bound water", bwater )
     fr.insert(4, "Layer SNR", SNR )
+    if CalcDOI:
+        fr.insert(5, "Resolution", DOIMetric )
 
     fr.to_csv("akvoInversion.csv")    
     #fr.to_excel("akvoInversion.xlsx")    

akvo/tressel/logbarrier.py

 
     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) ]
+        #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[-1]
         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
         #print("Returning L curve result")
         return x, ibreak, np.sqrt(phid/len(b)), PHIM, PHID/len(b), np.argmax(kappa), Wd, Phim_base, alphastar
     else:
-        print("")
+        print("Returning max iteration result")
         return x, ibreak, np.sqrt(phid/len(b))
 
 

akvo/tressel/nonlinearinv.py

     bounds = np.zeros((len(x0),2))
     bounds[:,0] = x0*0.75
     bounds[:,1] = x0*1.25
+    #bounds[:,0] = 0
+    #bounds[:,1] = np.max(x0)*1.25
     return so.minimize(PHI, x0, args, 'L-BFGS-B', bounds=bounds)     # Works well 
     #return so.minimize(PHI, x0, args, 'Powell', bounds=bounds)       # Slow but works 
     #return so.minimize(PHI, x0, args, 'trust-constr', bounds=bounds) # very Slow