source: trunk/GSASIIscriptable.py @ 3090

Last change on this file since 3090 was 3090, checked in by odonnell, 4 years ago

more refinement specifier documentation

File size: 90.0 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
18Look at :class:`G2Project` to start.
19
20=====================
21Refinement parameters
22=====================
23
24There are three classes of refinement parameters:
25
26    * Histogram. Turned on and off through :func:`~G2PwdrData.set_refinements`
27      and :func:`~G2PwdrData.clear_refinements`
28    * Phase. Turned on and off through :func:`~G2Phase.set_refinements`
29      and :func:`~G2Phase.clear_refinements`
30    * Histogram-and-phase (HAP). Turned on and off through
31      :func:`~G2Phase.set_HAP_refinements` and :func:`~G2Phase.clear_HAP_refinements`
32
33
34============================
35Refinement specifiers format
36============================
37
38.. _Histogram_parameters_table:
39
40--------------------
41Histogram parameters
42--------------------
43
44This table describes the dictionaries supplied to :func:`~G2PwdrData.set_refinements`
45and :func:`~G2PwdrData.clear_refinements`.
46
47.. tabularcolumns:: |l|p|
48
49===================== ====================  ====================
50key                   subkey                explanation
51===================== ====================  ====================
52Limits                                      The 2-theta range of values to consider. Can
53                                            be either a dictionary of 'low' and 'high',
54                                            or a list of 2 items [low, high]
55\                     low                   Sets the low limit
56\                     high                  Sets the high limit
57Sample Parameters                           Should be provided as a **list** of subkeys
58                                            to set or clear, e.g. ['DisplaceX', 'Scale']
59\                     Absorption
60\                     Contrast
61\                     DisplaceX             Sample displacement along the X direction
62\                     DisplaceY             Sample displacement along the Y direction
63\                     Scale                 Histogram Scale factor
64Background                                  Sample background. If value is a boolean,
65                                            the background's 'refine' parameter is set
66                                            to the given boolean. Usually should be a
67                                            dictionary with any of the following keys:
68\                     type                  The background model, e.g. 'chebyschev'
69\                     refine                The value of the refine flag, boolean
70\                     no. coeffs            Number of coefficients to use, integer
71\                     coeffs                List of floats, literal values for background
72\                     FixedPoints           List of (2-theta, intensity) values for fixed points
73\                     fit fixed points      If True, triggers a fit to the fixed points to be calculated. It is calculated when this key is detected, regardless of calls to refine.
74Instrument Parameters                       As in Sample Paramters, Should be provided as a **list** of subkeys to
75                                            set or clear, e.g. ['X', 'Y', 'Zero', 'SH/L']
76\                     U, V, W               All separate keys. Gaussian peak profile terms
77\                     X, Y                  Separate keys. Lorentzian peak profile terms
78\                     Zero                  Zero shift
79\                     SH/L
80\                     Polariz.              Polarization parameter
81\                     Lam                   Lambda, the incident wavelength
82===================== ====================  ====================
83
84
85.. _Phase_parameters_table:
86
87----------------
88Phase parameters
89----------------
90
91This table describes the dictionaries supplied to :func:`~G2Phase.set_refinements`
92and :func:`~G2Phase.clear_refinements`.
93
94===================== ====================
95key                   explanation
96===================== ====================
97Cell                  Whether or not to refine the unit cell.
98Atoms                 Dictionary of atoms and refinement flags.
99                      Each key should be an atom label, e.g.
100                      'O3', 'Mn5', and each value should be
101                      a string defining what values to refine.
102                      Values can be any combination of 'F'
103                      for fractional occupancy, 'X' for position,
104                      and 'U' for Debye-Waller factor
105LeBail                Enables LeBail intensity extraction.
106===================== ====================
107
108.. _HAP_parameters_table:
109
110------------------------------
111Histogram-and-phase parameters
112------------------------------
113
114This table describes the dictionaries supplied to :func:`~G2Phase.set_HAP_refinements`
115and :func:`~G2Phase.clear_HAP_refinements`.
116
117===================== ===================== ====================
118key                   subkey                explanation
119===================== ===================== ====================
120Babinet                                     Should be a **list** of the following
121                                            subkeys. If not, assumes both
122                                            BabA and BabU
123\                     BabA
124\                     BabU
125Extinction                                  Should be boolean, whether or not to
126                                            refine.
127HStrain                                     Should be boolean, whether or not to
128                                            refine.
129Mustrain
130\                     type                  Mustrain model. One of 'isotropic',
131                                            'uniaxial', or 'generalized'
132\                     direction             For uniaxial. A list of three integers,
133                                            the [hkl] direction of the axis.
134\                     refine                Usually boolean, whether or not to refine.
135                                            When in doubt, set it to true.
136                                            For uniaxial model, can specify list
137                                            of 'axial' or 'equatorial'. If boolean
138                                            given sets both axial and equatorial.
139Pref.Ori.                                   Boolean, whether to refine
140Show                                        Boolean, whether to refine
141Size                                        Not implemented
142Use                                         Boolean, whether to refine
143Scale                                       Boolean, whether to refine
144===================== ===================== ====================
145
146============================
147Scriptable API
148============================
149"""
150from __future__ import division, print_function # needed?
151import argparse
152import os.path as ospath
153import datetime as dt
154import sys
155import cPickle
156import imp
157import copy
158import os
159import random as ran
160
161import numpy.ma as ma
162import scipy.interpolate as si
163import numpy as np
164import scipy as sp
165
166import GSASIIpath
167GSASIIpath.SetBinaryPath(False) # would rather have this in __name__ == '__main__' stanza
168import GSASIIobj as G2obj
169import GSASIIpwd as G2pwd
170import GSASIIstrMain as G2strMain
171import GSASIIspc as G2spc
172import GSASIIElem as G2elem
173
174# Delay imports to not slow down small scripts
175G2fil = None
176PwdrDataReaders = []
177PhaseReaders = []
178
179
180def LoadG2fil():
181    """Delay importing this module, it is slow"""
182    global G2fil
183    global PwdrDataReaders
184    global PhaseReaders
185    if G2fil is None:
186        import GSASIIfiles
187        G2fil = GSASIIfiles
188        PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
189        PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
190
191
192def LoadDictFromProjFile(ProjFile):
193    '''Read a GSAS-II project file and load items to dictionary
194    :param str ProjFile: GSAS-II project (name.gpx) full file name
195    :returns: Project,nameList, where
196
197      * Project (dict) is a representation of gpx file following the GSAS-II tree struture
198        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
199        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
200      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
201
202    Example for fap.gpx::
203
204      Project = {                 #NB:dict order is not tree order
205        u'Phases':{'data':None,'fap':{phase dict}},
206        u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
207        u'Rigid bodies':{'data': {rigid body dict}},
208        u'Covariance':{'data':{covariance data dict}},
209        u'Controls':{'data':{controls data dict}},
210        u'Notebook':{'data':[notebook list]},
211        u'Restraints':{'data':{restraint data dict}},
212        u'Constraints':{'data':{constraint data dict}}]
213        }
214      nameList = [                #NB: reproduces tree order
215        [u'Notebook',],
216        [u'Controls',],
217        [u'Covariance',],
218        [u'Constraints',],
219        [u'Restraints',],
220        [u'Rigid bodies',],
221        [u'PWDR FAP.XRA Bank 1',
222             u'Comments',
223             u'Limits',
224             u'Background',
225             u'Instrument Parameters',
226             u'Sample Parameters',
227             u'Peak List',
228             u'Index Peak List',
229             u'Unit Cells List',
230             u'Reflection Lists'],
231        [u'Phases', u'fap']
232        ]
233    '''
234    # Let IOError be thrown if file does not exist
235    # if not ospath.exists(ProjFile):
236    #     print ('\n*** Error attempt to open project file that does not exist:\n   '+
237    #         str(ProjFile))
238    #     return
239    file = open(ProjFile,'rb')
240    # print('loading from file: {}'.format(ProjFile))
241    Project = {}
242    nameList = []
243    try:
244        while True:
245            try:
246                data = cPickle.load(file)
247            except EOFError:
248                break
249            datum = data[0]
250            Project[datum[0]] = {'data':datum[1]}
251            nameList.append([datum[0],])
252            for datus in data[1:]:
253                Project[datum[0]][datus[0]] = datus[1]
254                nameList[-1].append(datus[0])
255        file.close()
256        # print('project load successful')
257    except:
258        raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file")
259    finally:
260        file.close()
261    return Project,nameList
262
263def SaveDictToProjFile(Project,nameList,ProjFile):
264    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
265
266    :param dict Project: representation of gpx file following the GSAS-II
267        tree structure as described for LoadDictFromProjFile
268    :param list nameList: names of main tree entries & subentries used to reconstruct project file
269    :param str ProjFile: full file name for output project.gpx file (including extension)
270    '''
271    file = open(ProjFile,'wb')
272    try:
273        for name in nameList:
274            data = []
275            item = Project[name[0]]
276            data.append([name[0],item['data']])
277            for item2 in name[1:]:
278                data.append([item2,item[item2]])
279            cPickle.dump(data,file,1)
280    finally:
281        file.close()
282
283def ImportPowder(reader,filename):
284    '''Use a reader to import a powder diffraction data file
285
286    :param str reader: a scriptable reader
287    :param str filename: full name of powder data file; can be "multi-Bank" data
288
289    :returns list rdlist: list of reader objects containing powder data, one for each
290        "Bank" of data encountered in file. Items in reader object of interest are:
291
292          * rd.comments: list of str: comments found on powder file
293          * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist
294          * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in  GSASII data tree.
295    '''
296    rdfile,rdpath,descr = imp.find_module(reader)
297    rdclass = imp.load_module(reader,rdfile,rdpath,descr)
298    rd = rdclass.GSAS_ReaderClass()
299    if not rd.scriptable:
300        print(u'**** ERROR: '+reader+u' is not a scriptable reader')
301        return None
302    fl = open(filename,'rb')
303    rdlist = []
304    if rd.ContentsValidator(fl):
305        fl.seek(0)
306        repeat = True
307        rdbuffer = {} # create temporary storage for file reader
308        block = 0
309        while repeat: # loop if the reader asks for another pass on the file
310            block += 1
311            repeat = False
312            rd.objname = ospath.basename(filename)
313            flag = rd.Reader(filename,fl,None,buffer=rdbuffer,blocknum=block,)
314            if flag:
315                rdlist.append(copy.deepcopy(rd)) # save the result before it is written over
316                if rd.repeat:
317                    repeat = True
318        return rdlist
319    print(rd.errors)
320    return None
321
322def SetDefaultDData(dType,histoName,NShkl=0,NDij=0):
323    '''Create an initial Histogram dictionary
324
325    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
326    '''
327    if dType in ['SXC','SNC']:
328        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
329            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
330            'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955,
331            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}],
332            'Flack':[0.0,False]}
333    elif dType == 'SNT':
334        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
335            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
336            'Extinction':['Lorentzian','None', {
337            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]}
338    elif 'P' in dType:
339        return {'Histogram':histoName,'Show':False,'Scale':[1.0,False],
340            'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1],
341            'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1],
342                [1.,1.,1.,0.,0.,0.],6*[False,]],
343            'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1],
344                NShkl*[0.01,],NShkl*[False,]],
345            'HStrain':[NDij*[0.0,],NDij*[False,]],
346            'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}}
347
348
349def PreSetup(data):
350    '''Create part of an initial (empty) phase dictionary
351
352    from GSASIIphsGUI.py, near end of UpdatePhaseData
353
354    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
355    '''
356    if 'RBModels' not in data:
357        data['RBModels'] = {}
358    if 'MCSA' not in data:
359        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
360    if 'dict' in str(type(data['MCSA']['Results'])):
361        data['MCSA']['Results'] = []
362    if 'Modulated' not in data['General']:
363        data['General']['Modulated'] = False
364    if 'modulated' in data['General']['Type']:
365        data['General']['Modulated'] = True
366        data['General']['Type'] = 'nuclear'
367
368
369def SetupGeneral(data, dirname):
370    """Helps initialize phase data.
371
372    From GSASIIphsGui.py, function of the same name. Minor changes for imports etc.
373
374    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
375    """
376    mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True,
377                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
378    generalData = data['General']
379    atomData = data['Atoms']
380    generalData['AtomTypes'] = []
381    generalData['Isotopes'] = {}
382
383    if 'Isotope' not in generalData:
384        generalData['Isotope'] = {}
385    if 'Data plot type' not in generalData:
386        generalData['Data plot type'] = 'Mustrain'
387    if 'POhkl' not in generalData:
388        generalData['POhkl'] = [0,0,1]
389    if 'Map' not in generalData:
390        generalData['Map'] = mapDefault.copy()
391    if 'Flip' not in generalData:
392        generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
393            'k-factor':0.1,'k-Max':20.,}
394    if 'testHKL' not in generalData['Flip']:
395        generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]]
396    if 'doPawley' not in generalData:
397        generalData['doPawley'] = False     #ToDo: change to ''
398    if 'Pawley dmin' not in generalData:
399        generalData['Pawley dmin'] = 1.0
400    if 'Pawley neg wt' not in generalData:
401        generalData['Pawley neg wt'] = 0.0
402    if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
403        'MCSA controls' not in generalData:
404        generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50],
405        'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0,
406        'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True}
407    if 'AtomPtrs' not in generalData:
408        generalData['AtomPtrs'] = [3,1,7,9]
409        if generalData['Type'] == 'macromolecular':
410            generalData['AtomPtrs'] = [6,4,10,12]
411        elif generalData['Type'] == 'magnetic':
412            generalData['AtomPtrs'] = [3,1,10,12]
413    if generalData['Type'] in ['modulated',]:
414        generalData['Modulated'] = True
415        generalData['Type'] = 'nuclear'
416        if 'Super' not in generalData:
417            generalData['Super'] = 1
418            generalData['SuperVec'] = [[0,0,.1],False,4]
419            generalData['SSGData'] = {}
420        if '4DmapData' not in generalData:
421            generalData['4DmapData'] = mapDefault.copy()
422            generalData['4DmapData'].update({'MapType':'Fobs'})
423    if 'Modulated' not in generalData:
424        generalData['Modulated'] = False
425    if 'HydIds' not in generalData:
426        generalData['HydIds'] = {}
427    cx,ct,cs,cia = generalData['AtomPtrs']
428    generalData['NoAtoms'] = {}
429    generalData['BondRadii'] = []
430    generalData['AngleRadii'] = []
431    generalData['vdWRadii'] = []
432    generalData['AtomMass'] = []
433    generalData['Color'] = []
434    if generalData['Type'] == 'magnetic':
435        generalData['MagDmin'] = generalData.get('MagDmin',1.0)
436        landeg = generalData.get('Lande g',[])
437    generalData['Mydir'] = dirname
438    badList = {}
439    for iat,atom in enumerate(atomData):
440        atom[ct] = atom[ct].lower().capitalize()              #force to standard form
441        if generalData['AtomTypes'].count(atom[ct]):
442            generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1])
443        elif atom[ct] != 'UNK':
444            Info = G2elem.GetAtomInfo(atom[ct])
445            if not Info:
446                if atom[ct] not in badList:
447                    badList[atom[ct]] = 0
448                badList[atom[ct]] += 1
449                atom[ct] = 'UNK'
450                continue
451            atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
452            generalData['AtomTypes'].append(atom[ct])
453            generalData['Z'] = Info['Z']
454            generalData['Isotopes'][atom[ct]] = Info['Isotopes']
455            generalData['BondRadii'].append(Info['Drad'])
456            generalData['AngleRadii'].append(Info['Arad'])
457            generalData['vdWRadii'].append(Info['Vdrad'])
458            if atom[ct] in generalData['Isotope']:
459                if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
460                    isotope = generalData['Isotopes'][atom[ct]].keys()[-1]
461                    generalData['Isotope'][atom[ct]] = isotope
462                generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass'])
463            else:
464                generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
465                if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
466                    isotope = generalData['Isotopes'][atom[ct]].keys()[-1]
467                    generalData['Isotope'][atom[ct]] = isotope
468                generalData['AtomMass'].append(Info['Mass'])
469            generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1])
470            generalData['Color'].append(Info['Color'])
471            if generalData['Type'] == 'magnetic':
472                if len(landeg) < len(generalData['AtomTypes']):
473                    landeg.append(2.0)
474    if generalData['Type'] == 'magnetic':
475        generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])]
476
477    if badList:
478        msg = 'Warning: element symbol(s) not found:'
479        for key in badList:
480            msg += '\n\t' + key
481            if badList[key] > 1:
482                msg += ' (' + str(badList[key]) + ' times)'
483        raise Exception("Phase error:\n" + msg)
484        # wx.MessageBox(msg,caption='Element symbol error')
485    F000X = 0.
486    F000N = 0.
487    for i,elem in enumerate(generalData['AtomTypes']):
488        F000X += generalData['NoAtoms'][elem]*generalData['Z']
489        isotope = generalData['Isotope'][elem]
490        F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0]
491    generalData['F000X'] = F000X
492    generalData['F000N'] = F000N
493    import GSASIImath as G2mth
494    generalData['Mass'] = G2mth.getMass(generalData)
495
496
497def make_empty_project(author=None, filename=None):
498    """Creates an dictionary in the style of GSASIIscriptable, for an empty
499    project.
500
501    If no author name or filename is supplied, 'no name' and
502    <current dir>/test_output.gpx are used , respectively.
503
504    Returns: project dictionary, name list
505
506    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
507    """
508    if not filename:
509        filename = os.path.join(os.getcwd(), 'test_output.gpx')
510    else:
511        filename = os.path.abspath(filename)
512    gsasii_version = str(GSASIIpath.GetVersionNumber())
513    LoadG2fil()
514    import matplotlib as mpl
515    python_library_versions = G2fil.get_python_versions([mpl, np, sp])
516
517    controls_data = dict(G2obj.DefaultControls)
518    controls_data['LastSavedAs'] = unicode(filename)
519    controls_data['LastSavedUsing'] = gsasii_version
520    controls_data['PythonVersions'] = python_library_versions
521    if author:
522        controls_data['Author'] = author
523
524    output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [],
525                                       'Global': []}},
526              'Controls': {'data': controls_data},
527              u'Covariance': {'data': {}},
528              u'Notebook': {'data': ['']},
529              u'Restraints': {'data': {}},
530              u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []},
531                                'Residue': {'AtInfo': {}},
532                                'Vector':  {'AtInfo': {}}}}}
533
534    names = [[u'Notebook'], [u'Controls'], [u'Covariance'],
535             [u'Constraints'], [u'Restraints'], [u'Rigid bodies']]
536
537    return output, names
538
539
540class G2ImportException(Exception):
541    pass
542
543
544def import_generic(filename, readerlist):
545    """Attempt to import a filename, using a list of reader objects.
546
547    Returns the first reader object which worked."""
548    # Translated from OnImportGeneric method in GSASII.py
549    primaryReaders, secondaryReaders = [], []
550    for reader in readerlist:
551        flag = reader.ExtensionValidator(filename)
552        if flag is None:
553            secondaryReaders.append(reader)
554        elif flag:
555            primaryReaders.append(reader)
556    if not secondaryReaders and not primaryReaders:
557        raise G2ImportException("Could not read file: ", filename)
558
559    with open(filename, 'Ur') as fp:
560        rd_list = []
561
562        for rd in primaryReaders + secondaryReaders:
563            # Initialize reader
564            rd.selections = []
565            rd.dnames = []
566            rd.ReInitialize()
567            # Rewind file
568            fp.seek(0)
569            rd.errors = ""
570            if not rd.ContentsValidator(fp):
571                # Report error
572                pass
573            if len(rd.selections) > 1:
574                # Select data?
575                # GSASII.py:543
576                raise G2ImportException("Not sure what data to select")
577
578            block = 0
579            rdbuffer = {}
580            fp.seek(0)
581            repeat = True
582            while repeat:
583                repeat = False
584                block += 1
585                rd.objname = os.path.basename(filename)
586                flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block)
587                if flag:
588                    # Omitting image loading special cases
589                    rd.readfilename = filename
590                    rd_list.append(copy.deepcopy(rd))
591                    repeat = rd.repeat
592                else:
593                    raise G2ImportException("{}.Reader() returned:".format(rd),
594                                            flag)
595
596            if rd_list:
597                if rd.warnings:
598                    print("Read warning by", rd.formatName, "reader:",
599                          rd.warnings, file=sys.stderr)
600                return rd_list
601    raise G2ImportException("No reader could read file: " + filename)
602
603
604def load_iprms(instfile, reader):
605    """Loads instrument parameters from a file, and edits the
606    given reader.
607
608    Returns a 2-tuple of (Iparm1, Iparm2) parameters
609    """
610    LoadG2fil()
611    ext = os.path.splitext(instfile)[1]
612
613    if ext.lower() == '.instprm':
614        # New GSAS File, load appropriately
615        with open(instfile) as f:
616            lines = f.readlines()
617        bank = reader.powderentry[2]
618        numbanks = reader.numbanks
619        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
620
621        reader.instfile = instfile
622        reader.instmsg = 'GSAS-II file' + instfile
623        return iparms
624    elif ext.lower() not in ('.prm', '.inst', '.ins'):
625        raise ValueError('Expected .prm file, found: ', instfile)
626
627    # It's an old GSAS file, load appropriately
628    Iparm = {}
629    with open(instfile, 'Ur') as fp:
630        for line in fp:
631            if '#' in line:
632                continue
633            Iparm[line[:12]] = line[12:-1]
634    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
635    if ibanks == 1:
636        reader.instbank = 1
637        reader.powderentry[2] = 1
638        reader.instfile = instfile
639        reader.instmsg = instfile + ' bank ' + str(reader.instbank)
640        return G2fil.SetPowderInstParms(Iparm, reader)
641    # TODO handle >1 banks
642    raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm")
643
644def load_pwd_from_reader(reader, instprm, existingnames=[]):
645    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
646
647    :returns: (name, tree) - 2-tuple of the histogram name (str), and data
648
649    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
650    """
651    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
652    HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames))
653
654    try:
655        Iparm1, Iparm2 = instprm
656    except ValueError:
657        Iparm1, Iparm2 = load_iprms(instprm, reader)
658
659    Ymin = np.min(reader.powderdata[1])
660    Ymax = np.max(reader.powderdata[1])
661    valuesdict = {'wtFactor': 1.0,
662                  'Dummy': False,
663                  'ranId': ran.randint(0, sys.maxint),
664                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
665                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
666                  'qPlot': False, 'dPlot': False, 'sqrtPlot': False,
667                  'Yminmax': [Ymin, Ymax]}
668    reader.Sample['ranId'] = valuesdict['ranId']
669
670    # Ending keys:
671    # [u'Reflection Lists',
672    #  u'Limits',
673    #  'data',
674    #  u'Index Peak List',
675    #  u'Comments',
676    #  u'Unit Cells List',
677    #  u'Sample Parameters',
678    #  u'Peak List',
679    #  u'Background',
680    #  u'Instrument Parameters']
681    Tmin = np.min(reader.powderdata[0])
682    Tmax = np.max(reader.powderdata[0])
683
684    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
685                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
686
687    output_dict = {u'Reflection Lists': {},
688                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
689                   u'data': [valuesdict, reader.powderdata, HistName],
690                   u'Index Peak List': [[], []],
691                   u'Comments': reader.comments,
692                   u'Unit Cells List': [],
693                   u'Sample Parameters': reader.Sample,
694                   u'Peak List': {'peaks': [], 'sigDict': {}},
695                   u'Background': reader.pwdparms.get('Background', default_background),
696                   u'Instrument Parameters': [Iparm1, Iparm2],
697                   }
698
699    names = [u'Comments',
700             u'Limits',
701             u'Background',
702             u'Instrument Parameters',
703             u'Sample Parameters',
704             u'Peak List',
705             u'Index Peak List',
706             u'Unit Cells List',
707             u'Reflection Lists']
708
709    # TODO controls?? GSASII.py:1664-7
710
711    return HistName, [HistName] + names, output_dict
712
713
714def _deep_copy_into(from_, into):
715    """Helper function for reloading .gpx file. See G2Project.reload()
716
717    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
718    """
719    if isinstance(from_, dict) and isinstance(into, dict):
720        combined_keys = set(from_.keys()).union(into.keys())
721        for key in combined_keys:
722            if key in from_ and key in into:
723                both_dicts = (isinstance(from_[key], dict)
724                              and isinstance(into[key], dict))
725                both_lists = (isinstance(from_[key], list)
726                              and isinstance(into[key], list))
727                if both_dicts or both_lists:
728                    _deep_copy_into(from_[key], into[key])
729                else:
730                    into[key] = from_[key]
731            elif key in from_:
732                into[key] = from_[key]
733            else:  # key in into
734                del into[key]
735    elif isinstance(from_, list) and isinstance(into, list):
736        if len(from_) == len(into):
737            for i in xrange(len(from_)):
738                both_dicts = (isinstance(from_[i], dict)
739                              and isinstance(into[i], dict))
740                both_lists = (isinstance(from_[i], list)
741                              and isinstance(into[i], list))
742                if both_dicts or both_lists:
743                    _deep_copy_into(from_[i], into[i])
744                else:
745                    into[i] = from_[i]
746        else:
747            into[:] = from_
748
749
750class G2ObjectWrapper(object):
751    """Base class for all GSAS-II object wrappers.
752
753    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
754    of overrides are implemented so that the wrapper behaves like a dictionary.
755
756    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
757    """
758    def __init__(self, datadict):
759        self.data = datadict
760
761    def __getitem__(self, key):
762        return self.data[key]
763
764    def __setitem__(self, key, value):
765        self.data[key] = value
766
767    def __contains__(self, key):
768        return key in self.data
769
770    def get(self, k, d=None):
771        return self.data.get(k, d)
772
773    def keys(self):
774        return self.data.keys()
775
776    def values(self):
777        return self.data.values()
778
779    def items(self):
780        return self.data.items()
781
782
783class G2Project(G2ObjectWrapper):
784    """Represents an entire GSAS-II project."""
785    def __init__(self, gpxfile=None, author=None, filename=None):
786        """Loads a GSAS-II project from a specified filename.
787
788        :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
789            creates an empty project.
790        :param str author: Author's name. Optional.
791        :param str filename: The filename the project should be saved to in
792            the future. If both filename and gpxfile are present, the project is
793            loaded from the gpxfile, then set to  save to filename in the future"""
794        if gpxfile is None:
795            filename = os.path.abspath(os.path.expanduser(filename))
796            self.filename = filename
797            self.data, self.names = make_empty_project(author=author, filename=filename)
798        elif isinstance(gpxfile, str):
799            # TODO set author, filename
800            self.filename = os.path.abspath(os.path.expanduser(gpxfile))
801            self.data, self.names = LoadDictFromProjFile(gpxfile)
802            self.update_ids()
803        else:
804            raise ValueError("Not sure what to do with gpxfile")
805
806    @classmethod
807    def from_dict_and_names(cls, gpxdict, names, filename=None):
808        """Creates a :class:`G2Project` directly from
809        a dictionary and a list of names. If in doubt, do not use this.
810
811        :returns: a :class:`G2Project`
812        """
813        out = cls()
814        if filename:
815            filename = os.path.abspath(os.path.expanduser(filename))
816            out.filename = filename
817            gpxdict['Controls']['data']['LastSavedAs'] = filename
818        else:
819            try:
820                out.filename = gpxdict['Controls']['data']['LastSavedAs']
821            except KeyError:
822                out.filename = None
823        out.data = gpxdict
824        out.names = names
825
826    def save(self, filename=None):
827        """Saves the project, either to the current filename, or to a new file.
828
829        Updates self.filename if a new filename provided"""
830        # TODO update LastSavedUsing ?
831        if filename:
832            filename = os.path.abspath(os.path.expanduser(filename))
833            self.data['Controls']['data']['LastSavedAs'] = filename
834            self.filename = filename
835        elif not self.filename:
836            raise AttributeError("No file name to save to")
837        SaveDictToProjFile(self.data, self.names, self.filename)
838
839    def add_powder_histogram(self, datafile, iparams, phases=[]):
840        """Loads a powder data histogram into the project.
841
842        Automatically checks for an instrument parameter file, or one can be
843        provided.
844
845        :param str datafile: The powder data file to read, a filename.
846        :param str iparams: The instrument parameters file, a filename.
847        :param list phases: Phases to link to the new histogram
848
849        :returns: A :class:`G2PwdrData` object representing
850            the histogram
851        """
852        LoadG2fil()
853        datafile = os.path.abspath(os.path.expanduser(datafile))
854        iparams = os.path.abspath(os.path.expanduser(iparams))
855        pwdrreaders = import_generic(datafile, PwdrDataReaders)
856        histname, new_names, pwdrdata = load_pwd_from_reader(
857                                          pwdrreaders[0], iparams,
858                                          [h.name for h in self.histograms()])
859        if histname in self.data:
860            print("Warning - redefining histogram", histname)
861        elif self.names[-1][0] == 'Phases':
862            self.names.insert(-1, new_names)
863        else:
864            self.names.append(new_names)
865        self.data[histname] = pwdrdata
866        self.update_ids()
867
868        for phase in phases:
869            phase = self.phase(phase)
870            self.link_histogram_phase(histname, phase)
871
872        return self.histogram(histname)
873
874    def add_phase(self, phasefile, phasename=None, histograms=[]):
875        """Loads a phase into the project from a .cif file
876
877        :param str phasefile: The CIF file from which to import the phase.
878        :param str phasename: The name of the new phase, or None for the default
879        :param list histograms: The names of the histograms to associate with
880            this phase
881
882        :returns: A :class:`G2Phase` object representing the
883            new phase.
884        """
885        LoadG2fil()
886        histograms = [self.histogram(h).name for h in histograms]
887        phasefile = os.path.abspath(os.path.expanduser(phasefile))
888
889        # TODO handle multiple phases in a file
890        phasereaders = import_generic(phasefile, PhaseReaders)
891        phasereader = phasereaders[0]
892
893        phasename = phasename or phasereader.Phase['General']['Name']
894        phaseNameList = [p.name for p in self.phases()]
895        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
896        phasereader.Phase['General']['Name'] = phasename
897
898        if 'Phases' not in self.data:
899            self.data[u'Phases'] = { 'data': None }
900        assert phasename not in self.data['Phases'], "phase names should be unique"
901        self.data['Phases'][phasename] = phasereader.Phase
902
903        if phasereader.Constraints:
904            Constraints = self.data['Constraints']
905            for i in phasereader.Constraints:
906                if isinstance(i, dict):
907                    if '_Explain' not in Constraints:
908                        Constraints['_Explain'] = {}
909                    Constraints['_Explain'].update(i)
910                else:
911                    Constraints['Phase'].append(i)
912
913        data = self.data['Phases'][phasename]
914        generalData = data['General']
915        SGData = generalData['SGData']
916        NShkl = len(G2spc.MustrainNames(SGData))
917        NDij = len(G2spc.HStrainNames(SGData))
918        Super = generalData.get('Super', 0)
919        if Super:
920            SuperVec = np.array(generalData['SuperVec'][0])
921        else:
922            SuperVec = []
923        UseList = data['Histograms']
924
925        for hist in histograms:
926            self.link_histogram_phase(hist, phasename)
927
928        for obj in self.names:
929            if obj[0] == 'Phases':
930                phasenames = obj
931                break
932        else:
933            phasenames = [u'Phases']
934            self.names.append(phasenames)
935        phasenames.append(unicode(phasename))
936
937        # TODO should it be self.filename, not phasefile?
938        SetupGeneral(data, os.path.dirname(phasefile))
939        self.index_ids()
940
941        self.update_ids()
942        return self.phase(phasename)
943
944    def link_histogram_phase(self, histogram, phase):
945        """Associates a given histogram and phase.
946
947        .. seealso::
948
949            :func:`~G2Project.histogram`
950            :func:`~G2Project.phase`"""
951        hist = self.histogram(histogram)
952        phase = self.phase(phase)
953
954        generalData = phase['General']
955
956        if hist.name.startswith('HKLF '):
957            raise NotImplementedError("HKLF not yet supported")
958        elif hist.name.startswith('PWDR '):
959            hist['Reflection Lists'][generalData['Name']] = {}
960            UseList = phase['Histograms']
961            SGData = generalData['SGData']
962            NShkl = len(G2spc.MustrainNames(SGData))
963            NDij = len(G2spc.HStrainNames(SGData))
964            UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij)
965            UseList[hist.name]['hId'] = hist.id
966            for key, val in [('Use', True), ('LeBail', False),
967                             ('newLeBail', True),
968                             ('Babinet', {'BabA': [0.0, False],
969                                          'BabU': [0.0, False]})]:
970                if key not in UseList[hist.name]:
971                    UseList[hist.name][key] = val
972        else:
973            raise RuntimeError("Unexpected histogram" + hist.name)
974
975
976    def reload(self):
977        """Reload self from self.filename"""
978        data, names = LoadDictFromProjFile(self.filename)
979        self.names = names
980        # Need to deep copy the new data file data into the current tree,
981        # so that any existing G2Phase, or G2PwdrData objects will still be
982        # valid
983        _deep_copy_into(from_=data, into=self.data)
984
985    def refine(self, newfile=None, printFile=None, makeBack=False):
986        # TODO migrate to RefineCore
987        # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
988        #      calcControls,pawleyLookup,ifPrint,printFile,dlg)
989        # index_ids will automatically save the project
990        self.index_ids()
991        # TODO G2strMain does not properly use printFile
992        G2strMain.Refine(self.filename, makeBack=makeBack)
993        # Reload yourself
994        self.reload()
995
996    def histogram(self, histname):
997        """Returns the histogram named histname, or None if it does not exist.
998
999        :param histname: The name of the histogram (str), or ranId or index.
1000        :returns: A :class:`G2PwdrData` object, or None if
1001            the histogram does not exist
1002
1003        .. seealso::
1004            :func:`~GSASIIscriptable.G2Project.histograms`
1005            :func:`~GSASIIscriptable.G2Project.phase`
1006            :func:`~GSASIIscriptable.G2Project.phases`
1007            """
1008        if isinstance(histname, G2PwdrData):
1009            if histname.proj == self:
1010                return histname
1011        if histname in self.data:
1012            return G2PwdrData(self.data[histname], self)
1013        try:
1014            # see if histname is an id or ranId
1015            histname = int(histname)
1016        except ValueError:
1017            return
1018
1019        for histogram in self.histograms():
1020            if histogram.id == histname or histogram.ranId == histname:
1021                return histogram
1022
1023    def histograms(self):
1024        """Return a list of all histograms, as
1025        :class:`G2PwdrData` objects
1026
1027        .. seealso::
1028            :func:`~GSASIIscriptable.G2Project.histograms`
1029            :func:`~GSASIIscriptable.G2Project.phase`
1030            :func:`~GSASIIscriptable.G2Project.phases`
1031            """
1032        output = []
1033        for obj in self.names:
1034            if len(obj) > 1 and obj[0] != u'Phases':
1035                output.append(self.histogram(obj[0]))
1036        return output
1037
1038    def phase(self, phasename):
1039        """
1040        Gives an object representing the specified phase in this project.
1041
1042        :param str phasename: The name of the desired phase. Either the name
1043            (str), the phase's ranId, or the phase's index
1044        :returns: A :class:`G2Phase` object
1045        :raises: KeyError
1046
1047        .. seealso::
1048            :func:`~GSASIIscriptable.G2Project.histograms`
1049            :func:`~GSASIIscriptable.G2Project.phase`
1050            :func:`~GSASIIscriptable.G2Project.phases`
1051            """
1052        phases = self.data['Phases']
1053        if phasename in phases:
1054            return G2Phase(phases[phasename], phasename, self)
1055
1056        try:
1057            # phasename should be phase index or ranId
1058            phasename = int(phasename)
1059        except ValueError:
1060            return
1061
1062        for phase in self.phases():
1063            if phase.id == phasename or phase.ranId == phasename:
1064                return phase
1065
1066    def phases(self):
1067        """
1068        Returns a list of all the phases in the project.
1069
1070        :returns: A list of :class:`G2Phase` objects
1071
1072        .. seealso::
1073            :func:`~GSASIIscriptable.G2Project.histogram`
1074            :func:`~GSASIIscriptable.G2Project.histograms`
1075            :func:`~GSASIIscriptable.G2Project.phase`
1076            """
1077        for obj in self.names:
1078            if obj[0] == 'Phases':
1079                return [self.phase(p) for p in obj[1:]]
1080        return []
1081
1082    def update_ids(self):
1083        """Makes sure all phases and histograms have proper hId and pId"""
1084        # Translated from GetUsedHistogramsAndPhasesfromTree,
1085        #   GSASIIdataGUI.py:4107
1086        for i, h in enumerate(self.histograms()):
1087            h.id = i
1088        for i, p in enumerate(self.phases()):
1089            p.id = i
1090
1091    def do_refinements(self, refinements, histogram='all', phase='all',
1092                       outputnames=None, makeBack=False):
1093        """Conducts a series of refinements. Wrapper around iter_refinements
1094
1095        :param list refinements: A list of dictionaries defining refinements
1096        :param str histogram: Name of histogram for refinements to be applied
1097            to, or 'all'
1098        :param str phase: Name of phase for refinements to be applied to, or
1099            'all'
1100        """
1101        for proj in self.iter_refinements(refinements, histogram, phase,
1102                                          outputnames, makeBack):
1103            pass
1104        return self
1105
1106    def iter_refinements(self, refinements, histogram='all', phase='all',
1107                         outputnames=None, makeBack=False):
1108        """Conducts a series of refinements, iteratively. Stops after every
1109        refinement and yields this project, to allow error checking or
1110        logging of intermediate results.
1111
1112        >>> def checked_refinements(proj):
1113        ...     for p in proj.iter_refinements(refs):
1114        ...         # Track intermediate results
1115        ...         log(p.histogram('0').residuals)
1116        ...         log(p.phase('0').get_cell())
1117        ...         # Check if parameter diverged, nonsense answer, or whatever
1118        ...         if is_something_wrong(p):
1119        ...             raise Exception("I need a human!")
1120
1121        :param list refinements: A list of dictionaries defining refinements
1122        :param str histogram: Name of histogram for refinements to be applied
1123            to, or 'all'
1124        :param str phase: Name of phase for refinements to be applied to, or
1125            'all'
1126        """
1127        if outputnames:
1128            if len(refinements) != len(outputnames):
1129                raise ValueError("Should have same number of outuputs to"
1130                                 "refinements")
1131        else:
1132            outputnames = [None for r in refinements]
1133
1134        for output, refinement in zip(outputnames, refinements):
1135            self.set_refinement(refinement, histogram)
1136            # Handle 'once' args - refinements that are disabled after this
1137            # refinement
1138            if 'once' in refinement:
1139                temp = {'set': refinement['once']}
1140                self.set_refinement(temp, histogram, phase)
1141
1142            if output:
1143                self.save(output)
1144
1145            self.refine(makeBack=makeBack)
1146            yield self
1147
1148            # Handle 'once' args - refinements that are disabled after this
1149            # refinement
1150            if 'once' in refinement:
1151                temp = {'clear': refinement['once']}
1152                self.set_refinement(temp, histogram, phase)
1153
1154    def set_refinement(self, refinement, histogram='all', phase='all'):
1155        """Apply specified refinements to a given histogram(s) or phase(s).
1156
1157        Refinement parameters are categorize in three groups:
1158
1159        1. Histogram parameters
1160        2. Phase parameters
1161        3. Histogram-and-Phase (HAP) parameters
1162
1163        :param dict refinement: The refinements to be conducted
1164        :param histogram: Either a name of a histogram (str), a list of
1165            histogram names, or 'all' (default)
1166        :param phase: Either a name of a phase (str), a list of phase names, or
1167            'all' (default)
1168
1169        .. seealso::
1170            :func:`~G2PwdrData.set_refinements`
1171            :func:`~G2PwdrData.clear_refinements`
1172            :func:`~G2Phase.set_refinements`
1173            :func:`~G2Phase.clear_refinements`
1174            :func:`~G2Phase.set_HAP_refinements`
1175            :func:`~G2Phase.clear_HAP_refinements`"""
1176
1177        if histogram == 'all':
1178            hists = self.histograms()
1179        elif isinstance(histogram, str) or isinstance(histogram, int):
1180            hists = [self.histogram(histogram)]
1181        else:
1182            hists = [self.histogram(name) for name in histogram]
1183
1184        if phase == 'all':
1185            phases = self.phases()
1186        elif isinstance(phase, str) or isinstance(phase, int):
1187            phases = [self.phase(phase)]
1188        else:
1189            phases = [self.phase(name) for name in phase]
1190
1191        # TODO: HAP parameters:
1192        #   Babinet
1193        #   Extinction
1194        #   HStrain
1195        #   Mustrain
1196        #   Pref. Ori
1197        #   Size
1198
1199        pwdr_set = {}
1200        phase_set = {}
1201        hap_set = {}
1202        for key, val in refinement.get('set', {}).items():
1203            # Apply refinement options
1204            if G2PwdrData.is_valid_refinement_key(key):
1205                pwdr_set[key] = val
1206            elif G2Phase.is_valid_refinement_key(key):
1207                phase_set[key] = val
1208            elif G2Phase.is_valid_HAP_refinement_key(key):
1209                hap_set[key] = val
1210            else:
1211                raise ValueError("Unknown refinement key", key)
1212
1213        for hist in hists:
1214            hist.set_refinements(pwdr_set)
1215        for phase in phases:
1216            phase.set_refinements(phase_set)
1217        for phase in phases:
1218            phase.set_HAP_refinements(hap_set, hists)
1219
1220        pwdr_clear = {}
1221        phase_clear = {}
1222        hap_clear = {}
1223        for key, val in refinement.get('clear', {}).items():
1224            # Clear refinement options
1225            if G2PwdrData.is_valid_refinement_key(key):
1226                pwdr_clear[key] = val
1227            elif G2Phase.is_valid_refinement_key(key):
1228                phase_clear[key] = val
1229            elif G2Phase.is_valid_HAP_refinement_key(key):
1230                hap_set[key] = val
1231            else:
1232                raise ValueError("Unknown refinement key", key)
1233
1234        for hist in hists:
1235            hist.clear_refinements(pwdr_clear)
1236        for phase in phases:
1237            phase.clear_refinements(phase_clear)
1238        for phase in phases:
1239            phase.clear_HAP_refinements(hap_clear, hists)
1240
1241    def index_ids(self):
1242        import GSASIIstrIO as G2strIO
1243        self.save()
1244        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1245
1246    def add_constraint_raw(self, cons_scope, constr):
1247        """Adds a constraint of type consType to the project.
1248        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1249
1250        WARNING it does not check the constraint is well-constructed"""
1251        constrs = self.data['Constraints']['data']
1252        if 'Global' not in constrs:
1253            constrs['Global'] = []
1254        constrs[cons_scope].append(constr)
1255
1256    def hold_many(self, vars, type):
1257        """Apply holds for all the variables in vars, for constraint of a given type.
1258
1259        type is passed directly to add_constraint_raw as consType
1260
1261        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
1262            string variable specifiers, or arguments for :meth:`make_var_obj`
1263        :param str type: A string constraint type specifier. See
1264            :class:`~GSASIIscriptable.G2Project.add_constraint_raw`
1265
1266        """
1267        for var in vars:
1268            if isinstance(var, str):
1269                var = self.make_var_obj(var)
1270            elif not isinstance(var, G2obj.G2VarObj):
1271                var = self.make_var_obj(*var)
1272            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1273
1274    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1275                     reloadIdx=True):
1276        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1277        or individual names of phase, histogram, varname, and atomId.
1278
1279        Automatically converts string phase, hist, or atom names into the ID required
1280        by G2VarObj."""
1281
1282        if reloadIdx:
1283            self.index_ids()
1284
1285        # If string representation, short circuit
1286        if hist is None and varname is None and atomId is None:
1287            if isinstance(phase, str) and ':' in phase:
1288                return G2obj.G2VarObj(phase)
1289
1290        # Get phase index
1291        phaseObj = None
1292        if isinstance(phase, G2Phase):
1293            phaseObj = phase
1294            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1295        elif phase in self.data['Phases']:
1296            phaseObj = self.phase(phase)
1297            phaseRanId = phaseObj.ranId
1298            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1299
1300        # Get histogram index
1301        if isinstance(hist, G2PwdrData):
1302            hist = G2obj.HistRanIdLookup[hist.ranId]
1303        elif hist in self.data:
1304            histRanId = self.histogram(hist).ranId
1305            hist = G2obj.HistRanIdLookup[histRanId]
1306
1307        # Get atom index (if any)
1308        if isinstance(atomId, G2AtomRecord):
1309            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1310        elif phaseObj:
1311            atomObj = phaseObj.atom(atomId)
1312            if atomObj:
1313                atomRanId = atomObj.ranId
1314                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1315
1316        return G2obj.G2VarObj(phase, hist, varname, atomId)
1317
1318
1319class G2AtomRecord(G2ObjectWrapper):
1320    """Wrapper for an atom record. Has convenient accessors via @property.
1321
1322
1323    Available accessors: label, type, refinement_flags, coordinates,
1324        occupancy, ranId, id, adp_flag, uiso
1325
1326    Example:
1327
1328    >>> atom = some_phase.atom("O3")
1329    >>> # We can access the underlying data format
1330    >>> atom.data
1331    ['O3', 'O-2', '', ... ]
1332    >>> # We can also use wrapper accessors
1333    >>> atom.coordinates
1334    (0.33, 0.15, 0.5)
1335    >>> atom.refinement_flags
1336    u'FX'
1337    >>> atom.ranId
1338    4615973324315876477
1339    >>> atom.occupancy
1340    1.0
1341    """
1342    def __init__(self, data, indices, proj):
1343        self.data = data
1344        self.cx, self.ct, self.cs, self.cia = indices
1345        self.proj = proj
1346
1347    @property
1348    def label(self):
1349        return self.data[self.ct-1]
1350
1351    @property
1352    def type(self):
1353        return self.data[self.ct]
1354
1355    @property
1356    def refinement_flags(self):
1357        return self.data[self.ct+1]
1358
1359    @refinement_flags.setter
1360    def refinement_flags(self, other):
1361        # Automatically check it is a valid refinement
1362        for c in other:
1363            if c not in ' FXU':
1364                raise ValueError("Invalid atom refinement: ", other)
1365        self.data[self.ct+1] = unicode(other)
1366
1367    @property
1368    def coordinates(self):
1369        return tuple(self.data[self.cx:self.cx+3])
1370
1371    @property
1372    def occupancy(self):
1373        return self.data[self.cx+3]
1374
1375    @occupancy.setter
1376    def occupancy(self, val):
1377        self.data[self.cx+3] = float(val)
1378
1379    @property
1380    def ranId(self):
1381        return self.data[self.cia+8]
1382
1383    @property
1384    def adp_flag(self):
1385        # Either 'I' or 'A'
1386        return self.data[self.cia]
1387
1388    @property
1389    def uiso(self):
1390        if self.adp_flag == 'I':
1391            return self.data[self.cia+1]
1392        else:
1393            return self.data[self.cia+2:self.cia+8]
1394
1395    @uiso.setter
1396    def uiso(self, value):
1397        if self.adp_flag == 'I':
1398            self.data[self.cia+1] = float(value)
1399        else:
1400            assert len(value) == 6
1401            self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
1402
1403
1404class G2PwdrData(G2ObjectWrapper):
1405    """Wraps a Powder Data Histogram."""
1406    def __init__(self, data, proj):
1407        self.data = data
1408        self.proj = proj
1409
1410    @staticmethod
1411    def is_valid_refinement_key(key):
1412        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1413                      'Instrument Parameters']
1414        return key in valid_keys
1415
1416    @property
1417    def name(self):
1418        return self.data['data'][-1]
1419
1420    @property
1421    def ranId(self):
1422        return self.data['data'][0]['ranId']
1423
1424    @property
1425    def residuals(self):
1426        data = self.data['data'][0]
1427        return {key: data[key]
1428                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1429
1430    @property
1431    def id(self):
1432        self.proj.update_ids()
1433        return self.data['data'][0]['hId']
1434
1435    @id.setter
1436    def id(self, val):
1437        self.data['data'][0]['hId'] = val
1438
1439    def fit_fixed_points(self):
1440        """Attempts to apply a background fit to the fixed points currently specified."""
1441        def SetInstParms(Inst):
1442            dataType = Inst['Type'][0]
1443            insVary = []
1444            insNames = []
1445            insVals = []
1446            for parm in Inst:
1447                insNames.append(parm)
1448                insVals.append(Inst[parm][1])
1449                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1450                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1451                        Inst[parm][2] = False
1452            instDict = dict(zip(insNames, insVals))
1453            instDict['X'] = max(instDict['X'], 0.01)
1454            instDict['Y'] = max(instDict['Y'], 0.01)
1455            if 'SH/L' in instDict:
1456                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1457            return dataType, instDict, insVary
1458
1459        bgrnd = self.data['Background']
1460
1461        # Need our fixed points in order
1462        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1463        X = [x for x, y in bgrnd[1]['FixedPoints']]
1464        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1465
1466        limits = self.data['Limits'][1]
1467        if X[0] > limits[0]:
1468            X = [limits[0]] + X
1469            Y = [Y[0]] + Y
1470        if X[-1] < limits[1]:
1471            X += [limits[1]]
1472            Y += [Y[-1]]
1473
1474        # Some simple lookups
1475        controls = self.proj['Controls']['data']
1476        inst, inst2 = self.data['Instrument Parameters']
1477        pwddata = self.data['data'][1]
1478
1479        # Construct the data for background fitting
1480        xBeg = np.searchsorted(pwddata[0], limits[0])
1481        xFin = np.searchsorted(pwddata[0], limits[1])
1482        xdata = pwddata[0][xBeg:xFin]
1483        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1484
1485        W = [1]*len(xdata)
1486        Z = [0]*len(xdata)
1487
1488        dataType, insDict, insVary = SetInstParms(inst)
1489        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1490
1491        # Do the fit
1492        data = np.array([xdata, ydata, W, Z, Z, Z])
1493        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1494                        prevVaryList=bakVary, controls=controls)
1495
1496        # Post-fit
1497        parmDict = {}
1498        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1499        parmDict.update(bakDict)
1500        parmDict.update(insDict)
1501        pwddata[3][xBeg:xFin] *= 0
1502        pwddata[5][xBeg:xFin] *= 0
1503        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1504
1505        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1506        # TODO update background
1507        self.proj.save()
1508
1509    def y_calc(self):
1510        return self.data['data'][1][3]
1511
1512    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1513        try:
1514            import matplotlib.pyplot as plt
1515            data = self.data['data'][1]
1516            if Yobs:
1517                plt.plot(data[0], data[1], label='Yobs')
1518            if Ycalc:
1519                plt.plot(data[0], data[3], label='Ycalc')
1520            if Background:
1521                plt.plot(data[0], data[4], label='Background')
1522            if Residual:
1523                plt.plot(data[0], data[5], label="Residual")
1524        except ImportError:
1525            pass
1526
1527    def set_refinements(self, refs):
1528        """Sets the refinement parameter 'key' to the specification 'value'
1529
1530        :param dict refs: A dictionary of the parameters to be set. See
1531                          :ref:`Histogram_parameters_table` for a description of
1532                          what these dictionaries should be.
1533
1534        :returns: None
1535
1536        """
1537        do_fit_fixed_points = False
1538        for key, value in refs.items():
1539            if key == 'Limits':
1540                old_limits = self.data['Limits'][1]
1541                new_limits = value
1542                if isinstance(new_limits, dict):
1543                    if 'low' in new_limits:
1544                        old_limits[0] = new_limits['low']
1545                    if 'high' in new_limits:
1546                        old_limits[1] = new_limits['high']
1547                else:
1548                    old_limits[0], old_limits[1] = new_limits
1549            elif key == 'Sample Parameters':
1550                sample = self.data['Sample Parameters']
1551                for sparam in value:
1552                    if sparam not in sample:
1553                        raise ValueError("Unknown refinement parameter, "
1554                                         + str(sparam))
1555                    sample[sparam][1] = True
1556            elif key == 'Background':
1557                bkg, peaks = self.data['Background']
1558
1559                # If True or False, just set the refine parameter
1560                if value in (True, False):
1561                    bkg[1] = value
1562                    return
1563
1564                if 'type' in value:
1565                    bkg[0] = value['type']
1566                if 'refine' in value:
1567                    bkg[1] = value['refine']
1568                if 'no. coeffs' in value:
1569                    cur_coeffs = bkg[2]
1570                    n_coeffs = value['no. coeffs']
1571                    if n_coeffs > cur_coeffs:
1572                        for x in range(n_coeffs - cur_coeffs):
1573                            bkg.append(0.0)
1574                    else:
1575                        for _ in range(cur_coeffs - n_coeffs):
1576                            bkg.pop()
1577                    bkg[2] = n_coeffs
1578                if 'coeffs' in value:
1579                    bkg[3:] = value['coeffs']
1580                if 'FixedPoints' in value:
1581                    peaks['FixedPoints'] = [(float(a), float(b))
1582                                            for a, b in value['FixedPoints']]
1583                if value.get('fit fixed points', False):
1584                    do_fit_fixed_points = True
1585
1586            elif key == 'Instrument Parameters':
1587                instrument, secondary = self.data['Instrument Parameters']
1588                for iparam in value:
1589                    try:
1590                        instrument[iparam][2] = True
1591                    except IndexError:
1592                        raise ValueError("Invalid key:", iparam)
1593            else:
1594                raise ValueError("Unknown key:", key)
1595        # Fit fixed points after the fact - ensure they are after fixed points
1596        # are added
1597        if do_fit_fixed_points:
1598            # Background won't be fit if refinement flag not set
1599            orig = self.data['Background'][0][1]
1600            self.data['Background'][0][1] = True
1601            self.fit_fixed_points()
1602            # Restore the previous value
1603            self.data['Background'][0][1] = orig
1604
1605    def clear_refinements(self, refs):
1606        """Clears the refinement parameter 'key' and its associated value.
1607
1608        :param dict refs: A dictionary of parameters to clear."""
1609        for key, value in refs.items():
1610            if key == 'Limits':
1611                old_limits, cur_limits = self.data['Limits']
1612                cur_limits[0], cur_limits[1] = old_limits
1613            elif key == 'Sample Parameters':
1614                sample = self.data['Sample Parameters']
1615                for sparam in value:
1616                    sample[sparam][1] = False
1617            elif key == 'Background':
1618                bkg, peaks = self.data['Background']
1619
1620                # If True or False, just set the refine parameter
1621                if value in (True, False):
1622                    bkg[1] = False
1623                    return
1624
1625                bkg[1] = False
1626                if 'FixedPoints' in value:
1627                    if 'FixedPoints' in peaks:
1628                        del peaks['FixedPoints']
1629            elif key == 'Instrument Parameters':
1630                instrument, secondary = self.data['Instrument Parameters']
1631                for iparam in value:
1632                    instrument[iparam][2] = False
1633            else:
1634                raise ValueError("Unknown key:", key)
1635
1636
1637class G2Phase(G2ObjectWrapper):
1638    """A wrapper object around a given phase.
1639
1640    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1641    """
1642    def __init__(self, data, name, proj):
1643        self.data = data
1644        self.name = name
1645        self.proj = proj
1646
1647    @staticmethod
1648    def is_valid_refinement_key(key):
1649        valid_keys = ["Cell", "Atoms", "LeBail"]
1650        return key in valid_keys
1651
1652    @staticmethod
1653    def is_valid_HAP_refinement_key(key):
1654        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1655                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1656        return key in valid_keys
1657
1658    def atom(self, atomlabel):
1659        """Returns the atom specified by atomlabel, or None if it does not
1660        exist.
1661
1662        :param str atomlabel: The name of the atom (e.g. "O2")
1663        :returns: A :class:`G2AtomRecord` object
1664            representing the atom.
1665        """
1666        # Consult GSASIIobj.py for the meaning of this
1667        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1668        ptrs = [cx, ct, cs, cia]
1669        atoms = self.data['Atoms']
1670        for atom in atoms:
1671            if atom[ct-1] == atomlabel:
1672                return G2AtomRecord(atom, ptrs, self.proj)
1673
1674    def atoms(self):
1675        """Returns a list of atoms present in the phase.
1676
1677        :returns: A list of :class:`G2AtomRecord` objects.
1678
1679        .. seealso::
1680            :func:`~GSASIIscriptable.G2Phase.atom`
1681            :class:`G2AtomRecord`
1682        """
1683        ptrs = self.data['General']['AtomPtrs']
1684        output = []
1685        atoms = self.data['Atoms']
1686        for atom in atoms:
1687            output.append(G2AtomRecord(atom, ptrs, self.proj))
1688        return output
1689
1690    def histograms(self):
1691        output = []
1692        for hname in self.data.get('Histograms', {}).keys():
1693            output.append(self.proj.histogram(hname))
1694        return output
1695
1696    @property
1697    def ranId(self):
1698        return self.data['ranId']
1699
1700    @property
1701    def id(self):
1702        return self.data['pId']
1703
1704    @id.setter
1705    def id(self, val):
1706        self.data['pId'] = val
1707
1708    def get_cell(self):
1709        """Returns a dictionary of the cell parameters, with keys:
1710            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
1711
1712        :returns: a dict
1713
1714        .. seealso::
1715           :func:`~G2Phase.get_cell_and_esd`
1716
1717        """
1718        cell = self.data['General']['Cell']
1719        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
1720                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
1721                'volume': cell[7]}
1722
1723    def get_cell_and_esd(self):
1724        """
1725        Returns a pair of dictionaries, the first representing the unit cell, the second
1726        representing the estimated standard deviations of the unit cell.
1727
1728        :returns: a tuple of two dictionaries
1729
1730        .. seealso::
1731           :func:`~G2Phase.get_cell`
1732
1733        """
1734        # translated from GSASIIstrIO.ExportBaseclass.GetCell
1735        import GSASIIstrIO as G2stIO
1736        import GSASIIlattice as G2lat
1737        import GSASIImapvars as G2mv
1738        try:
1739            pfx = str(self.id) + '::'
1740            sgdata = self['General']['SGData']
1741            covDict = self.proj['Covariance']['data']
1742
1743            parmDict = dict(zip(covDict.get('varyList',[]),
1744                                covDict.get('variables',[])))
1745            sigDict = dict(zip(covDict.get('varyList',[]),
1746                               covDict.get('sig',[])))
1747
1748            if covDict.get('covMatrix') is not None:
1749                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1750                                                  covDict['varyList'],
1751                                                  parmDict))
1752
1753            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
1754            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
1755            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
1756            cellDict, cellSigDict = {}, {}
1757            for i, key in enumerate(['length_a', 'length_b', 'length_c',
1758                                     'angle_alpha', 'angle_beta', 'angle_gamma',
1759                                     'volume']):
1760                cellDict[key] = cellList[i]
1761                cellSigDict[key] = cellSig[i]
1762            return cellDict, cellSigDict
1763        except KeyError:
1764            cell = self.get_cell()
1765            return cell, {key: 0.0 for key in cell}
1766
1767    def export_CIF(self, outputname, quickmode=True):
1768        """Write this phase to a .cif file named outputname
1769
1770        :param str outputname: The name of the .cif file to write to
1771        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
1772        # This code is all taken from exports/G2export_CIF.py
1773        # Functions copied have the same names
1774        import GSASIImath as G2mth
1775        import GSASIImapvars as G2mv
1776        from exports import G2export_CIF as cif
1777
1778        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
1779        CIFname = os.path.splitext(self.proj.filename)[0]
1780        CIFname = os.path.split(CIFname)[1]
1781        CIFname = ''.join([c if ord(c) < 128 else ''
1782                           for c in CIFname.replace(' ', '_')])
1783        try:
1784            author = self.proj['Controls']['data'].get('Author','').strip()
1785        except KeyError:
1786            pass
1787        oneblock = True
1788
1789        covDict = self.proj['Covariance']['data']
1790        parmDict = dict(zip(covDict.get('varyList',[]),
1791                            covDict.get('variables',[])))
1792        sigDict = dict(zip(covDict.get('varyList',[]),
1793                           covDict.get('sig',[])))
1794
1795        if covDict.get('covMatrix') is not None:
1796            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1797                                              covDict['varyList'],
1798                                              parmDict))
1799
1800        with open(outputname, 'w') as fp:
1801            fp.write(' \n' + 70*'#' + '\n')
1802            cif.WriteCIFitem(fp, 'data_' + CIFname)
1803            # from exports.G2export_CIF.WritePhaseInfo
1804            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
1805            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
1806            # TODO get esds
1807            cellDict = self.get_cell()
1808            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
1809            names = ['length_a','length_b','length_c',
1810                     'angle_alpha','angle_beta ','angle_gamma',
1811                     'volume']
1812            for key, val in cellDict.items():
1813                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
1814
1815            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
1816                         self.data['General']['SGData']['SGSys'])
1817
1818            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
1819            # regularize capitalization and remove trailing H/R
1820            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
1821            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
1822
1823            # generate symmetry operations including centering and center of symmetry
1824            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
1825                self.data['General']['SGData'])
1826            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
1827            for i, op in enumerate(SymOpList,start=1):
1828                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
1829
1830            # TODO skipped histograms, exports/G2export_CIF.py:880
1831
1832            # report atom params
1833            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
1834                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
1835                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
1836            else:
1837                raise Exception("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
1838            # report cell contents
1839            cif.WriteComposition(fp, self.data, self.name, parmDict)
1840            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
1841                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
1842                raise NotImplementedError("only quickmode currently supported")
1843            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
1844                cif.WriteCIFitem(fp,'\n# Difference density results')
1845                MinMax = self.data['General']['Map']['minmax']
1846                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
1847                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
1848
1849
1850    def set_refinements(self, refs):
1851        """Sets the refinement parameter 'key' to the specification 'value'
1852
1853        :param dict refs: A dictionary of the parameters to be set. See
1854                          :ref:`Phase_parameters_table` for a description of
1855                          this dictionary.
1856
1857        :returns: None"""
1858        for key, value in refs.items():
1859            if key == "Cell":
1860                self.data['General']['Cell'][0] = value
1861
1862            elif key == "Atoms":
1863                for atomlabel, atomrefinement in value.items():
1864                    if atomlabel == 'all':
1865                        for atom in self.atoms():
1866                            atom.refinement_flags = atomrefinement
1867                    else:
1868                        atom = self.atom(atomlabel)
1869                        if atom is None:
1870                            raise ValueError("No such atom: " + atomlabel)
1871                        atom.refinement_flags = atomrefinement
1872
1873            elif key == "LeBail":
1874                hists = self.data['Histograms']
1875                for hname, hoptions in hists.items():
1876                    if 'LeBail' not in hoptions:
1877                        hoptions['newLeBail'] = bool(True)
1878                    hoptions['LeBail'] = bool(value)
1879            else:
1880                raise ValueError("Unknown key:", key)
1881
1882    def clear_refinements(self, refs):
1883        """Clears a given set of parameters.
1884
1885        :param dict refs: The parameters to clear"""
1886        for key, value in refs.items():
1887            if key == "Cell":
1888                self.data['General']['Cell'][0] = False
1889            elif key == "Atoms":
1890                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1891
1892                for atomlabel in value:
1893                    atom = self.atom(atomlabel)
1894                    # Set refinement to none
1895                    atom.refinement_flags = ' '
1896            elif key == "LeBail":
1897                hists = self.data['Histograms']
1898                for hname, hoptions in hists.items():
1899                    if 'LeBail' not in hoptions:
1900                        hoptions['newLeBail'] = True
1901                    hoptions['LeBail'] = False
1902            else:
1903                raise ValueError("Unknown key:", key)
1904
1905    def set_HAP_refinements(self, refs, histograms='all'):
1906        """Sets the given HAP refinement parameters between this phase and
1907        the given histograms
1908
1909        :param dict refs: A dictionary of the parameters to be set. See
1910                          :ref:`HAP_parameters_table` for a description of this
1911                          dictionary.
1912        :param histograms: Either 'all' (default) or a list of the histograms
1913            whose HAP parameters will be set with this phase. Histogram and phase
1914            must already be associated
1915
1916        :returns: None
1917        """
1918        if histograms == 'all':
1919            histograms = self.data['Histograms'].values()
1920        else:
1921            histograms = [self.data['Histograms'][h.name] for h in histograms]
1922
1923        for key, val in refs.items():
1924            for h in histograms:
1925                if key == 'Babinet':
1926                    try:
1927                        sets = list(val)
1928                    except ValueError:
1929                        sets = ['BabA', 'BabU']
1930                    for param in sets:
1931                        if param not in ['BabA', 'BabU']:
1932                            raise ValueError("Not sure what to do with" + param)
1933                        for hist in histograms:
1934                            hist['Babinet'][param][1] = True
1935                elif key == 'Extinction':
1936                    for h in histograms:
1937                        h['Extinction'][1] = bool(val)
1938                elif key == 'HStrain':
1939                    for h in histograms:
1940                        hist['HStrain'][1] = [bool(val) for p in hist['Hstrain'][0]]
1941                elif key == 'Mustrain':
1942                    for h in histograms:
1943                        mustrain = h['Mustrain']
1944                        newType = None
1945                        if isinstance(val, (unicode, str)):
1946                            if val in ['isotropic', 'uniaxial', 'generalized']:
1947                                newType = val
1948                            else:
1949                                raise ValueError("Not a Mustrain type: " + val)
1950                        elif isinstance(val, dict):
1951                            newType = val.get('type', None)
1952                            direction = val.get('direction', None)
1953
1954                        if newType:
1955                            mustrain[0] = newType
1956                            if newType == 'isotropic':
1957                                mustrain[2][0] = True
1958                                mustrain[5] = [False for p in mustrain[4]]
1959                            elif newType == 'uniaxial':
1960                                if 'refine' in val:
1961                                    types = val['refine']
1962                                    if isinstance(types, (unicode, str)):
1963                                        types = [types]
1964                                    elif isinstance(types, bool):
1965                                        mustrain[2][0] = types
1966                                        mustrain[2][1] = types
1967                                        types = []
1968                                    else:
1969                                        raise ValueError("Not sure what to do with: "
1970                                                         + str(types))
1971                                else:
1972                                    types = []
1973
1974                                for unitype in types:
1975                                    if unitype == 'equatorial':
1976                                        mustrain[2][0] = True
1977                                    elif unitype == 'axial':
1978                                        mustrain[2][1] = True
1979                                    else:
1980                                        msg = 'Invalid uniaxial mustrain type'
1981                                        raise ValueError(msg + ': ' + unitype)
1982                            else:  # newtype == 'generalized'
1983                                mustrain[2] = [False for p in mustrain[1]]
1984
1985                        if direction:
1986                            direction = [int(n) for n in direction]
1987                            if len(direction) != 3:
1988                                raise ValueError("Expected hkl, found", direction)
1989                            mustrain[3] = direction
1990                elif key == 'Pref.Ori.':
1991                    for h in histograms:
1992                        h['Pref.Ori.'][2] = bool(val)
1993                elif key == 'Show':
1994                    for h in histograms:
1995                        h['Show'] = bool(val)
1996                elif key == 'Size':
1997                    # TODO
1998                    raise NotImplementedError()
1999                elif key == 'Use':
2000                    for h in histograms:
2001                        h['Use'] = bool(val)
2002                elif key == 'Scale':
2003                    for h in histograms:
2004                        h['Scale'][1] = bool(val)
2005
2006    def clear_HAP_refinements(self, refs, histograms='all'):
2007        """Clears the given HAP refinement parameters between this phase and
2008        the given histograms
2009
2010        :param dict refs: A dictionary of the parameters to be cleared.
2011        :param histograms: Either 'all' (default) or a list of the histograms
2012            whose HAP parameters will be cleared with this phase. Histogram and
2013            phase must already be associated
2014
2015        :returns: None
2016        """
2017        if histograms == 'all':
2018            histograms = self.data['Histograms'].values()
2019        else:
2020            histograms = [self.data['Histograms'][h.name] for h in histograms]
2021
2022        for key, val in refs.items():
2023            for h in histograms:
2024                if key == 'Babinet':
2025                    try:
2026                        sets = list(val)
2027                    except ValueError:
2028                        sets = ['BabA', 'BabU']
2029                    for param in sets:
2030                        if param not in ['BabA', 'BabU']:
2031                            raise ValueError("Not sure what to do with" + param)
2032                        for hist in histograms:
2033                            hist['Babinet'][param][1] = False
2034                elif key == 'Extinction':
2035                    for h in histograms:
2036                        h['Extinction'][1] = False
2037                elif key == 'HStrain':
2038                    for h in histograms:
2039                        hist['HStrain'][1] = [False for p in hist['Hstrain'][0]]
2040                elif key == 'Mustrain':
2041                    for h in histograms:
2042                        mustrain = h['Mustrain']
2043                        mustrain[2] = [False for p in mustrain[1]]
2044                        mustrain[5] = [False for p in mustrain[4]]
2045                elif key == 'Pref.Ori.':
2046                    for h in histograms:
2047                        h['Pref.Ori.'][2] = False
2048                elif key == 'Show':
2049                    for h in histograms:
2050                        h['Show'] = False
2051                elif key == 'Size':
2052                    raise NotImplementedError()
2053                elif key == 'Use':
2054                    for h in histograms:
2055                        h['Use'] = False
2056                elif key == 'Scale':
2057                    for h in histograms:
2058                        h['Scale'][1] = False
2059
2060
2061##########################
2062# Command Line Interface #
2063##########################
2064# Each of these takes an argparse.Namespace object as their argument,
2065# representing the parsed command-line arguments for the relevant subcommand.
2066# The argument specification for each is in the subcommands dictionary (see
2067# below)
2068
2069
2070def create(args):
2071    """The create subcommand."""
2072    proj = G2Project(filename=args.filename)
2073
2074    hist_objs = []
2075    for h in args.histograms:
2076        hist_objs.append(proj.add_powder_histogram(h, args.iparams))
2077
2078    for p in args.phases:
2079        proj.add_phase(p, histograms=hist_objs)
2080
2081    proj.save()
2082
2083
2084def dump(args):
2085    """The dump subcommand"""
2086    if not args.histograms and not args.phases:
2087        args.raw = True
2088    if args.raw:
2089        import IPython.lib.pretty as pretty
2090
2091    for fname in args.files:
2092        if args.raw:
2093            proj, nameList = LoadDictFromProjFile(fname)
2094            print("file:", fname)
2095            print("names:", nameList)
2096            for key, val in proj.items():
2097                print(key, ":")
2098                pretty.pprint(val)
2099        else:
2100            proj = G2Project(fname)
2101            if args.histograms:
2102                hists = proj.histograms()
2103                for h in hists:
2104                    print(fname, "hist", h.id, h.name)
2105            if args.phases:
2106                phase_list = proj.phases()
2107                for p in phase_list:
2108                    print(fname, "phase", p.id, p.name)
2109
2110
2111def IPyBrowse(args):
2112    """Load a .gpx file and then open a IPython shell to browse it
2113    """
2114    for fname in args.files:
2115        proj, nameList = LoadDictFromProjFile(fname)
2116        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
2117        GSASIIpath.IPyBreak_base(msg)
2118        break
2119
2120
2121def refine(args):
2122    """The refine subcommand"""
2123    proj = G2Project(args.gpxfile)
2124    if args.refinements is None:
2125        proj.refine()
2126    else:
2127        import json
2128        with open(args.refinements) as refs:
2129            refs = json.load(refs)
2130        proj.do_refinements(refs['refinements'])
2131
2132
2133def seqrefine(args):
2134    """The seqrefine subcommand"""
2135    raise NotImplementedError("seqrefine is not yet implemented")
2136
2137
2138def export(args):
2139    """The export subcommand"""
2140    # Export phase as CIF to args.exportfile
2141    proj = G2Project(args.gpxfile)
2142    phase = proj.phase(args.phase)
2143    phase.export_CIF(args.exportfile)
2144
2145
2146def _args_kwargs(*args, **kwargs):
2147    return args, kwargs
2148
2149# A dictionary of the name of each subcommand, and a tuple
2150# of its associated function and a list of its arguments
2151# The arguments are passed directly to the add_argument() method
2152# of an argparse.ArgumentParser
2153subcommands = {"create":
2154               (create, [_args_kwargs('filename',
2155                                      help='the project file to create. should end in .gpx'),
2156
2157                         _args_kwargs('-i', '--iparams',
2158                                      help='instrument parameter file'),
2159
2160                         _args_kwargs('-d', '--histograms',
2161                                      nargs='+',
2162                                      help='list of datafiles to add as histograms'),
2163
2164                         _args_kwargs('-p', '--phases',
2165                                      nargs='+',
2166                                      help='list of phases to add. phases are '
2167                                           'automatically associated with all '
2168                                           'histograms given.')]),
2169
2170               "dump": (dump, [_args_kwargs('-d', '--histograms',
2171                                     action='store_true',
2172                                     help='list histograms in files, overrides --raw'),
2173
2174                               _args_kwargs('-p', '--phases',
2175                                            action='store_true',
2176                                            help='list phases in files, overrides --raw'),
2177
2178                               _args_kwargs('-r', '--raw',
2179                                      action='store_true', help='dump raw file contents, default'),
2180
2181                               _args_kwargs('files', nargs='+')]),
2182
2183               "refine":
2184               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
2185                         _args_kwargs('refinements',
2186                                      help='json file of refinements to apply. if not present'
2187                                           ' refines file as-is',
2188                                      default=None,
2189                                      nargs='?')]),
2190
2191               "seqrefine": (seqrefine, []),
2192               "export": (export, [_args_kwargs('gpxfile',
2193                                                help='the project file from which to export'),
2194                                   _args_kwargs('phase', help='identifier of phase to export'),
2195                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
2196               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
2197                                                   help='list of files to browse')])}
2198
2199
2200def main():
2201    '''The command line interface for GSASIIscriptable.
2202
2203    Executes one of the following subcommands:
2204
2205        * :func:`create`
2206        * :func:`dump`
2207        * :func:`refine`
2208        * :func:`seqrefine`
2209        * :func:`export`
2210        * :func:`browse`
2211
2212    .. seealso::
2213        :func:`create`
2214        :func:`dump`
2215        :func:`refine`
2216        :func:`seqrefine`
2217        :func:`export`
2218        :func:`browse`
2219    '''
2220    parser = argparse.ArgumentParser()
2221    subs = parser.add_subparsers()
2222
2223    # Create all of the specified subparsers
2224    for name, (func, args) in subcommands.items():
2225        new_parser = subs.add_parser(name)
2226        for listargs, kwds in args:
2227            new_parser.add_argument(*listargs, **kwds)
2228        new_parser.set_defaults(func=func)
2229
2230    # Parse and trigger subcommand
2231    result = parser.parse_args()
2232    result.func(result)
2233
2234    # argv = sys.argv
2235    # if len(argv) > 1 and argv[1] in subcommands:
2236    #     subcommands[argv[1]](*argv[2:])
2237    # elif len(argv) == 1 or argv[1] in ('help', '--help', '-h'):
2238    #     # TODO print usage
2239    #     subcommand_names = ' | '.join(sorted(subcommands.keys()))
2240    #     print("USAGE: {} [ {} ] ...".format(argv[0], subcommand_names))
2241    # else:
2242    #     print("Unknown subcommand: {}".format(argv[1]))
2243    #     print("Available subcommands:")
2244    #     for name in sorted(subcommands.keys()):
2245    #         print("\t{}".format(name))
2246    #     sys.exit(-1)
2247    # sys.exit(0)
2248
2249if __name__ == '__main__':
2250    main()
2251
2252
2253    # from refinements.py
2254#      USAGE = """USAGE: {} datafile instparams phasefile projectname refinements
2255
2256# datafile:    Input powder data
2257# intparams:   Corresponding instrument parameter file
2258# phasefile:   Phase to refine against data
2259# projectname: Project file to be created, should end in .gpx
2260# refinements: JSON file of refinements to be executed
2261# """
2262#     try:
2263#         datafile, instprm, phasefile, projectname, refinements = sys.argv[1:]
2264#     except ValueError:
2265#         print(USAGE.format(sys.argv[0]))
2266#         sys.exit(-1)
2267
2268#     try:
2269#         with open(refinements) as f:
2270#             refinements = json.load(f)
2271#     except IOError:
2272#         print("No such refinements file: {}".format(refinements))
2273
2274#     print("Creating project file \"{}\"...".format(projectname))
2275#     proj = G2Project(filename=projectname)
2276#     # Add the histogram
2277#     hist = proj.add_powder_histogram(datafile, instprm)
2278#     # Add the phase, and associate it with the histogram
2279#     proj.add_phase(phasefile, histograms=[hist.name])
2280
2281#     proj.do_refinements(refinements['refinements'])
2282#     proj.save()
2283
2284
2285    # from gpx_dumper
2286    # import IPython.lib.pretty as pretty
2287    # proj, nameList = LoadDictFromProjFile(sys.argv[1])
2288    # print("names:", nameList)
2289    # for key, val in proj.items():
2290    #     print(key, ":", sep='')
2291    #     pretty.pprint(val)
Note: See TracBrowser for help on using the repository browser.