Changeset 4789 for trunk/GSASIImath.py


Ignore:
Timestamp:
Feb 1, 2021 1:36:50 PM (10 months ago)
Author:
toby
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r4770 r4789  
    126126       
    127127def setHcorr(info,Amat,xtol,problem=False):
     128    '''Find & report high correlation terms in covariance matrix
     129    '''
    128130    Adiag = np.sqrt(np.diag(Amat))
    129131    if np.any(np.abs(Adiag) < 1.e-14): raise G2NormException  # test for any hard singularities
     
    141143    info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
    142144    return Bmat,Nzeros
    143    
     145
     146def setSVDwarn(info,Amat,Nzeros,indices):
     147    '''Find & report terms causing SVN zeros
     148    '''
     149    if Nzeros == 0: return
     150    d = np.abs(np.diag(nl.qr(Amat)[1]))
     151    svdsing = list(np.where(d < 1.e-14)[0])
     152    if len(svdsing) < Nzeros: # try to get the Nzeros worst terms
     153        svdsing = list(np.where(d <= sorted(d)[Nzeros-1])[0])
     154
     155
     156    if not len(svdsing): # make sure at least the worst term is shown
     157        svdsing = [np.argmin(d)]
     158    info['SVDsing'] = [indices[i] for i in svdsing]
     159
    144160def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
    145161    '''
     
    187203    ifConverged = False
    188204    deltaChi2 = -10.
    189     x0 = np.array(x0, ndmin=1, dtype=np.float64)      #might be redundant?
     205    x0 = np.array(x0, ndmin=1, dtype=np.float64) # make sure that x0 float 1-D
     206    # array (in case any parameters were set to int)
    190207    n = len(x0)
    191208    if type(args) != type(()):
     
    193210       
    194211    icycle = 0
     212    AmatAll = None # changed only if cycles > 0
    195213    lamMax = lam = 0  #  start w/o Marquardt and add it in if needed
    196214    nfev = 0
     
    221239        nfev += 1
    222240        chisq0 = np.sum(M**2)
    223         Yvec,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables
     241        YvecAll,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables
     242        Yvec = copy.copy(YvecAll)
    224243        Amat = copy.copy(AmatAll)
    225244        Xvec = copy.copy(XvecAll)
     
    257276                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
    258277                info['psing'] = [i for i in range(n) if i not in indices]
    259 
    260278                info['msg'] = 'SVD inversion failure\n'
    261279                return [x0,None,info]
    262             if Nzeros:     # remove bad vars from the matrix
    263                 d = np.abs(np.diag(nl.qr(Amatlam)[1]))
    264                 psing = list(np.where(d < 1.e-14)[0])
    265                 if not len(psing): # make sure at least the worst term is removed
    266                     psing = [np.argmin(d)]
    267                 G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.
    268                     format(Nzeros,psing), mode='warn')
    269                 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
    270                 Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
    271                 Ainv,nz = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
    272                 if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop'.format(nz),
    273                     mode='warn')
    274280            Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
    275281            XvecAll[indices] = Xvec         # expand
     
    295301                    G2fil.G2Print(Msg.msg,mode='warn')
    296302                info['msg'] = Msg.msg + '\n\n'
     303                setSVDwarn(info,Amatlam,Nzeros,indices)
    297304                return [x0,None,info]
    298305            nfev += 1
     
    311318                    try:         # report highly correlated parameters from full Hessian, if we can
    312319                        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,
    313                          'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}
     320                            'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll,
     321                            'chisq0':chisq00, 'Ouch#4':True}
    314322                        info['psing'] = [i for i in range(n) if i not in indices]
    315323                        Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True)
    316324                        info['SVD0'] = Nzeros
     325                        setSVDwarn(info,Amatlam,Nzeros,indices)
    317326                        return [x0,Bmat,info]
    318327                    except Exception as Msg:
     
    331340        if Print and n-len(indices):
    332341            G2fil.G2Print(
    333                 'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d'%
    334                 (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices)))
     342                'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d, SVD=%d'%
     343                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices),Nzeros))
    335344        elif Print:
    336345            G2fil.G2Print(
    337                 'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g'%
    338                 (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
     346                'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g, SVD=%d'%
     347                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,Nzeros))
    339348        Histograms = args[0][0]
    340349        if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle)   # update plot
     
    355364        info = {'num cyc':icycle,'fvec':M2,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-2}
    356365        info['msg'] = Msg.msg + '\n'
     366        setSVDwarn(info,Amatlam,Nzeros,indices)
    357367        return [x0,None,info]
    358368    chisqf = np.sum(M**2) # ending chi**2
    359     Yvec,Amat = Hess(x0,*args)    # we could save some time and use the last Hessian from the last refinement cycle
    360369    psing_prev = [i for i in range(n) if i not in indices] # save dropped vars
     370    if AmatAll is None: # Save some time and use Hessian from the last refinement cycle
     371        Yvec,Amat = Hess(x0,*args)
     372    else:
     373        Yvec = copy.copy(YvecAll)
     374        Amat = copy.copy(AmatAll)
    361375    indices = range(n)
    362376    info = {}
     
    366380            'Converged':ifConverged, 'DelChi2':deltaChi2, 'chisq0':chisq00})
    367381        if icycle > 0: info.update({'Xvec':XvecAll})
     382        setSVDwarn(info,Amatlam,Nzeros,indices)
     383        # expand Bmat by filling with zeros if columns have been dropped
     384        if len(psing_prev):
     385            ins = [j-i for i,j in enumerate(psing_prev)]
     386            Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0)
    368387        return [x0,Bmat,info]
    369388    except nl.LinAlgError:
     
    392411        info['psing'] = [i for i in range(n) if i not in indices]
    393412        return [x0,None,info]
    394     if Nzeros: # sometimes SVD zeros show up when lam=0
    395         d = np.abs(np.diag(nl.qr(Amat)[1]))
    396         psing = list(np.where(d < 1.e-14)[0])
    397         if len(psing) < Nzeros:
    398             psing = list(np.where(d <= sorted(d)[Nzeros-1])[0])
    399         if Print:
    400             G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+
    401                 'Likely problem with variable(s) #%s'%(Nzeros,psing), mode='warn')
    402         Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
    403413    # expand Bmat by filling with zeros if columns have been dropped
    404414    psing = [i for i in range(n) if i not in indices]
     
    409419        'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
    410420    info['psing'] = [i for i in range(n) if i not in indices]
     421    setSVDwarn(info,Amat,Nzeros,indices)
    411422    return [x0,Bmat,info]
    412423           
Note: See TracChangeset for help on using the changeset viewer.