Changeset 4693 for trunk


Ignore:
Timestamp:
Dec 31, 2020 4:11:48 PM (2 years ago)
Author:
toby
Message:

updates to math_new after testing; improve output

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r4690 r4693  
    51805180                text += '\nMax shift/sigma={:.3f}\n'.format(Rvals['Max shft/sig'])
    51815181            if 'msg' in Rvals: text += '\n' + Rvals['msg'] + '\n'
    5182             text += '\nLoad new result?'
    51835182            if lamMax >= 10.:
    51845183                text += '\nWARNING: Steepest descents dominates;'+   \
    51855184                ' minimum may not have been reached or result may be false minimum.'+  \
    5186                 ' You should reconsider which parameters you refine'
     5185                ' You should reconsider which parameters you refine. Check covariance matrix.\n'
     5186            text += '\nLoad new result?'
    51875187            dlg2 = wx.MessageDialog(self,text,'Refinement results, Rw =%.3f'%(Rw),wx.OK|wx.CANCEL)
    51885188            dlg2.CenterOnParent()

  • trunk/GSASIImath_new.py

    r4690 r4693  
    106106    return res,nzero
    107107
    108 def dropTerms(hessian, xvector, yvector, indices, bad):
     108def dropTerms(bad, hessian, indices, *vectors):
    109109    '''Remove the 'bad' terms from the Hessian and vector
    110110
     111    :param tuple bad: a list of variable (row/column) numbers that should be
     112      removed from the hessian and vector. Example: (0,3) removes the 1st and
     113      4th column/row
    111114    :param np.array hessian: a square matrix of length n x n
    112     :param np.array xvector: the least-squares model values, length n
    113     :param np.array yvector: the least-square vector of length n
    114115    :param np.array indices: the indices of the least-squares vector of length n
    115116       referenced to the initial variable list; as this routine is called
    116117       multiple times, more terms may be removed from this list
    117     :param tuple bad: a list of variable (row/column) numbers that should be
    118       removed from the hessian and vector. Example: (0,3) removes the 1st and
    119       4th column/row
    120     :returns: hessian, xvector, yvector, indices where the lengths are
     118    :param additional-args: various least-squares model values, length n
     119
     120    :returns: hessian, indices, vector0, vector1,...  where the lengths are
    121121      now n' x n' and n', with n' = n - len(bad)
    122122    '''
    123     return (np.delete(np.delete(hessian,bad,1),bad,0),
    124                 np.delete(xvector,bad),
    125                 np.delete(yvector,bad),
    126                 np.delete(indices,bad))
    127 
     123    out = [np.delete(np.delete(hessian,bad,1),bad,0),np.delete(indices,bad)]
     124    for v in vectors:
     125        out.append(np.delete(v,bad))
     126    return out
     127
     128       
     129       
    128130def setHcorr(info,Amat,xtol,problem=False):
    129131    Adiag = np.sqrt(np.diag(Amat))
    … …  
    141143        m = max(0.95,0.99*np.amax(AcovUp))
    142144    info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))]
    143     return Bmat
     145    return Bmat,Nzeros
    144146   
    145147def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
    … …  
    194196       
    195197    icycle = 0
    196     lam = 0  #  Could start w/o Marquardt and add it in if needed
    197     lam = 10.**lamda
    198     lamMax = lam
     198    lamMax = lam = 0  #  start w/o Marquardt and add it in if needed
    199199    nfev = 0
    200200    if Print:
    … …  
    215215    if n == 0:
    216216        info = {'num cyc':0,'fvec':M2,'nfev':1,'lamMax':0,'psing':[],'SVD0':0}
    217         info['msg'] = 'no variables, no refinement\n'
     217        info['msg'] = 'no variables: skipping refinement\n'
    218218        return [x0,None,info]
    219219    while icycle < maxcyc:
    … …  
    226226        Xvec = copy.copy(XvecAll)
    227227        # we could remove vars that were dropped in previous cycles here (use indices), but for
    228         # now let's reset them each cycle in case the singularities go awat
     228        # now let's reset them each cycle in case the singularities go away
    229229        indices = range(n)
    230230        Adiag = np.sqrt(np.diag(Amat))
    231         psing = np.where(np.abs(Adiag) < 1.e-14)[0]         # remove any hard singularities
     231        psing = np.where(np.abs(Adiag) < 1.e-14)[0]       # find any hard singularities
    232232        if len(psing):
    233             G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format(psing), mode='warn')
    234             Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
    235             Adiag = np.sqrt(np.diag(Amat))
     233            G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format(
     234                psing), mode='warn')
     235            Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,
     236                                        Amat, indices, Xvec, Yvec, Adiag)
    236237        Anorm = np.outer(Adiag,Adiag)        # normalize matrix & vector
    237238        Yvec /= Adiag
    … …  
    240241        maxdrop = 3 # max number of drop var attempts
    241242        loops = 0
    242         while True: #---------------------- loop as we increase lamda and possibly drop terms from hessian, x & y vectors
     243        while True: #--------- loop as we increase lamda and possibly drop terms
     244            lamMax = max(lamMax,lam)
    243245            Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
    244246            try:
    245                 Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
     247                Ainv,Nzeros = pinv(Amatlam,xtol)    # Moore-Penrose SVD inversion
    246248            except nl.LinAlgError:
    247249                loops += 1
    … …  
    252254                G2fil.G2Print('ouch #2 bad SVD inversion; dropping terms for for variable(s) #{}'.
    253255                                  format(psing), mode='warn')
    254                 Amatlam, Xvec, Yvec, indices = dropTerms(Amatlam, Xvec, Yvec, indices, psing)
    255                 if loops < maxdrop: continue
     256                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,
     257                                        Amat, indices, Xvec, Yvec, Adiag)
     258                if loops < maxdrop: continue # try again, same lam but fewer vars
    256259                G2fil.G2Print('giving up with ouch #2', mode='error')
    257                 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':Nzeros}
     260                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
     261                info['psing'] = [i for i in range(n) if i not in indices]
     262
    258263                info['msg'] = 'SVD inversion failure\n'
    259264                return [x0,None,info]
    … …  
    265270                G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.
    266271                                  format(Nzeros,psing), mode='warn')
    267                 Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
     272                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,
     273                                        Amat, indices, Xvec, Yvec, Adiag)
    268274                Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
    269                 Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
    270             Adiag = np.sqrt(np.diag(Amat))
     275                Ainv,nz = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
     276                if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop',
     277                                             mode='warn')
    271278            Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
    272279            XvecAll[indices] = Xvec         # expand
    … …  
    275282            except Exception as Msg:
    276283                if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     284                G2fil.G2Print(Msg.msg,mode='warn')
    277285                loops += 1
    278286                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
    279                 psing = list(np.where(d < 1.e-14)[0])
    280                 if len(psing) < Nzeros:
    281                     psing = list(np.where(d <= sorted(d)[max(Nzeros,1)-1])[0])
    282                 G2fil.G2Print(Msg.msg,mode='warn')
    283                 # dropping variables does not seem to help
    284                 if len(psing):
    285                     G2fil.G2Print('ouch #3 unable to evaluate objective function;'+
    286                                     '\n\tdropping terms for variable(s) #{}'.format(psing),
    287                                       mode='warn')
    288                     Amatlam, Xvec, Yvec, indices = dropTerms(Amatlam, Xvec, Yvec, indices, psing)
    289                     if loops < maxdrop: continue
    290                     G2fil.G2Print('giving up with ouch #3', mode='error')
    291                 else: #nothing to improve
    292                     G2fil.G2Print('ouch #3 unable to evaluate objective function. Giving up',
    293                                   mode='warn')
     287                G2fil.G2Print('ouch #3 unable to evaluate objective function;',
     288                                      mode='error')
    294289                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
    295290                info['psing'] = [i for i in range(n) if i not in indices]
    296                 try:         # report highly correlated parameters from full Hessian, if we can
     291                try:         # try to report highly correlated parameters from full Hessian
    297292                    setHcorr(info,AmatAll,xtol,problem=True)
    298293                except nl.LinAlgError:
    … …  
    304299                    G2fil.G2Print('Error determining highly correlated vars',mode='warn')
    305300                    G2fil.G2Print(Msg.msg,mode='warn')
    306                 info['msg'] = Msg.msg + '\n'
     301                info['msg'] = Msg.msg + '\n\n'
    307302                return [x0,None,info]
    308303            nfev += 1
    309304            chisq1 = np.sum(M2**2)
    310305            if chisq1 > chisq0*(1.+chitol):     #TODO put Alan Coehlo's criteria for lambda here?
    311                 lam *= 10.
    312                 if lam == 0: lam = 10.**lamda  #  this is here if we start w/o Marquardt term
     306                if lam == 0:
     307                    lam = 10.**lamda  #  set to initial Marquardt term
     308                else:
     309                    lam *= 10.
    313310                if Print:
    314                     G2fil.G2Print(('divergence: chi^2 %.5g on %d obs. w/%d SVD zeros\n'+
     311                    G2fil.G2Print(('divergence: chi^2 %.5g on %d obs. (%d SVD zeros)\n'+
    315312                                    '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam))
    316313                if lam > 10.:
    … …  
    318315                                  (chisq1,chisq0,lam), mode='warn')
    319316                    try:         # report highly correlated parameters from full Hessian, if we can
    320                         info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
     317                        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,
    321318                         'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}
    322319                        info['psing'] = [i for i in range(n) if i not in indices]
    323                         Bmat = setHcorr(info,AmatAll,xtol,problem=True)
     320                        Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True)
     321                        info['SVD0'] = Nzeros
    324322                        return [x0,Bmat,info]
    325323                    except Exception as Msg:
    … …  
    335333                break
    336334        # complete current cycle
    337         lamMax = max(lamMax,lam)
    338335        deltaChi2 = (chisq0-chisq1)/chisq0
    339336        if Print and n-len(indices):
    … …  
    366363    chisqf = np.sum(M**2) # ending chi**2
    367364    Yvec,Amat = Hess(x0,*args)    # we could save some time and use the last Hessian from the last refinement cycle
    368     psing_prev = [i for i in range(n) if i not in indices]
     365    psing_prev = [i for i in range(n) if i not in indices] # save dropped vars
    369366    indices = range(n)
    370367    info = {}
    371368    try:         # report highly correlated parameters from full Hessian, if we can
    372         Bmat = setHcorr(info,Amat,xtol,problem=False)
     369        Bmat,Nzeros = setHcorr(info,Amat,xtol,problem=False)
    373370        info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,'psing':psing_prev,
    374371                         'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
    … …  
    382379        G2fil.G2Print('Error determining highly correlated vars',mode='warn')
    383380        G2fil.G2Print(Msg.msg,mode='warn')
    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)
     381    # matrix above inversion failed, drop previously removed variables & try again
     382    Amat, indices, Yvec = dropTerms(psing_prev, Amat, indices, Yvec)
    387383    Adiag = np.sqrt(np.diag(Amat))
    388384    Anorm = np.outer(Adiag,Adiag)
    … …  
    394390        G2fil.G2Print('ouch #5 linear algebra error in making final v-cov matrix', mode='error')
    395391        psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e-14)[0])
    396         if not len(psing): # make sure at least the worst is returned
     392        if not len(psing): # make sure at least the worst term is flagged
    397393            psing = [np.argmin(d)]
    398         Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
     394        Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
    399395        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':-1,'Xvec':None, 'chisq0':chisqf}
    400396        info['psing'] = [i for i in range(n) if i not in indices]
    401397        return [x0,None,info]
    402     if Nzeros: # this is unexpected, since previous refinements should have removed the zeros
     398    if Nzeros: # sometimes SVD zeros show up when lam=0
    403399        d = np.abs(np.diag(nl.qr(Amat)[1]))
    404400        psing = list(np.where(d < 1.e-14)[0])
    … …  
    406402            psing = list(np.where(d <= sorted(d)[Nzeros-1])[0])
    407403        if Print:
    408             G2fil.G2Print('Found %d SVD zeros w/o Lambda. Change problem with variable(s) #%s'%
     404            G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+
     405                              'Likely problem with variable(s) #%s'%
    409406                              (Nzeros,psing), mode='warn')
    410         Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing)
    411     # expand Bmat by filling with zeros
     407        Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
     408    # expand Bmat by filling with zeros if columns have been dropped
     409    psing = [i for i in range(n) if i not in indices]
    412410    if len(psing):
    413         Bmat = np.insert(np.insert(Bmat,psing,0,1),psing,0,0)
     411        ins = [j-i for i,j in enumerate(psing)]
     412        Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0)
    414413    info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
    415414                         'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})

  • trunk/GSASIIstrMain.py

    r4691 r4693  
    6060    if len(psing):
    6161        if msg: msg += '\n'
    62         msg += 'Parameters dropped due to singularities:'
    63     for i,val in enumerate(psing):
    64         if i == 0:
    65             msg += '\n\t{}'.format(varyList[val])
    66         else:
    67             msg += ', {}'.format(varyList[val])
     62        msg += '{} Parameters dropped due to singularities:'.format(len(psing))
     63        for i,val in enumerate(psing):
     64            if i == 0:
     65                msg += '\n\t{}'.format(varyList[val])
     66            else:
     67                msg += ', {}'.format(varyList[val])
     68        G2fil.G2Print(msg, mode='warn')
    6869    #report on highly correlated variables
    6970    Hcorr = result[2].get('Hcorr',[])
    70     for i,(v1,v2,corr) in enumerate(Hcorr):
     71    if len(Hcorr) > 0:
    7172        if msg: msg += '\n'
    72         if i == 0:
    73             msg += 'Note highly correlated parameters:\n'
     73        msg += 'Note highly correlated parameters:'
     74    elif SVD0 > 0:
     75        if msg: msg += '\n'
     76        msg += 'Check covariance matrix for parameter correlation'
     77    for v1,v2,corr in Hcorr:
    7478        if corr > .95:
    7579            stars = '**'
    7680        else:
    7781            stars = '   '
    78         msg += ' {} {} and {} (@{:.2f}%)'.format(
     82        msg += '\n {} {} and {} (@{:.2f}%)'.format(
    7983            stars,varyList[v1],varyList[v2],100.*corr)
    8084    if msg:
Note: See TracChangeset for help on using the changeset viewer.