Changeset 2854


Ignore:
Timestamp:
Jun 3, 2017 12:23:04 PM (6 years ago)
Author:
vondreele
Message:

implement Hessian SVD refinement option (no Levenberg-Marquardt). Works for final refinements but not good for initial refinements (far from minimum). Warning in Controls status bar

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r2850 r2854  
    28312831        LSSizer = wx.FlexGridSizer(cols=4,vgap=5,hgap=5)
    28322832        LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Refinement derivatives: '),0,WACV)
    2833         Choice=['analytic Jacobian','numeric','analytic Hessian']   #TODO +'SVD refine' - what flags will it need?
     2833        Choice=['analytic Jacobian','numeric','analytic Hessian','Hessian SVD']   #TODO +'SVD refine' - what flags will it need?
    28342834        derivSel = wx.ComboBox(parent=G2frame.dataDisplay,value=data['deriv type'],choices=Choice,
    28352835            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    … …  
    28472847            maxCyc.Bind(wx.EVT_COMBOBOX, OnMaxCycles)
    28482848            LSSizer.Add(maxCyc,0,WACV)
    2849             LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Initial lambda = 10**'),0,WACV)
    2850             MarqChoice = ['-3','-2','-1','0','1','2','3','4']
    2851             marqLam = wx.ComboBox(parent=G2frame.dataDisplay,value=str(data['Marquardt']),choices=MarqChoice,
    2852                 style=wx.CB_READONLY|wx.CB_DROPDOWN)
    2853             marqLam.Bind(wx.EVT_COMBOBOX,OnMarqLam)
    2854             LSSizer.Add(marqLam,0,WACV)
     2849            if 'SVD' not in data['deriv type']:
     2850                LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Initial lambda = 10**'),0,WACV)
     2851                MarqChoice = ['-3','-2','-1','0','1','2','3','4']
     2852                marqLam = wx.ComboBox(parent=G2frame.dataDisplay,value=str(data['Marquardt']),choices=MarqChoice,
     2853                    style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2854                marqLam.Bind(wx.EVT_COMBOBOX,OnMarqLam)
     2855                LSSizer.Add(marqLam,0,WACV)
    28552856            LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' SVD zero tolerance:'),0,WACV)
    28562857            LSSizer.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'SVDtol',nDig=(10,1,'g'),min=1.e-9,max=.01),0,WACV)
    … …  
    28972898    if not G2frame.dataFrame.GetStatusBar():
    28982899        Status = G2frame.dataFrame.CreateStatusBar()
     2900    else:
     2901        Status = G2frame.dataFrame.GetStatusBar()
     2902    if 'SVD' in data['deriv type']:
     2903        Status.SetStatusText('Hessian SVD not recommended for initial refinements; use analytic Hessian or Jacobian')
     2904    else:
    28992905        Status.SetStatusText('')
    29002906    G2frame.dataFrame.SetLabel('Controls')
    … …  
    46124618                Nvars = len(data['varyList'])
    46134619                Rvals = data['Rvals']
    4614                 text = '\nFinal residuals: \nwR = %.3f%% \nchi**2 = %.1f \nGOF = %.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
    4615                 text += '\nNobs = %d \nNvals = %d'%(Rvals['Nobs'],Nvars)
     4620                text = '\nFinal residuals: \nwR = %.3f%% \nchi**2 = %.1f \nGOF = %.2f ' \
     4621                    %(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
     4622                text += '\nNobs = %d \nNvals = %d \nSVD zeros = %d'%(Rvals['Nobs'],Nvars,Rvals.get('SVD0',0.))
    46164623                if 'lamMax' in Rvals:
    46174624                    text += '\nlog10 MaxLambda = %.1f'%(np.log10(Rvals['lamMax']))

  • trunk/GSASIImath.py

    r2853 r2854  
    146146         * 'lamMax':
    147147         * 'psing':
     148         * 'SVD0':
    148149           
    149150    """
    … …  
    162163    nfev = 0
    163164    if Print:
    164         print ' Hessian SVD refinement on %d variables:'%(n)
     165        print ' Hessian Levenburg-Marquardt SVD refinement on %d variables:'%(n)
    165166    Lam = np.zeros((n,n))
    166167    while icycle < maxcyc:
    … …  
    170171        chisq0 = np.sum(M**2)
    171172        Yvec,Amat = Hess(x0,*args)
    172 #        print 'unscaled SVD zeros',pinv(Amat)[1]
    173173        Adiag = np.sqrt(np.diag(Amat))
    174174        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
    175175        if np.any(psing):                #hard singularity in matrix
    176             return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
     176            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
    177177        Anorm = np.outer(Adiag,Adiag)
    178178        Yvec /= Adiag
    … …  
    185185            try:
    186186                Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
    187 #                Ainv = nl.inv(Amatlam); Nzeros = 0
    188187            except nl.LinAlgError:
    189188                print 'ouch #1 bad SVD inversion; change parameterization'
    190189                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
    191                 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
     190                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
    192191            Xvec = np.inner(Ainv,Yvec)      #solve
    193192            Xvec /= Adiag
    … …  
    228227#        Bmat = nl.inv(Amatlam); Nzeros = 0
    229228        Bmat = Bmat/Anorm
    230         return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
     229        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
    231230    except nl.LinAlgError:
    232231        print 'ouch #2 linear algebra error in making v-cov matrix'
    … …  
    234233        if maxcyc:
    235234            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    236         return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]         
     235        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]         
     236           
     237def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
     238   
     239    """
     240    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
     241    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`   
     242    where :math:`x = arg min(\sum_{y=0}^{N_{obs}} (func(y)^2,axis=0))`
     243
     244    :param function func: callable method or function
     245        should take at least one (possibly length N vector) argument and
     246        returns M floating point numbers.
     247    :param np.ndarray x0: The starting estimate for the minimization of length N
     248    :param function Hess: callable method or function
     249        A required function or method to compute the weighted vector and Hessian for func.
     250        It must be a symmetric NxN array
     251    :param tuple args: Any extra arguments to func are placed in this tuple.
     252    :param float ftol: Relative error desired in the sum of squares.
     253    :param float xtol: Relative tolerance of zeros in the SVD solution in nl.pinv.
     254    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
     255        until other limits are met (ftol, xtol)
     256    :param bool Print: True for printing results (residuals & times) by cycle
     257
     258    :returns: (x,cov_x,infodict) where
     259
     260      * x : ndarray
     261        The solution (or the result of the last iteration for an unsuccessful
     262        call).
     263      * cov_x : ndarray
     264        Uses the fjac and ipvt optional outputs to construct an
     265        estimate of the jacobian around the solution.  ``None`` if a
     266        singular matrix encountered (indicates very flat curvature in
     267        some direction).  This matrix must be multiplied by the
     268        residual standard deviation to get the covariance of the
     269        parameter estimates -- see curve_fit.
     270      * infodict : dict
     271        a dictionary of optional outputs with the keys:
     272
     273         * 'fvec' : the function evaluated at the output
     274         * 'num cyc':
     275         * 'nfev':
     276         * 'lamMax':0.
     277         * 'psing':
     278         * 'SVD0':
     279           
     280    """
     281               
     282    ifConverged = False
     283    deltaChi2 = -10.
     284    x0 = np.array(x0, ndmin=1)      #might be redundant?
     285    n = len(x0)
     286    if type(args) != type(()):
     287        args = (args,)
     288       
     289    icycle = 0
     290    nfev = 0
     291    if Print:
     292        print ' Hessian SVD refinement on %d variables:'%(n)
     293    while icycle < maxcyc:
     294        time0 = time.time()
     295        M = func(x0,*args)
     296        nfev += 1
     297        chisq0 = np.sum(M**2)
     298        Yvec,Amat = Hess(x0,*args)
     299        Adiag = np.sqrt(np.diag(Amat))
     300        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
     301        if np.any(psing):                #hard singularity in matrix
     302            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
     303        Anorm = np.outer(Adiag,Adiag)
     304        Yvec /= Adiag
     305        Amat /= Anorm
     306        if Print:
     307            print 'initial chi^2 %.5g'%(chisq0)
     308        try:
     309            Ainv,Nzeros = pinv(Amat,xtol)    #do Moore-Penrose inversion (via SVD)
     310        except nl.LinAlgError:
     311            print 'ouch #1 bad SVD inversion; change parameterization'
     312            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
     313            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
     314        Xvec = np.inner(Ainv,Yvec)      #solve
     315        Xvec /= Adiag
     316        M2 = func(x0+Xvec,*args)
     317        nfev += 1
     318        chisq1 = np.sum(M2**2)
     319        deltaChi2 = (chisq0-chisq1)/chisq0
     320        if Print:
     321            print ' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Delta: %.3g'%(
     322                icycle,time.time()-time0,chisq1,deltaChi2)
     323        if deltaChi2 < ftol:
     324            ifConverged = True
     325            if Print: print "converged"
     326            break
     327        icycle += 1
     328    M = func(x0,*args)
     329    nfev += 1
     330    Yvec,Amat = Hess(x0,*args)
     331    Adiag = np.sqrt(np.diag(Amat))
     332    Anorm = np.outer(Adiag,Adiag)
     333    Amat = Amat/Anorm       
     334    try:
     335        Bmat,Nzero = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
     336        print 'Found %d SVD zeros'%(Nzero)
     337#        Bmat = nl.inv(Amatlam); Nzeros = 0
     338        Bmat = Bmat/Anorm
     339        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':[],
     340            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
     341    except nl.LinAlgError:
     342        print 'ouch #2 linear algebra error in making v-cov matrix'
     343        psing = []
     344        if maxcyc:
     345            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
     346        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]         
    237347           
    238348def getVCov(varyNames,varyList,covMatrix):

  • trunk/GSASIIstrMain.py

    r2847 r2854  
    6666                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
    6767            ncyc = int(result[2]['nfev']/2)
    68         elif 'Hessian' in Controls['deriv type']:
     68        elif 'analytic Hessian' in Controls['deriv type']:
    6969            Lamda = Controls.get('Marquardt',-3)
    7070            maxCyc = Controls['max cyc']
    … …  
    7474            Rvals['lamMax'] = result[2]['lamMax']
    7575            Controls['Marquardt'] = -3  #reset to default
     76        elif 'Hessian SVD' in Controls['deriv type']:
     77            maxCyc = Controls['max cyc']
     78            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
     79                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
     80            ncyc = result[2]['num cyc']+1
    7681        else:           #'numeric'
    7782            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
    … …  
    8388#        for item in table: print item,table[item]               #useful debug - are things shifting?
    8489        runtime = time.time()-begin
     90        Rvals['SVD0'] = result[2].get('SVD0',0)
    8591        Rvals['converged'] = result[2].get('Converged')
    8692        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
    … …  
    105111#            for item in table: print item,table[item]               #useful debug - are things shifting?
    106112            break                   #refinement succeeded - finish up!
    107         except TypeError:          #result[1] is None on singular matrix
     113        except TypeError:          #result[1] is None on singular matrix or LinAlgError
    108114            IfOK = False
    109115            if not len(varyList):
Note: See TracChangeset for help on using the changeset viewer.