Changeset 4690 for trunk


Ignore:
Timestamp:
Dec 30, 2020 5:52:57 PM (3 years ago)
Author:
toby
Message:

major updates to trial version of HessianLSQ

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r4689 r4690  
    51835183            if lamMax >= 10.:
    51845184                text += '\nWARNING: Steepest descents dominates;'+   \
    5185                 ' minimum may not have been reached\nor result may be false minimum.'+  \
     5185                ' minimum may not have been reached or result may be false minimum.'+  \
    51865186                ' You should reconsider which parameters you refine'
    51875187            dlg2 = wx.MessageDialog(self,text,'Refinement results, Rw =%.3f'%(Rw),wx.OK|wx.CANCEL)

  • trunk/GSASIImath_new.py

    r4683 r4690  
    5656##### Hessian least-squares Levenberg-Marquardt routine
    5757################################################################################
    58 
     58class G2NormException(Exception): pass
     59   
    5960def pinv(a, rcond=1e-15 ):
    6061    '''
    … …  
    124125                np.delete(yvector,bad),
    125126                np.delete(indices,bad))
     127
     128def setHcorr(info,Amat,xtol,problem=False):
     129    Adiag = np.sqrt(np.diag(Amat))
     130    if np.any(np.abs(Adiag) < 1.e-14): raise G2NormException  # test for any hard singularities
     131    Anorm = np.outer(Adiag,Adiag)
     132    Amat = Amat/Anorm
     133    Bmat,Nzeros = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
     134    Bmat = Bmat/Anorm
     135    sig = np.sqrt(np.diag(Bmat))
     136    xvar = np.outer(sig,np.ones_like(sig))
     137    AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal
     138    if Nzeros or problem: # something is wrong, so report what is found
     139        m = min(0.99,0.99*np.amax(AcovUp))
     140    else:
     141        m = max(0.95,0.99*np.amax(AcovUp))
     142    info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
     143    return Bmat
    126144   
    127145def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
    … …  
    195213    if Print:
    196214        G2fil.G2Print('initial chi^2 %.5g with %d obs.'%(chisq00,Nobs))
     215    if n == 0:
     216        info = {'num cyc':0,'fvec':M2,'nfev':1,'lamMax':0,'psing':[],'SVD0':0}
     217        info['msg'] = 'no variables, no refinement\n'
     218        return [x0,None,info]
    197219    while icycle < maxcyc:
    198220        M = M2
    … …  
    208230        Adiag = np.sqrt(np.diag(Amat))
    209231        psing = np.where(np.abs(Adiag) < 1.e-14)[0]         # remove any hard singularities
    210         if psing:
     232        if len(psing):
    211233            G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format(psing), mode='warn')
    212234            Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
    213235            Adiag = np.sqrt(np.diag(Amat))
    214         Anorm = np.outer(Adiag[indices],Adiag[indices])        # normalize matrix & vector
    215         Yvec /= Adiag[indices]
     236        Anorm = np.outer(Adiag,Adiag)        # normalize matrix & vector
     237        Yvec /= Adiag
    216238        Amat /= Anorm
    217239        chitol = ftol
    … …  
    226248                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
    227249                psing = list(np.where(d < 1.e-14)[0])
    228                 if not psing: # make sure at least the worst term is removed
     250                if not len(psing): # make sure at least the worst term is removed
    229251                    psing = [np.argmin(d)]
    230252                G2fil.G2Print('ouch #2 bad SVD inversion; dropping terms for for variable(s) #{}'.
    231253                                  format(psing), mode='warn')
    232254                Amatlam, Xvec, Yvec, indices = dropTerms(Amatlam, Xvec, Yvec, indices, psing)
    233                 Anorm = np.outer(Adiag[indices],Adiag[indices])
    234255                if loops < maxdrop: continue
    235256                G2fil.G2Print('giving up with ouch #2', mode='error')
    … …  
    240261                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
    241262                psing = list(np.where(d < 1.e-14)[0])
    242                 if not psing: # make sure at least the worst term is removed
     263                if not len(psing): # make sure at least the worst term is removed
    243264                    psing = [np.argmin(d)]
    244265                G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.
    … …  
    246267                Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
    247268                Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
    248                 Anorm = np.outer(Adiag[indices],Adiag[indices])
    249269                Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
    250             Xvec = np.inner(Ainv,Yvec)/Adiag[indices]      #solve for LS terms
     270            Adiag = np.sqrt(np.diag(Amat))
     271            Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
    251272            XvecAll[indices] = Xvec         # expand
    252273            try:
    … …  
    261282                G2fil.G2Print(Msg.msg,mode='warn')
    262283                # dropping variables does not seem to help
    263                 if psing:
     284                if len(psing):
    264285                    G2fil.G2Print('ouch #3 unable to evaluate objective function;'+
    265286                                    '\n\tdropping terms for variable(s) #{}'.format(psing),
    … …  
    274295                info['psing'] = [i for i in range(n) if i not in indices]
    275296                try:         # report highly correlated parameters from full Hessian, if we can
    276                     Adiag = np.sqrt(np.diag(AmatAll))
    277                     Anorm = np.outer(Adiag,Adiag)
    278                     if np.where(np.abs(Adiag) < 1.e-14)[0]: raise Exception     # test for any hard singularities
    279                     Amat = AmatAll/Anorm
    280                     Bmat,Nzeros = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
    281                     Bmat = Bmat/Anorm
    282                     sig = np.sqrt(np.diag(Bmat))
    283                     xvar = np.outer(sig,np.ones_like(sig))
    284                     AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal
    285                     m = min(0.99,0.99*np.amax(AcovUp))
    286                     info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
    287                 except:
    288                     pass
     297                    setHcorr(info,AmatAll,xtol,problem=True)
     298                except nl.LinAlgError:
     299                    G2fil.G2Print('Warning: Hessian too ill-conditioned to get full covariance matrix')
     300                except G2NormException:
     301                    G2fil.G2Print('Warning: Hessian normalization problem')
     302                except Exception as Msg:
     303                    if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     304                    G2fil.G2Print('Error determining highly correlated vars',mode='warn')
     305                    G2fil.G2Print(Msg.msg,mode='warn')
    289306                info['msg'] = Msg.msg + '\n'
    290307                return [x0,None,info]
    … …  
    298315                                    '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam))
    299316                if lam > 10.:
    300                     G2fil.G2Print('ouch #4 stuck: chisq-new %g.4 > chisq0 %g.4 with lambda %g.1'%
     317                    G2fil.G2Print('ouch #4 stuck: chisq-new %.4g > chisq0 %.4g with lambda %.1g'%
    301318                                  (chisq1,chisq0,lam), mode='warn')
     319                    try:         # report highly correlated parameters from full Hessian, if we can
     320                        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
     321                         'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}
     322                        info['psing'] = [i for i in range(n) if i not in indices]
     323                        Bmat = setHcorr(info,AmatAll,xtol,problem=True)
     324                        return [x0,Bmat,info]
     325                    except Exception as Msg:
     326                        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     327                        G2fil.G2Print('Error determining highly correlated vars',mode='warn')
     328                        G2fil.G2Print(Msg.msg,mode='warn')
    302329                    maxcyc = -1 # no more cycles
    303330                    break
    304331                chitol *= 2
    305332            else:   # refinement succeeded
    306                 x0[indices] += Xvec
     333                x0 += XvecAll
    307334                lam /= 10. # drop lam on next cycle
    308335                break
    … …  
    343370    info = {}
    344371    try:         # report highly correlated parameters from full Hessian, if we can
    345         Adiag = np.sqrt(np.diag(Amat))
    346         Anorm = np.outer(Adiag,Adiag)
    347         if np.where(np.abs(Adiag) < 1.e-14)[0]: raise Exception     # test for any hard singularities
    348         AmatN = Amat/Anorm
    349         Bmat,Nzeros = pinv(AmatN,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
    350         Bmat = Bmat/Anorm
    351         sig = np.sqrt(np.diag(Bmat))
    352         xvar = np.outer(sig,np.ones_like(sig))
    353         AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal
    354         m = min(0.99,0.99*np.amax(AcovUp))
    355         info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
     372        Bmat = setHcorr(info,Amat,xtol,problem=False)
    356373        info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,'psing':psing_prev,
    357374                         'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
    … …  
    359376    except nl.LinAlgError:
    360377        G2fil.G2Print('Warning: Hessian too ill-conditioned to get full covariance matrix')
     378    except G2NormException:
     379        G2fil.G2Print('Warning: Hessian normalization problem')
    361380    except Exception as Msg:
    362381        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     382        G2fil.G2Print('Error determining highly correlated vars',mode='warn')
    363383        G2fil.G2Print(Msg.msg,mode='warn')
    364         G2fil.G2Print('Warning: Hessian normalization problem')
    365     # matrix inversion failed drop previously removed variables
    366     Xvec = np.zeros(n) # dummy variable
    367     Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing_prev)
     384    # matrix above inversion failed, drop previously removed variables
     385    dummy = np.zeros(n)
     386    Amat, Xvec, Yvec, indices = dropTerms(Amat, dummy, Yvec, indices, psing_prev)
    368387    Adiag = np.sqrt(np.diag(Amat))
    369388    Anorm = np.outer(Adiag,Adiag)
    … …  
    375394        G2fil.G2Print('ouch #5 linear algebra error in making final v-cov matrix', mode='error')
    376395        psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e-14)[0])
    377         if not psing: # make sure at least the worst is returned
     396        if not len(psing): # make sure at least the worst is returned
    378397            psing = [np.argmin(d)]
    379398        Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
    … …  
    391410        Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
    392411    # expand Bmat by filling with zeros
    393     if psing:
     412    if len(psing):
    394413        Bmat = np.insert(np.insert(Bmat,psing,0,1),psing,0,0)
    395414    info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
Note: See TracChangeset for help on using the changeset viewer.