Changeset 2846
- Timestamp:
- May 30, 2017 4:40:30 PM (6 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIgrid.py
r2842 r2846 421 421 def Draw(self): 422 422 423 # def OnRadius(event):424 # event.Skip()425 # try:426 # val = float(radius.GetValue())427 # if val < 0.5:428 # raise ValueError429 # self.Sphere = val430 # except ValueError:431 # pass432 # radius.SetValue('%.3f'%(self.Sphere))433 #434 423 def OnAtomType(event): 435 424 Obj = event.GetEventObject() … … 446 435 topSizer.Add(wx.StaticText(self.panel,label=' Sphere centered at atoms: '),0,WACV) 447 436 cx,ct,cs = self.Drawing['atomPtrs'][:3] 448 # print self.Drawing.keys()449 437 for id in self.indx: 450 438 atom = self.Drawing['Atoms'][id] … … 459 447 sphereSizer = wx.BoxSizer(wx.HORIZONTAL) 460 448 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)) 465 450 sphereSizer.Add(radius,0,WACV) 466 451 mainSizer.Add(sphereSizer,0,WACV) … … 617 602 self.panel.Destroy() 618 603 self.panel = wx.Panel(self) 619 # Ind = {}620 604 mainSizer = wx.BoxSizer(wx.VERTICAL) 621 605 MatSizer = wx.BoxSizer(wx.HORIZONTAL) … … 639 623 for iy,line in enumerate(self.Trans): 640 624 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)) 642 626 Trmat.Add(item) 643 627 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)) 645 629 Trmat.Add(vec) 646 630 transSizer.Add(Trmat) … … 811 795 def Draw(self): 812 796 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 #838 797 def OnExpand(event): 839 798 self.Expand = expand.GetValue() … … 859 818 self.panel.Destroy() 860 819 self.panel = wx.Panel(self) 861 # Ind = {}862 820 mainSizer = wx.BoxSizer(wx.VERTICAL) 863 821 MatSizer = wx.BoxSizer(wx.HORIZONTAL) … … 868 826 for iy,line in enumerate(self.Trans): 869 827 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)) 876 829 Trmat.Add(item) 877 830 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)) 884 832 Trmat.Add(vec) 885 833 transSizer.Add(Trmat) … … 1112 1060 def Draw(self): 1113 1061 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 #1120 1062 def OnCent(event): 1121 1063 Obj = event.GetEventObject() … … 1143 1085 self.panel.Destroy() 1144 1086 self.panel = wx.Panel(self) 1145 # Ind = {}1146 1087 mainSizer = wx.BoxSizer(wx.VERTICAL) 1147 1088 MatSizer = wx.BoxSizer(wx.HORIZONTAL) … … 1161 1102 for iy,line in enumerate(self.Trans): 1162 1103 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)) 1169 1105 Trmat.Add(item) 1170 1106 transSizer.Add(Trmat) … … 1411 1347 for id,item in enumerate(self.data['AtomTypes']): 1412 1348 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)) 1414 1350 radiiSizer.Add(bRadii,0,WACV) 1415 1351 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)) 1417 1353 radiiSizer.Add(aRadii,0,WACV) 1418 1354 mainSizer.Add(radiiSizer,0,wx.EXPAND) … … 1422 1358 for i,name in enumerate(Names): 1423 1359 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)) 1425 1361 factorSizer.Add(bondFact) 1426 1362 mainSizer.Add(factorSizer,0,wx.EXPAND) … … 2776 2712 data = {} 2777 2713 data['deriv type'] = 'analytic Hessian' 2778 data['min dM/M'] = 0.00 012714 data['min dM/M'] = 0.001 2779 2715 data['shift factor'] = 1. 2780 2716 data['max cyc'] = 3 2781 2717 data['F**2'] = False 2718 if 'SVDtol' not in data: 2719 data['SVDtol'] = 1.e-6 2782 2720 if 'shift factor' not in data: 2783 2721 data['shift factor'] = 1. … … 2901 2839 LSSizer.Add(derivSel,0,WACV) 2902 2840 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) 2905 2842 if 'Hessian' in data['deriv type']: 2906 2843 LSSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Max cycles: '),0,WACV) … … 2916 2853 marqLam.Bind(wx.EVT_COMBOBOX,OnMarqLam) 2917 2854 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) 2918 2857 else: #TODO what for SVD refine? 2919 2858 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.) 2921 2860 LSSizer.Add(Factr,0,WACV) 2922 2861 if G2frame.Sngl: … … 2934 2873 LSSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=usrRej[item][0]),0,WACV) 2935 2874 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]) 2937 2876 LSSizer.Add(usrrej,0,WACV) 2938 2877 return LSSizer … … 4390 4329 G2plt.PlotDeltSig(G2frame,kind) 4391 4330 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'] = val4399 # wtval.SetValue('%.3f'%(val))4400 #4401 4331 # def OnCompression(event): 4402 4332 # data[0] = int(comp.GetValue()) … … 4451 4381 wtSizer = wx.BoxSizer(wx.HORIZONTAL) 4452 4382 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) 4457 4384 wtSizer.Add(wtval,0,WACV) 4458 4385 # if kind == 'PWDR': #possible future compression feature; NB above patch as well -
trunk/GSASIImath.py
r2840 r2846 47 47 ################################################################################ 48 48 49 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,lamda=-3,Print=False): 49 def 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 106 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False): 50 107 51 108 """ … … 63 120 :param tuple args: Any extra arguments to func are placed in this tuple. 64 121 :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. 66 123 :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine 67 124 until other limits are met (ftol, xtol) … … 105 162 nfev = 0 106 163 if Print: 107 print ' Hessian refinement on %d variables:'%(n)164 print ' Hessian SVD refinement on %d variables:'%(n) 108 165 Lam = np.zeros((n,n)) 109 166 while icycle < maxcyc: … … 126 183 Amatlam = Amat*(One+Lam) 127 184 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 129 187 except nl.LinAlgError: 130 188 print 'ouch #1' … … 163 221 Lam = np.eye(Amat.shape[0])*lam 164 222 Amatlam = Amat/Anorm 165 # Amatlam = Amat*(One+Lam) #scale Amat to Marquardt array?166 223 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 169 227 return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}] 170 228 except nl.LinAlgError: … … 173 231 if maxcyc: 174 232 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 177 235 def getVCov(varyNames,varyList,covMatrix): 178 236 '''obtain variance-covariance terms for a set of variables. NB: the varyList -
trunk/GSASIIobj.py
r2845 r2846 874 874 DefaultControls = { 875 875 'deriv type':'analytic Hessian', 876 'min dM/M':0.00 01,'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, 877 877 'UsrReject':{'minF/sig':0,'MinExt':0.01,'MaxDF/F':100.,'MaxD':500.,'MinD':0.05}, 878 878 'Copy2Next':False,'Reverse Seq':False,'HatomFix':False, -
trunk/GSASIIpwdGUI.py
r2820 r2846 560 560 dlg = wx.ProgressDialog('Sequential peak fit','Data set name = '+names[0],len(names), 561 561 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.00 01,}562 controls = {'deriv type':'analytic','min dM/M':0.001,} 563 563 print 'Peak Fitting with '+controls['deriv type']+' derivatives:' 564 564 oneCycle = False … … 618 618 controls = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.root, 'Controls')) 619 619 if not controls: 620 controls = {'deriv type':'analytic','min dM/M':0.00 01,} #fill in defaults if needed620 controls = {'deriv type':'analytic','min dM/M':0.001,} #fill in defaults if needed 621 621 print 'Peak Fitting with '+controls['deriv type']+' derivatives:' 622 622 PatternId = G2frame.PatternId -
trunk/GSASIIstrMain.py
r2824 r2846 59 59 #fl.close() 60 60 Ftol = Controls['min dM/M'] 61 Xtol = Controls['SVDtol'] 61 62 Factor = Controls['shift factor'] 62 63 if 'Jacobian' in Controls['deriv type']: … … 68 69 Lamda = Controls.get('Marquardt',-3) 69 70 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, 71 72 args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)) 72 73 ncyc = result[2]['num cyc']+1
Note: See TracChangeset
for help on using the changeset viewer.