source: branch/2frame/GSASIIscriptable.py @ 2984

Last change on this file since 2984 was 2984, checked in by odonnell, 5 years ago

added .plot() method to histograms, changed to (set|clear)_refinements

File size: 61.6 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2017-04-12 15:12:45 -0500 (Wed, 12 Apr 2017) $
5# $Author: vondreele $
6# $Revision: 2777 $
7# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIIscriptable.py $
8# $Id: GSASIIIO.py 2777 2017-04-12 20:12:45Z vondreele $
9########### SVN repository information ###################
10"""
11*GSASIIscriptable: Scripting Tools*
12-----------------------------------
13
14Routines for reading, writing, modifying and creating GSAS-II project (.gpx) files.
15
16Supports a command line interface as well.
17"""
18from __future__ import division, print_function # needed?
19import os.path as ospath
20import sys
21import cPickle
22import imp
23import copy
24import platform
25import os
26import random as ran
27import json
28
29import numpy.ma as ma
30import scipy.interpolate as si
31import numpy as np
32import scipy as sp
33import matplotlib as mpl
34
35import GSASIIpath
36GSASIIpath.SetBinaryPath() # would rather have this in __name__ == '__main__' stanza
37import GSASIIIO as G2IO
38import GSASIIfiles as G2fil
39import GSASIIobj as G2obj
40import GSASIIpwd as G2pwd
41import GSASIIstrIO as G2strIO
42import GSASIIspc as G2spc
43import GSASIIstrMain as G2strMain
44import GSASIImath as G2mth
45import GSASIIElem as G2elem
46
47try:
48    import matplotlib.pyplot as plt
49except ImportError:
50    plt = None
51
52def LoadDictFromProjFile(ProjFile):
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
58        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
59        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
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},
67        u'Rigid bodies':{'data': {rigid body dict}},
68        u'Covariance':{'data':{covariance data dict}},
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
75        [u'Notebook',],
76        [u'Controls',],
77        [u'Covariance',],
78        [u'Constraints',],
79        [u'Restraints',],
80        [u'Rigid bodies',],
81        [u'PWDR FAP.XRA Bank 1',
82             u'Comments',
83             u'Limits',
84             u'Background',
85             u'Instrument Parameters',
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        ]
93    '''
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
99    file = open(ProjFile,'rb')
100    # print('loading from file: {}'.format(ProjFile))
101    Project = {}
102    nameList = []
103    try:
104        while True:
105            try:
106                data = cPickle.load(file)
107            except EOFError:
108                break
109            datum = data[0]
110            Project[datum[0]] = {'data':datum[1]}
111            nameList.append([datum[0],])
112            for datus in data[1:]:
113                Project[datum[0]][datus[0]] = datus[1]
114                nameList[-1].append(datus[0])
115        file.close()
116        # print('project load successful')
117    except:
118        raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file")
119    finally:
120        file.close()
121    return Project,nameList
122
123def SaveDictToProjFile(Project,nameList,ProjFile):
124    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
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)
130    '''
131    file = open(ProjFile,'wb')
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
145def ImportPowder(reader,filename):
146    '''Use a reader to import a powder diffraction data file
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.
157    '''
158    rdfile,rdpath,descr = imp.find_module(reader)
159    rdclass = imp.load_module(reader,rdfile,rdpath,descr)
160    rd = rdclass.GSAS_ReaderClass()
161    if not rd.scriptable:
162        print(u'**** ERROR: '+reader+u' is not a scriptable reader')
163        return None
164    fl = open(filename,'rb')
165    rdlist = []
166    if rd.ContentsValidator(fl):
167        fl.seek(0)
168        repeat = True
169        rdbuffer = {} # create temporary storage for file reader
170        block = 0
171        while repeat: # loop if the reader asks for another pass on the file
172            block += 1
173            repeat = False
174            rd.objname = ospath.basename(filename)
175            flag = rd.Reader(filename,fl,None,buffer=rdbuffer,blocknum=block,)
176            if flag:
177                rdlist.append(copy.deepcopy(rd)) # save the result before it is written over
178                if rd.repeat:
179                    repeat = True
180        return rdlist
181    print(rd.errors)
182    return None
183
184def 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
211def 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
230def 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
355def 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
395class G2ImportException(Exception):
396    pass
397
398
399def 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
459def 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
498def 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
567def _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
602class 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
635class 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
1081class 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
1144class 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    def fit_fixed_points(self):
1165        """Attempts to apply a background fit to the fixed points currently specified."""
1166
1167        def SetInstParms(Inst):
1168            dataType = Inst['Type'][0]
1169            insVary = []
1170            insNames = []
1171            insVals = []
1172            for parm in Inst:
1173                insNames.append(parm)
1174                insVals.append(Inst[parm][1])
1175                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1176                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1177                        Inst[parm][2] = False
1178            instDict = dict(zip(insNames, insVals))
1179            instDict['X'] = max(instDict['X'], 0.01)
1180            instDict['Y'] = max(instDict['Y'], 0.01)
1181            if 'SH/L' in instDict:
1182                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1183            return dataType, instDict, insVary
1184
1185        bgrnd = self['Background']
1186
1187        # Need our fixed points in order
1188        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1189        X = [x for x, y in bgrnd[1]['FixedPoints']]
1190        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1191
1192        limits = self['Limits'][1]
1193        if X[0] > limits[0]:
1194            X = [limits[0]] + X
1195            Y = [Y[0]] + Y
1196        if X[-1] < limits[1]:
1197            X += [limits[1]]
1198            Y += [Y[-1]]
1199
1200        # Some simple lookups
1201        inst, inst2 = self['Instrument Parameters']
1202        pwddata = self['data'][1]
1203
1204        # Construct the data for background fitting
1205        xBeg = np.searchsorted(pwddata[0], limits[0])
1206        xFin = np.searchsorted(pwddata[0], limits[1])
1207        xdata = pwddata[0][xBeg:xFin]
1208        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1209
1210        W = [1]*len(xdata)
1211        Z = [0]*len(xdata)
1212
1213        dataType, insDict, insVary = SetInstParms(inst)
1214        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1215
1216        # Do the fit
1217        # TODO get controls from root
1218        data = np.array([xdata, ydata, W, Z, Z, Z])
1219        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1220                        bakVary, controls={})
1221
1222        # Post-fit
1223        parmDict = {}
1224        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1225        parmDict.update(bakDict)
1226        parmDict.update(insDict)
1227        pwddata[3][xBeg:xFin] *= 0
1228        pwddata[5][xBeg:xFin] *= 0
1229        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1230
1231        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1232        # TODO update background
1233        self.proj.save()
1234
1235    def y_calc(self):
1236        return self['data'][1][3]
1237
1238    def plot(self, Yobs=True, Ycalc=True, Background=True):
1239        if plt:
1240            data = self['data'][1]
1241            if Yobs:
1242                plt.plot(data[0], data[1], label='Yobs')
1243            if Ycalc:
1244                plt.plot(data[0], data[3], label='Ycalc')
1245            if Background:
1246                plt.plot(data[0], data[4], label='Background')
1247
1248    def set_refinements(self, refs):
1249        """Sets the refinement parameter 'key' to the specification 'value'
1250
1251        :param dict refs: A dictionary of the parameters to be set.
1252
1253        :returns: None"""
1254        do_fit_fixed_points = False
1255        for key, value in refs.items():
1256            if key == 'Limits':
1257                old_limits = self.data['Limits'][1]
1258                new_limits = value
1259                if isinstance(new_limits, dict):
1260                    if 'low' in new_limits:
1261                        old_limits[0] = new_limits['low']
1262                    if 'high' in new_limits:
1263                        old_limits[1] = new_limits['high']
1264                else:
1265                    old_limits[0], old_limits[1] = new_limits
1266            elif key == 'Sample Parameters':
1267                sample = self.data['Sample Parameters']
1268                for sparam in value:
1269                    sample[sparam][1] = True
1270            elif key == 'Background':
1271                bkg, peaks = self.data['Background']
1272
1273                # If True or False, just set the refine parameter
1274                if value in (True, False):
1275                    bkg[1] = value
1276                    return
1277
1278                if 'type' in value:
1279                    bkg[0] = value['type']
1280                if 'refine' in value:
1281                    bkg[1] = value['refine']
1282                if 'no. coeffs' in value:
1283                    cur_coeffs = bkg[2]
1284                    n_coeffs = value['no. coeffs']
1285                    if n_coeffs > cur_coeffs:
1286                        for x in range(n_coeffs - cur_coeffs):
1287                            bkg.append(0.0)
1288                    else:
1289                        for _ in range(cur_coeffs - n_coeffs):
1290                            bkg.pop()
1291                    bkg[2] = n_coeffs
1292                if 'coeffs' in value:
1293                    bkg[3:] = value['coeffs']
1294                if 'FixedPoints' in value:
1295                    peaks['FixedPoints'] = [(float(a), float(b))
1296                                            for a, b in value['FixedPoints']]
1297                if value.get('fit fixed points', False):
1298                    do_fit_fixed_points = True
1299
1300            elif key == 'Instrument Parameters':
1301                instrument, secondary = self.data['Instrument Parameters']
1302                for iparam in value:
1303                    try:
1304                        instrument[iparam][2] = True
1305                    except IndexError:
1306                        raise ValueError("Invalid key:", iparam)
1307            else:
1308                raise ValueError("Unknown key:", key)
1309        # Fit fixed points after the fact - ensure they are after fixed points
1310        # are added
1311        if do_fit_fixed_points:
1312            self.fit_fixed_points()
1313
1314    def clear_refinements(self, refs):
1315        """Clears the refinement parameter 'key' and its associated value."""
1316        for key, value in refs.items():
1317            if key == 'Limits':
1318                old_limits, cur_limits = self.data['Limits']
1319                cur_limits[0], cur_limits[1] = old_limits
1320            elif key == 'Sample Parameters':
1321                sample = self.data['Sample Parameters']
1322                for sparam in value:
1323                    sample[sparam][1] = False
1324            elif key == 'Background':
1325                bkg, peaks = self.data['Background']
1326
1327                # If True or False, just set the refine parameter
1328                if value in (True, False):
1329                    bkg[1] = False
1330                    return
1331
1332                bkg[1] = False
1333                if 'FixedPoints' in value:
1334                    if 'FixedPoints' in peaks:
1335                        del peaks['FixedPoints']
1336            elif key == 'Instrument Parameters':
1337                instrument, secondary = self.data['Instrument Parameters']
1338                for iparam in value:
1339                    instrument[iparam][2] = False
1340            else:
1341                raise ValueError("Unknown key:", key)
1342
1343
1344class G2Phase(G2ObjectWrapper):
1345    """A wrapper object around a given phase.
1346    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
1347    """
1348    def __init__(self, data, name, proj):
1349        self.data = data
1350        self.name = name
1351        self.proj = proj
1352
1353    @staticmethod
1354    def is_valid_refinement_key(key):
1355        valid_keys = ["Cell", "Atoms", "LeBail"]
1356        return key in valid_keys
1357
1358    def atom(self, atomlabel):
1359        """Returns the atom specified by atomlabel, or None if it does not
1360        exist.
1361
1362        :param str atomlabel: The name of the atom (e.g. "O2")
1363
1364        :returns: A :class:`G2AtomRecord` object
1365        representing the atom."""
1366        # Consult GSASIIobj.py for the meaning of this
1367        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1368        ptrs = [cx, ct, cs, cia]
1369        atoms = self.data['Atoms']
1370        for atom in atoms:
1371            if atom[ct-1] == atomlabel:
1372                return G2AtomRecord(atom, ptrs, self.proj)
1373
1374    def atoms(self):
1375        """Returns a list of atoms present in the phase.
1376
1377        :returns: A list of :class:`G2AtomRecord` objects. See
1378        also :func:`~GSASIIscriptable.G2Phase.atom`"""
1379        ptrs = self.data['General']['AtomPtrs']
1380        output = []
1381        atoms = self.data['Atoms']
1382        for atom in atoms:
1383            output.append(G2AtomRecord(atom, ptrs, self.proj))
1384        return output
1385
1386    def histograms(self):
1387        output = []
1388        for hname in self.data.get('Histograms', {}).keys():
1389            output.append(self.proj.histogram(hname))
1390        return output
1391
1392    @property
1393    def ranId(self):
1394        return self.data['ranId']
1395
1396    def cell_dict(self):
1397        """Returns a dictionary of the cell parameters, with keys:
1398            'a', 'b', 'c', 'alpha', 'beta', 'gamma', 'vol'
1399
1400        :returns: a dict"""
1401        cell = self['General']['Cell']
1402        return {'a': cell[1], 'b': cell[2], 'c': cell[3],
1403                'alpha': cell[4], 'beta': cell[5], 'gamma': cell[6],
1404                'vol': cell[7]}
1405
1406    def set_refinements(self, refs):
1407        for key, value in refs.items():
1408            if key == "Cell":
1409                self.data['General']['Cell'][0] = True
1410            elif key == "Atoms":
1411                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1412
1413                for atomlabel, atomrefinement in value.items():
1414                    if atomlabel == 'all':
1415                        for atom in self.atoms():
1416                            atom.refinement_flags = atomrefinement
1417                    else:
1418                        atom = self.atom(atomlabel)
1419                        atom.refinement_flags = atomrefinement
1420            elif key == "LeBail":
1421                hists = self.data['Histograms']
1422                for hname, hoptions in hists.items():
1423                    if 'LeBail' not in hoptions:
1424                        hoptions['newLeBail'] = True
1425                    hoptions['LeBail'] = bool(value)
1426            else:
1427                raise ValueError("Unknown key:", key)
1428
1429    def clear_refinements(self, refs):
1430        for key, value in refs.items():
1431            if key == "Cell":
1432                self.data['General']['Cell'][0] = False
1433            elif key == "Atoms":
1434                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1435
1436                for atomlabel in value:
1437                    atom = self.atom(atomlabel)
1438                    # Set refinement to none
1439                    atom.refinement_flags = ' '
1440            elif key == "LeBail":
1441                hists = self.data['Histograms']
1442                for hname, hoptions in hists.items():
1443                    if 'LeBail' not in hoptions:
1444                        hoptions['newLeBail'] = True
1445                    hoptions['LeBail'] = False
1446            else:
1447                raise ValueError("Unknown key:", key)
1448
1449
1450##########################
1451# Command Line Interface #
1452##########################
1453
1454
1455def create(*args):
1456    """The create subcommand.
1457
1458    Should be passed all the command-line arguments after `create`"""
1459    proj = G2Project(filename=args[0])
1460
1461    isPhase = False
1462    isPowderData = False
1463    isInstPrms = False
1464    instPrms = None
1465
1466    # TODO how to associate phase with histogram?
1467    for arg in args[1:]:
1468        if arg == '--phases':
1469            isPhase = True
1470            isPowderData = False
1471            isInstPrms = False
1472        elif arg == '--powder':
1473            isPhase = False
1474            isPowderData = True
1475            isInstPrms = False
1476        # Instrument parameters must be specified before
1477        # any powder data files are passed
1478        elif arg == '--iparams':
1479            isPhase = False
1480            isPowderData = False
1481            isInstPrms = True
1482        elif isPhase:
1483            proj.add_phase(arg)
1484        elif isPowderData:
1485            proj.add_powder_histogram(arg, instPrms)
1486        elif isInstPrms:
1487            instPrms = arg
1488            isInstPrms = False
1489        else:
1490            print("Not sure what to do with: {}".format(arg))
1491
1492    proj.save()
1493
1494
1495def dump(*args):
1496    """The dump subcommand"""
1497    raw = True
1498    histograms = False
1499    phases = False
1500    filenames = []
1501    for arg in args:
1502        if arg in ['-h', '--histograms']:
1503            histograms = True
1504            raw = False
1505        elif arg in ['-p', '--phases']:
1506            phases = True
1507            raw = False
1508        elif arg in ['-r', '--raw']:
1509            raw = True
1510        else:
1511            filenames.append(arg)
1512    if raw:
1513        import IPython.lib.pretty as pretty
1514    for fname in filenames:
1515        if raw:
1516            proj, nameList = LoadDictFromProjFile(fname)
1517            print("file:", fname)
1518            print("names:", nameList)
1519            for key, val in proj.items():
1520                print(key, ":")
1521                pretty.pprint(val)
1522        else:
1523            proj = G2Project(fname)
1524            if histograms:
1525                hists = proj.histograms()
1526                for h in hists:
1527                    print(fname, h.name)
1528            if phases:
1529                phases = proj.phases()
1530                for p in phases:
1531                    print(fname, p.name)
1532
1533
1534def IPyBrowse(*args):
1535    """Load a .gpx file and then open a IPython shell to browse it
1536    """
1537    filename = []
1538    for arg in args:
1539        fname = arg
1540        proj, nameList = LoadDictFromProjFile(fname)
1541        msg = "\nfile {} loaded into proj (dict) with names in nameList".format(fname)
1542        GSASIIpath.IPyBreak_base(msg)
1543        break
1544
1545
1546def refine(*args):
1547    """The refine subcommand"""
1548    proj = G2Project(args[0])
1549    proj.refine()
1550
1551
1552def seqrefine(*args):
1553    """The seqrefine subcommand"""
1554    raise NotImplementedError("seqrefine is not yet implemented")
1555
1556
1557def export(*args):
1558    """The export subcommand"""
1559    # Export CIF or Structure or ...
1560    raise NotImplementedError("export is not yet implemented")
1561
1562
1563subcommands = {"create": create,
1564               "dump": dump,
1565               "refine": refine,
1566               "seqrefine": seqrefine,
1567               "export": export,
1568               "browse": IPyBrowse}
1569
1570
1571def main():
1572    '''TODO: the command-line options need some thought
1573    '''
1574    argv = sys.argv
1575    if len(argv) > 1 and argv[1] in subcommands:
1576        subcommands[argv[1]](*argv[2:])
1577    elif len(argv) == 1 or argv[1] in ('help', '--help', '-h'):
1578        # TODO print usage
1579        subcommand_names = ' | '.join(sorted(subcommands.keys()))
1580        print("USAGE: {} [ {} ] ...".format(argv[0], subcommand_names))
1581    else:
1582        print("Unknown subcommand: {}".format(argv[1]))
1583        print("Available subcommands:")
1584        for name in sorted(subcommands.keys()):
1585            print("\t{}".format(name))
1586        sys.exit(-1)
1587    # arg = sys.argv
1588    # print(arg)
1589    # if len(arg) > 1:
1590    #     GPXfile = arg[1]
1591    #     if not ospath.exists(GPXfile):
1592    #         print(u'ERROR - '+GPXfile+u" doesn't exist!")
1593    #         exit()
1594    #     Project,nameList = LoadDictFromProjFile(GPXfile)
1595    #     SaveDictToProjFile(Project,nameList,'testout.gpx')
1596    # else:
1597    #     print('ERROR - missing filename')
1598    #     exit()
1599    # print("Done")
1600
1601PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
1602PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
1603if __name__ == '__main__':
1604    main()
1605
1606    # from gpx_manipulatons.py
1607    # try:
1608    #     filename, authorname = sys.argv[1:3]
1609    #     proj, names = make_empty_project(authorname, filename)
1610    #     SaveDictToProjFile(proj, names, os.path.abspath(filename))
1611    # except ValueError:
1612    #     print("Usage: {} <filename> <author>".format(sys.argv[0]))
1613    #     sys.exit(-1)
1614
1615
1616    # from refinements.py
1617#      USAGE = """USAGE: {} datafile instparams phasefile projectname refinements
1618
1619# datafile:    Input powder data
1620# intparams:   Corresponding instrument parameter file
1621# phasefile:   Phase to refine against data
1622# projectname: Project file to be created, should end in .gpx
1623# refinements: JSON file of refinements to be executed
1624# """
1625#     try:
1626#         datafile, instprm, phasefile, projectname, refinements = sys.argv[1:]
1627#     except ValueError:
1628#         print(USAGE.format(sys.argv[0]))
1629#         sys.exit(-1)
1630
1631#     try:
1632#         with open(refinements) as f:
1633#             refinements = json.load(f)
1634#     except IOError:
1635#         print("No such refinements file: {}".format(refinements))
1636
1637#     print("Creating project file \"{}\"...".format(projectname))
1638#     proj = G2Project(filename=projectname)
1639#     # Add the histogram
1640#     hist = proj.add_powder_histogram(datafile, instprm)
1641#     # Add the phase, and associate it with the histogram
1642#     proj.add_phase(phasefile, histograms=[hist.name])
1643
1644#     proj.do_refinements(refinements['refinements'])
1645#     proj.save()
1646
1647
1648    # from gpx_dumper
1649    # import IPython.lib.pretty as pretty
1650    # proj, nameList = LoadDictFromProjFile(sys.argv[1])
1651    # print("names:", nameList)
1652    # for key, val in proj.items():
1653    #     print(key, ":", sep='')
1654    #     pretty.pprint(val)
1655
Note: See TracBrowser for help on using the repository browser.