Changeset 399
 Timestamp:
 Oct 25, 2011 4:53:33 PM (11 years ago)
 Location:
 trunk
 Files:

 7 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASII.py
r397 r399 64 64 65 65 [wxID_FILECLOSE, wxID_FILEEXIT, wxID_FILEOPEN, wxID_FILESAVE, wxID_FILESAVEAS, 66 wxID_REFINE, wxID_SOLVE, wxID_MAKEPDFS, wxID_VIEWLSPARMS, 67 ] = [wx.NewId() for _init_coll_File_Items in range( 9)]66 wxID_REFINE, wxID_SOLVE, wxID_MAKEPDFS, wxID_VIEWLSPARMS, wxID_SEQREFINE, 67 ] = [wx.NewId() for _init_coll_File_Items in range(10)] 68 68 69 69 [wxID_PWDRREAD,wxID_SNGLREAD,wxID_ADDPHASE,wxID_DELETEPHASE, … … 153 153 self.Refine.Enable(False) 154 154 self.Bind(wx.EVT_MENU, self.OnRefine, id=wxID_REFINE) 155 self.SeqRefine = parent.Append(help='', id=wxID_SEQREFINE, kind=wx.ITEM_NORMAL, 156 text='Sequental refine') 157 self.SeqRefine.Enable(False) 158 self.Bind(wx.EVT_MENU, self.OnSeqRefine, id=wxID_SEQREFINE) 155 159 self.Solve = parent.Append(help='', id=wxID_SOLVE, kind=wx.ITEM_NORMAL, 156 160 text='Solve') … … 299 303 self.PatternTree.Expand(self.root) 300 304 self.Refine.Enable(True) 305 self.SeqRefine.Enable(True) 301 306 self.Solve.Enable(True) 302 307 … … 1077 1082 if data: 1078 1083 self.Refine.Enable(True) 1084 self.SeqRefine.Enable(True) 1079 1085 self.Solve.Enable(True) #not right but something needed here 1080 1086 item, cookie = self.PatternTree.GetNextChild(self.root, cookie) … … 1517 1523 finally: 1518 1524 dlg.Destroy() 1525 1526 def OnSeqRefine(self,event): 1527 self.OnFileSave(event) 1528 Id = G2gd.GetPatternTreeItemId(self,self.root,'Sequental results') 1529 if not Id: 1530 Id = self.PatternTree.AppendItem(self.root,text='Sequental results') 1531 self.PatternTree.SetItemPyData(Id,{}) 1532 dlg = wx.ProgressDialog('Residual for histogram 0','Powder profile Rwp =',101.0, 1533 style = wx.PD_ELAPSED_TIMEwx.PD_AUTO_HIDEwx.PD_CAN_ABORT) 1534 screenSize = wx.ClientDisplayRect() 1535 Size = dlg.GetSize() 1536 dlg.SetPosition(wx.Point(screenSize[2]Size[0]305,screenSize[1]+5)) 1537 try: 1538 G2str.SeqRefine(self.GSASprojectfile,dlg) 1539 finally: 1540 dlg.Destroy() 1541 dlg = wx.MessageDialog(self,'Load new result?','Refinement results',wx.OKwx.CANCEL) 1542 try: 1543 if dlg.ShowModal() == wx.ID_OK: 1544 Id = 0 1545 self.PatternTree.DeleteChildren(self.root) 1546 if self.HKL: self.HKL = [] 1547 if self.G2plotNB.plotList: 1548 self.G2plotNB.clear() 1549 G2IO.ProjFileOpen(self) 1550 item, cookie = self.PatternTree.GetFirstChild(self.root) 1551 while item and not Id: 1552 name = self.PatternTree.GetItemText(item) 1553 if name[:4] in ['PWDR','HKLF']: 1554 Id = item 1555 item, cookie = self.PatternTree.GetNextChild(self.root, cookie) 1556 if Id: 1557 self.PatternTree.SelectItem(Id) 1558 finally: 1559 dlg.Destroy() 1519 1560 1520 1561 def OnSolve(self,event): 
trunk/GSASIIIO.py
r397 r399 801 801 break 802 802 datum = data[0] 803 print 'load: ',datum[0]804 803 805 804 Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0]) … … 809 808 self.PatternTree.SetItemPyData(Id,datum[1]) 810 809 for datus in data[1:]: 811 print ' load: ',datus[0]812 810 sub = self.PatternTree.AppendItem(Id,datus[0]) 813 811 self.PatternTree.SetItemPyData(sub,datus[1]) … … 833 831 data = [] 834 832 name = self.PatternTree.GetItemText(item) 835 print 'save: ',name836 833 data.append([name,self.PatternTree.GetItemPyData(item)]) 837 834 item2, cookie2 = self.PatternTree.GetFirstChild(item) 838 835 while item2: 839 836 name = self.PatternTree.GetItemText(item2) 840 print ' save: ',name841 837 data.append([name,self.PatternTree.GetItemPyData(item2)]) 842 838 item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2) 
trunk/GSASIIgrid.py
r397 r399 11 11 import time 12 12 import cPickle 13 import numpy as np 13 14 import GSASIIpath 14 15 import GSASIIplot as G2plt … … 461 462 return self.rowLabels[row] 462 463 464 def GetColValues(self, col): 465 data = [] 466 for row in range(self.GetNumberRows()): 467 data.append(self.GetValue(row, col)) 468 return data 469 463 470 def GetRowValues(self, row): 464 471 data = [] … … 580 587 data['shift factor'] = value 581 588 Factr.SetValue('%.5f'%(value)) 589 590 def OnSelectData(event): 591 choices = ['All',]+GetPatternTreeDataNames(self,['PWDR',]) 592 if 'Seq Data' in data: 593 sel = [] 594 for item in data['Seq Data']: 595 sel.append(choices.index(item)) 596 names = [] 597 dlg = wx.MultiChoiceDialog(self,'Select data:','Sequential refinement',choices) 598 dlg.SetSelections(sel) 599 if dlg.ShowModal() == wx.ID_OK: 600 sel = dlg.GetSelections() 601 for i in sel: names.append(choices[i]) 602 if 'All' in names: 603 names = choices[1:] 604 605 dlg.Destroy() 606 data['Seq Data'] = names 607 reverseSel.Enable(True) 608 609 def OnReverse(event): 610 data['Reverse Seq'] = reverseSel.GetValue() 582 611 583 612 if self.dataDisplay: … … 610 639 Factr.Bind(wx.EVT_KILL_FOCUS,OnFactor) 611 640 LSSizer.Add(Factr,0,wx.ALIGN_CENTER_VERTICAL) 612 613 641 mainSizer.Add(LSSizer) 614 642 mainSizer.Add((5,5),0) 643 644 SeqSizer = wx.BoxSizer(wx.HORIZONTAL) 645 SeqSizer.Add(wx.StaticText(self.dataDisplay,label=' Sequential Refinement Powder Data: '),0,wx.ALIGN_CENTER_VERTICAL) 646 selSeqData = wx.Button(self.dataDisplay,1,label=' Select data') 647 selSeqData.Bind(wx.EVT_BUTTON,OnSelectData) 648 SeqSizer.Add(selSeqData,0,wx.ALIGN_CENTER_VERTICAL) 649 SeqSizer.Add((5,0),0) 650 reverseSel = wx.CheckBox(self.dataDisplay,1,label=' Reverse order?') 651 reverseSel.Bind(wx.EVT_CHECKBOX,OnReverse) 652 if 'Seq Data' not in data: 653 reverseSel.Enable(False) 654 if 'Reverse Seq' in data: 655 reverseSel.SetValue(data['Reverse Seq']) 656 SeqSizer.Add(reverseSel,0,wx.ALIGN_CENTER_VERTICAL) 657 658 659 mainSizer.Add(SeqSizer) 660 mainSizer.Add((5,5),0) 661 615 662 mainSizer.Add(wx.StaticText(self.dataDisplay,label=' Density Map Controls:'),0,wx.ALIGN_CENTER_VERTICAL) 616 663 … … 633 680 self.dataDisplay.AppendText(line+'\n') 634 681 682 def UpdateSeqResults(self,data): 683 if not data: 684 print 'No sequential refinement results' 685 return 686 histNames = data['histNames'] 687 688 def ColSelect(event): 689 cols = self.dataDisplay.GetSelectedCols() 690 plotData = [] 691 plotNames = [] 692 for col in cols: 693 plotData.append(self.SeqTable.GetColValues(col)) 694 plotNames.append(self.SeqTable.GetColLabelValue(col)) 695 plotData = np.array(plotData) 696 G2plt.PlotSeq(self,plotData,plotNames) 697 698 if self.dataDisplay: 699 self.dataDisplay.Destroy() 700 self.dataFrame.SetMenuBar(self.dataFrame.BlankMenu) 701 self.dataFrame.SetLabel('Sequental refinement results') 702 self.dataFrame.CreateStatusBar() 703 colLabels = data['varyList'] 704 Types = len(data['varyList'])*[wg.GRID_VALUE_FLOAT,] 705 seqList = [list(data[name]['variables']) for name in histNames] 706 self.SeqTable = Table(seqList,colLabels=colLabels,rowLabels=histNames,types=Types) 707 self.dataDisplay = GSGrid(parent=self.dataFrame) 708 self.dataDisplay.SetTable(self.SeqTable, True) 709 self.dataDisplay.EnableEditing(False) 710 self.dataDisplay.Bind(wg.EVT_GRID_LABEL_LEFT_DCLICK, ColSelect) 711 self.dataDisplay.SetRowLabelSize(8*len(histNames[0])) #pretty arbitrary 8 712 self.dataDisplay.SetMargins(0,0) 713 self.dataDisplay.AutoSizeColumns(True) 714 self.dataFrame.setSizePosLeft([700,350]) 715 635 716 def UpdateConstraints(self,data): 636 717 Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree() … … 933 1014 self.dataDisplay.SetSize(mainSizer.Fit(self.dataFrame)) 934 1015 self.dataFrame.setSizePosLeft(mainSizer.Fit(self.dataFrame)) 935 1016 1017 def GetPatternTreeDataNames(self,dataTypes): 1018 names = [] 1019 item, cookie = self.PatternTree.GetFirstChild(self.root) 1020 while item: 1021 name = self.PatternTree.GetItemText(item) 1022 if name[:4] in dataTypes: 1023 names.append(name) 1024 item, cookie = self.PatternTree.GetNextChild(self.root, cookie) 1025 return names 936 1026 937 1027 def GetPatternTreeItemId(self, parentId, itemText): … … 997 1087 self.PatternTree.SetItemPyData(item,data) 998 1088 self.Refine.Enable(True) 1089 self.SeqRefine.Enable(True) 999 1090 UpdateControls(self,data) 1091 elif self.PatternTree.GetItemText(item) == 'Sequental results': 1092 data = self.PatternTree.GetItemPyData(item) 1093 UpdateSeqResults(self,data) 1000 1094 elif self.PatternTree.GetItemText(item) == 'Covariance': 1001 1095 data = self.PatternTree.GetItemPyData(item) 
trunk/GSASIIimgGUI.py
r398 r399 214 214 data['rings'] = [] 215 215 data['ellipses'] = [] 216 self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMCLEARCALIB,enable=False) 217 self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=False) 216 # self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=False) 218 217 G2plt.PlotExposedImage(self,event=event) 219 218 220 219 def OnCalibrate(event): 221 self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMCLEARCALIB,enable=True)222 220 self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=True) 223 221 self.dataFrame.GetStatusBar().SetStatusText('Select > 4 points on 1st used ring; LB to pick, RB on point to delete else RB to finish') … … 399 397 self.dataFrame.Bind(wx.EVT_MENU, OnClearCalib, id=G2gd.wxID_IMCLEARCALIB) 400 398 if not data['rings']: 401 self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMCLEARCALIB,enable=False)402 399 self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=False) 403 400 self.dataFrame.Bind(wx.EVT_MENU, OnIntegrate, id=G2gd.wxID_IMINTEGRATE) 
trunk/GSASIIphsGUI.py
r397 r399 2093 2093 #how about HKLF data? This is only for PWDR data 2094 2094 Obj = event.GetEventObject() 2095 Obj.SetValue(False)2096 2095 hist = Indx[Obj.GetId()] 2097 2096 sourceDict = UseList[hist] … … 2338 2337 showData.Bind(wx.EVT_CHECKBOX, OnShowData) 2339 2338 showSizer.Add(showData,0,wx.ALIGN_CENTER_VERTICAL) 2340 copyData = wx. CheckBox(dataDisplay,1,label=' Copy?')2339 copyData = wx.Button(dataDisplay,1,label=' Copy?') 2341 2340 Indx[copyData.GetId()] = item 2342 copyData.Bind(wx.EVT_ CHECKBOX,OnCopyData)2341 copyData.Bind(wx.EVT_BUTTON,OnCopyData) 2343 2342 showSizer.Add(copyData,wx.ALIGN_CENTER_VERTICAL) 2344 2343 mainSizer.Add(showSizer,0,wx.ALIGN_CENTER_VERTICAL) 
trunk/GSASIIplot.py
r396 r399 1433 1433 Plot.set_ylabel('Variable name') 1434 1434 Page.canvas.draw() 1435 1435 1436 def PlotSeq(self,SeqData,SeqNames): 1437 1438 try: 1439 plotNum = self.G2plotNB.plotList.index('Sequential refinement') 1440 Page = self.G2plotNB.nb.GetPage(plotNum) 1441 Page.figure.clf() 1442 Plot = Page.figure.gca() 1443 if not Page.IsShown(): 1444 Page.Show() 1445 except ValueError: 1446 Plot = self.G2plotNB.addMpl('Sequential refinement').gca() 1447 plotNum = self.G2plotNB.plotList.index('Sequential refinement') 1448 Page = self.G2plotNB.nb.GetPage(plotNum) 1449 1450 Page.SetFocus() 1451 self.G2plotNB.status.SetFields(['','']) 1452 if len(SeqData): 1453 X = np.arange(0,len(SeqData[0]),1) 1454 for Y,name in zip(SeqData,SeqNames): 1455 Plot.plot(X,Y,label=name) 1456 Plot.legend(loc='best') 1457 Page.canvas.draw() 1436 1458 1437 1459 def PlotExposedImage(self,newPlot=False,event=None): … … 1811 1833 Plot.plot([arcxI[1],arcxO[1]],[arcyI[1],arcyO[1]],picker=3) 1812 1834 for i in range(Nazm): 1813 cake = LRAzim[0]+i*delAzm 1835 cake = LRAzim[0]+i*delAzmAzmthOff 1814 1836 ind = np.searchsorted(Azm,cake) 1815 1837 Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5)) 
trunk/GSASIIstruct.py
r397 r399 109 109 return datus[1] 110 110 111 def GetHistogramNames(GPXfile ):112 '''Returns a list of histogram names found in GSASII gpx file111 def GetHistogramNames(GPXfile,hType): 112 """ Returns a list of histogram names found in GSASII gpx file 113 113 input: 114 114 GPXfile = .gpx full file name 115 hType = list ['PWDR','HKLF'] 115 116 return: 116 117 HistogramNames = list of histogram names (types = PWDR & HKLF) 117 '''118 """ 118 119 file = open(GPXfile,'rb') 119 120 HistogramNames = [] … … 124 125 break 125 126 datum = data[0] 126 if datum[0][:4] in ['PWDR','HKLF']:127 if datum[0][:4] in hType: 127 128 HistogramNames.append(datum[0]) 128 129 file.close() … … 139 140 ''' 140 141 phaseNames = GetPhaseNames(GPXfile) 142 histoList = GetHistogramNames(GPXfile,['PWDR','HKLF']) 141 143 phaseData = {} 142 144 for name in phaseNames: … … 144 146 Histograms = {} 145 147 Phases = {} 146 pId = 0147 hId = 0148 148 for phase in phaseData: 149 149 Phase = phaseData[phase] 150 150 if Phase['Histograms']: 151 151 if phase not in Phases: 152 pId = phaseNames.index(phase) 152 153 Phase['pId'] = pId 153 pId += 1154 154 Phases[phase] = Phase 155 155 for hist in Phase['Histograms']: … … 160 160 Histograms[hist] = GetHKLFdata(GPXfile,hist) 161 161 #future restraint, etc. histograms here 162 hId = histoList.index(hist) 162 163 Histograms[hist]['hId'] = hId 163 hId += 1164 164 return Histograms,Phases 165 165 166 def GPXBackup(GPXfile): 167 import distutils.file_util as dfu 168 GPXpath,GPXname = ospath.split(GPXfile) 169 Name = ospath.splitext(GPXname)[0] 170 files = os.listdir(GPXpath) 171 last = 0 172 for name in files: 173 name = name.split('.') 174 if len(name) == 3 and name[0] == Name and 'bak' in name[1]: 175 last = max(last,int(name[1].strip('bak'))+1) 176 GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx') 177 dfu.copy_file(GPXfile,GPXback) 178 return GPXback 179 166 180 def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData): 167 181 ''' Updates gpxfile from all histograms that are found in any phase … … 174 188 ''' 175 189 176 def GPXBackup(GPXfile):177 import distutils.file_util as dfu178 GPXpath,GPXname = ospath.split(GPXfile)179 Name = ospath.splitext(GPXname)[0]180 files = os.listdir(GPXpath)181 last = 0182 for name in files:183 name = name.split('.')184 if len(name) == 3 and name[0] == Name and 'bak' in name[1]:185 last = max(last,int(name[1].strip('bak'))+1)186 GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')187 dfu.copy_file(GPXfile,GPXback)188 return GPXback189 190 190 GPXback = GPXBackup(GPXfile) 191 191 print '\n',135*'' … … 200 200 break 201 201 datum = data[0] 202 print 'read: ',datum[0]202 # print 'read: ',datum[0] 203 203 if datum[0] == 'Phases': 204 204 for iphase in range(len(data)): … … 210 210 try: 211 211 histogram = Histograms[datum[0]] 212 print 'found ',datum[0]212 # print 'found ',datum[0] 213 213 data[0][1][1] = histogram['Data'] 214 214 for datus in data[1:]: 215 print ' read: ',datus[0]215 # print ' read: ',datus[0] 216 216 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: 217 217 datus[1] = histogram[datus[0]] … … 223 223 outfile.close() 224 224 print 'refinement save successful' 225 226 def SetSeqResult(GPXfile,Histograms,SeqResult): 227 GPXback = GPXBackup(GPXfile) 228 print '\n',135*'' 229 print 'Read from file:',GPXback 230 print 'Save to file :',GPXfile 231 infile = open(GPXback,'rb') 232 outfile = open(GPXfile,'wb') 233 while True: 234 try: 235 data = cPickle.load(infile) 236 except EOFError: 237 break 238 datum = data[0] 239 if datum[0] == 'Sequental results': 240 data[0][1] = SeqResult 241 try: 242 histogram = Histograms[datum[0]] 243 data[0][1][1] = histogram['Data'] 244 for datus in data[1:]: 245 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: 246 datus[1] = histogram[datus[0]] 247 except KeyError: 248 pass 249 250 cPickle.dump(data,outfile,1) 251 infile.close() 252 outfile.close() 253 print 'refinement save successful' 254 225 255 226 256 def GetPWDRdata(GPXfile,PWDRname): … … 878 908 pId = Phases[phase]['pId'] 879 909 for histogram in HistoPhase: 910 try: 911 Histogram = Histograms[histogram] 912 except KeyError: 913 #skip if histogram not included e.g. in a sequential refinement 914 continue 880 915 hapData = HistoPhase[histogram] 881 Histogram = Histograms[histogram]882 916 hId = Histogram['hId'] 883 917 limits = Histogram['Limits'][1] … … 968 1002 return hapVary,hapDict,controlDict 969 1003 970 def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms ):1004 def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True): 971 1005 972 1006 def PrintSizeAndSig(hapData,sizeSig): … … 1064 1098 pId = Phases[phase]['pId'] 1065 1099 for histogram in HistoPhase: 1100 try: 1101 Histogram = Histograms[histogram] 1102 except KeyError: 1103 #skip if histogram not included e.g. in a sequential refinement 1104 continue 1066 1105 print '\n Phase: ',phase,' in histogram: ',histogram 1067 1106 print 130*'' 1068 1107 hapData = HistoPhase[histogram] 1069 Histogram = Histograms[histogram]1070 1108 hId = Histogram['hId'] 1071 1109 pfx = str(pId)+':'+str(hId)+':' … … 1087 1125 if pfx+item in sigDict: 1088 1126 PhFrExtPOSig[item] = sigDict[pfx+item] 1089 if 'Scale' in PhFrExtPOSig: 1090 print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale']) 1091 if 'Extinction' in PhFrExtPOSig: 1092 print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction']) 1093 if hapData['Pref.Ori.'][0] == 'MD': 1094 if 'MD' in PhFrExtPOSig: 1095 print ' MarchDollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD']) 1096 else: 1097 PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig) 1127 if Print: 1128 if 'Scale' in PhFrExtPOSig: 1129 print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale']) 1130 if 'Extinction' in PhFrExtPOSig: 1131 print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction']) 1132 if hapData['Pref.Ori.'][0] == 'MD': 1133 if 'MD' in PhFrExtPOSig: 1134 print ' MarchDollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD']) 1135 else: 1136 PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig) 1098 1137 SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]], 1099 1138 'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]], … … 1123 1162 if pfx+name in sigDict: 1124 1163 SizeMuStrSig['HStrain'][name] = sigDict[pfx+name] 1125 PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size']) 1126 PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData) 1127 PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData) 1164 if Print: 1165 PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size']) 1166 PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData) 1167 PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData) 1128 1168 1129 1169 def GetHistogramData(Histograms,Print=True): … … 1267 1307 return histVary,histDict,controlDict 1268 1308 1269 def SetHistogramData(parmDict,sigDict,Histograms ):1309 def SetHistogramData(parmDict,sigDict,Histograms,Print=True): 1270 1310 1271 1311 def SetBackgroundParms(pfx,Background,parmDict,sigDict): … … 1376 1416 print '\n Histogram: ',histogram,' histogram Id: ',hId 1377 1417 print 135*'' 1378 print ' Instrument type: ',Sample['Type']1379 1418 print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs']) 1380 PrintSampleParmsSig(Sample,sampSig) 1381 PrintInstParmsSig(Inst,instSig) 1382 PrintBackgroundSig(Background,backSig) 1419 if Print: 1420 print ' Instrument type: ',Sample['Type'] 1421 PrintSampleParmsSig(Sample,sampSig) 1422 PrintInstParmsSig(Inst,instSig) 1423 PrintBackgroundSig(Background,backSig) 1383 1424 1384 1425 def GetAtomFXU(pfx,FFtables,calcControls,parmDict): … … 2176 2217 except IndexError: 2177 2218 pass 2178 return dMdv 2219 return dMdv 2220 2221 def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): 2222 parmdict.update(zip(varylist,values)) 2223 G2mv.Dict2Map(parmdict) 2224 Histograms,Phases = HistoPhases 2225 dMdv = np.empty(0) 2226 for histogram in Histograms: 2227 if 'PWDR' in histogram[:4]: 2228 Histogram = Histograms[histogram] 2229 hId = Histogram['hId'] 2230 hfx = ':%d:'%(hId) 2231 Limits = calcControls[hfx+'Limits'] 2232 x,y,w,yc,yb,yd = Histogram['Data'] 2233 xB = np.searchsorted(x,Limits[0]) 2234 xF = np.searchsorted(x,Limits[1]) 2235 dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], 2236 varylist,Histogram,Phases,calcControls,pawleyLookup) 2237 if len(dMdv): 2238 dMdv = np.concatenate((dMdv.T,dMdvh.T)).T 2239 else: 2240 dMdv = dMdvh 2241 return dMdv 2242 2243 def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): 2244 parmdict.update(zip(varylist,values)) 2245 Values2Dict(parmdict, varylist, values) 2246 G2mv.Dict2Map(parmdict) 2247 Histograms,Phases = HistoPhases 2248 M = np.empty(0) 2249 sumwYo = 0 2250 Nobs = 0 2251 for histogram in Histograms: 2252 if 'PWDR' in histogram[:4]: 2253 Histogram = Histograms[histogram] 2254 hId = Histogram['hId'] 2255 hfx = ':%d:'%(hId) 2256 Limits = calcControls[hfx+'Limits'] 2257 x,y,w,yc,yb,yd = Histogram['Data'] 2258 yc *= 0.0 #zero full calcd profiles 2259 yb *= 0.0 2260 yd *= 0.0 2261 xB = np.searchsorted(x,Limits[0]) 2262 xF = np.searchsorted(x,Limits[1]) 2263 Histogram['Nobs'] = xFxB 2264 Nobs += Histogram['Nobs'] 2265 Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) 2266 sumwYo += Histogram['sumwYo'] 2267 yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF], 2268 varylist,Histogram,Phases,calcControls,pawleyLookup) 2269 yc[xB:xF] += yb[xB:xF] 2270 yd[xB:xF] = yc[xB:xF]y[xB:xF] #ycyo then all dydv have no '' needed 2271 Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) 2272 wdy = np.sqrt(w[xB:xF])*(yd[xB:xF]) 2273 Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) 2274 M = np.concatenate((M,wdy)) 2275 Histograms['sumwYo'] = sumwYo 2276 Histograms['Nobs'] = Nobs 2277 Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.) 2278 if dlg: 2279 GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile wRp =',Rwp,'%'))[0] 2280 if not GoOn: 2281 parmDict['saved values'] = values 2282 raise Exception #Abort!! 2283 return M 2284 2179 2285 2180 2286 def Refine(GPXfile,dlg): … … 2182 2288 import pytexture as ptx 2183 2289 ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics 2184 2185 def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):2186 parmdict.update(zip(varylist,values))2187 G2mv.Dict2Map(parmDict)2188 Histograms,Phases = HistoPhases2189 dMdv = np.empty(0)2190 for histogram in Histograms:2191 if 'PWDR' in histogram[:4]:2192 Histogram = Histograms[histogram]2193 hId = Histogram['hId']2194 hfx = ':%d:'%(hId)2195 Limits = calcControls[hfx+'Limits']2196 x,y,w,yc,yb,yd = Histogram['Data']2197 xB = np.searchsorted(x,Limits[0])2198 xF = np.searchsorted(x,Limits[1])2199 dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],2200 varylist,Histogram,Phases,calcControls,pawleyLookup)2201 if len(dMdv):2202 dMdv = np.concatenate((dMdv.T,dMdvh.T)).T2203 else:2204 dMdv = dMdvh2205 return dMdv2206 2207 def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):2208 parmdict.update(zip(varylist,values))2209 Values2Dict(parmdict, varylist, values)2210 G2mv.Dict2Map(parmDict)2211 Histograms,Phases = HistoPhases2212 M = np.empty(0)2213 sumwYo = 02214 Nobs = 02215 for histogram in Histograms:2216 if 'PWDR' in histogram[:4]:2217 Histogram = Histograms[histogram]2218 hId = Histogram['hId']2219 hfx = ':%d:'%(hId)2220 Limits = calcControls[hfx+'Limits']2221 x,y,w,yc,yb,yd = Histogram['Data']2222 yc *= 0.0 #zero full calcd profiles2223 yb *= 0.02224 yd *= 0.02225 xB = np.searchsorted(x,Limits[0])2226 xF = np.searchsorted(x,Limits[1])2227 Histogram['Nobs'] = xFxB2228 Nobs += Histogram['Nobs']2229 Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)2230 sumwYo += Histogram['sumwYo']2231 yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],2232 varylist,Histogram,Phases,calcControls,pawleyLookup)2233 yc[xB:xF] += yb[xB:xF]2234 yd[xB:xF] = yc[xB:xF]y[xB:xF] #ycyo then all dydv have no '' needed2235 Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))2236 wdy = np.sqrt(w[xB:xF])*(yd[xB:xF])2237 Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)2238 M = np.concatenate((M,wdy))2239 Histograms['sumwYo'] = sumwYo2240 Histograms['Nobs'] = Nobs2241 Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)2242 if dlg:2243 GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile wRp =',Rwp,'%'))[0]2244 if not GoOn:2245 parmDict['saved values'] = values2246 raise Exception #Abort!!2247 return M2248 2290 2249 2291 ShowBanner() … … 2351 2393 # file.close() 2352 2394 2395 def SeqRefine(GPXfile,dlg): 2396 import cPickle 2397 import pytexture as ptx 2398 ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics 2399 2400 ShowBanner() 2401 print ' Sequential Refinement' 2402 G2mv.InitVars() 2403 Controls = GetControls(GPXfile) 2404 ShowControls(Controls) 2405 Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) 2406 if not Phases: 2407 print ' *** ERROR  you have no histograms to refine! ***' 2408 print ' *** Refine aborted ***' 2409 raise Exception 2410 if not Histograms: 2411 print ' *** ERROR  you have no data to refine with! ***' 2412 print ' *** Refine aborted ***' 2413 raise Exception 2414 Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases,False) 2415 if 'Seq Data' in Controls: 2416 histNames = Controls['Seq Data'] 2417 else: 2418 histNames = GetHistogramNames(GPXfile,['PWDR',]) 2419 if 'Reverse Seq' in Controls: 2420 if Controls['Reverse Seq']: 2421 histNames.reverse() 2422 SeqResult = {'histNames':histNames} 2423 2424 for ihst,histogram in enumerate(histNames): 2425 print ihst,histogram 2426 ifPrint = False 2427 if dlg: 2428 dlg.SetTitle('Residual for histogram '+str(ihst)) 2429 calcControls = {} 2430 calcControls['Natoms'] = Natoms 2431 calcControls['FFtables'] = FFtables 2432 varyList = [] 2433 parmDict = {} 2434 Histo = {histogram:Histograms[histogram],} 2435 hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False) 2436 calcControls.update(controlDict) 2437 histVary,histDict,controlDict = GetHistogramData(Histo,False) 2438 calcControls.update(controlDict) 2439 varyList = phaseVary+hapVary+histVary 2440 if not ihst: 2441 saveVaryList = varyList[:] 2442 for i,item in enumerate(saveVaryList): 2443 items = item.split(':') 2444 if items[1]: 2445 items[1] = '' 2446 item = ':'.join(items) 2447 saveVaryList[i] = item 2448 SeqResult['varyList'] = saveVaryList 2449 else: 2450 newVaryList = varyList[:] 2451 for i,item in enumerate(newVaryList): 2452 items = item.split(':') 2453 if items[1]: 2454 items[1] = '' 2455 item = ':'.join(items) 2456 newVaryList[i] = item 2457 if newVaryList != SeqResult['varyList']: 2458 print newVaryList 2459 print SeqResult['varyList'] 2460 print '**** ERROR  variable list for this histogram does not match previous' 2461 raise Exception 2462 parmDict.update(phaseDict) 2463 parmDict.update(hapDict) 2464 parmDict.update(histDict) 2465 GetFprime(calcControls,Histo) 2466 constrDict,constrFlag,fixedList = G2mv.InputParse([]) #constraints go here? 2467 groups,parmlist = G2mv.GroupConstraints(constrDict) 2468 G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList) 2469 G2mv.Map2Dict(parmDict,varyList) 2470 2471 while True: 2472 begin = time.time() 2473 values = np.array(Dict2Values(parmDict, varyList)) 2474 Ftol = Controls['min dM/M'] 2475 Factor = Controls['shift factor'] 2476 if Controls['deriv type'] == 'analytic': 2477 result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True, 2478 ftol=Ftol,col_deriv=True,factor=Factor, 2479 args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 2480 ncyc = int(result[2]['nfev']/2) 2481 else: #'numeric' 2482 result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e8,factor=Factor, 2483 args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 2484 ncyc = int(result[2]['nfev']/len(varyList)) 2485 runtime = time.time()begin 2486 chisq = np.sum(result[2]['fvec']**2) 2487 Values2Dict(parmDict, varyList, result[0]) 2488 G2mv.Dict2Map(parmDict) 2489 ApplyXYZshifts(parmDict,varyList) 2490 2491 Rwp = np.sqrt(chisq/Histo['sumwYo'])*100. #to % 2492 GOF = chisq/(Histo['Nobs']len(varyList)) 2493 print '\n Refinement results for histogram: v'+histogram 2494 print 135*'' 2495 print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList) 2496 print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc) 2497 print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF) 2498 try: 2499 covMatrix = result[1]*GOF 2500 sig = np.sqrt(np.diag(covMatrix)) 2501 if np.any(np.isnan(sig)): 2502 print '*** Least squares aborted  some invalid esds possible ***' 2503 ifPrint = True 2504 break #refinement succeeded  finish up! 2505 except TypeError: #result[1] is None on singular matrix 2506 print '**** Refinement failed  singular matrix ****' 2507 Ipvt = result[2]['ipvt'] 2508 for i,ipvt in enumerate(Ipvt): 2509 if not np.sum(result[2]['fjac'],axis=1)[i]: 2510 print 'Removing parameter: ',varyList[ipvt1] 2511 del(varyList[ipvt1]) 2512 break 2513 2514 GetFobsSq(Histo,Phases,parmDict,calcControls) 2515 sigDict = dict(zip(varyList,sig)) 2516 covData = {'variables':result[0],'covMatrix':covMatrix} 2517 SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint) 2518 SetHistogramData(parmDict,sigDict,Histo,ifPrint) 2519 SeqResult[histogram] = covData 2520 SetSeqResult(GPXfile,Histograms,SeqResult) 2521 2522 2353 2523 def main(): 2354 2524 arg = sys.argv
Note: See TracChangeset
for help on using the changeset viewer.