Changeset 3000 for trunk/GSASIIscriptable.py
- Timestamp:
- Aug 11, 2017 5:34:54 PM (5 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk
-
Property
svn:mergeinfo
set to
/branch/2frame merged eligible
-
Property
svn:mergeinfo
set to
-
trunk/GSASIIscriptable.py
r2836 r3000 1 #!/usr/bin/env python 1 2 # -*- coding: utf-8 -*- 2 3 ########### SVN repository information ################### … … 8 9 ########### SVN repository information ################### 9 10 """ 10 11 *GSASIIscriptable: Scripting Tools* 12 ----------------------------------- 13 14 Routines for reading, writing, modifying and creating GSAS-II project (.gpx) files. 15 16 Supports a command line interface as well. 11 17 """ 18 from __future__ import division, print_function # needed? 12 19 import os.path as ospath 13 20 import sys … … 15 22 import imp 16 23 import copy 24 import platform 25 import os 26 import random as ran 27 import json 28 29 import numpy.ma as ma 30 import scipy.interpolate as si 31 import numpy as np 32 import scipy as sp 33 import matplotlib as mpl 34 17 35 import GSASIIpath 36 GSASIIpath.SetBinaryPath(False) # would rather have this in __name__ == '__main__' stanza 37 import GSASIIIO as G2IO 38 import GSASIIfiles as G2fil 18 39 import GSASIIobj as G2obj 40 import GSASIIpwd as G2pwd 41 import GSASIIstrIO as G2strIO 42 import GSASIIspc as G2spc 43 import GSASIIstrMain as G2strMain 44 import GSASIImath as G2mth 45 import GSASIIElem as G2elem 46 47 try: 48 import matplotlib.pyplot as plt 49 except ImportError: 50 plt = None 19 51 20 52 def LoadDictFromProjFile(ProjFile): 21 ''''Read a GSAS-II project file and load items to dictionary 22 :param: str ProjFile: GSAS-II project (name.gpx) full file name 23 :returns: dict Project: representation of gpx file following the GSAS-II tree struture 53 '''Read a GSAS-II project file and load items to dictionary 54 :param str ProjFile: GSAS-II project (name.gpx) full file name 55 :returns: Project,nameList, where 56 57 * Project (dict) is a representation of gpx file following the GSAS-II tree struture 24 58 for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict 25 59 data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)} 26 :returns: list nameList: names of main tree entries & subentries used to reconstruct project file 27 Example for fap.gpx: 28 Project = { #NB:dict order is not tree order 29 u'Phases':{'data':None,'fap':{phase dict}}, 30 u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc}, 60 * nameList (list) has names of main tree entries & subentries used to reconstruct project file 61 62 Example for fap.gpx:: 63 64 Project = { #NB:dict order is not tree order 65 u'Phases':{'data':None,'fap':{phase dict}}, 66 u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc}, 31 67 u'Rigid bodies':{'data': {rigid body dict}}, 32 68 u'Covariance':{'data':{covariance data dict}}, 33 u'Controls':{'data':{controls data dict}}, 34 u'Notebook':{'data':[notebook list]}, 35 u'Restraints':{'data':{restraint data dict}}, 36 u'Constraints':{'data':{constraint data dict}}]} 37 nameList = [ #NB: reproduces tree order 69 u'Controls':{'data':{controls data dict}}, 70 u'Notebook':{'data':[notebook list]}, 71 u'Restraints':{'data':{restraint data dict}}, 72 u'Constraints':{'data':{constraint data dict}}] 73 } 74 nameList = [ #NB: reproduces tree order 38 75 [u'Notebook',], 39 [u'Controls',], 40 [u'Covariance',], 76 [u'Controls',], 77 [u'Covariance',], 41 78 [u'Constraints',], 42 [u'Restraints',], 43 [u'Rigid bodies',], 79 [u'Restraints',], 80 [u'Rigid bodies',], 44 81 [u'PWDR FAP.XRA Bank 1', 45 u'Comments', 46 u'Limits', 47 u'Background', 82 u'Comments', 83 u'Limits', 84 u'Background', 48 85 u'Instrument Parameters', 49 u'Sample Parameters', 50 u'Peak List', 51 u'Index Peak List', 52 u'Unit Cells List', 53 u'Reflection Lists'], 54 [u'Phases', 55 u'fap']]86 u'Sample Parameters', 87 u'Peak List', 88 u'Index Peak List', 89 u'Unit Cells List', 90 u'Reflection Lists'], 91 [u'Phases', u'fap'] 92 ] 56 93 ''' 57 if not ospath.exists(ProjFile): 58 print ('\n*** Error attempt to open project file that does not exist:\n '+ 59 str(ProjFile)) 60 return 94 # Let IOError be thrown if file does not exist 95 # if not ospath.exists(ProjFile): 96 # print ('\n*** Error attempt to open project file that does not exist:\n '+ 97 # str(ProjFile)) 98 # return 61 99 file = open(ProjFile,'rb') 62 print 'loading from file: ',ProjFile100 # print('loading from file: {}'.format(ProjFile)) 63 101 Project = {} 64 102 nameList = [] … … 76 114 nameList[-1].append(datus[0]) 77 115 file.close() 78 print('project load successful')116 # print('project load successful') 79 117 except: 80 print("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file") 81 return None 118 raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file") 119 finally: 120 file.close() 82 121 return Project,nameList 83 122 84 123 def SaveDictToProjFile(Project,nameList,ProjFile): 85 124 '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile 86 param: dict Project: representation of gpx file following the GSAS-II 87 tree struture as described for LoadDictFromProjFile 88 param: list nameList: names of main tree entries & subentries used to reconstruct project file 89 param: str ProjFile: full file name for output project.gpx file (including extension) 125 126 :param dict Project: representation of gpx file following the GSAS-II 127 tree structure as described for LoadDictFromProjFile 128 :param list nameList: names of main tree entries & subentries used to reconstruct project file 129 :param str ProjFile: full file name for output project.gpx file (including extension) 90 130 ''' 91 131 file = open(ProjFile,'wb') 92 print 'save to file: ',ProjFile 93 for name in nameList: 94 data = [] 95 item = Project[name[0]] 96 data.append([name[0],item['data']]) 97 for item2 in name[1:]: 98 data.append([item2,item[item2]]) 99 cPickle.dump(data,file,1) 100 file.close() 101 print('project save successful') 102 132 # print('save to file: {}'.format(ProjFile)) 133 try: 134 for name in nameList: 135 data = [] 136 item = Project[name[0]] 137 data.append([name[0],item['data']]) 138 for item2 in name[1:]: 139 data.append([item2,item[item2]]) 140 cPickle.dump(data,file,1) 141 finally: 142 file.close() 143 # print('project save successful') 144 103 145 def ImportPowder(reader,filename): 104 146 '''Use a reader to import a powder diffraction data file 105 param: str reader: a scriptable reader 106 param: str filename: full name of powder data file; can be "multi-Bank" data107 returns: list rdlist: list of rrader objects containing powder data, one for each108 "Bank" of data encountered in file 109 items in reader object of interest are:110 rd.comments: list of str: comments found on powder file111 rd.dnames: list of str: data nammes suitable for use in GSASII data tree 112 NB: duplicated in all rd entries in rdlist113 rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed114 for a PWDR entry in GSASII data tree.147 148 :param str reader: a scriptable reader 149 :param str filename: full name of powder data file; can be "multi-Bank" data 150 151 :returns list rdlist: list of reader objects containing powder data, one for each 152 "Bank" of data encountered in file. Items in reader object of interest are: 153 154 * rd.comments: list of str: comments found on powder file 155 * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist 156 * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in GSASII data tree. 115 157 ''' 116 158 rdfile,rdpath,descr = imp.find_module(reader) 117 159 rdclass = imp.load_module(reader,rdfile,rdpath,descr) 118 rd = rdclass.GSAS_ReaderClass() 160 rd = rdclass.GSAS_ReaderClass() 119 161 if not rd.scriptable: 120 print '**** ERROR: ',reader,' is not a scriptable reader'162 print(u'**** ERROR: '+reader+u' is not a scriptable reader') 121 163 return None 122 164 fl = open(filename,'rb') … … 137 179 repeat = True 138 180 return rdlist 139 print rd.errors181 print(rd.errors) 140 182 return None 141 183 184 def SetDefaultDData(dType,histoName,NShkl=0,NDij=0): 185 '''Create an initial Histogram dictionary 186 187 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 188 ''' 189 if dType in ['SXC','SNC']: 190 return {'Histogram':histoName,'Show':False,'Scale':[1.0,True], 191 'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}, 192 'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955, 193 'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}], 194 'Flack':[0.0,False]} 195 elif dType == 'SNT': 196 return {'Histogram':histoName,'Show':False,'Scale':[1.0,True], 197 'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}, 198 'Extinction':['Lorentzian','None', { 199 'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]} 200 elif 'P' in dType: 201 return {'Histogram':histoName,'Show':False,'Scale':[1.0,False], 202 'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1], 203 'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1], 204 [1.,1.,1.,0.,0.,0.],6*[False,]], 205 'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1], 206 NShkl*[0.01,],NShkl*[False,]], 207 'HStrain':[NDij*[0.0,],NDij*[False,]], 208 'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}} 209 210 211 def PreSetup(data): 212 '''Create part of an initial (empty) phase dictionary 213 214 from GSASIIphsGUI.py, near end of UpdatePhaseData 215 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 216 ''' 217 if 'RBModels' not in data: 218 data['RBModels'] = {} 219 if 'MCSA' not in data: 220 data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}} 221 if 'dict' in str(type(data['MCSA']['Results'])): 222 data['MCSA']['Results'] = [] 223 if 'Modulated' not in data['General']: 224 data['General']['Modulated'] = False 225 if 'modulated' in data['General']['Type']: 226 data['General']['Modulated'] = True 227 data['General']['Type'] = 'nuclear' 228 229 230 def SetupGeneral(data, dirname): 231 """Helps initialize phase data. 232 233 From GSASIIphsGui.py, function of the same name. Minor changes for imports etc. 234 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 235 """ 236 mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True, 237 'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False} 238 generalData = data['General'] 239 atomData = data['Atoms'] 240 generalData['AtomTypes'] = [] 241 generalData['Isotopes'] = {} 242 243 if 'Isotope' not in generalData: 244 generalData['Isotope'] = {} 245 if 'Data plot type' not in generalData: 246 generalData['Data plot type'] = 'Mustrain' 247 if 'POhkl' not in generalData: 248 generalData['POhkl'] = [0,0,1] 249 if 'Map' not in generalData: 250 generalData['Map'] = mapDefault.copy() 251 if 'Flip' not in generalData: 252 generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None', 253 'k-factor':0.1,'k-Max':20.,} 254 if 'testHKL' not in generalData['Flip']: 255 generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]] 256 if 'doPawley' not in generalData: 257 generalData['doPawley'] = False #ToDo: change to '' 258 if 'Pawley dmin' not in generalData: 259 generalData['Pawley dmin'] = 1.0 260 if 'Pawley neg wt' not in generalData: 261 generalData['Pawley neg wt'] = 0.0 262 if 'Algolrithm' in generalData.get('MCSA controls',{}) or \ 263 'MCSA controls' not in generalData: 264 generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50], 265 'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0, 266 'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True} 267 if 'AtomPtrs' not in generalData: 268 generalData['AtomPtrs'] = [3,1,7,9] 269 if generalData['Type'] == 'macromolecular': 270 generalData['AtomPtrs'] = [6,4,10,12] 271 elif generalData['Type'] == 'magnetic': 272 generalData['AtomPtrs'] = [3,1,10,12] 273 if generalData['Type'] in ['modulated',]: 274 generalData['Modulated'] = True 275 generalData['Type'] = 'nuclear' 276 if 'Super' not in generalData: 277 generalData['Super'] = 1 278 generalData['SuperVec'] = [[0,0,.1],False,4] 279 generalData['SSGData'] = {} 280 if '4DmapData' not in generalData: 281 generalData['4DmapData'] = mapDefault.copy() 282 generalData['4DmapData'].update({'MapType':'Fobs'}) 283 if 'Modulated' not in generalData: 284 generalData['Modulated'] = False 285 if 'HydIds' not in generalData: 286 generalData['HydIds'] = {} 287 cx,ct,cs,cia = generalData['AtomPtrs'] 288 generalData['NoAtoms'] = {} 289 generalData['BondRadii'] = [] 290 generalData['AngleRadii'] = [] 291 generalData['vdWRadii'] = [] 292 generalData['AtomMass'] = [] 293 generalData['Color'] = [] 294 if generalData['Type'] == 'magnetic': 295 generalData['MagDmin'] = generalData.get('MagDmin',1.0) 296 landeg = generalData.get('Lande g',[]) 297 generalData['Mydir'] = dirname 298 badList = {} 299 for iat,atom in enumerate(atomData): 300 atom[ct] = atom[ct].lower().capitalize() #force to standard form 301 if generalData['AtomTypes'].count(atom[ct]): 302 generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1]) 303 elif atom[ct] != 'UNK': 304 Info = G2elem.GetAtomInfo(atom[ct]) 305 if not Info: 306 if atom[ct] not in badList: 307 badList[atom[ct]] = 0 308 badList[atom[ct]] += 1 309 atom[ct] = 'UNK' 310 continue 311 atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo 312 generalData['AtomTypes'].append(atom[ct]) 313 generalData['Z'] = Info['Z'] 314 generalData['Isotopes'][atom[ct]] = Info['Isotopes'] 315 generalData['BondRadii'].append(Info['Drad']) 316 generalData['AngleRadii'].append(Info['Arad']) 317 generalData['vdWRadii'].append(Info['Vdrad']) 318 if atom[ct] in generalData['Isotope']: 319 if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]: 320 isotope = generalData['Isotopes'][atom[ct]].keys()[-1] 321 generalData['Isotope'][atom[ct]] = isotope 322 generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass']) 323 else: 324 generalData['Isotope'][atom[ct]] = 'Nat. Abund.' 325 if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]: 326 isotope = generalData['Isotopes'][atom[ct]].keys()[-1] 327 generalData['Isotope'][atom[ct]] = isotope 328 generalData['AtomMass'].append(Info['Mass']) 329 generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1]) 330 generalData['Color'].append(Info['Color']) 331 if generalData['Type'] == 'magnetic': 332 if len(landeg) < len(generalData['AtomTypes']): 333 landeg.append(2.0) 334 if generalData['Type'] == 'magnetic': 335 generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])] 336 337 if badList: 338 msg = 'Warning: element symbol(s) not found:' 339 for key in badList: 340 msg += '\n\t' + key 341 if badList[key] > 1: 342 msg += ' (' + str(badList[key]) + ' times)' 343 raise Exception("Phase error:\n" + msg) 344 # wx.MessageBox(msg,caption='Element symbol error') 345 F000X = 0. 346 F000N = 0. 347 for i,elem in enumerate(generalData['AtomTypes']): 348 F000X += generalData['NoAtoms'][elem]*generalData['Z'] 349 isotope = generalData['Isotope'][elem] 350 F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0] 351 generalData['F000X'] = F000X 352 generalData['F000N'] = F000N 353 generalData['Mass'] = G2mth.getMass(generalData) 354 355 def make_empty_project(author=None, filename=None): 356 """Creates an dictionary in the style of GSASIIscriptable, for an empty 357 project. 358 359 If no author name or filename is supplied, 'no name' and 360 <current dir>/test_output.gpx are used , respectively. 361 362 Returns: project dictionary, name list 363 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 364 """ 365 if not filename: 366 filename = os.path.join(os.getcwd(), 'test_output.gpx') 367 else: 368 filename = os.path.abspath(filename) 369 gsasii_version = str(GSASIIpath.GetVersionNumber()) 370 python_library_versions = G2fil.get_python_versions([mpl, np, sp]) 371 372 controls_data = dict(G2obj.DefaultControls) 373 controls_data['LastSavedAs'] = unicode(filename) 374 controls_data['LastSavedUsing'] = gsasii_version 375 controls_data['PythonVersions'] = python_library_versions 376 if author: 377 controls_data['Author'] = author 378 379 output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [], 380 'Global': []}}, 381 'Controls': {'data': controls_data}, 382 u'Covariance': {'data': {}}, 383 u'Notebook': {'data': ['']}, 384 u'Restraints': {'data': {}}, 385 u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []}, 386 'Residue': {'AtInfo': {}}, 387 'Vector': {'AtInfo': {}}}}} 388 389 names = [[u'Notebook'], [u'Controls'], [u'Covariance'], 390 [u'Constraints'], [u'Restraints'], [u'Rigid bodies']] 391 392 return output, names 393 394 395 class G2ImportException(Exception): 396 pass 397 398 399 def import_generic(filename, readerlist): 400 """Attempt to import a filename, using a list of reader objects. 401 402 Returns the first reader object which worked.""" 403 # Translated from OnImportGeneric method in GSASII.py 404 primaryReaders, secondaryReaders = [], [] 405 for reader in readerlist: 406 flag = reader.ExtensionValidator(filename) 407 if flag is None: 408 secondaryReaders.append(reader) 409 elif flag: 410 primaryReaders.append(reader) 411 if not secondaryReaders and not primaryReaders: 412 raise G2ImportException("Could not read file: ", filename) 413 414 with open(filename, 'Ur') as fp: 415 rd_list = [] 416 417 for rd in primaryReaders + secondaryReaders: 418 # Initialize reader 419 rd.selections = [] 420 rd.dnames = [] 421 rd.ReInitialize() 422 # Rewind file 423 fp.seek(0) 424 rd.errors = "" 425 if not rd.ContentsValidator(fp): 426 # Report error 427 pass 428 if len(rd.selections) > 1: 429 # Select data? 430 # GSASII.py:543 431 raise G2ImportException("Not sure what data to select") 432 433 block = 0 434 rdbuffer = {} 435 fp.seek(0) 436 repeat = True 437 while repeat: 438 repeat = False 439 block += 1 440 rd.objname = os.path.basename(filename) 441 flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block) 442 if flag: 443 # Omitting image loading special cases 444 rd.readfilename = filename 445 rd_list.append(copy.deepcopy(rd)) 446 repeat = rd.repeat 447 else: 448 raise G2ImportException("{}.Reader() returned:".format(rd), 449 flag) 450 451 if rd_list: 452 if rd.warnings: 453 print("Read warning by", rd.formatName, "reader:", 454 rd.warnings, file=sys.stderr) 455 return rd_list 456 raise G2ImportException("No reader could read file: " + filename) 457 458 459 def load_iprms(instfile, reader): 460 """Loads instrument parameters from a file, and edits the 461 given reader. 462 463 Returns a 2-tuple of (Iparm1, Iparm2) parameters 464 """ 465 ext = os.path.splitext(instfile)[1] 466 467 if ext.lower() == '.instprm': 468 # New GSAS File, load appropriately 469 with open(instfile) as f: 470 lines = f.readlines() 471 bank = reader.powderentry[2] 472 numbanks = reader.numbanks 473 iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader) 474 475 reader.instfile = instfile 476 reader.instmsg = 'GSAS-II file' + instfile 477 return iparms 478 elif ext.lower() not in ('.prm', '.inst', '.ins'): 479 raise ValueError('Expected .prm file, found: ', instfile) 480 481 # It's an old GSAS file, load appropriately 482 Iparm = {} 483 with open(instfile, 'Ur') as fp: 484 for line in fp: 485 if '#' in line: 486 continue 487 Iparm[line[:12]] = line[12:-1] 488 ibanks = int(Iparm.get('INS BANK ', '1').strip()) 489 if ibanks == 1: 490 reader.instbank = 1 491 reader.powderentry[2] = 1 492 reader.instfile = instfile 493 reader.instmsg = instfile + ' bank ' + str(reader.instbank) 494 return G2fil.SetPowderInstParms(Iparm, reader) 495 # TODO handle >1 banks 496 raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm") 497 498 def load_pwd_from_reader(reader, instprm, existingnames=[]): 499 """Loads powder data from a reader object, and assembles it into a GSASII data tree. 500 501 Returns: (name, tree) - 2-tuple of the histogram name (str), and data 502 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 503 """ 504 HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_') 505 HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames)) 506 507 try: 508 Iparm1, Iparm2 = instprm 509 except ValueError: 510 Iparm1, Iparm2 = load_iprms(instprm, reader) 511 512 Ymin = np.min(reader.powderdata[1]) 513 Ymax = np.max(reader.powderdata[1]) 514 valuesdict = {'wtFactor': 1.0, 515 'Dummy': False, 516 'ranId': ran.randint(0, sys.maxint), 517 'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax, 518 'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax, 519 'qPlot': False, 'dPlot': False, 'sqrtPlot': False, 520 'Yminmax': [Ymin, Ymax]} 521 reader.Sample['ranId'] = valuesdict['ranId'] 522 523 # Ending keys: 524 # [u'Reflection Lists', 525 # u'Limits', 526 # 'data', 527 # u'Index Peak List', 528 # u'Comments', 529 # u'Unit Cells List', 530 # u'Sample Parameters', 531 # u'Peak List', 532 # u'Background', 533 # u'Instrument Parameters'] 534 Tmin = np.min(reader.powderdata[0]) 535 Tmax = np.max(reader.powderdata[0]) 536 537 default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0], 538 {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}] 539 540 output_dict = {u'Reflection Lists': {}, 541 u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]), 542 u'data': [valuesdict, reader.powderdata, HistName], 543 u'Index Peak List': [[], []], 544 u'Comments': reader.comments, 545 u'Unit Cells List': [], 546 u'Sample Parameters': reader.Sample, 547 u'Peak List': {'peaks': [], 'sigDict': {}}, 548 u'Background': reader.pwdparms.get('Background', default_background), 549 u'Instrument Parameters': [Iparm1, Iparm2], 550 } 551 552 names = [u'Comments', 553 u'Limits', 554 u'Background', 555 u'Instrument Parameters', 556 u'Sample Parameters', 557 u'Peak List', 558 u'Index Peak List', 559 u'Unit Cells List', 560 u'Reflection Lists'] 561 562 # TODO controls?? GSASII.py:1664-7 563 564 return HistName, [HistName] + names, output_dict 565 566 567 def _deep_copy_into(from_, into): 568 """Helper function for reloading .gpx file. See G2Project.reload() 569 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 570 """ 571 if isinstance(from_, dict) and isinstance(into, dict): 572 combined_keys = set(from_.keys()).union(into.keys()) 573 for key in combined_keys: 574 if key in from_ and key in into: 575 both_dicts = (isinstance(from_[key], dict) 576 and isinstance(into[key], dict)) 577 both_lists = (isinstance(from_[key], list) 578 and isinstance(into[key], list)) 579 if both_dicts or both_lists: 580 _deep_copy_into(from_[key], into[key]) 581 else: 582 into[key] = from_[key] 583 elif key in from_: 584 into[key] = from_[key] 585 else: # key in into 586 del into[key] 587 elif isinstance(from_, list) and isinstance(into, list): 588 if len(from_) == len(into): 589 for i in xrange(len(from_)): 590 both_dicts = (isinstance(from_[i], dict) 591 and isinstance(into[i], dict)) 592 both_lists = (isinstance(from_[i], list) 593 and isinstance(into[i], list)) 594 if both_dicts or both_lists: 595 _deep_copy_into(from_[i], into[i]) 596 else: 597 into[i] = from_[i] 598 else: 599 into[:] = from_ 600 601 602 class G2ObjectWrapper(object): 603 """Base class for all GSAS-II object wrappers. 604 605 The underlying GSAS-II format can be accessed as `wrapper.data`. A number 606 of overrides are implemented so that the wrapper behaves like a dictionary. 607 608 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 609 """ 610 def __init__(self, datadict): 611 self.data = datadict 612 613 def __getitem__(self, key): 614 return self.data[key] 615 616 def __setitem__(self, key, value): 617 self.data[key] = value 618 619 def __contains__(self, key): 620 return key in self.data 621 622 def get(self, k, d=None): 623 return self.data.get(k, d) 624 625 def keys(self): 626 return self.data.keys() 627 628 def values(self): 629 return self.data.values() 630 631 def items(self): 632 return self.data.items() 633 634 635 class G2Project(G2ObjectWrapper): 636 def __init__(self, gpxfile=None, author=None, filename=None): 637 """Loads a GSAS-II project from a specified filename. 638 639 :param str gpxfile: Existing .gpx file to be loaded. If nonexistent, 640 creates an empty project. 641 :param str author: Author's name. Optional. 642 :param str filename: The filename the project should be saved to in 643 the future. If both filename and gpxfile are present, the project is 644 loaded from the gpxfile, then set to save to filename in the future""" 645 if gpxfile is None: 646 filename = os.path.abspath(os.path.expanduser(filename)) 647 self.filename = filename 648 self.data, self.names = make_empty_project(author=author, filename=filename) 649 elif isinstance(gpxfile, str): 650 # TODO set author, filename 651 self.filename = os.path.abspath(os.path.expanduser(gpxfile)) 652 self.data, self.names = LoadDictFromProjFile(gpxfile) 653 self.index_ids() 654 else: 655 raise ValueError("Not sure what to do with gpxfile") 656 657 @classmethod 658 def from_dict_and_names(cls, gpxdict, names, filename=None): 659 """Creates a :class:`G2Project` directly from 660 a dictionary and a list of names. If in doubt, do not use this. 661 662 :returns: a :class:`G2Project` 663 """ 664 out = cls() 665 if filename: 666 filename = os.path.abspath(os.path.expanduser(filename)) 667 out.filename = filename 668 gpxdict['Controls']['data']['LastSavedAs'] = filename 669 else: 670 try: 671 out.filename = gpxdict['Controls']['data']['LastSavedAs'] 672 except KeyError: 673 out.filename = None 674 out.data = gpxdict 675 out.names = names 676 677 def save(self, filename=None): 678 """Saves the project, either to the current filename, or to a new file. 679 680 Updates self.filename if a new filename provided""" 681 # TODO update LastSavedUsing ? 682 if filename: 683 filename = os.path.abspath(os.path.expanduser(filename)) 684 self.data['Controls']['data']['LastSavedAs'] = filename 685 self.filename = filename 686 elif not self.filename: 687 raise AttributeError("No file name to save to") 688 SaveDictToProjFile(self.data, self.names, self.filename) 689 690 def add_powder_histogram(self, datafile, iparams): 691 """Loads a powder data histogram into the project. 692 693 Automatically checks for an instrument parameter file, or one can be 694 provided. 695 696 :param str datafile: The powder data file to read, a filename. 697 :param str iparams: The instrument parameters file, a filename. 698 699 :returns: A :class:`G2PwdrData` object representing 700 the histogram""" 701 datafile = os.path.abspath(os.path.expanduser(datafile)) 702 iparams = os.path.abspath(os.path.expanduser(iparams)) 703 pwdrreaders = import_generic(datafile, PwdrDataReaders) 704 histname, new_names, pwdrdata = load_pwd_from_reader( 705 pwdrreaders[0], iparams, 706 list(self.histogram_names())) 707 if histname in self.data: 708 print("Warning - redefining histogram", histname) 709 else: 710 if self.names[-1][0] == 'Phases': 711 self.names.insert(-1, new_names) 712 else: 713 self.names.append(new_names) 714 self.data[histname] = pwdrdata 715 return self.histogram(histname) 716 717 def add_phase(self, phasefile, phasename=None, histograms=[]): 718 """Loads a phase into the project from a .cif file 719 720 :param str phasefile: The CIF file from which to import the phase. 721 :param str phasename: The name of the new phase, or None for the default 722 :param list histograms: The names of the histograms to associate with 723 this phase 724 725 :returns: A :class:`G2Phase` object representing the 726 new phase.""" 727 phasefile = os.path.abspath(os.path.expanduser(phasefile)) 728 729 # TODO handle multiple phases in a file 730 phasereaders = import_generic(phasefile, PhaseReaders) 731 phasereader = phasereaders[0] 732 733 phasename = phasename or phasereader.Phase['General']['Name'] 734 phaseNameList = list(self.phase_names()) 735 phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList) 736 phasereader.Phase['General']['Name'] = phasename 737 738 if 'Phases' not in self.data: 739 self.data[u'Phases'] = { 'data': None } 740 assert phasename not in self.data['Phases'], "phase names should be unique" 741 self.data['Phases'][phasename] = phasereader.Phase 742 743 if phasereader.Constraints: 744 Constraints = self.data['Constraints'] 745 for i in phasereader.Constraints: 746 if isinstance(i, dict): 747 if '_Explain' not in Constraints: 748 Constraints['_Explain'] = {} 749 Constraints['_Explain'].update(i) 750 else: 751 Constraints['Phase'].append(i) 752 753 data = self.data['Phases'][phasename] 754 generalData = data['General'] 755 SGData = generalData['SGData'] 756 NShkl = len(G2spc.MustrainNames(SGData)) 757 NDij = len(G2spc.HStrainNames(SGData)) 758 Super = generalData.get('Super', 0) 759 if Super: 760 SuperVec = np.array(generalData['SuperVec'][0]) 761 else: 762 SuperVec = [] 763 UseList = data['Histograms'] 764 765 for histname in histograms: 766 if histname.startswith('HKLF '): 767 raise NotImplementedError("Does not support HKLF yet") 768 elif histname.startswith('PWDR '): 769 hist = self.histogram(histname) 770 hist['Reflection Lists'][generalData['Name']] = {} 771 UseList[histname] = SetDefaultDData('PWDR', histname, NShkl=NShkl, NDij=NDij) 772 for key, val in [('Use', True), ('LeBail', False), 773 ('newLeBail', True), 774 ('Babinet', {'BabA': [0.0, False], 775 'BabU': [0.0, False]})]: 776 if key not in UseList[histname]: 777 UseList[histname][key] = val 778 else: 779 raise NotImplementedError("Unexpected histogram" + histname) 780 781 for obj in self.names: 782 if obj[0] == 'Phases': 783 phasenames = obj 784 break 785 else: 786 phasenames = [u'Phases'] 787 self.names.append(phasenames) 788 phasenames.append(unicode(phasename)) 789 790 # TODO should it be self.filename, not phasefile? 791 SetupGeneral(data, os.path.dirname(phasefile)) 792 self.index_ids() 793 794 return self.phase(phasename) 795 796 def reload(self): 797 """Reload self from self.filename""" 798 data, names = LoadDictFromProjFile(self.filename) 799 self.names = names 800 # Need to deep copy the new data file data into the current tree, 801 # so that any existing G2Phase, or G2PwdrData objects will still be 802 # valid 803 _deep_copy_into(from_=data, into=self.data) 804 805 def refine(self, newfile=None, printFile=None): 806 # index_ids will automatically save the project 807 self.index_ids() 808 # TODO G2strMain does not properly use printFile 809 G2strMain.Refine(self.filename) 810 # Reload yourself 811 self.reload() 812 813 def histogram_names(self): 814 """Gives a list of the names of each histogram in the project. 815 816 :returns: list of strings 817 818 .. seealso:: 819 :func:`~GSASIIscriptable.G2Project.histogram` 820 :func:`~GSASIIscriptable.G2Project.histograms` 821 :func:`~GSASIIscriptable.G2Project.phase` 822 :func:`~GSASIIscriptable.G2Project.phases` 823 :func:`~GSASIIscriptable.G2Project.phase_names` 824 """ 825 output = [] 826 for obj in self.names: 827 if len(obj) > 1 and obj[0] != u'Phases': 828 output.append(obj[0]) 829 return output 830 831 def histogram(self, histname): 832 """Returns the histogram named histname, or None if it does not exist. 833 834 :returns: A :class:`G2PwdrData` object, or None if 835 the histogram does not exist 836 .. seealso:: 837 :func:`~GSASIIscriptable.G2Project.histogram_names` 838 :func:`~GSASIIscriptable.G2Project.histograms` 839 :func:`~GSASIIscriptable.G2Project.phase` 840 :func:`~GSASIIscriptable.G2Project.phases` 841 :func:`~GSASIIscriptable.G2Project.phase_names` 842 """ 843 if histname in self.data: 844 return G2PwdrData(self.data[histname], self) 845 return None 846 847 def histograms(self): 848 """Return a list of all histograms, as 849 :class:`G2PwdrData` objects 850 .. seealso:: 851 :func:`~GSASIIscriptable.G2Project.histogram_names` 852 :func:`~GSASIIscriptable.G2Project.histograms` 853 :func:`~GSASIIscriptable.G2Project.phase` 854 :func:`~GSASIIscriptable.G2Project.phases` 855 :func:`~GSASIIscriptable.G2Project.phase_names` 856 """ 857 return [self.histogram(name) for name in self.histogram_names()] 858 859 def phase_names(self): 860 """Gives an list of the names of each phase in the project. 861 862 :returns: A list of strings 863 .. seealso:: 864 :func:`~GSASIIscriptable.G2Project.histogram` 865 :func:`~GSASIIscriptable.G2Project.histograms` 866 :func:`~GSASIIscriptable.G2Project.histogram_names` 867 :func:`~GSASIIscriptable.G2Project.phase` 868 :func:`~GSASIIscriptable.G2Project.phases` 869 """ 870 for obj in self.names: 871 if obj[0] == 'Phases': 872 return obj[1:] 873 return [] 874 875 def phase(self, phasename): 876 """ 877 Gives an object representing the specified phase in this project. 878 879 :param str phasename: The name of the desired phase 880 :returns: A :class:`G2Phase` object 881 :raises: KeyError 882 .. seealso:: 883 :func:`~GSASIIscriptable.G2Project.histogram_names` 884 :func:`~GSASIIscriptable.G2Project.histograms` 885 :func:`~GSASIIscriptable.G2Project.phase` 886 :func:`~GSASIIscriptable.G2Project.phases` 887 :func:`~GSASIIscriptable.G2Project.phase_names` 888 """ 889 phases = self.data['Phases'] 890 if phasename in phases: 891 return G2Phase(phases[phasename], phasename, self) 892 893 def phases(self): 894 """ 895 Returns a list of all the phases in the project. 896 897 :returns: A :class:`G2Phase` 898 .. seealso:: 899 :func:`~GSASIIscriptable.G2Project.histogram` 900 :func:`~GSASIIscriptable.G2Project.histograms` 901 :func:`~GSASIIscriptable.G2Project.histogram_names` 902 :func:`~GSASIIscriptable.G2Project.phase` 903 :func:`~GSASIIscriptable.G2Project.phase_names` 904 """ 905 return [self.phase(name) for name in self.phase_names()] 906 907 def do_refinements(self, refinements, histogram='all', phase='all', 908 outputnames=None): 909 """Conducts a series of refinements. 910 911 :param list refinements: A list of dictionaries defining refinements 912 :param str histogram: Name of histogram for refinements to be applied 913 to, or 'all' 914 :param str phase: Name of phase for refinements to be applied to, or 915 'all' 916 """ 917 if outputnames: 918 if len(refinements) != len(outputnames): 919 raise ValueError("Should have same number of outuputs to" 920 "refinements") 921 else: 922 outputnames = [None for r in refinements] 923 924 for output, refinement in zip(outputnames, refinements): 925 self.set_refinement(refinement, histogram) 926 # Handle 'once' args - refinements that are disabled after this 927 # refinement 928 if 'once' in refinement: 929 temp = {'set': refinement['once']} 930 self.set_refinement(temp, histogram, phase) 931 932 if output: 933 self.save(output) 934 935 self.refine() # newFile=output) 936 937 # Handle 'once' args - refinements that are disabled after this 938 # refinement 939 if 'once' in refinement: 940 temp = {'clear': refinement['once']} 941 self.set_refinement(temp, histogram, phase) 942 943 def set_refinement(self, refinement, histogram='all', phase='all'): 944 """Apply specified refinements to a given histogram(s) or phase(s). 945 946 :param dict refinement: The refinements to be conducted 947 :param histogram: Either a name of a histogram (str), a list of 948 histogram names, or 'all' (default) 949 :param phase: Either a name of a phase (str), a list of phase names, or 950 'all' (default)""" 951 if histogram == 'all': 952 hists = [self.histogram(name) 953 for name in self.histogram_names()] 954 elif isinstance(histogram, str): 955 hists = [self.histogram(histogram)] 956 else: 957 hists = [self.histogram(name) for name in histogram] 958 959 if phase == 'all': 960 phases = [self.phase(name) for name in self.phase_names()] 961 elif isinstance(phase, str): 962 phases = [self.phase(phase)] 963 else: 964 phases = [self.phase(name) for name in phase] 965 966 967 # TODO: HAP parameters: 968 # Babinet 969 # Extinction 970 # HStrain 971 # Mustrain 972 # Pref. Ori 973 # Size 974 975 pwdr_set = {} 976 phase_set = {} 977 for key, val in refinement.get('set', {}).items(): 978 # Apply refinement options 979 if G2PwdrData.is_valid_refinement_key(key): 980 pwdr_set[key] = val 981 elif G2Phase.is_valid_refinement_key(key): 982 phase_set[key] = val 983 else: 984 print("Unknown refinement key:", key) 985 986 for hist in hists: 987 hist.set_refinements(pwdr_set) 988 for phase in phases: 989 phase.set_refinements(phase_set) 990 991 pwdr_clear = {} 992 phase_clear = {} 993 for key, val in refinement.get('clear', {}).items(): 994 # Clear refinement options 995 if G2PwdrData.is_valid_refinement_key(key): 996 pwdr_clear[key] = val 997 elif G2Phase.is_valid_refinement_key(key): 998 phase_clear[key] = val 999 else: 1000 print("Unknown refinement key:", key) 1001 1002 for hist in hists: 1003 hist.clear_refinements(pwdr_clear) 1004 for phase in phases: 1005 phase.clear_refinements(phase_clear) 1006 1007 def index_ids(self): 1008 self.save() 1009 return G2strIO.GetUsedHistogramsAndPhases(self.filename) 1010 1011 def add_constraint_raw(self, cons_scope, constr): 1012 """Adds a constraint of type consType to the project. 1013 cons_scope should be one of "Hist", "Phase", "HAP", or "Global". 1014 1015 WARNING it does not check the constraint is well-constructed""" 1016 constrs = self['Constraints']['data'] 1017 if 'Global' not in constrs: 1018 constrs['Global'] = [] 1019 constrs[cons_scope].append(constr) 1020 1021 def hold_many(self, vars, type): 1022 """Apply holds for all the variables in vars, for constraint of a given type. 1023 1024 type is passed directly to add_constraint_raw as consType 1025 1026 :param list vars: A list of variables to hold. Either G2obj.G2VarObj objects, 1027 string variable specifiers, or arguments for :meth:`make_var_obj` 1028 :param str type: A string constraint type specifier. See 1029 :class:`~GSASIIscriptable.G2Project.add_constraint_raw`""" 1030 for var in vars: 1031 if isinstance(var, str): 1032 var = self.make_var_obj(var) 1033 elif not isinstance(var, G2obj.G2VarObj): 1034 var = self.make_var_obj(*var) 1035 self.add_constraint_raw(type, [[1.0, var], None, None, 'h']) 1036 1037 def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None, 1038 reloadIdx=True): 1039 """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a") 1040 or individual names of phase, histogram, varname, and atomId. 1041 1042 Automatically converts string phase, hist, or atom names into the ID required 1043 by G2VarObj.""" 1044 if reloadIdx: 1045 self.index_ids() 1046 1047 # If string representation, short circuit 1048 if hist is None and varname is None and atomId is None: 1049 if isinstance(phase, str) and ':' in phase: 1050 return G2obj.G2VarObj(phase) 1051 1052 # Get phase index 1053 phaseObj = None 1054 if isinstance(phase, G2Phase): 1055 phaseObj = phase 1056 phase = G2obj.PhaseRanIdLookup[phase.ranId] 1057 elif phase in self['Phases']: 1058 phaseObj = self.phase(phase) 1059 phaseRanId = phaseObj.ranId 1060 phase = G2obj.PhaseRanIdLookup[phaseRanId] 1061 1062 # Get histogram index 1063 if isinstance(hist, G2PwdrData): 1064 hist = G2obj.HistRanIdLookup[hist.ranId] 1065 elif hist in self.data: 1066 histRanId = self.histogram(hist).ranId 1067 hist = G2obj.HistRanIdLookup[histRanId] 1068 1069 # Get atom index (if any) 1070 if isinstance(atomId, G2AtomRecord): 1071 atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId] 1072 elif phaseObj: 1073 atomObj = phaseObj.atom(atomId) 1074 if atomObj: 1075 atomRanId = atomObj.ranId 1076 atomId = G2obj.AtomRanIdLookup[phase][atomRanId] 1077 1078 return G2obj.G2VarObj(phase, hist, varname, atomId) 1079 1080 1081 class G2AtomRecord(G2ObjectWrapper): 1082 """Wrapper for an atom record. Has convenient accessors via @property. 1083 1084 Example: 1085 1086 >>> atom = some_phase.atom("O3") 1087 >>> # We can access the underlying data format 1088 >>> atom.data 1089 ['O3', 'O-2', '', ... ] 1090 >>> # We can also use wrapper accessors 1091 >>> atom.coordinates 1092 (0.33, 0.15, 0.5) 1093 >>> atom.refinement_flags 1094 u'FX' 1095 >>> atom.ranId 1096 4615973324315876477 1097 """ 1098 def __init__(self, data, indices, proj): 1099 self.data = data 1100 self.cx, self.ct, self.cs, self.cia = indices 1101 self.proj = proj 1102 1103 @property 1104 def label(self): 1105 return self.data[self.ct-1] 1106 1107 @property 1108 def type(self): 1109 return self.data[self.ct] 1110 1111 @property 1112 def refinement_flags(self): 1113 return self.data[self.ct+1] 1114 1115 @refinement_flags.setter 1116 def refinement_flags(self, other): 1117 # Automatically check it is a valid refinement 1118 for c in other: 1119 if c not in ' FXU': 1120 raise ValueError("Invalid atom refinement: ", other) 1121 self.data[self.ct+1] = unicode(other) 1122 1123 @property 1124 def coordinates(self): 1125 return tuple(self.data[self.cx:self.cx+3]) 1126 1127 @property 1128 def ranId(self): 1129 return self.data[self.cia+8] 1130 1131 @property 1132 def adp_flag(self): 1133 # Either 'I' or 'A' 1134 return self.data[self.cia] 1135 1136 @property 1137 def uiso(self): 1138 if self.adp_flag == 'I': 1139 return self.data[self.cia+1] 1140 else: 1141 return self.data[self.cia+2:self.cia+8] 1142 1143 1144 class G2PwdrData(G2ObjectWrapper): 1145 """Wraps a Poweder Data Histogram.""" 1146 def __init__(self, data, proj): 1147 self.data = data 1148 self.proj = proj 1149 1150 @staticmethod 1151 def is_valid_refinement_key(key): 1152 valid_keys = ['Limits', 'Sample Parameters', 'Background', 1153 'Instrument Parameters'] 1154 return key in valid_keys 1155 1156 @property 1157 def name(self): 1158 return self['data'][-1] 1159 1160 @property 1161 def ranId(self): 1162 return self['data'][0]['ranId'] 1163 1164 @property 1165 def residuals(self): 1166 data = self['data'][0] 1167 return {key: data[key] 1168 for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']} 1169 1170 def fit_fixed_points(self): 1171 """Attempts to apply a background fit to the fixed points currently specified.""" 1172 1173 def SetInstParms(Inst): 1174 dataType = Inst['Type'][0] 1175 insVary = [] 1176 insNames = [] 1177 insVals = [] 1178 for parm in Inst: 1179 insNames.append(parm) 1180 insVals.append(Inst[parm][1]) 1181 if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha', 1182 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]: 1183 Inst[parm][2] = False 1184 instDict = dict(zip(insNames, insVals)) 1185 instDict['X'] = max(instDict['X'], 0.01) 1186 instDict['Y'] = max(instDict['Y'], 0.01) 1187 if 'SH/L' in instDict: 1188 instDict['SH/L'] = max(instDict['SH/L'], 0.002) 1189 return dataType, instDict, insVary 1190 1191 bgrnd = self['Background'] 1192 1193 # Need our fixed points in order 1194 bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0]) 1195 X = [x for x, y in bgrnd[1]['FixedPoints']] 1196 Y = [y for x, y in bgrnd[1]['FixedPoints']] 1197 1198 limits = self['Limits'][1] 1199 if X[0] > limits[0]: 1200 X = [limits[0]] + X 1201 Y = [Y[0]] + Y 1202 if X[-1] < limits[1]: 1203 X += [limits[1]] 1204 Y += [Y[-1]] 1205 1206 # Some simple lookups 1207 controls = self.proj['Controls']['data'] 1208 inst, inst2 = self['Instrument Parameters'] 1209 pwddata = self['data'][1] 1210 1211 # Construct the data for background fitting 1212 xBeg = np.searchsorted(pwddata[0], limits[0]) 1213 xFin = np.searchsorted(pwddata[0], limits[1]) 1214 xdata = pwddata[0][xBeg:xFin] 1215 ydata = si.interp1d(X,Y)(ma.getdata(xdata)) 1216 1217 W = [1]*len(xdata) 1218 Z = [0]*len(xdata) 1219 1220 dataType, insDict, insVary = SetInstParms(inst) 1221 bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd) 1222 1223 # Do the fit 1224 data = np.array([xdata, ydata, W, Z, Z, Z]) 1225 G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data, 1226 prevVaryList=bakVary, controls=controls) 1227 1228 # Post-fit 1229 parmDict = {} 1230 bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd) 1231 parmDict.update(bakDict) 1232 parmDict.update(insDict) 1233 pwddata[3][xBeg:xFin] *= 0 1234 pwddata[5][xBeg:xFin] *= 0 1235 pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0] 1236 1237 # TODO adjust pwddata? GSASIIpwdGUI.py:1041 1238 # TODO update background 1239 self.proj.save() 1240 1241 def y_calc(self): 1242 return self['data'][1][3] 1243 1244 def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True): 1245 if plt: 1246 data = self['data'][1] 1247 if Yobs: 1248 plt.plot(data[0], data[1], label='Yobs') 1249 if Ycalc: 1250 plt.plot(data[0], data[3], label='Ycalc') 1251 if Background: 1252 plt.plot(data[0], data[4], label='Background') 1253 if Residual: 1254 plt.plot(data[0], data[5], label="Residual") 1255 1256 def set_refinements(self, refs): 1257 """Sets the refinement parameter 'key' to the specification 'value' 1258 1259 :param dict refs: A dictionary of the parameters to be set. 1260 1261 :returns: None""" 1262 do_fit_fixed_points = False 1263 for key, value in refs.items(): 1264 if key == 'Limits': 1265 old_limits = self.data['Limits'][1] 1266 new_limits = value 1267 if isinstance(new_limits, dict): 1268 if 'low' in new_limits: 1269 old_limits[0] = new_limits['low'] 1270 if 'high' in new_limits: 1271 old_limits[1] = new_limits['high'] 1272 else: 1273 old_limits[0], old_limits[1] = new_limits 1274 elif key == 'Sample Parameters': 1275 sample = self.data['Sample Parameters'] 1276 for sparam in value: 1277 sample[sparam][1] = True 1278 elif key == 'Background': 1279 bkg, peaks = self.data['Background'] 1280 1281 # If True or False, just set the refine parameter 1282 if value in (True, False): 1283 bkg[1] = value 1284 return 1285 1286 if 'type' in value: 1287 bkg[0] = value['type'] 1288 if 'refine' in value: 1289 bkg[1] = value['refine'] 1290 if 'no. coeffs' in value: 1291 cur_coeffs = bkg[2] 1292 n_coeffs = value['no. coeffs'] 1293 if n_coeffs > cur_coeffs: 1294 for x in range(n_coeffs - cur_coeffs): 1295 bkg.append(0.0) 1296 else: 1297 for _ in range(cur_coeffs - n_coeffs): 1298 bkg.pop() 1299 bkg[2] = n_coeffs 1300 if 'coeffs' in value: 1301 bkg[3:] = value['coeffs'] 1302 if 'FixedPoints' in value: 1303 peaks['FixedPoints'] = [(float(a), float(b)) 1304 for a, b in value['FixedPoints']] 1305 if value.get('fit fixed points', False): 1306 do_fit_fixed_points = True 1307 1308 elif key == 'Instrument Parameters': 1309 instrument, secondary = self.data['Instrument Parameters'] 1310 for iparam in value: 1311 try: 1312 instrument[iparam][2] = True 1313 except IndexError: 1314 raise ValueError("Invalid key:", iparam) 1315 else: 1316 raise ValueError("Unknown key:", key) 1317 # Fit fixed points after the fact - ensure they are after fixed points 1318 # are added 1319 if do_fit_fixed_points: 1320 # Background won't be fit if refinement flag not set 1321 orig = self['Background'][0][1] 1322 self['Background'][0][1] = True 1323 self.fit_fixed_points() 1324 # Restore the previous value 1325 self['Background'][0][1] = orig 1326 1327 def clear_refinements(self, refs): 1328 """Clears the refinement parameter 'key' and its associated value.""" 1329 for key, value in refs.items(): 1330 if key == 'Limits': 1331 old_limits, cur_limits = self.data['Limits'] 1332 cur_limits[0], cur_limits[1] = old_limits 1333 elif key == 'Sample Parameters': 1334 sample = self.data['Sample Parameters'] 1335 for sparam in value: 1336 sample[sparam][1] = False 1337 elif key == 'Background': 1338 bkg, peaks = self.data['Background'] 1339 1340 # If True or False, just set the refine parameter 1341 if value in (True, False): 1342 bkg[1] = False 1343 return 1344 1345 bkg[1] = False 1346 if 'FixedPoints' in value: 1347 if 'FixedPoints' in peaks: 1348 del peaks['FixedPoints'] 1349 elif key == 'Instrument Parameters': 1350 instrument, secondary = self.data['Instrument Parameters'] 1351 for iparam in value: 1352 instrument[iparam][2] = False 1353 else: 1354 raise ValueError("Unknown key:", key) 1355 1356 1357 class G2Phase(G2ObjectWrapper): 1358 """A wrapper object around a given phase. 1359 Author: Jackson O'Donnell jacksonhodonnell@gmail.com 1360 """ 1361 def __init__(self, data, name, proj): 1362 self.data = data 1363 self.name = name 1364 self.proj = proj 1365 1366 @staticmethod 1367 def is_valid_refinement_key(key): 1368 valid_keys = ["Cell", "Atoms", "LeBail"] 1369 return key in valid_keys 1370 1371 @staticmethod 1372 def is_valid_HAP_refinement_key(key): 1373 valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain", 1374 "Pref.Ori.", "Show", "Size", "Use", "Scale"] 1375 return key in valid_keys 1376 1377 def atom(self, atomlabel): 1378 """Returns the atom specified by atomlabel, or None if it does not 1379 exist. 1380 1381 :param str atomlabel: The name of the atom (e.g. "O2") 1382 1383 :returns: A :class:`G2AtomRecord` object 1384 representing the atom.""" 1385 # Consult GSASIIobj.py for the meaning of this 1386 cx, ct, cs, cia = self.data['General']['AtomPtrs'] 1387 ptrs = [cx, ct, cs, cia] 1388 atoms = self.data['Atoms'] 1389 for atom in atoms: 1390 if atom[ct-1] == atomlabel: 1391 return G2AtomRecord(atom, ptrs, self.proj) 1392 1393 def atoms(self): 1394 """Returns a list of atoms present in the phase. 1395 1396 :returns: A list of :class:`G2AtomRecord` objects. See 1397 also :func:`~GSASIIscriptable.G2Phase.atom`""" 1398 ptrs = self.data['General']['AtomPtrs'] 1399 output = [] 1400 atoms = self.data['Atoms'] 1401 for atom in atoms: 1402 output.append(G2AtomRecord(atom, ptrs, self.proj)) 1403 return output 1404 1405 def histograms(self): 1406 output = [] 1407 for hname in self.data.get('Histograms', {}).keys(): 1408 output.append(self.proj.histogram(hname)) 1409 return output 1410 1411 @property 1412 def ranId(self): 1413 return self.data['ranId'] 1414 1415 def cell_dict(self): 1416 """Returns a dictionary of the cell parameters, with keys: 1417 'a', 'b', 'c', 'alpha', 'beta', 'gamma', 'vol' 1418 1419 :returns: a dict""" 1420 cell = self['General']['Cell'] 1421 return {'a': cell[1], 'b': cell[2], 'c': cell[3], 1422 'alpha': cell[4], 'beta': cell[5], 'gamma': cell[6], 1423 'vol': cell[7]} 1424 1425 def set_refinements(self, refs): 1426 for key, value in refs.items(): 1427 if key == "Cell": 1428 self.data['General']['Cell'][0] = True 1429 elif key == "Atoms": 1430 cx, ct, cs, cia = self.data['General']['AtomPtrs'] 1431 1432 for atomlabel, atomrefinement in value.items(): 1433 if atomlabel == 'all': 1434 for atom in self.atoms(): 1435 atom.refinement_flags = atomrefinement 1436 else: 1437 atom = self.atom(atomlabel) 1438 atom.refinement_flags = atomrefinement 1439 elif key == "LeBail": 1440 hists = self.data['Histograms'] 1441 for hname, hoptions in hists.items(): 1442 if 'LeBail' not in hoptions: 1443 hoptions['newLeBail'] = True 1444 hoptions['LeBail'] = bool(value) 1445 else: 1446 raise ValueError("Unknown key:", key) 1447 1448 def clear_refinements(self, refs): 1449 for key, value in refs.items(): 1450 if key == "Cell": 1451 self.data['General']['Cell'][0] = False 1452 elif key == "Atoms": 1453 cx, ct, cs, cia = self.data['General']['AtomPtrs'] 1454 1455 for atomlabel in value: 1456 atom = self.atom(atomlabel) 1457 # Set refinement to none 1458 atom.refinement_flags = ' ' 1459 elif key == "LeBail": 1460 hists = self.data['Histograms'] 1461 for hname, hoptions in hists.items(): 1462 if 'LeBail' not in hoptions: 1463 hoptions['newLeBail'] = True 1464 hoptions['LeBail'] = False 1465 else: 1466 raise ValueError("Unknown key:", key) 1467 1468 1469 ########################## 1470 # Command Line Interface # 1471 ########################## 1472 1473 1474 def create(*args): 1475 """The create subcommand. 1476 1477 Should be passed all the command-line arguments after `create`""" 1478 proj = G2Project(filename=args[0]) 1479 1480 isPhase = False 1481 isPowderData = False 1482 isInstPrms = False 1483 instPrms = None 1484 1485 # TODO how to associate phase with histogram? 1486 for arg in args[1:]: 1487 if arg == '--phases': 1488 isPhase = True 1489 isPowderData = False 1490 isInstPrms = False 1491 elif arg == '--powder': 1492 isPhase = False 1493 isPowderData = True 1494 isInstPrms = False 1495 # Instrument parameters must be specified before 1496 # any powder data files are passed 1497 elif arg == '--iparams': 1498 isPhase = False 1499 isPowderData = False 1500 isInstPrms = True 1501 elif isPhase: 1502 proj.add_phase(arg) 1503 elif isPowderData: 1504 proj.add_powder_histogram(arg, instPrms) 1505 elif isInstPrms: 1506 instPrms = arg 1507 isInstPrms = False 1508 else: 1509 print("Not sure what to do with: {}".format(arg)) 1510 1511 proj.save() 1512 1513 1514 def dump(*args): 1515 """The dump subcommand""" 1516 raw = True 1517 histograms = False 1518 phases = False 1519 filenames = [] 1520 for arg in args: 1521 if arg in ['-h', '--histograms']: 1522 histograms = True 1523 raw = False 1524 elif arg in ['-p', '--phases']: 1525 phases = True 1526 raw = False 1527 elif arg in ['-r', '--raw']: 1528 raw = True 1529 else: 1530 filenames.append(arg) 1531 if raw: 1532 import IPython.lib.pretty as pretty 1533 for fname in filenames: 1534 if raw: 1535 proj, nameList = LoadDictFromProjFile(fname) 1536 print("file:", fname) 1537 print("names:", nameList) 1538 for key, val in proj.items(): 1539 print(key, ":") 1540 pretty.pprint(val) 1541 else: 1542 proj = G2Project(fname) 1543 if histograms: 1544 hists = proj.histograms() 1545 for h in hists: 1546 print(fname, h.name) 1547 if phases: 1548 phases = proj.phases() 1549 for p in phases: 1550 print(fname, p.name) 1551 1552 1553 def IPyBrowse(*args): 1554 """Load a .gpx file and then open a IPython shell to browse it 1555 """ 1556 filename = [] 1557 for arg in args: 1558 fname = arg 1559 proj, nameList = LoadDictFromProjFile(fname) 1560 msg = "\nfile {} loaded into proj (dict) with names in nameList".format(fname) 1561 GSASIIpath.IPyBreak_base(msg) 1562 break 1563 1564 1565 def refine(*args): 1566 """The refine subcommand""" 1567 proj = G2Project(args[0]) 1568 proj.refine() 1569 1570 1571 def seqrefine(*args): 1572 """The seqrefine subcommand""" 1573 raise NotImplementedError("seqrefine is not yet implemented") 1574 1575 1576 def export(*args): 1577 """The export subcommand""" 1578 # Export CIF or Structure or ... 1579 raise NotImplementedError("export is not yet implemented") 1580 1581 1582 subcommands = {"create": create, 1583 "dump": dump, 1584 "refine": refine, 1585 "seqrefine": seqrefine, 1586 "export": export, 1587 "browse": IPyBrowse} 1588 142 1589 143 1590 def main(): 144 'Needs a doc string' 145 arg = sys.argv 146 print arg 147 if len(arg) > 1: 148 GPXfile = arg[1] 149 if not ospath.exists(GPXfile): 150 print 'ERROR - ',GPXfile," doesn't exist!" 151 exit() 152 Project,nameList = LoadDictFromProjFile(GPXfile) 153 SaveDictToProjFile(Project,nameList,'testout.gpx') 1591 '''TODO: the command-line options need some thought 1592 ''' 1593 argv = sys.argv 1594 if len(argv) > 1 and argv[1] in subcommands: 1595 subcommands[argv[1]](*argv[2:]) 1596 elif len(argv) == 1 or argv[1] in ('help', '--help', '-h'): 1597 # TODO print usage 1598 subcommand_names = ' | '.join(sorted(subcommands.keys())) 1599 print("USAGE: {} [ {} ] ...".format(argv[0], subcommand_names)) 154 1600 else: 155 print 'ERROR - missing filename' 156 exit() 157 print "Done" 158 1601 print("Unknown subcommand: {}".format(argv[1])) 1602 print("Available subcommands:") 1603 for name in sorted(subcommands.keys()): 1604 print("\t{}".format(name)) 1605 sys.exit(-1) 1606 # arg = sys.argv 1607 # print(arg) 1608 # if len(arg) > 1: 1609 # GPXfile = arg[1] 1610 # if not ospath.exists(GPXfile): 1611 # print(u'ERROR - '+GPXfile+u" doesn't exist!") 1612 # exit() 1613 # Project,nameList = LoadDictFromProjFile(GPXfile) 1614 # SaveDictToProjFile(Project,nameList,'testout.gpx') 1615 # else: 1616 # print('ERROR - missing filename') 1617 # exit() 1618 # print("Done") 1619 1620 PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data") 1621 PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase") 159 1622 if __name__ == '__main__': 160 1623 main() 1624 1625 # from gpx_manipulatons.py 1626 # try: 1627 # filename, authorname = sys.argv[1:3] 1628 # proj, names = make_empty_project(authorname, filename) 1629 # SaveDictToProjFile(proj, names, os.path.abspath(filename)) 1630 # except ValueError: 1631 # print("Usage: {} <filename> <author>".format(sys.argv[0])) 1632 # sys.exit(-1) 1633 1634 1635 # from refinements.py 1636 # USAGE = """USAGE: {} datafile instparams phasefile projectname refinements 1637 1638 # datafile: Input powder data 1639 # intparams: Corresponding instrument parameter file 1640 # phasefile: Phase to refine against data 1641 # projectname: Project file to be created, should end in .gpx 1642 # refinements: JSON file of refinements to be executed 1643 # """ 1644 # try: 1645 # datafile, instprm, phasefile, projectname, refinements = sys.argv[1:] 1646 # except ValueError: 1647 # print(USAGE.format(sys.argv[0])) 1648 # sys.exit(-1) 1649 1650 # try: 1651 # with open(refinements) as f: 1652 # refinements = json.load(f) 1653 # except IOError: 1654 # print("No such refinements file: {}".format(refinements)) 1655 1656 # print("Creating project file \"{}\"...".format(projectname)) 1657 # proj = G2Project(filename=projectname) 1658 # # Add the histogram 1659 # hist = proj.add_powder_histogram(datafile, instprm) 1660 # # Add the phase, and associate it with the histogram 1661 # proj.add_phase(phasefile, histograms=[hist.name]) 1662 1663 # proj.do_refinements(refinements['refinements']) 1664 # proj.save() 1665 1666 1667 # from gpx_dumper 1668 # import IPython.lib.pretty as pretty 1669 # proj, nameList = LoadDictFromProjFile(sys.argv[1]) 1670 # print("names:", nameList) 1671 # for key, val in proj.items(): 1672 # print(key, ":", sep='') 1673 # pretty.pprint(val) 1674
Note: See TracChangeset
for help on using the changeset viewer.