Changeset 4789


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

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

Location:
trunk
Files:
1 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r4780 r4789  
    51945194                text += '\nMax shift/sigma={:.3f}\n'.format(Rvals['Max shft/sig'])
    51955195            if 'msg' in Rvals: text += '\n' + Rvals['msg'] + '\n'
     5196            if 'Aborted' in Rvals:
     5197                text += '\ERROR: Minimizer failed to reduce chi**2!\n'
    51965198            if lamMax >= 10.:
    51975199                text += '\nWARNING: Steepest descents dominates;'+   \

  • 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           

  • trunk/GSASIIphsGUI.py

    r4783 r4789  
    1069910699            Sytsym,Mult = G2spc.SytSym(rbObj['Orig'][0],data['General']['SGData'])[:2]
    1070010700            sytsymtxt = wx.StaticText(RigidBodies,label='Origin site symmetry: %s, multiplicity: %d '%(Sytsym,Mult))
    10701             mainSizer.Add(sytsymtxt,0,WACV)
     10701            mainSizer.Add(sytsymtxt)
    1070210702            OriSizer1 = wx.FlexGridSizer(0,5,5,5)
    1070310703            if len(atomData):

  • trunk/GSASIIplot.py

    r4788 r4789  
    131131    import mpl_toolkits.mplot3d.axes3d as mp3d
    132132    from scipy.ndimage.interpolation import map_coordinates
    133 except ImportError:
     133except (ImportError, ValueError):
    134134    pass
    135135import GSASIIpath

  • trunk/GSASIIpwd.py

    r4672 r4789  
    22# -*- coding: utf-8 -*-
    33'''
    4 *GSASII powder calculation module*
    5 ==================================
     4*GSASII powder calculation module (GSASIIpwd)*
     5==============================================
    66
    77'''
    … …  
    21442144                return -M           #abort!!
    21452145        return M
    2146        
     2146
     2147    # beginning of DoPeakFit
    21472148    if controls:
    21482149        Ftol = controls['min dM/M']

  • trunk/GSASIIstrMain.py

    r4783 r4789  
    188188            ncyc = result[2]['num cyc']+1
    189189            Rvals['lamMax'] = result[2]['lamMax']
     190            if 'Ouch#4' in  result[2]:
     191                Rvals['Aborted'] = True
    190192            if 'msg' in result[2]:
    191193                Rvals['msg'] = result[2]['msg']

  • trunk/docs/source/index.rst

    r4335 r4789  
    4141  G2tools.rst
    4242
     43.. only:: html
     44
     45          `General Index <./genindex.html>`_
     46
     47          `Module Index <./py-modindex.html>`_
     48

  • trunk/fprime.py

    r4672 r4789  
    183183        self.SpinButton = wx.SpinButton(id=wxID_SPINBUTTON, parent=panel,
    184184              size=wx.Size(25,24), style=wx.SP_VERTICAL | wx.SP_ARROW_KEYS)
    185         slideSizer.Add(self.SpinButton,0,wx.ALIGN_RIGHT)
     185        slideSizer.Add(self.SpinButton)
    186186        self.SpinButton.SetRange(-1,1)
    187187        self.SpinButton.SetValue(0)
Note: See TracChangeset for help on using the changeset viewer.