Changeset 620 for trunk/GSASIIstruct.py
 Timestamp:
 May 17, 2012 7:52:59 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIstruct.py
r612 r620 424 424 HKLFdata = single crystal data list of reflections: for each reflection: 425 425 HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase] 426 need this [h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}] 426 427 ''' 427 428 fl = open(GPXfile,'rb') … … 495 496 refI = int(varyI.split(':')[1]) 496 497 ih,ik,il = PawleyRef[refI][:3] 498 dspI = PawleyRef[refI][4] 497 499 for varyJ in pawleyVary[i+1:]: 498 500 refJ = int(varyJ.split(':')[1]) 499 501 jh,jk,jl = PawleyRef[refJ][:3] 502 dspJ = PawleyRef[refJ][4] 500 503 if SGLaue in ['4/m','4/mmm']: 501 504 isum = ih**2+ik**2 … … 519 522 jsum = jh**2+jk**2+jl**2 520 523 if isum == jsum: 521 eqvDict[varyI].append(varyJ) 524 eqvDict[varyI].append(varyJ) 525 # if abs(dspIdspJ) < 3.e4: 526 # eqvDict[varyI].append(varyJ) 522 527 for item in pawleyVary: 523 528 if eqvDict[item]: … … 1746 1751 ################################################################################ 1747 1752 1748 def getPenalties(): 1749 return [] 1750 1751 def penaltyFxn(penalties): 1752 return [] 1753 1754 def penaltyDeriv(penalties): 1755 return [] 1753 def penaltyFxn(parmDict,varyList): 1754 pFxn = np.zeros(len(varyList)) 1755 for i,item in enumerate(varyList): 1756 if 'PWLref' in item and parmDict[item] < 0.: 1757 pFxn[i] = parmDict[item] 1758 return pFxn 1759 1760 def penaltyDeriv(parmDict,varyList): 1761 pDerv = np.zeros(len(varyList)) 1762 for i,item in enumerate(varyList): 1763 if 'PWLref' in item and parmDict[item] < 0.: 1764 pDerv[i] += 1./parmDict[item] 1765 return pDerv 1756 1766 1757 1767 ################################################################################ … … 2742 2752 return dMdv 2743 2753 2744 def dervRefine(values,HistoPhases,p enalties,parmdict,varylist,calcControls,pawleyLookup,dlg):2754 def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): 2745 2755 parmdict.update(zip(varylist,values)) 2746 2756 G2mv.Dict2Map(parmdict,varylist) … … 2765 2775 else: 2766 2776 dMdv = dMdvh 2767 if penalties: 2768 dmdv = np.concatenate((dmdv.T,penaltyDeriv(penalties).T)).T2777 # dpdv = penaltyDeriv(parmdict,varylist) 2778 # dMdv = np.concatenate((dMdv.T,np.outer(dpdv,dpdv))).T 2769 2779 return dMdv 2770 2780 2771 #def ComputePowderHessian(args): 2772 # Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args 2773 # hId = Histogram['hId'] 2774 # hfx = ':%d:'%(hId) 2775 # Limits = calcControls[hfx+'Limits'] 2776 # x,y,w,yc,yb,yd = Histogram['Data'] 2777 # dy = yyc 2778 # xB = np.searchsorted(x,Limits[0]) 2779 # xF = np.searchsorted(x,Limits[1]) 2780 # dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv( 2781 # parmdict,x[xB:xF], 2782 # varylist,Histogram,Phases,calcControls,pawleyLookup) 2783 # Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) 2784 # Hess = np.inner(dMdvh,dMdvh) 2785 # return Vec,Hess 2786 # 2787 def HessRefine(values,HistoPhases,penalties,parmdict,varylist,calcControls,pawleyLookup,dlg): 2781 def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): 2788 2782 parmdict.update(zip(varylist,values)) 2789 2783 G2mv.Dict2Map(parmdict,varylist) … … 2813 2807 Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) 2814 2808 Hess = np.inner(dMdvh,dMdvh) 2809 # dpdv = penaltyDeriv(parmdict,varylist) 2810 # Vec += dpdv*penaltyFxn(parmdict,varylist) 2811 # Hess += np.outer(dpdv,dpdv) 2815 2812 return Vec,Hess 2816 2813 2817 #def ComputePowderProfile(args): 2818 # Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args 2819 # hId = Histogram['hId'] 2820 # hfx = ':%d:'%(hId) 2821 # x,y,w,yc,yb,yd = Histogram['Data'] 2822 # Limits = calcControls[hfx+'Limits'] 2823 # xB = np.searchsorted(x,Limits[0]) 2824 # xF = np.searchsorted(x,Limits[1]) 2825 # yc,yb = getPowderProfile(parmdict,x[xB:xF],varylist,Histogram,Phases, 2826 # calcControls,pawleyLookup) 2827 # return xB,xF,yc,yb,Histogram['Reflection Lists'] 2828 # 2829 def errRefine(values,HistoPhases,penalties,parmdict,varylist,calcControls,pawleyLookup,dlg): 2814 def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): 2830 2815 parmdict.update(zip(varylist,values)) 2831 2816 Values2Dict(parmdict, varylist, values) … … 2871 2856 parmDict['saved values'] = values 2872 2857 raise Exception #Abort!! 2873 if penalties: 2874 M = np.concatenate((M,penaltyFxn(penalties))) 2858 # M = np.concatenate((M,penaltyFxn(parmdict,varylist))) 2875 2859 return M 2876 2860 … … 2935 2919 begin = time.time() 2936 2920 values = np.array(Dict2Values(parmDict, varyList)) 2937 penalties = getPenalties()2938 2921 Ftol = Controls['min dM/M'] 2939 2922 Factor = Controls['shift factor'] … … 2942 2925 result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True, 2943 2926 ftol=Ftol,col_deriv=True,factor=Factor, 2944 args=([Histograms,Phases],p enalties,parmDict,varyList,calcControls,pawleyLookup,dlg))2927 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 2945 2928 ncyc = int(result[2]['nfev']/2) 2946 2929 elif 'Hessian' in Controls['deriv type']: 2947 2930 result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc, 2948 args=([Histograms,Phases],p enalties,parmDict,varyList,calcControls,pawleyLookup,dlg))2931 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 2949 2932 ncyc = result[2]['num cyc']+1 2950 2933 Rvals['lamMax'] = result[2]['lamMax'] 2951 2934 else: #'numeric' 2952 2935 result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e8,factor=Factor, 2953 args=([Histograms,Phases],p enalties,parmDict,varyList,calcControls,pawleyLookup,dlg))2936 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 2954 2937 ncyc = int(result[2]['nfev']/len(varyList)) 2955 2938 # table = dict(zip(varyList,zip(values,result[0],(result[0]values)))) … … 3121 3104 begin = time.time() 3122 3105 values = np.array(Dict2Values(parmDict,varyList)) 3123 penalties = getPenalties()3124 3106 Ftol = Controls['min dM/M'] 3125 3107 Factor = Controls['shift factor'] … … 3129 3111 result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True, 3130 3112 ftol=Ftol,col_deriv=True,factor=Factor, 3131 args=([Histo,Phases],p enalties,parmDict,varyList,calcControls,pawleyLookup,dlg))3113 args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 3132 3114 ncyc = int(result[2]['nfev']/2) 3133 3115 elif 'Hessian' in Controls['deriv type']: 3134 3116 result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc, 3135 args=([Histo,Phases],p enalties,parmDict,varyList,calcControls,pawleyLookup,dlg))3117 args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 3136 3118 ncyc = result[2]['num cyc']+1 3137 3119 else: #'numeric' 3138 3120 result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e8,factor=Factor, 3139 args=([Histo,Phases],p enalties,parmDict,varyList,calcControls,pawleyLookup,dlg))3121 args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 3140 3122 ncyc = int(result[2]['nfev']/len(varyList)) 3141 3123 … … 3226 3208 Factor = DisAglCtls['Factors'] 3227 3209 Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii']))) 3228 # Units = np.array([ #is there a nicer way to make this?3229 # [1,1,1],[1,1,0],[1,1,1],[1,0,1],[1,0,0],[1,0,1],[1,1,1],[1,1,0],[1,1,1],3230 # [0,1,1],[0,1,0],[0,1,1],[0,0,1],[0,0,0],[0,0,1],[0,1,1],[0,1,0],[0,1,1],3231 # [1,1,1],[1,1,0],[1,1,1],[1,0,1],[1,0,0],[1,0,1],[1,1,1],[1,1,0],[1,1,1]])3232 3210 indices = (1,0,1) 3233 3211 Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
Note: See TracChangeset
for help on using the changeset viewer.