source: trunk/GSASIIscriptable.py @ 3000

Last change on this file since 3000 was 3000, checked in by toby, 8 years ago

make the two frame version the trunk as we hit version 3000
File size: 62.3 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(False) # 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    @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
1357class 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
1474def 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
1514def 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
1553def 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
1565def refine(*args):
1566    """The refine subcommand"""
1567    proj = G2Project(args[0])
1568    proj.refine()
1569
1570
1571def seqrefine(*args):
1572    """The seqrefine subcommand"""
1573    raise NotImplementedError("seqrefine is not yet implemented")
1574
1575
1576def export(*args):
1577    """The export subcommand"""
1578    # Export CIF or Structure or ...
1579    raise NotImplementedError("export is not yet implemented")
1580
1581
1582subcommands = {"create": create,
1583               "dump": dump,
1584               "refine": refine,
1585               "seqrefine": seqrefine,
1586               "export": export,
1587               "browse": IPyBrowse}
1588
1589
1590def main():
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))
1600    else:
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
1620PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
1621PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
1622if __name__ == '__main__':
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 TracBrowser for help on using the repository browser.