Changeset 2846


Ignore:
Timestamp:
May 30, 2017 4:40:30 PM (4 years ago)
Author:
vondreele
Message:

implement SVD inversion (Moore Penrose invert routine) in Hessian LSQ
clean up commented out code replaced by ValidatedTextCtrl?
set default dM/M to 0.001 (from 0.0001) - more reasonable convergence criterion

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r2842 r2846  
    421421    def Draw(self):
    422422       
    423 #        def OnRadius(event):
    424 #            event.Skip()
    425 #            try:
    426 #                val = float(radius.GetValue())
    427 #                if val < 0.5:
    428 #                    raise ValueError
    429 #                self.Sphere = val
    430 #            except ValueError:
    431 #                pass
    432 #            radius.SetValue('%.3f'%(self.Sphere))
    433 #           
    434423        def OnAtomType(event):
    435424            Obj = event.GetEventObject()
    … …  
    446435            topSizer.Add(wx.StaticText(self.panel,label=' Sphere centered at atoms: '),0,WACV)
    447436            cx,ct,cs = self.Drawing['atomPtrs'][:3]
    448 #            print self.Drawing.keys()
    449437            for id in self.indx:
    450438                atom = self.Drawing['Atoms'][id]
    … …  
    459447        sphereSizer = wx.BoxSizer(wx.HORIZONTAL)
    460448        sphereSizer.Add(wx.StaticText(self.panel,label=' Sphere radius: '),0,WACV)
    461         radius = G2G.ValidatedTxtCtrl(self.panel,self.Sphere,nDig=(10,3),typeHint=float,size=(65,25))
    462 #        radius = wx.TextCtrl(self.panel,value='%.3f'%(self.Sphere),style=wx.TE_PROCESS_ENTER)
    463 #        radius.Bind(wx.EVT_TEXT_ENTER,OnRadius)
    464 #        radius.Bind(wx.EVT_KILL_FOCUS,OnRadius)
     449        radius = G2G.ValidatedTxtCtrl(self.panel,self.Sphere,nDig=(10,3),size=(65,25))
    465450        sphereSizer.Add(radius,0,WACV)
    466451        mainSizer.Add(sphereSizer,0,WACV)
    … …  
    617602        self.panel.Destroy()
    618603        self.panel = wx.Panel(self)
    619 #        Ind = {}
    620604        mainSizer = wx.BoxSizer(wx.VERTICAL)
    621605        MatSizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    639623        for iy,line in enumerate(self.Trans):
    640624            for ix,val in enumerate(line):
    641                 item = G2G.ValidatedTxtCtrl(self.panel,self.Trans[iy],ix,nDig=(10,3),typeHint=float,size=(65,25))
     625                item = G2G.ValidatedTxtCtrl(self.panel,self.Trans[iy],ix,nDig=(10,3),size=(65,25))
    642626                Trmat.Add(item)
    643627            Trmat.Add((25,0),0)
    644             vec = G2G.ValidatedTxtCtrl(self.panel,self.Vec,iy,nDig=(10,3),typeHint=float,size=(65,25))
     628            vec = G2G.ValidatedTxtCtrl(self.panel,self.Vec,iy,nDig=(10,3),size=(65,25))
    645629            Trmat.Add(vec)
    646630        transSizer.Add(Trmat)
    … …  
    811795    def Draw(self):
    812796
    813 #        def OnMatValue(event):
    814 #            event.Skip()
    815 #            Obj = event.GetEventObject()
    816 #            ix,iy = Ind[Obj.GetId()]
    817 #            val = Obj.GetValue()
    818 #            if '/' in val:
    819 #                vals = val.split('/')
    820 #                self.Trans[iy,ix] = float(vals[0])/float(vals[1])
    821 #            else:   
    822 #                self.Trans[iy,ix] = float(Obj.GetValue())
    823 #            Obj.SetValue('%5.3f'%(self.Trans[iy,ix]))
    824 #           
    825 #           
    826 #        def OnVecValue(event):
    827 #            event.Skip()
    828 #            Obj = event.GetEventObject()
    829 #            iy = Ind[Obj.GetId()]
    830 #            val = Obj.GetValue()
    831 #            if '/' in val:
    832 #                vals = val.split('/')
    833 #                self.Vec[iy] = float(vals[0])/float(vals[1])
    834 #            else:   
    835 #                self.Vec[iy] = float(Obj.GetValue())
    836 #            Obj.SetValue('%5.3f'%(self.Vec[iy]))
    837 #           
    838797        def OnExpand(event):
    839798            self.Expand = expand.GetValue()
    … …  
    859818        self.panel.Destroy()
    860819        self.panel = wx.Panel(self)
    861 #        Ind = {}
    862820        mainSizer = wx.BoxSizer(wx.VERTICAL)
    863821        MatSizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    868826        for iy,line in enumerate(self.Trans):
    869827            for ix,val in enumerate(line):
    870                 item = G2G.ValidatedTxtCtrl(self.panel,self.Trans[iy],ix,nDig=(10,3),typeHint=float,size=(65,25))
    871 #                item = wx.TextCtrl(self.panel,value='%5.3f'%(val),
    872 #                    size=(50,25),style=wx.TE_PROCESS_ENTER)
    873 #                Ind[item.GetId()] = [ix,iy]
    874 #                item.Bind(wx.EVT_TEXT_ENTER,OnMatValue)
    875 #                item.Bind(wx.EVT_KILL_FOCUS,OnMatValue)
     828                item = G2G.ValidatedTxtCtrl(self.panel,self.Trans[iy],ix,nDig=(10,3),size=(65,25))
    876829                Trmat.Add(item)
    877830            Trmat.Add((25,0),0)
    878             vec = G2G.ValidatedTxtCtrl(self.panel,self.Vec,iy,nDig=(10,3),typeHint=float,size=(65,25))
    879 #            vec = wx.TextCtrl(self.panel,value='%5.3f'%(self.Vec[iy]),
    880 #                    size=(50,25),style=wx.TE_PROCESS_ENTER)
    881 #            Ind[vec.GetId()] = [iy]       
    882 #            vec.Bind(wx.EVT_TEXT_ENTER,OnVecValue)
    883 #            vec.Bind(wx.EVT_KILL_FOCUS,OnVecValue)
     831            vec = G2G.ValidatedTxtCtrl(self.panel,self.Vec,iy,nDig=(10,3),size=(65,25))
    884832            Trmat.Add(vec)
    885833        transSizer.Add(Trmat)
    … …  
    11121060    def Draw(self):
    11131061               
    1114 #        def OnMatValue(event):
    1115 #            event.Skip()
    1116 #            Obj = event.GetEventObject()
    1117 #            ix,iy = Ind[Obj.GetId()]
    1118 #            self.Trans[ix,iy] = float(Obj.GetValue())
    1119 #               
    11201062        def OnCent(event):
    11211063            Obj = event.GetEventObject()
    … …  
    11431085        self.panel.Destroy()
    11441086        self.panel = wx.Panel(self)
    1145 #        Ind = {}
    11461087        mainSizer = wx.BoxSizer(wx.VERTICAL)
    11471088        MatSizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    11611102        for iy,line in enumerate(self.Trans):
    11621103            for ix,val in enumerate(line):
    1163                 item = G2G.ValidatedTxtCtrl(self.panel,self.Trans[iy],ix,nDig=(10,3),typeHint=float,size=(65,25))
    1164 #                item = wx.TextCtrl(self.panel,value='%5.3f'%(val),
    1165 #                    size=(50,25),style=wx.TE_PROCESS_ENTER)
    1166 #                Ind[item.GetId()] = [ix,iy]
    1167 #                item.Bind(wx.EVT_TEXT_ENTER,OnMatValue)
    1168 #                item.Bind(wx.EVT_KILL_FOCUS,OnMatValue)
     1104                item = G2G.ValidatedTxtCtrl(self.panel,self.Trans[iy],ix,nDig=(10,3),size=(65,25))
    11691105                Trmat.Add(item)
    11701106        transSizer.Add(Trmat)
    … …  
    14111347        for id,item in enumerate(self.data['AtomTypes']):
    14121348            radiiSizer.Add(wx.StaticText(self.panel,-1,' '+item),0,WACV)
    1413             bRadii = G2G.ValidatedTxtCtrl(self.panel,data['BondRadii'],id,nDig=(10,3),typeHint=float)
     1349            bRadii = G2G.ValidatedTxtCtrl(self.panel,data['BondRadii'],id,nDig=(10,3))
    14141350            radiiSizer.Add(bRadii,0,WACV)
    14151351            if self.Angle:
    1416                 aRadii = G2G.ValidatedTxtCtrl(self.panel,data['AngleRadii'],id,nDig=(10,3),typeHint=float)
     1352                aRadii = G2G.ValidatedTxtCtrl(self.panel,data['AngleRadii'],id,nDig=(10,3))
    14171353                radiiSizer.Add(aRadii,0,WACV)
    14181354        mainSizer.Add(radiiSizer,0,wx.EXPAND)
    … …  
    14221358            for i,name in enumerate(Names):
    14231359                factorSizer.Add(wx.StaticText(self.panel,-1,name+' search factor'),0,WACV)
    1424                 bondFact = G2G.ValidatedTxtCtrl(self.panel,data['Factors'],i,nDig=(10,3),typeHint=float)
     1360                bondFact = G2G.ValidatedTxtCtrl(self.panel,data['Factors'],i,nDig=(10,3))
    14251361                factorSizer.Add(bondFact)
    14261362            mainSizer.Add(factorSizer,0,wx.EXPAND)
    … …  
    27762712        data = {}
    27772713        data['deriv type'] = 'analytic Hessian'
    2778         data['min dM/M'] = 0.0001
     2714        data['min dM/M'] = 0.001
    27792715        data['shift factor'] = 1.
    27802716        data['max cyc'] = 3       
    27812717        data['F**2'] = False
     2718    if 'SVDtol' not in data:
     2719        data['SVDtol'] = 1.e-6
    27822720    if 'shift factor' not in data:
    27832721        data['shift factor'] = 1.
    … …  
    29012839        LSSizer.Add(derivSel,0,WACV)
    29022840        LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Min delta-M/M: '),0,WACV)
    2903         Cnvrg = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'min dM/M',nDig=(10,2,'g'),min=1.e-9,max=1.,typeHint=float)
    2904         LSSizer.Add(Cnvrg,0,WACV)
     2841        LSSizer.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'min dM/M',nDig=(10,2,'g'),min=1.e-9,max=1.),0,WACV)
    29052842        if 'Hessian' in data['deriv type']:
    29062843            LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Max cycles: '),0,WACV)
    … …  
    29162853            marqLam.Bind(wx.EVT_COMBOBOX,OnMarqLam)
    29172854            LSSizer.Add(marqLam,0,WACV)
     2855            LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' SVD zero tolerance:'),0,WACV)
     2856            LSSizer.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'SVDtol',nDig=(10,1,'g'),min=1.e-9,max=.01),0,WACV)
    29182857        else:       #TODO what for SVD refine?
    29192858            LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Initial shift factor: '),0,WACV)
    2920             Factr = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'shift factor',nDig=(10,5),min=1.e-5,max=100.,typeHint=float)
     2859            Factr = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'shift factor',nDig=(10,5),min=1.e-5,max=100.)
    29212860            LSSizer.Add(Factr,0,WACV)
    29222861        if G2frame.Sngl:
    … …  
    29342873                LSSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=usrRej[item][0]),0,WACV)
    29352874                usrrej = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,userReject,item,nDig=(10,2),
    2936                     min=usrRej[item][1][0],max=usrRej[item][1][1],typeHint=float)
     2875                    min=usrRej[item][1][0],max=usrRej[item][1][1])
    29372876                LSSizer.Add(usrrej,0,WACV)
    29382877        return LSSizer
    … …  
    43904329        G2plt.PlotDeltSig(G2frame,kind)
    43914330       
    4392 #    def OnWtFactor(event):
    4393 #        event.Skip()
    4394 #        try:
    4395 #            val = float(wtval.GetValue())
    4396 #        except ValueError:
    4397 #            val = data[0]['wtFactor']
    4398 #        data[0]['wtFactor'] = val
    4399 #        wtval.SetValue('%.3f'%(val))
    4400 #       
    44014331#    def OnCompression(event):
    44024332#        data[0] = int(comp.GetValue())
    … …  
    44514381    wtSizer = wx.BoxSizer(wx.HORIZONTAL)
    44524382    wtSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' Weight factor: '),0,WACV)
    4453     wtval = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data[0],'wtFactor',nDig=(10,3),min=1.e-9,typeHint=float)
    4454 #    wtval = wx.TextCtrl(G2frame.dataDisplay,-1,'%.3f'%(data[0]['wtFactor']),style=wx.TE_PROCESS_ENTER)
    4455 #    wtval.Bind(wx.EVT_TEXT_ENTER,OnWtFactor)
    4456 #    wtval.Bind(wx.EVT_KILL_FOCUS,OnWtFactor)
     4383    wtval = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data[0],'wtFactor',nDig=(10,3),min=1.e-9)
    44574384    wtSizer.Add(wtval,0,WACV)
    44584385#    if kind == 'PWDR':         #possible future compression feature; NB above patch as well

  • trunk/GSASIImath.py

    r2840 r2846  
    4747################################################################################
    4848
    49 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,lamda=-3,Print=False):
     49def pinv(a, rcond=1e-15 ):
     50    """
     51    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
     52    Modified from numpy.linalg.pinv; assumes a is Hessian & returns no. zeros found
     53    Calculate the generalized inverse of a matrix using its
     54    singular-value decomposition (SVD) and including all
     55    *large* singular values.
     56
     57    Parameters
     58    ----------
     59    a : (M, M) array_like - here assumed to be LS Hessian
     60      Matrix to be pseudo-inverted.
     61    rcond : float
     62      Cutoff for small singular values.
     63      Singular values smaller (in modulus) than
     64      `rcond` * largest_singular_value (again, in modulus)
     65      are set to zero.
     66
     67    Returns
     68    -------
     69    B : (M, M) ndarray
     70      The pseudo-inverse of `a`
     71
     72    Raises
     73    ------
     74    LinAlgError
     75      If the SVD computation does not converge.
     76
     77    Notes
     78    -----
     79    The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
     80    defined as: "the matrix that 'solves' [the least-squares problem]
     81    :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
     82    :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
     83
     84    It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
     85    value decomposition of A, then
     86    :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
     87    orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
     88    of A's so-called singular values, (followed, typically, by
     89    zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
     90    consisting of the reciprocals of A's singular values
     91    (again, followed by zeros). [1]_
     92
     93    References
     94    ----------
     95    .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
     96           FL, Academic Press, Inc., 1980, pp. 139-142.
     97
     98    """
     99    u, s, vt = nl.svd(a, 0)
     100    cutoff = rcond*np.maximum.reduce(s)
     101    s = np.where(s>cutoff,1./s,0.)
     102    nzero = s.shape[0]-np.count_nonzero(s)
     103    res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))
     104    return res,nzero
     105
     106def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
    50107   
    51108    """
    … …  
    63120    :param tuple args: Any extra arguments to func are placed in this tuple.
    64121    :param float ftol: Relative error desired in the sum of squares.
    65     :param float xtol: Relative error desired in the approximate solution.
     122    :param float xtol: Relative tolerance of zeros in the SVD solution in nl.pinv.
    66123    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
    67124        until other limits are met (ftol, xtol)
    … …  
    105162    nfev = 0
    106163    if Print:
    107         print ' Hessian refinement on %d variables:'%(n)
     164        print ' Hessian SVD refinement on %d variables:'%(n)
    108165    Lam = np.zeros((n,n))
    109166    while icycle < maxcyc:
    … …  
    126183            Amatlam = Amat*(One+Lam)
    127184            try:
    128                 Xvec = nl.solve(Amatlam,Yvec)
     185                Ainv = pinv(Amatlam,xtol)[0]    #do Moore-Penrose inversion (via SVD)
     186                Xvec = np.inner(Ainv,Yvec)      #solve
    129187            except nl.LinAlgError:
    130188                print 'ouch #1'
    … …  
    163221    Lam = np.eye(Amat.shape[0])*lam
    164222    Amatlam = Amat/Anorm       
    165 #    Amatlam = Amat*(One+Lam)                #scale Amat to Marquardt array?       
    166223    try:
    167         Bmat = nl.inv(Amatlam)/Anorm
    168 #        Bmat = Bmat*(One+Lam)               #rescale Bmat to Marquardt array?
     224        Bmat,Nzero = pinv(Amatlam,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
     225        print 'Found %d SVD zeros'%(Nzero)
     226        Bmat = Bmat/Anorm
    169227        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
    170228    except nl.LinAlgError:
    … …  
    173231        if maxcyc:
    174232            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    175         return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
    176    
     233        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]         
     234            
    177235def getVCov(varyNames,varyList,covMatrix):
    178236    '''obtain variance-covariance terms for a set of variables. NB: the varyList

  • trunk/GSASIIobj.py

    r2845 r2846  
    874874DefaultControls = {
    875875    'deriv type':'analytic Hessian',
    876     'min dM/M':0.0001,'shift factor':1.,'max cyc':3,'F**2':False,
     876    'min dM/M':0.001,'shift factor':1.,'max cyc':3,'F**2':False,'SVDtol':1.e-6,
    877877    'UsrReject':{'minF/sig':0,'MinExt':0.01,'MaxDF/F':100.,'MaxD':500.,'MinD':0.05},
    878878    'Copy2Next':False,'Reverse Seq':False,'HatomFix':False,

  • trunk/GSASIIpwdGUI.py

    r2820 r2846  
    560560        dlg = wx.ProgressDialog('Sequential peak fit','Data set name = '+names[0],len(names),
    561561            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
    562         controls = {'deriv type':'analytic','min dM/M':0.0001,}
     562        controls = {'deriv type':'analytic','min dM/M':0.001,}
    563563        print 'Peak Fitting with '+controls['deriv type']+' derivatives:'
    564564        oneCycle = False
    … …  
    618618        controls = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.root, 'Controls'))
    619619        if not controls:
    620             controls = {'deriv type':'analytic','min dM/M':0.0001,}     #fill in defaults if needed
     620            controls = {'deriv type':'analytic','min dM/M':0.001,}     #fill in defaults if needed
    621621        print 'Peak Fitting with '+controls['deriv type']+' derivatives:'
    622622        PatternId = G2frame.PatternId

  • trunk/GSASIIstrMain.py

    r2824 r2846  
    5959        #fl.close()
    6060        Ftol = Controls['min dM/M']
     61        Xtol = Controls['SVDtol']
    6162        Factor = Controls['shift factor']
    6263        if 'Jacobian' in Controls['deriv type']:           
    … …  
    6869            Lamda = Controls.get('Marquardt',-3)
    6970            maxCyc = Controls['max cyc']
    70             result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
     71            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
    7172                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
    7273            ncyc = result[2]['num cyc']+1
Note: See TracChangeset for help on using the changeset viewer.