Ignore:
Timestamp:
Aug 11, 2017 5:34:54 PM (4 years ago)
Author:
toby
Message:

make the two frame version the trunk as we hit version 3000

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/GSASIIscriptable.py

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