Changeset 1775
- Timestamp:
- Apr 2, 2015 3:43:31 PM (8 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIstrIO.py
r1773 r1775 2071 2071 controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4] 2072 2072 controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5]) 2073 controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]] 2073 controlDict[pfx+'SHhkl'] = [] 2074 if hapData['Pref.Ori.'][6][0] != '': 2075 controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]] 2074 2076 controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7] 2075 2077 for item in hapData['Pref.Ori.'][5]: -
trunk/GSASIIstrMath.py
r1770 r1775 304 304 ################################################################################ 305 305 306 def penaltyFxn(HistoPhases, parmDict,varyList):306 def penaltyFxn(HistoPhases,calcControls,parmDict,varyList): 307 307 'Needs a doc string' 308 308 Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases … … 327 327 names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'], 328 328 ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'], 329 ['ChemComp','Sites'],['Texture','HKLs'] ]329 ['ChemComp','Sites'],['Texture','HKLs'],] 330 330 for name,rest in names: 331 331 pWsum[name] = 0. … … 389 389 pVals.append(Z1[ind[0]][ind[1]]) 390 390 pWt.append(wt/esd1**2) 391 pWsum[name] += wt*( (obs-calc)/esd)**2391 pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd)**2 392 392 if ifesd2: 393 393 Z2 = 1.-Z … … 397 397 pWt.append(wt/esd2**2) 398 398 pWsum[name] += wt*((obs-calc)/esd)**2 399 399 400 name = 'SH-Pref.Ori.' 401 pWsum[name] = 0.0 402 for hist in Phases[phase]['Histograms']: 403 if hist in Histograms: 404 hId = Histograms[hist]['hId'] 405 phfx = '%d:%d:'%(pId,hId) 406 if calcControls[phfx+'poType'] == 'SH': 407 toler = calcControls[phfx+'SHtoler'] 408 wt = 1./toler**2 409 HKLs = calcControls[phfx+'SHhkl'] 410 SHcofNames = Phases[phase]['Histograms'][hist]['Pref.Ori.'][5].keys() 411 SHcof = dict(zip(SHcofNames,[parmDict[phfx+cof] for cof in SHcofNames])) 412 for i,PH in enumerate(HKLs): 413 phi,beta = G2lat.CrsAng(PH,cell,SGData) 414 SH3Coef = {} 415 for item in SHcof: 416 L,N = eval(item.strip('C')) 417 SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item] 418 ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData) 419 X = np.linspace(0,90.0,26) 420 Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0) 421 IndY = ma.nonzero(Y) 422 for ind in IndY[0]: 423 pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind])) 424 pVals.append(Y[ind]) 425 pWt.append(wt) 426 pWsum[name] += wt*(-Y[ind])**2 400 427 pWsum['PWLref'] = 0. 401 428 for item in varyList: … … 411 438 return pNames,pVals,pWt,pWsum 412 439 413 def penaltyDeriv(pNames,pVal,HistoPhases, parmDict,varyList):440 def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList): 414 441 'Needs a doc string' 415 442 Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases … … 443 470 if 'PWL' in pName: 444 471 pDerv[varyList.index(pName)][ip] += 1. 472 continue 473 elif 'SH-' in pName: 445 474 continue 446 475 id = int(pnames[2]) … … 504 533 except ValueError: 505 534 pass 535 536 # for ip,pName in enumerate(pNames): 537 # pnames = pNames.split(':') 538 # if 'SH-' in pName and pId == int(pnames[0]): 539 # name = pnames[2] 540 # L,N = eval(name.strip('C')) 541 # hId = int(pnames[1]) 542 # phfx = '%d:%d:'%(pId,hId) 543 # hklId = int(pnames[3]) 544 # gam = float(pnames[4]) 545 # HKLs = calcControls[phfx+'SHhkl'] 546 # phi,beta = G2lat.CrsAng(HKLs[hklId],cell,SGData) 547 # SHcofNames = Phases[phase]['Histograms'][hist]['Pref.Ori.'][5].keys() 548 # SHcof = dict(zip(SHcofNames,[parmDict[phfx+cof] for cof in SHcofNames])) 549 # 550 # raise Exception 551 506 552 return pDerv 507 553 … … 1939 1985 refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict) #corrected reflection position 1940 1986 Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.)) #Lorentz correction 1941 # refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift1942 1987 refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam 1943 1988 refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) … … 2459 2504 dMdv = np.sqrt(wtFactor)*dMdvh 2460 2505 2461 pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases, parmDict,varylist)2506 pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist) 2462 2507 if np.any(pVals): 2463 dpdv = penaltyDeriv(pNames,pVals,HistoPhases, parmDict,varylist)2508 dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist) 2464 2509 dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T 2465 2510 … … 2534 2579 else: 2535 2580 continue #skip non-histogram entries 2536 pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases, parmDict,varylist)2581 pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist) 2537 2582 if np.any(pVals): 2538 dpdv = penaltyDeriv(pNames,pVals,HistoPhases, parmDict,varylist)2583 dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist) 2539 2584 Vec += np.sum(dpdv*pWt*pVals,axis=1) 2540 2585 Hess += np.inner(dpdv*pWt,dpdv) … … 2686 2731 dlg.Destroy() 2687 2732 raise Exception #Abort!! 2688 pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases, parmDict,varylist)2733 pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist) 2689 2734 if len(pVals): 2690 2735 pSum = np.sum(pWt*pVals**2) 2691 2736 for name in pWsum: 2692 if pWsum :2737 if pWsum[name]: 2693 2738 print ' Penalty function for %8s = %12.5g'%(name,pWsum[name]) 2694 2739 print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
Note: See TracChangeset
for help on using the changeset viewer.