Changeset 4789 for trunk/GSASIImath.py

Timestamp:

Feb 1, 2021 1:36:50 PM (2 years ago)

Author:

toby

Message:

switch to new HessianLSQ; Ouch#4 msg in GUI; misc wx4.1 fixes; docs improvements; switch to new HessianLSQ

File:

• trunk/GSASIImath.py (modified) (13 diffs)

trunk/GSASIImath.py

@@ -126 +126 @@
         
 def setHcorr(info,Amat,xtol,problem=False):
+    '''Find & report high correlation terms in covariance matrix
+    '''
     Adiag = np.sqrt(np.diag(Amat))
     if np.any(np.abs(Adiag) < 1.e-14): raise G2NormException  # test for any hard singularities
@@ -141 +143 @@
     info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
     return Bmat,Nzeros
-    
+
+def setSVDwarn(info,Amat,Nzeros,indices):
+    '''Find & report terms causing SVN zeros
+    '''
+    if Nzeros == 0: return
+    d = np.abs(np.diag(nl.qr(Amat)[1]))
+    svdsing = list(np.where(d < 1.e-14)[0])
+    if len(svdsing) < Nzeros: # try to get the Nzeros worst terms
+        svdsing = list(np.where(d <= sorted(d)[Nzeros-1])[0])
+
+
+    if not len(svdsing): # make sure at least the worst term is shown
+        svdsing = [np.argmin(d)]
+    info['SVDsing'] = [indices[i] for i in svdsing]
+
 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
     '''
@@ -187 +203 @@
     ifConverged = False
     deltaChi2 = -10.
-    x0 = np.array(x0, ndmin=1, dtype=np.float64)      #might be redundant?
+    x0 = np.array(x0, ndmin=1, dtype=np.float64) # make sure that x0 float 1-D
+    # array (in case any parameters were set to int)
     n = len(x0)
     if type(args) != type(()):
@@ -193 +210 @@
         
     icycle = 0
+    AmatAll = None # changed only if cycles > 0
     lamMax = lam = 0  #  start w/o Marquardt and add it in if needed
     nfev = 0
@@ -221 +239 @@
         nfev += 1
         chisq0 = np.sum(M**2)
-        Yvec,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables
+        YvecAll,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables
+        Yvec = copy.copy(YvecAll)
         Amat = copy.copy(AmatAll)
         Xvec = copy.copy(XvecAll)
@@ -257 +276 @@
                 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
                 info['psing'] = [i for i in range(n) if i not in indices]
-
                 info['msg'] = 'SVD inversion failure\n'
                 return [x0,None,info]
-            if Nzeros:     # remove bad vars from the matrix
-                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
-                psing = list(np.where(d < 1.e-14)[0])
-                if not len(psing): # make sure at least the worst term is removed
-                    psing = [np.argmin(d)]
-                G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.
-                    format(Nzeros,psing), mode='warn')
-                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
-                Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
-                Ainv,nz = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
-                if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop'.format(nz),
-                    mode='warn')
             Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
             XvecAll[indices] = Xvec         # expand
@@ -295 +301 @@
                     G2fil.G2Print(Msg.msg,mode='warn')
                 info['msg'] = Msg.msg + '\n\n'
+                setSVDwarn(info,Amatlam,Nzeros,indices)
                 return [x0,None,info]
             nfev += 1
@@ -311 +318 @@
                     try:         # report highly correlated parameters from full Hessian, if we can
                         info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,
-                         'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}
+                            'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll,
+                            'chisq0':chisq00, 'Ouch#4':True}
                         info['psing'] = [i for i in range(n) if i not in indices]
                         Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True)
                         info['SVD0'] = Nzeros
+                        setSVDwarn(info,Amatlam,Nzeros,indices)
                         return [x0,Bmat,info]
                     except Exception as Msg:
@@ -331 +340 @@
         if Print and n-len(indices):
             G2fil.G2Print(
-                'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d'%
-                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices)))
+                'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d, SVD=%d'%
+                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices),Nzeros))
         elif Print:
             G2fil.G2Print(
-                'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g'%
-                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
+                'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g, SVD=%d'%
+                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,Nzeros))
         Histograms = args[0][0]
         if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle)   # update plot
@@ -355 +364 @@
         info = {'num cyc':icycle,'fvec':M2,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-2}
         info['msg'] = Msg.msg + '\n'
+        setSVDwarn(info,Amatlam,Nzeros,indices)
         return [x0,None,info]
     chisqf = np.sum(M**2) # ending chi**2
-    Yvec,Amat = Hess(x0,*args)    # we could save some time and use the last Hessian from the last refinement cycle
     psing_prev = [i for i in range(n) if i not in indices] # save dropped vars
+    if AmatAll is None: # Save some time and use Hessian from the last refinement cycle
+        Yvec,Amat = Hess(x0,*args)
+    else:
+        Yvec = copy.copy(YvecAll)
+        Amat = copy.copy(AmatAll)
     indices = range(n)
     info = {}
@@ -366 +380 @@
             'Converged':ifConverged, 'DelChi2':deltaChi2, 'chisq0':chisq00})
         if icycle > 0: info.update({'Xvec':XvecAll})
+        setSVDwarn(info,Amatlam,Nzeros,indices)
+        # expand Bmat by filling with zeros if columns have been dropped
+        if len(psing_prev):
+            ins = [j-i for i,j in enumerate(psing_prev)]
+            Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0)
         return [x0,Bmat,info]
     except nl.LinAlgError:
@@ -392 +411 @@
         info['psing'] = [i for i in range(n) if i not in indices]
         return [x0,None,info]
-    if Nzeros: # sometimes SVD zeros show up when lam=0
-        d = np.abs(np.diag(nl.qr(Amat)[1]))
-        psing = list(np.where(d < 1.e-14)[0])
-        if len(psing) < Nzeros:
-            psing = list(np.where(d <= sorted(d)[Nzeros-1])[0])
-        if Print:
-            G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+
-                'Likely problem with variable(s) #%s'%(Nzeros,psing), mode='warn')
-        Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
     # expand Bmat by filling with zeros if columns have been dropped
     psing = [i for i in range(n) if i not in indices]
@@ -409 +419 @@
         'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
     info['psing'] = [i for i in range(n) if i not in indices]
+    setSVDwarn(info,Amat,Nzeros,indices)
     return [x0,Bmat,info]