Changeset 975


Ignore:
Timestamp:
Jun 29, 2013 11:39:57 AM (8 years ago)
Author:
vondreele
Message:

a few mods to HessianLSQ

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r974 r975  
    4242atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    4343
    44 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
     44def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
    4545   
    4646    """
     
    6666    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
    6767        until other limits are met (ftol, xtol)
     68    :param bool Print: True for printing results (residuals & times) by cycle
    6869
    6970    :returns: (x,cov_x,infodict) where
     
    100101    lamMax = lam
    101102    nfev = 0
     103    if Print:
     104        print ' Hessian refinement on %d variables:'%(n)
    102105    while icycle < maxcyc:
    103         lamMax = max(lamMax,lam)
     106        time0 = time.time()
    104107        M = func(x0,*args)
    105108        nfev += 1
     
    112115        Anorm = np.outer(Adiag,Adiag)
    113116        Yvec /= Adiag
    114         Amat /= Anorm       
     117        Amat /= Anorm
    115118        while True:
    116119            Lam = np.eye(Amat.shape[0])*lam
     
    135138                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
    136139                break
     140        lamMax = max(lamMax,lam)
     141        if Print:
     142            print ' Cycle: %d, Time: %.2fs, Chi**2: %.3g, Lambda: %.3g'%(icycle,time.time()-time0,chisq1,lam)
    137143        if (chisq0-chisq1)/chisq0 < ftol:
    138144            break
     
    141147    nfev += 1
    142148    Yvec,Amat = Hess(x0,*args)
     149    Lam = np.eye(Amat.shape[0])*lam
     150    Amatlam = Amat*(One+Lam)
    143151    try:
    144         Bmat = nl.inv(Amat)
     152        Bmat = nl.inv(Amatlam)
    145153        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
    146154    except nl.LinAlgError:
  • trunk/GSASIIstrMain.py

    r960 r975  
    119119            ncyc = int(result[2]['nfev']/2)
    120120        elif 'Hessian' in Controls['deriv type']:
    121             result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,
     121            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,Print=True,
    122122                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
    123123            ncyc = result[2]['num cyc']+1
Note: See TracChangeset for help on using the changeset viewer.