Changeset 3136 for trunk/GSASIIpwd.py
- Timestamp:
- Oct 23, 2017 11:39:16 AM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified trunk/GSASIIpwd.py ¶
r3070 r3136 13 13 # $Id$ 14 14 ########### SVN repository information ################### 15 from __future__ import division, print_function 15 16 import sys 16 17 import math … … 40 41 import pydiffax as pyx 41 42 except ImportError: 42 print 'pydiffax is not available for this platform - under develpment'43 print ('pydiffax is not available for this platform - under develpment') 43 44 44 45 … … 434 435 parmDict,varyList = MakeParms(peaks) 435 436 if not len(varyList): 436 print ' Nothing varied'437 print (' Nothing varied') 437 438 return newpeaks,None,None,None,None,None 438 439 … … 529 530 nmax = max(nmin+1,nmax) 530 531 for p in range(nmin,nmax): 531 if max( factorize(p).keys()) < thresh:532 if max(list(factorize(p).keys())) < thresh: 532 533 plist.append(p) 533 534 return plist … … 822 823 break 823 824 except ValueError: 824 print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'825 print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****') 825 826 break 826 827 return yb,sumBk … … 949 950 break 950 951 except ValueError: 951 print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'952 print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****') 952 953 break 953 954 return dydb,dyddb,dydpk … … 1432 1433 1433 1434 def InstPrint(Inst,sigDict): 1434 print 'Instrument Parameters:'1435 print ('Instrument Parameters:') 1435 1436 if 'C' in Inst['Type'][0]: 1436 1437 ptfmt = "%12.6f" … … 1448 1449 else: 1449 1450 sigstr += 12*' ' 1450 print ptlbls1451 print ptstr1452 print sigstr1451 print (ptlbls) 1452 print (ptstr) 1453 print (sigstr) 1453 1454 1454 1455 def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList): … … 1473 1474 varyList = insVary 1474 1475 if not len(varyList): 1475 print '**** ERROR - nothing to refine! ****'1476 print ('**** ERROR - nothing to refine! ****') 1476 1477 return False 1477 1478 while True: … … 1485 1486 Values2Dict(parmDict, varyList, result[0]) 1486 1487 GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2 1487 print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)1488 print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)1489 print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)1488 print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList))) 1489 print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)) 1490 print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)) 1490 1491 try: 1491 1492 sig = np.sqrt(np.diag(result[1])*GOF) 1492 1493 if np.any(np.isnan(sig)): 1493 print '*** Least squares aborted - some invalid esds possible ***'1494 print ('*** Least squares aborted - some invalid esds possible ***') 1494 1495 break #refinement succeeded - finish up! 1495 1496 except ValueError: #result[1] is None on singular matrix 1496 print '**** Refinement failed - singular matrix ****'1497 print ('**** Refinement failed - singular matrix ****') 1497 1498 1498 1499 sigDict = dict(zip(varyList,sig)) … … 1559 1560 1560 1561 def BackgroundPrint(Background,sigDict): 1561 print 'Background coefficients for',Background[0][0],'function'1562 print ('Background coefficients for '+Background[0][0]+' function') 1562 1563 ptfmt = "%12.5f" 1563 1564 ptstr = 'value: ' … … 1572 1573 sigstr += " "*12 1573 1574 if len(ptstr) > 75: 1574 print ptstr1575 if Background[0][1]: print sigstr1575 print (ptstr) 1576 if Background[0][1]: print (sigstr) 1576 1577 ptstr = 'value: ' 1577 1578 sigstr = 'esd : ' 1578 1579 if len(ptstr) > 8: 1579 print ptstr1580 if Background[0][1]: print sigstr1580 print (ptstr) 1581 if Background[0][1]: print (sigstr) 1581 1582 1582 1583 if Background[1]['nDebye']: 1583 1584 parms = ['DebyeA;','DebyeR;','DebyeU;'] 1584 print 'Debye diffuse scattering coefficients'1585 print ('Debye diffuse scattering coefficients') 1585 1586 ptfmt = "%12.5f" 1586 print ' term DebyeA esd DebyeR esd DebyeU esd'1587 print (' term DebyeA esd DebyeR esd DebyeU esd') 1587 1588 for term in range(Background[1]['nDebye']): 1588 1589 line = ' term %d'%(term) … … 1593 1594 else: 1594 1595 line += " "*12 1595 print line1596 print (line) 1596 1597 if Background[1]['nPeaks']: 1597 print 'Coefficients for Background Peaks'1598 print ('Coefficients for Background Peaks') 1598 1599 ptfmt = "%15.3f" 1599 1600 for j,pl in enumerate(Background[1]['peaksList']): … … 1610 1611 else: 1611 1612 sigstr += " "*15 1612 print names1613 print ptstr1614 print sigstr1613 print (names) 1614 print (ptstr) 1615 print (sigstr) 1615 1616 1616 1617 def SetInstParms(Inst): … … 1658 1659 1659 1660 def InstPrint(Inst,sigDict): 1660 print 'Instrument Parameters:'1661 print ('Instrument Parameters:') 1661 1662 ptfmt = "%12.6f" 1662 1663 ptlbls = 'names :' … … 1672 1673 else: 1673 1674 sigstr += 12*' ' 1674 print ptlbls1675 print ptstr1676 print sigstr1675 print (ptlbls) 1676 print (ptstr) 1677 print (sigstr) 1677 1678 1678 1679 def SetPeaksParms(dataType,Peaks): … … 1722 1723 1723 1724 def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW): 1724 print 'Peak coefficients:'1725 print ('Peak coefficients:') 1725 1726 if 'C' in dataType: 1726 1727 names = ['pos','int','sig','gam'] … … 1734 1735 head += name.center(10)+'esd'.center(10) 1735 1736 head += 'bins'.center(8) 1736 print head1737 print (head) 1737 1738 if 'C' in dataType: 1738 1739 ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"} … … 1753 1754 ptstr += 10*' ' 1754 1755 ptstr += '%9.2f'%(ptsperFW[i]) 1755 print '%s'%(('Peak'+str(i+1)).center(8)),ptstr1756 print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr) 1756 1757 1757 1758 def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg): … … 1809 1810 Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100. #to % 1810 1811 Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2 1811 print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)1812 print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList))) 1812 1813 if ncyc: 1813 print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)1814 print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])1814 print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)) 1815 print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])) 1815 1816 sig = [0]*len(varyList) 1816 1817 if len(varyList) == 0: break # if nothing was refined … … 1818 1819 sig = np.sqrt(np.diag(result[1])*Rvals['GOF']) 1819 1820 if np.any(np.isnan(sig)): 1820 print '*** Least squares aborted - some invalid esds possible ***'1821 print ('*** Least squares aborted - some invalid esds possible ***') 1821 1822 break #refinement succeeded - finish up! 1822 1823 except ValueError: #result[1] is None on singular matrix 1823 print '**** Refinement failed - singular matrix ****'1824 print ('**** Refinement failed - singular matrix ****') 1824 1825 Ipvt = result[2]['ipvt'] 1825 1826 for i,ipvt in enumerate(Ipvt): 1826 1827 if not np.sum(result[2]['fjac'],axis=1)[i]: 1827 print 'Removing parameter: ',varyList[ipvt-1]1828 print ('Removing parameter: '+varyList[ipvt-1]) 1828 1829 badVary.append(varyList[ipvt-1]) 1829 1830 del(varyList[ipvt-1]) … … 1851 1852 if len(binsperFWHM): 1852 1853 if min(binsperFWHM) < 1.: 1853 print '*** Warning: calculated peak widths are too narrow to refine profile coefficients ***'1854 print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***') 1854 1855 if 'T' in Inst['Type'][0]: 1855 print ' Manually increase sig-0, 1, or 2 in Instrument Parameters'1856 print (' Manually increase sig-0, 1, or 2 in Instrument Parameters') 1856 1857 else: 1857 print ' Manually increase W in Instrument Parameters'1858 print (' Manually increase W in Instrument Parameters') 1858 1859 elif min(binsperFWHM) < 4.: 1859 print '*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***'1860 print '*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM))1860 print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***') 1861 print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM))) 1861 1862 return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary 1862 1863 … … 1921 1922 1922 1923 def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data): 1923 print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])1924 print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])) 1924 1925 1925 1926 class RandomDisplacementBounds(object): … … 1980 1981 data['Scale'][0] = parmDict['Scale'] 1981 1982 line += ' esd: %.4g'%(sigDict['Scale']) 1982 print line1983 print (line) 1983 1984 line = ' Flat background: %15.4g'%(parmDict['FltBack']) 1984 1985 if 'FltBack' in varyList: 1985 1986 data['FltBack'][0] = parmDict['FltBack'] 1986 1987 line += ' esd: %15.3g'%(sigDict['FltBack']) 1987 print line1988 print (line) 1988 1989 for ilay,layer in enumerate(data['Layers']): 1989 1990 name = layer['Name'] 1990 print ' Parameters for layer: %d %s'%(ilay,name)1991 print (' Parameters for layer: %d %s'%(ilay,name)) 1991 1992 cid = str(ilay)+';' 1992 1993 line = ' ' … … 2002 2003 if cid+parm in varyList: 2003 2004 line += ' esd: %.3g'%(sigDict[cid+parm]) 2004 print line2005 print line22005 print (line) 2006 print (line2) 2006 2007 2007 2008 def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList): … … 2082 2083 take_step = RandomDisplacementBounds(xyrng[0], xyrng[1]) 2083 2084 T0 = estimateT0(take_step) 2084 print ' Estimated temperature: %.3g'%(T0)2085 print (' Estimated temperature: %.3g'%(T0)) 2085 2086 result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac, 2086 2087 interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds, … … 2102 2103 ncalc = result[3] 2103 2104 covM = [] 2104 print ' MC/SA final temperature: %.4g'%(result[2])2105 print (' MC/SA final temperature: %.4g'%(result[2])) 2105 2106 elif data['Minimizer'] == 'L-BFGS-B': 2106 2107 result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds, #ftol=Ftol, … … 2147 2148 covMatrix = [] 2148 2149 sigDict = dict(zip(varyList,sig)) 2149 print ' Results of reflectometry data modelling fit:'2150 print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)2151 print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])2150 print (' Results of reflectometry data modelling fit:') 2151 print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList))) 2152 print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])) 2152 2153 SetModelParms() 2153 2154 return True,result,varyList,sig,Rvals,covMatrix,parmDict,'' 2154 2155 except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans 2155 print Msg2156 print (Msg) 2156 2157 return False,0,0,0,0,0,0,Msg 2157 2158 … … 2336 2337 2337 2338 def makeRefdFFT(Limits,Profile): 2338 print 'make fft'2339 print ('make fft') 2339 2340 Q,Io = Profile[:2] 2340 2341 Qmin = Limits[1][0] … … 2401 2402 if 'bin' in name: 2402 2403 DIFFaX = name+'/DIFFaX.exe' 2403 print ' Execute ',DIFFaX2404 print (' Execute '+DIFFaX) 2404 2405 break 2405 2406 # make form factor file that DIFFaX wants - atom types are GSASII style 2406 2407 sf = open('data.sfc','w') 2407 2408 sf.write('GSASII special form factor file for DIFFaX\n\n') 2408 atTypes = Layers['AtInfo'].keys()2409 atTypes = list(Layers['AtInfo'].keys()) 2409 2410 if 'H' not in atTypes: 2410 2411 atTypes.insert(0,'H') … … 2520 2521 sumPx += p 2521 2522 if sumPx != 1.0: #this has to be picky since DIFFaX is. 2522 print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'2523 print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx) 2523 2524 df.close() 2524 2525 os.remove('data.sfc') … … 2531 2532 subp.call(DIFFaX) 2532 2533 except OSError: 2533 print ' DIFFax.exe is not available for this platform - under development'2534 print ' DIFFaX time = %.2fs'%(time.time()-time0)2534 print (' DIFFax.exe is not available for this platform - under development') 2535 print (' DIFFaX time = %.2fs'%(time.time()-time0)) 2535 2536 if os.path.exists('GSASII-DIFFaX.spc'): 2536 2537 Xpat = np.loadtxt('GSASII-DIFFaX.spc').T … … 2628 2629 Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7] 2629 2630 # atoms in layers 2630 atTypes = Layers['AtInfo'].keys()2631 atTypes = list(Layers['AtInfo'].keys()) 2631 2632 AtomXOU = [] 2632 2633 AtomTp = [] … … 2697 2698 time0 = time.time() 2698 2699 pyx.pygetspc(controls,Nspec,spec) 2699 print ' GETSPC time = %.2fs'%(time.time()-time0)2700 print (' GETSPC time = %.2fs'%(time.time()-time0)) 2700 2701 time0 = time.time() 2701 2702 U = ateln2*inst['U'][1]/10000. … … 2728 2729 profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0) 2729 2730 profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin] 2730 print ' Broadening time = %.2fs'%(time.time()-time0)2731 print (' Broadening time = %.2fs'%(time.time()-time0)) 2731 2732 2732 2733 def CalcStackingSADP(Layers,debug): … … 2763 2764 iB += Nblk 2764 2765 Layers['Sadp']['Img'] = Sapd 2765 print ' GETSAD time = %.2fs'%(time.time()-time0)2766 print (' GETSAD time = %.2fs'%(time.time()-time0)) 2766 2767 # GSASIIpath.IPyBreak() 2767 2768 … … 2820 2821 for i in range(100): 2821 2822 getPeakProfile(parmDict1,xdata,varyList,bakType) 2822 print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time02823 print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0) 2823 2824 2824 2825 def test2(name,delt): … … 2858 2859 for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]: 2859 2860 test3(name,shft) 2860 print "OK"2861 print ("OK") 2861 2862 plotter.StartEventLoop()
Note: See TracChangeset
for help on using the changeset viewer.