Changeset 4780 for trunk


Ignore:
Timestamp:
Jan 17, 2021 3:12:29 PM (2 years ago)
Author:
toby
Message:

separate hard & SVD singularities, new code in GSASIImath_new.py not GSASIImath.py

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r4764 r4780  
    52165216            else:
    52175217                msg = ''
    5218             # if len(Rvals['psing']) == 1:
    5219             #     h = 'This variable appears'
    5220             #     # that = 'that variable'
    5221             # else:
    5222             #     h = 'These variables appear'
    5223             #     # that = 'those variables'
    5224             # msg += h + ' to cause a singular matrix:\n\t'
    5225             # for i,var in enumerate(Rvals['psing']):
    5226             #     if i: msg += ', '
    5227             #     msg += varyList[var]
    5228             # # msg += '\n\nRefine again with '+that+' frozen?'
    52295218            result = wx.ID_NO
    52305219            try:
    5231                 # dlg = wx.MessageDialog(self, msg,'Refine again?',
    5232                 #     wx.YES_NO | wx.ICON_QUESTION)
    52335220                dlg = wx.MessageDialog(self, msg,'Note singularities',
    52345221                    wx.OK)
    … …  
    52385225            finally:
    52395226                dlg.Destroy()
    5240             # if result != wx.ID_YES: return
    5241             # Controls = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,self.root, 'Controls'))
    5242             # if 'parmFrozen' not in Controls:
    5243             #     Controls['parmFrozen'] = {}
    5244             # if 'FrozenList' not in Controls['parmFrozen']:
    5245             #     Controls['parmFrozen']['FrozenList'] = []
    5246             # Controls['parmFrozen']['FrozenList'] += [
    5247             #     G2obj.G2VarObj(varyList[i]) for i in Rvals['plist']]
    5248             # wx.CallAfter(self.OnRefine,event)
    52495227        else:
    52505228            self.ErrorDialog('Refinement error',Rvals['msg'])

  • trunk/GSASIImath_new.py

    r4759 r4780  
    22#GSASIImath - major mathematics routines
    33########### SVN repository information ###################
    4 # $Date: 2021-01-04 10:41:46 -0600 (Mon, 04 Jan 2021) $
     4# $Date: 2021-01-14 14:10:10 -0600 (Thu, 14 Jan 2021) $
    55# $Author: vondreele $
    6 # $Revision: 4709 $
     6# $Revision: 4767 $
    77# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIImath.py $
    8 # $Id: GSASIImath.py 4709 2021-01-04 16:41:46Z vondreele $
     8# $Id: GSASIImath.py 4767 2021-01-14 20:10:10Z vondreele $
    99########### SVN repository information ###################
    1010'''
    … …  
    2424import copy
    2525import GSASIIpath
    26 GSASIIpath.SetVersionNumber("$Revision: 4709 $")
     26GSASIIpath.SetVersionNumber("$Revision: 4767 $")
    2727import GSASIIElem as G2el
    2828import GSASIIlattice as G2lat
    … …  
    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)      #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
    … …  
    213231        info = {'num cyc':0,'fvec':M2,'nfev':1,'lamMax':0,'psing':[],'SVD0':0}
    214232        info['msg'] = 'no variables: skipping refinement\n'
    215         return [x0,None,info]
     233        info.update({'Converged':True, 'DelChi2':0, 'Xvec':None, 'chisq0':chisq00})
     234        return [x0,np.array([]),info]
     235    indices = range(n)
    216236    while icycle < maxcyc:
    217237        M = M2
    … …  
    219239        nfev += 1
    220240        chisq0 = np.sum(M**2)
    221         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)
    222243        Amat = copy.copy(AmatAll)
    223244        Xvec = copy.copy(XvecAll)
    … …  
    255276                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
    256277                info['psing'] = [i for i in range(n) if i not in indices]
    257 
    258278                info['msg'] = 'SVD inversion failure\n'
    259279                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')
    272280            Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
    273281            XvecAll[indices] = Xvec         # expand
    … …  
    293301                    G2fil.G2Print(Msg.msg,mode='warn')
    294302                info['msg'] = Msg.msg + '\n\n'
     303                setSVDwarn(info,Amatlam,Nzeros,indices)
    295304                return [x0,None,info]
    296305            nfev += 1
    … …  
    313322                        Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True)
    314323                        info['SVD0'] = Nzeros
     324                        setSVDwarn(info,Amatlam,Nzeros,indices)
    315325                        return [x0,Bmat,info]
    316326                    except Exception as Msg:
    … …  
    329339        if Print and n-len(indices):
    330340            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)))
     341                'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d, SVD=%d'%
     342                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices),Nzeros))
    333343        elif Print:
    334344            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))
     345                'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g, SVD=%d'%
     346                (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,Nzeros))
    337347        Histograms = args[0][0]
    338348        if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle)   # update plot
    … …  
    353363        info = {'num cyc':icycle,'fvec':M2,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-2}
    354364        info['msg'] = Msg.msg + '\n'
     365        setSVDwarn(info,Amatlam,Nzeros,indices)
    355366        return [x0,None,info]
    356367    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
    358368    psing_prev = [i for i in range(n) if i not in indices] # save dropped vars
     369    if AmatAll is None: # Save some time and use Hessian from the last refinement cycle
     370        Yvec,Amat = Hess(x0,*args)
     371    else:
     372        Yvec = copy.copy(YvecAll)
     373        Amat = copy.copy(AmatAll)
    359374    indices = range(n)
    360375    info = {}
    … …  
    362377        Bmat,Nzeros = setHcorr(info,Amat,xtol,problem=False)
    363378        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})
     379            'Converged':ifConverged, 'DelChi2':deltaChi2, 'chisq0':chisq00})
     380        if icycle > 0: info.update({'Xvec':XvecAll})
     381        setSVDwarn(info,Amatlam,Nzeros,indices)
     382        # expand Bmat by filling with zeros if columns have been dropped
     383        if len(psing_prev):
     384            ins = [j-i for i,j in enumerate(psing_prev)]
     385            Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0)
    365386        return [x0,Bmat,info]
    366387    except nl.LinAlgError:
    … …  
    389410        info['psing'] = [i for i in range(n) if i not in indices]
    390411        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)
    400412    # expand Bmat by filling with zeros if columns have been dropped
    401413    psing = [i for i in range(n) if i not in indices]
    … …  
    406418        'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
    407419    info['psing'] = [i for i in range(n) if i not in indices]
     420    setSVDwarn(info,Amat,Nzeros,indices)
    408421    return [x0,Bmat,info]
    409422           
    … …  
    554567##### Atom manipulations
    555568################################################################################
     569
     570def getAtomPtrs(data,draw=False):
     571    ''' get atom data pointers cx,ct,cs,cia in Atoms or Draw Atoms lists
     572    NB:may not match column numbers in displayed table
     573   
     574        param: dict: data phase data structure
     575        draw: boolean True if Draw Atoms list pointers are required
     576        return: cx,ct,cs,cia pointers to atom xyz, type, site sym, uiso/aniso flag
     577    '''
     578    if draw:
     579        return data['Drawing']['atomPtrs']
     580    else:
     581        return data['General']['AtomPtrs']
    556582
    557583def FindMolecule(ind,generalData,atomData):                    #uses numpy & masks - very fast even for proteins!
    … …  
    749775    generalData = data['General']
    750776    SGData = generalData['SGData']
    751     cx,ct,cs,cia = generalData['AtomPtrs']
     777    cx,ct,cs,cia = getAtomPtrs(data)
    752778    drawingData = data['Drawing']
    753     dcx,dct,dcs,dci = drawingData['atomPtrs']
     779    dcx,dct,dcs,dci = getAtomPtrs(data,True)
    754780    atoms = data['Atoms']
    755781    drawAtoms = drawingData['Atoms']
    … …  
    782808def FindNeighbors(phase,FrstName,AtNames,notName=''):
    783809    General = phase['General']
    784     cx,ct,cs,cia = General['AtomPtrs']
     810    cx,ct,cs,cia = getAtomPtrs(phase)
    785811    Atoms = phase['Atoms']
    786812    atNames = [atom[ct-1] for atom in Atoms]
    … …  
    870896def FindAllNeighbors(phase,FrstName,AtNames,notName='',Orig=None,Short=False):
    871897    General = phase['General']
    872     cx,ct,cs,cia = General['AtomPtrs']
     898    cx,ct,cs,cia = getAtomPtrs(phase)
    873899    Atoms = phase['Atoms']
    874900    atNames = [atom[ct-1] for atom in Atoms]
    … …  
    18821908    SGData = generalData['SGData']
    18831909    SSGData = generalData['SSGData']
    1884     cx,ct,cs,cia = generalData['AtomPtrs']
     1910    cx,ct,cs,cia = getAtomPtrs(data)
    18851911    drawingData = data['Drawing']
    18861912    modul = generalData['SuperVec'][0]
    1887     dcx,dct,dcs,dci = drawingData['atomPtrs']
     1913    dcx,dct,dcs,dci = getAtomPtrs(data,True)
    18881914    atoms = data['Atoms']
    18891915    drawAtoms = drawingData['Atoms']
    … …  
    29402966    General = Phase['General']
    29412967    Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
    2942     cx,ct,cs,cia = General['AtomPtrs']
     2968    cx,ct,cs,cia = getAtomPtrs(Phase)
    29432969    Atoms = Phase['Atoms']
    29442970    cartAtoms = []
    … …  
    41444170            Unique.append(icntr)
    41454171    return Unique
     4172
     4173def AtomsCollect(data,Ind,Sel):
     4174    '''Finds the symmetry set of atoms for those selected. Selects
     4175    the one closest to the selected part of the unit cell.
     4176    Works on the contents of data['Map Peaks']. Called from OnPeaksUnique in
     4177    GSASIIphsGUI.py,
     4178
     4179    :param data: the phase data structure
     4180    :param list Ind: list of selected peak indices
     4181    :param int Sel: selected part of unit to find atoms closest to
     4182
     4183    :returns: the list of symmetry unique peaks from among those given in Ind
     4184    '''       
     4185    cx,ct,cs,ci = getAtomPtrs(data)
     4186    cent = np.ones(3)*.5     
     4187    generalData = data['General']
     4188    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     4189    SGData = generalData['SGData']
     4190    Atoms = copy.deepcopy(data['Atoms'])
     4191    Indx = [True for ind in Ind]
     4192    # scan through peaks, finding all peaks equivalent to peak ind
     4193    for ind in Ind:
     4194        if Indx[ind]:
     4195            xyz = Atoms[ind][cx:cx+3]
     4196            uij = Atoms[ind][ci+2:ci+8]
     4197            if Atoms[ind][ci] == 'A':
     4198                Equiv = list(G2spc.GenAtom(xyz,SGData,Uij=uij))
     4199                Uijs = np.array([x[1] for x in Equiv])
     4200            else:
     4201                Equiv = G2spc.GenAtom(xyz,SGData)
     4202            xyzs = np.array([x[0] for x in Equiv])
     4203            dzeros = np.sqrt(np.sum(np.inner(Amat,xyzs)**2,axis=0))
     4204            dcent = np.sqrt(np.sum(np.inner(Amat,xyzs-cent)**2,axis=0))
     4205            xyzs = np.hstack((xyzs,dzeros[:,nxs],dcent[:,nxs]))
     4206            cind = np.argmin(xyzs.T[Sel-1])
     4207            Atoms[ind][cx:cx+3] = xyzs[cind][:3]
     4208            if Atoms[ind][ci] == 'A':
     4209                Atoms[ind][ci+2:ci+8] = Uijs[cind]
     4210    return Atoms
    41464211
    41474212################################################################################
    … …  
    50585123    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
    50595124    fixAtoms = data['Atoms']                       #if any
    5060     cx,ct,cs = generalData['AtomPtrs'][:3]
     5125    cx,ct,cs = getAtomPtrs(data)[:3]
    50615126    aTypes = set([])
    50625127    parmDict = {'Bmat':Bmat,'Gmat':Gmat}

  • trunk/GSASIIstrMain.py

    r4761 r4780  
    5555    #report on SVD 0's and highly correlated variables
    5656    msg = ''
    57     SVD0 = result[2].get('SVD0')
    58     if SVD0 > 0:
    59         msg += 'Warning: There were {} singularities in the Hessian'.format(SVD0)
    60     # process singular variables
     57    # process singular variables; all vars go to console, first 10 to
     58    # dialog window
    6159    psing = result[2].get('psing',[])
    6260    if len(psing):
    6361        if msg: msg += '\n'
    64         m = '{} Parameters dropped due to singularities:'.format(len(psing))
     62        m = 'Error: {} Parameter(s) dropped:'.format(len(psing))
    6563        msg += m
    6664        G2fil.G2Print(m, mode='warn')
    6765        m = ''
    6866        for i,val in enumerate(psing):
     67            if i == 0:
     68                msg += '\n{}'.format(varyList[val])
     69                m = '  {}'.format(varyList[val])
     70            else:
     71                if len(m) > 70:
     72                    G2fil.G2Print(m, mode='warn')
     73                    m = '  '
     74                else:
     75                    m += ', '
     76                m += '{}'.format(varyList[val])
     77                if i == 10:
     78                    msg += ', {}... see console for full list'.format(varyList[val])
     79                elif i > 10:
     80                    pass
     81                else:
     82                    msg += ', {}'.format(varyList[val])
     83        if m: G2fil.G2Print(m, mode='warn')
     84    SVD0 = result[2].get('SVD0')
     85    if SVD0 == 1:
     86        msg += 'Warning: Soft (SVD) singularity in the Hessian'
     87    elif SVD0 > 0:
     88        msg += 'Warning: {} soft (SVD) Hessian singularities'.format(SVD0)
     89    SVDsing = result[2].get('SVDsing',[])
     90    if len(SVDsing):
     91        if msg: msg += '\n'
     92        m = 'SVD problem(s) likely from:'
     93        msg += m
     94        G2fil.G2Print(m, mode='warn')
     95        m = ''
     96        for i,val in enumerate(SVDsing):
    6997            if i == 0:
    7098                msg += '\n{}'.format(varyList[val])
    … …  
    209237            if 'Hessian' in Controls['deriv type']:
    210238                num = len(varyList)-1
     239                # BHT -- I am not sure if this works correctly:
    211240                for i,val in enumerate(np.flipud(result[2]['psing'])):
    212241                    if val:
    … …  
    260289                if dlg: break # refining interactively
    261290                # non-interactive refinement
    262                 for val in sorted(psing,reverse=True):
    263                     G2fil.G2Print ('Removing parameter: '+varyList[val])
    264                     del(varyList[val])
    265                 if not psing: break    # removed variable(s), try again
     291                #for val in sorted(psing,reverse=True):
     292                #    G2fil.G2Print ('Removing parameter: '+varyList[val])
     293                #    del(varyList[val])
     294                #if not psing: break    # removed variable(s), try again
    266295            else:
    267296                G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
Note: See TracChangeset for help on using the changeset viewer.