Ignore:
Timestamp:
Jan 10, 2021 5:12:01 PM (10 months ago)
Author:
toby
Message:

wx4.1 & mpl3.3 fixes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath_new.py

    r4699 r4759  
    22#GSASIImath - major mathematics routines
    33########### SVN repository information ###################
    4 # $Date: 2020-12-12 13:30:31 -0600 (Sat, 12 Dec 2020) $
     4# $Date: 2021-01-04 10:41:46 -0600 (Mon, 04 Jan 2021) $
    55# $Author: vondreele $
    6 # $Revision: 4671 $
     6# $Revision: 4709 $
    77# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIImath.py $
    8 # $Id: GSASIImath.py 4671 2020-12-12 19:30:31Z vondreele $
     8# $Id: GSASIImath.py 4709 2021-01-04 16:41:46Z vondreele $
    99########### SVN repository information ###################
    1010'''
     
    2424import copy
    2525import GSASIIpath
    26 GSASIIpath.SetVersionNumber("$Revision: 4671 $")
     26GSASIIpath.SetVersionNumber("$Revision: 4709 $")
    2727import GSASIIElem as G2el
    2828import GSASIIlattice as G2lat
     
    102102    s = np.where(s>cutoff,1./s,0.)
    103103    nzero = s.shape[0]-np.count_nonzero(s)
    104 #    res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))
    105104    res = np.dot(vt.T,s[:,nxs]*u.T)
    106105    return res,nzero
     
    125124        out.append(np.delete(v,bad))
    126125    return out
    127 
    128        
    129126       
    130127def setHcorr(info,Amat,xtol,problem=False):
     
    233230            G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format(
    234231                psing), mode='warn')
    235             Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,
    236                                         Amat, indices, Xvec, Yvec, Adiag)
     232            Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
    237233        Anorm = np.outer(Adiag,Adiag)        # normalize matrix & vector
    238234        Yvec /= Adiag
     
    253249                    psing = [np.argmin(d)]
    254250                G2fil.G2Print('ouch #2 bad SVD inversion; dropping terms for for variable(s) #{}'.
    255                                   format(psing), mode='warn')
    256                 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,
    257                                         Amat, indices, Xvec, Yvec, Adiag)
     251                    format(psing), mode='warn')
     252                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
    258253                if loops < maxdrop: continue # try again, same lam but fewer vars
    259254                G2fil.G2Print('giving up with ouch #2', mode='error')
     
    269264                    psing = [np.argmin(d)]
    270265                G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.
    271                                   format(Nzeros,psing), mode='warn')
    272                 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,
    273                                         Amat, indices, Xvec, Yvec, Adiag)
     266                    format(Nzeros,psing), mode='warn')
     267                Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)
    274268                Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)
    275269                Ainv,nz = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
    276270                if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop'.format(nz),
    277                                              mode='warn')
     271                    mode='warn')
    278272            Xvec = np.inner(Ainv,Yvec)/Adiag      #solve for LS terms
    279273            XvecAll[indices] = Xvec         # expand
     
    285279                loops += 1
    286280                d = np.abs(np.diag(nl.qr(Amatlam)[1]))
    287                 G2fil.G2Print('ouch #3 unable to evaluate objective function;',
    288                                       mode='error')
     281                G2fil.G2Print('ouch #3 unable to evaluate objective function;',mode='error')
    289282                info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros}
    290283                info['psing'] = [i for i in range(n) if i not in indices]
     
    310303                if Print:
    311304                    G2fil.G2Print(('divergence: chi^2 %.5g on %d obs. (%d SVD zeros)\n'+
    312                                     '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam))
     305                        '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam))
    313306                if lam > 10.:
    314307                    G2fil.G2Print('ouch #4 stuck: chisq-new %.4g > chisq0 %.4g with lambda %.1g'%
    315                                   (chisq1,chisq0,lam), mode='warn')
     308                        (chisq1,chisq0,lam), mode='warn')
    316309                    try:         # report highly correlated parameters from full Hessian, if we can
    317310                        info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,
     
    369362        Bmat,Nzeros = setHcorr(info,Amat,xtol,problem=False)
    370363        info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,'psing':psing_prev,
    371                          'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
     364            'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
    372365        return [x0,Bmat,info]
    373366    except nl.LinAlgError:
     
    403396        if Print:
    404397            G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+
    405                               'Likely problem with variable(s) #%s'%
    406                               (Nzeros,psing), mode='warn')
     398                'Likely problem with variable(s) #%s'%(Nzeros,psing), mode='warn')
    407399        Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)
    408400    # expand Bmat by filling with zeros if columns have been dropped
     
    412404        Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0)
    413405    info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
    414                          'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
     406        'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00})
    415407    info['psing'] = [i for i in range(n) if i not in indices]
    416408    return [x0,Bmat,info]
     
    492484            Ainv,Nzeros = pinv(Amat,xtol)    #do Moore-Penrose inversion (via SVD)
    493485        except nl.LinAlgError:
    494             G2fil.G2Print('ouch #1 bad SVD inversion; change parameterization', mode='error')
     486            G2fil.G2Print('ouch #1 bad SVD inversion; change parameterization', mode='warn')
    495487            psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e-14)[0])
    496488            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
     
    518510    Amat = Amat/Anorm       
    519511    try:
    520         Bmat,Nzeros = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
    521         G2fil.G2Print('Found %d SVD zeros'%(Nzeros), mode='warn')
     512        Bmat,Nzero = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
     513        G2fil.G2Print('Found %d SVD zeros'%(Nzero), mode='warn')
     514#        Bmat = nl.inv(Amatlam); Nzeros = 0
    522515        Bmat = Bmat/Anorm
    523516        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':[],
    524             'SVD0':Nzeros,'Converged': ifConverged, 'DelChi2':deltaChi2,
     517            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,
    525518                             'chisq0':chisq00}]
    526519    except nl.LinAlgError:
     
    528521        psing = []
    529522        if maxcyc:
    530             psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e-14)[0])
     523            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    531524        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1,
    532525                             'chisq0':chisq00}]
Note: See TracChangeset for help on using the changeset viewer.