Changeset 4709


Ignore:
Timestamp:
Jan 4, 2021 10:41:46 AM (9 months ago)
Author:
vondreele
Message:

new (BT mostly) version of HessianLSQ & additions; fix nan problem in max shift/esd output
fix file/image name/replace handling for gain maps

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimgGUI.py

    r4698 r4709  
    185185            GainMap = np.where(GainMap < 800,1000,GainMap)
    186186            G2IO.PutG2Image(newimagefile,[],data,Npix,GainMap)
    187             Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text='IMG '+os.path.split(newimagefile)[1])
    188             G2frame.GPXtree.SetItemPyData(Id,[Npix,newimagefile])
    189             G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Comments'),[])
    190             G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Image Controls'),Data)
    191             G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Masks'),masks)
     187            GMname = 'IMG '+os.path.split(newimagefile)[1]
     188            Id = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,GMname)
     189            if not Id:           
     190                Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=GMname)
     191                G2frame.GPXtree.SetItemPyData(Id,[Npix,newimagefile])
     192                G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Comments'),[])
     193                G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Image Controls'),Data)
     194                G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Masks'),masks)
     195            else:
     196                G2frame.GPXtree.SetItemPyData(Id,[Npix,newimagefile])
     197                G2frame.GPXtree.SetItemPyData(G2gd.GetGPXtreeItemId(G2frame,Id,'Comments'),[])
     198                G2frame.GPXtree.SetItemPyData(G2gd.GetGPXtreeItemId(G2frame,Id,'Image Controls'),Data)
     199                G2frame.GPXtree.SetItemPyData(G2gd.GetGPXtreeItemId(G2frame,Id,'Masks'),masks)
    192200            G2frame.GPXtree.Expand(Id)
    193201            G2frame.GPXtree.SelectItem(Id)      #to show the gain map & put it in the list
  • trunk/GSASIImath.py

    r4688 r4709  
    5656##### Hessian least-squares Levenberg-Marquardt routine
    5757################################################################################
    58 
     58class G2NormException(Exception): pass
     59   
    5960def pinv(a, rcond=1e-15 ):
    6061    '''
     
    101102    s = np.where(s>cutoff,1./s,0.)
    102103    nzero = s.shape[0]-np.count_nonzero(s)
    103 #    res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))
    104104    res = np.dot(vt.T,s[:,nxs]*u.T)
    105105    return res,nzero
    106106
     107def dropTerms(bad, hessian, indices, *vectors):
     108    '''Remove the 'bad' terms from the Hessian and vector
     109
     110    :param tuple bad: a list of variable (row/column) numbers that should be
     111      removed from the hessian and vector. Example: (0,3) removes the 1st and
     112      4th column/row
     113    :param np.array hessian: a square matrix of length n x n
     114    :param np.array indices: the indices of the least-squares vector of length n
     115       referenced to the initial variable list; as this routine is called
     116       multiple times, more terms may be removed from this list
     117    :param additional-args: various least-squares model values, length n
     118
     119    :returns: hessian, indices, vector0, vector1,...  where the lengths are
     120      now n' x n' and n', with n' = n - len(bad)
     121    '''
     122    out = [np.delete(np.delete(hessian,bad,1),bad,0),np.delete(indices,bad)]
     123    for v in vectors:
     124        out.append(np.delete(v,bad))
     125    return out
     126       
     127def setHcorr(info,Amat,xtol,problem=False):
     128    Adiag = np.sqrt(np.diag(Amat))
     129    if np.any(np.abs(Adiag) < 1.e-14): raise G2NormException  # test for any hard singularities
     130    Anorm = np.outer(Adiag,Adiag)
     131    Amat = Amat/Anorm
     132    Bmat,Nzeros = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
     133    Bmat = Bmat/Anorm
     134    sig = np.sqrt(np.diag(Bmat))
     135    xvar = np.outer(sig,np.ones_like(sig))
     136    AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal
     137    if Nzeros or problem: # something is wrong, so report what is found
     138        m = min(0.99,0.99*np.amax(AcovUp))
     139    else:
     140        m = max(0.95,0.99*np.amax(AcovUp))
     141    info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
     142    return Bmat,Nzeros
     143   
    107144def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
    108145    '''
     
    138175        residual standard deviation to get the covariance of the
    139176        parameter estimates -- see curve_fit.
    140       * infodict : dict
    141         a dictionary of optional outputs with the keys:
     177      * infodict : dict, a dictionary of optional outputs with the keys:
    142178
    143179         * 'fvec' : the function evaluated at the output
    144180         * 'num cyc':
    145          * 'nfev':
     181         * 'nfev': number of objective function evaluation calls
    146182         * 'lamMax':
    147          * 'psing':
    148          * 'SVD0':
    149            
    150     '''
    151                
     183         * 'psing': list of variable variables that have been removed from the refinement
     184         * 'SVD0': -1 for singlar matrix, -2 for objective function exception, Nzeroes = # of SVD 0's
     185         * 'Hcorr': list entries (i,j,c) where i & j are of highly correlated variables & c is correlation coeff.
     186    '''
    152187    ifConverged = False
    153188    deltaChi2 = -10.
     
    158193       
    159194    icycle = 0
    160     One = np.ones((n,n))
    161     lam = 10.**lamda
    162     lamMax = lam
     195    lamMax = lam = 0  #  start w/o Marquardt and add it in if needed
    163196    nfev = 0
    164197    if Print:
    165198        G2fil.G2Print(' Hessian Levenberg-Marquardt SVD refinement on %d variables:'%(n))
    166     Lam = np.zeros((n,n))
    167     Xvec = np.zeros(len(x0))
    168     chisq00 = None
     199    XvecAll = np.zeros(n)
     200    try:
     201        M2 = func(x0,*args)
     202    except Exception as Msg:
     203        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     204        G2fil.G2Print('\nouch #0 unable to evaluate initial objective function\nCheck for an invalid parameter value',mode='error')
     205        G2fil.G2Print('Use Calculate/"View LS parms" to list all parameter values\n',mode='warn')
     206        G2fil.G2Print('Error message: '+Msg.msg,mode='warn')
     207        raise Exception('HessianLSQ -- ouch #0: look for invalid parameter (see console)')
     208    chisq00 = np.sum(M2**2) # starting Chi**2
     209    Nobs = len(M2)
     210    if Print:
     211        G2fil.G2Print('initial chi^2 %.5g with %d obs.'%(chisq00,Nobs))
     212    if n == 0:
     213        info = {'num cyc':0,'fvec':M2,'nfev':1,'lamMax':0,'psing':[],'SVD0':0}
     214        info['msg'] = 'no variables: skipping refinement\n'
     215        return [x0,None,info]
    169216    while icycle < maxcyc:
     217        M = M2
    170218        time0 = time.time()
    171         M = func(x0,*args)
    172         Nobs = len(M)
    173219        nfev += 1
    174220        chisq0 = np.sum(M**2)
    175         if chisq00 is None: chisq00 = chisq0
    176         Yvec,Amat = Hess(x0,*args)
     221        Yvec,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables
     222        Amat = copy.copy(AmatAll)
     223        Xvec = copy.copy(XvecAll)
     224        # we could remove vars that were dropped in previous cycles here (use indices), but for
     225        # now let's reset them each cycle in case the singularities go away
     226        indices = range(n)
    177227        Adiag = np.sqrt(np.diag(Amat))
    178         #psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
    179         psing = np.where(np.abs(Adiag) < 1.e-14)[0]
    180         if np.any(psing):                #hard singularity in matrix
    181             return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
    182         Anorm = np.outer(Adiag,Adiag)
     228        psing = np.where(np.abs(Adiag) < 1.e-14)[0]       # find any hard singularities
     229        if len(psing):
     230            G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format(
     231                psing), mode='warn')
     232            Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
     233        Anorm = np.outer(Adiag,Adiag)        # normalize matrix & vector
    183234        Yvec /= Adiag
    184235        Amat /= Anorm
    185         if Print:
    186             G2fil.G2Print('initial chi^2 %.5g on %d obs.'%(chisq0,Nobs))
    187236        chitol = ftol
    188         while True:
    189             Lam = np.eye(Amat.shape[0])*lam
    190             Amatlam = Amat*(One+Lam)
     237        maxdrop = 3 # max number of drop var attempts
     238        loops = 0
     239        while True: #--------- loop as we increase lamda and possibly drop terms
     240            lamMax = max(lamMax,lam)
     241            Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
    191242            try:
    192                 Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
     243                Ainv,Nzeros = pinv(Amatlam,xtol)    # Moore-Penrose SVD inversion
    193244            except nl.LinAlgError:
    194                 psing = list(np.where(np.abs(np.diag(nl.qr(Amatlam)[1])) < 1.e-14)[0])
    195                 G2fil.G2Print('ouch #1 bad SVD inversion; change parameterization for ',psing, mode='error')
    196                 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
    197             Xvec = np.inner(Ainv,Yvec)      #solve
    198             Xvec /= Adiag
    199             M2 = func(x0+Xvec,*args)
     245                loops += 1
     246                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
     247                psing = list(np.where(d < 1.e-14)[0])
     248                if not len(psing): # make sure at least the worst term is removed
     249                    psing = [np.argmin(d)]
     250                G2fil.G2Print('ouch #2 bad SVD inversion; dropping terms for for variable(s) #{}'.
     251                    format(psing), mode='warn')
     252                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
     253                if loops < maxdrop: continue # try again, same lam but fewer vars
     254                G2fil.G2Print('giving up with ouch #2', mode='error')
     255                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
     256                info['psing'] = [i for i in range(n) if i not in indices]
     257
     258                info['msg'] = 'SVD inversion failure\n'
     259                return [x0,None,info]
     260            if Nzeros:     # remove bad vars from the matrix
     261                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
     262                psing = list(np.where(d < 1.e-14)[0])
     263                if not len(psing): # make sure at least the worst term is removed
     264                    psing = [np.argmin(d)]
     265                G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.
     266                    format(Nzeros,psing), mode='warn')
     267                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
     268                Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
     269                Ainv,nz = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
     270                if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop'.format(nz),
     271                    mode='warn')
     272            Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
     273            XvecAll[indices] = Xvec         # expand
     274            try:
     275                M2 = func(x0+XvecAll,*args)
     276            except Exception as Msg:
     277                if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     278                G2fil.G2Print(Msg.msg,mode='warn')
     279                loops += 1
     280                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
     281                G2fil.G2Print('ouch #3 unable to evaluate objective function;',mode='error')
     282                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
     283                info['psing'] = [i for i in range(n) if i not in indices]
     284                try:         # try to report highly correlated parameters from full Hessian
     285                    setHcorr(info,AmatAll,xtol,problem=True)
     286                except nl.LinAlgError:
     287                    G2fil.G2Print('Warning: Hessian too ill-conditioned to get full covariance matrix')
     288                except G2NormException:
     289                    G2fil.G2Print('Warning: Hessian normalization problem')
     290                except Exception as Msg:
     291                    if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     292                    G2fil.G2Print('Error determining highly correlated vars',mode='warn')
     293                    G2fil.G2Print(Msg.msg,mode='warn')
     294                info['msg'] = Msg.msg + '\n\n'
     295                return [x0,None,info]
    200296            nfev += 1
    201297            chisq1 = np.sum(M2**2)
    202298            if chisq1 > chisq0*(1.+chitol):     #TODO put Alan Coehlo's criteria for lambda here?
    203                 lam *= 10.
     299                if lam == 0:
     300                    lam = 10.**lamda  #  set to initial Marquardt term
     301                else:
     302                    lam *= 10.
    204303                if Print:
    205                     G2fil.G2Print('new chi^2 %.5g on %d obs., %d SVD zeros ; matrix modification needed; lambda now %.1e'  \
    206                            %(chisq1,Nobs,Nzeros,lam))
    207             else:
    208                 x0 += Xvec
    209                 lam /= 10.
     304                    G2fil.G2Print(('divergence: chi^2 %.5g on %d obs. (%d SVD zeros)\n'+
     305                        '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam))
     306                if lam > 10.:
     307                    G2fil.G2Print('ouch #4 stuck: chisq-new %.4g > chisq0 %.4g with lambda %.1g'%
     308                        (chisq1,chisq0,lam), mode='warn')
     309                    try:         # report highly correlated parameters from full Hessian, if we can
     310                        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,
     311                         'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}
     312                        info['psing'] = [i for i in range(n) if i not in indices]
     313                        Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True)
     314                        info['SVD0'] = Nzeros
     315                        return [x0,Bmat,info]
     316                    except Exception as Msg:
     317                        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     318                        G2fil.G2Print('Error determining highly correlated vars',mode='warn')
     319                        G2fil.G2Print(Msg.msg,mode='warn')
     320                    maxcyc = -1 # no more cycles
     321                    break
     322                chitol *= 2
     323            else:   # refinement succeeded
     324                x0 += XvecAll
     325                lam /= 10. # drop lam on next cycle
    210326                break
    211             if lam > 10.:
    212                 G2fil.G2Print('ouch #3 chisq1 %g.4 stuck > chisq0 %g.4'%(chisq1,chisq0), mode='warn')
    213                 break
    214             chitol *= 2
    215         lamMax = max(lamMax,lam)
     327        # complete current cycle
    216328        deltaChi2 = (chisq0-chisq1)/chisq0
    217         if Print:
    218             G2fil.G2Print(' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g'%(
    219                 icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
     329        if Print and n-len(indices):
     330            G2fil.G2Print(
     331                'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d'%
     332                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices)))
     333        elif Print:
     334            G2fil.G2Print(
     335                'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g'%
     336                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
    220337        Histograms = args[0][0]
    221338        if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle)   # update plot
     
    225342            break
    226343        icycle += 1
    227     M = func(x0,*args)
     344    #----------------------- refinement complete, compute Covariance matrix w/o Levenberg-Marquardt
    228345    nfev += 1
    229     Yvec,Amat = Hess(x0,*args)
     346    try:
     347        M = func(x0,*args)
     348    except Exception as Msg:
     349        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     350        G2fil.G2Print(Msg.msg,mode='warn')
     351        G2fil.G2Print('ouch #5 final objective function re-eval failed',mode='error')
     352        psing = [i for i in range(n) if i not in indices]
     353        info = {'num cyc':icycle,'fvec':M2,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-2}
     354        info['msg'] = Msg.msg + '\n'
     355        return [x0,None,info]
     356    chisqf = np.sum(M**2) # ending chi**2
     357    Yvec,Amat = Hess(x0,*args)    # we could save some time and use the last Hessian from the last refinement cycle
     358    psing_prev = [i for i in range(n) if i not in indices] # save dropped vars
     359    indices = range(n)
     360    info = {}
     361    try:         # report highly correlated parameters from full Hessian, if we can
     362        Bmat,Nzeros = setHcorr(info,Amat,xtol,problem=False)
     363        info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,'psing':psing_prev,
     364            'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
     365        return [x0,Bmat,info]
     366    except nl.LinAlgError:
     367        G2fil.G2Print('Warning: Hessian too ill-conditioned to get full covariance matrix')
     368    except G2NormException:
     369        G2fil.G2Print('Warning: Hessian normalization problem')
     370    except Exception as Msg:
     371        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     372        G2fil.G2Print('Error determining highly correlated vars',mode='warn')
     373        G2fil.G2Print(Msg.msg,mode='warn')
     374    # matrix above inversion failed, drop previously removed variables & try again
     375    Amat, indices, Yvec = dropTerms(psing_prev, Amat, indices, Yvec)
    230376    Adiag = np.sqrt(np.diag(Amat))
    231377    Anorm = np.outer(Adiag,Adiag)
    232     Lam = np.eye(Amat.shape[0])*lam
    233     Amatlam = Amat/Anorm       
     378    Amat = Amat/Anorm
    234379    try:
    235         Bmat,Nzero = pinv(Amatlam,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
    236         if Print: G2fil.G2Print('Found %d SVD zeros'%(Nzero), mode='warn')
     380        Bmat,Nzeros = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
    237381        Bmat = Bmat/Anorm
    238         return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],
    239             'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,'Xvec':Xvec, 'chisq0':chisq00}]
    240     except nl.LinAlgError:
    241         G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error')
    242         psing = []
    243         if maxcyc:
    244             psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    245         return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1,'Xvec':None, 'chisq0':chisq00}]
     382    except nl.LinAlgError: # this is unexpected. How did we get this far with a singular matrix?
     383        G2fil.G2Print('ouch #5 linear algebra error in making final v-cov matrix', mode='error')
     384        psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e-14)[0])
     385        if not len(psing): # make sure at least the worst term is flagged
     386            psing = [np.argmin(d)]
     387        Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
     388        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':-1,'Xvec':None, 'chisq0':chisqf}
     389        info['psing'] = [i for i in range(n) if i not in indices]
     390        return [x0,None,info]
     391    if Nzeros: # sometimes SVD zeros show up when lam=0
     392        d = np.abs(np.diag(nl.qr(Amat)[1]))
     393        psing = list(np.where(d < 1.e-14)[0])
     394        if len(psing) < Nzeros:
     395            psing = list(np.where(d <= sorted(d)[Nzeros-1])[0])
     396        if Print:
     397            G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+
     398                'Likely problem with variable(s) #%s'%(Nzeros,psing), mode='warn')
     399        Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
     400    # expand Bmat by filling with zeros if columns have been dropped
     401    psing = [i for i in range(n) if i not in indices]
     402    if len(psing):
     403        ins = [j-i for i,j in enumerate(psing)]
     404        Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0)
     405    info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
     406        'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
     407    info['psing'] = [i for i in range(n) if i not in indices]
     408    return [x0,Bmat,info]
    246409           
    247410def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
  • trunk/GSASIIstrMain.py

    r4699 r4709  
    224224                Rvals['Max shft/sig'] = 0.0
    225225            else:
    226                 Rvals['Max shft/sig'] = np.max(Lastshft/sig)
     226                Rvals['Max shft/sig'] = np.max(np.nan_to_num(Lastshft/sig))
    227227            if np.any(np.isnan(sig)) or not sig.shape:
    228228                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***',mode='error')
Note: See TracChangeset for help on using the changeset viewer.