#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
"""
Classes and routines defined in :mod:`GSASIIscriptable` follow.
A script will create one or more :class:`G2Project` objects by reading
a GSAS-II project (.gpx) file or creating a new one and will then
perform actions such as adding a histogram (method :meth:`G2Project.add_powder_histogram`),
adding a phase (method :meth:`G2Project.add_phase`),
or setting parameters and performing a refinement
(method :meth:`G2Project.do_refinements`).
To change settings within histograms, images and phases, one usually needs to use
methods inside :class:`G2PwdrData`, :class:`G2Image` or :class:`G2Phase`.
"""
# Note that documentation for GSASIIscriptable.py has been moved
# to file docs/source/GSASIIscriptable.rst
#============================================================================
# Notes for adding a new object type
# 1) add a new object class (e.g. G2PDF)
# 2) add the wrapper into G2Project (e.g. _pdfs, pdf, pdfs)
# 3) add a new method to add the object into a project (G2Project.add_PDF)
# Document: (in ../docs/source/GSASIIscriptable.rst)
# 4) add to documentation in section :class:`G2Project`
# 5) add a new documentation section for the new class
#============================================================================
from __future__ import division, print_function
import argparse
import os.path as ospath
import sys
import platform
import pickle
import copy
import os
import random as ran
import numpy as np
import numpy.ma as ma
import GSASIIpath
GSASIIpath.SetBinaryPath(True) # for now, this is needed before some of these modules can be imported
import GSASIIobj as G2obj
import GSASIIpwd as G2pwd
import GSASIIstrMain as G2strMain
import GSASIIstrIO as G2stIO
import GSASIIspc as G2spc
import GSASIIElem as G2elem
import GSASIIfiles as G2fil
import GSASIIimage as G2img
import GSASIIlattice as G2lat
import GSASIImapvars as G2mv
# Delay imports loading to not slow down small scripts that don't need them
Readers = {'Pwdr':[], 'Phase':[], 'Image':[]}
'''Readers by reader type'''
exportersByExtension = {}
'''Specifies the list of extensions that are supported for Powder data export'''
npsind = lambda x: np.sin(x*np.pi/180.)
[docs]
def SetPrintLevel(level):
'''Set the level of output from calls to :func:`GSASIIfiles.G2Print`,
which should be used in place of print() where possible. This is a
wrapper for :func:`GSASIIfiles.G2SetPrintLevel` so that this routine is
documented here.
:param str level: a string used to set the print level, which may be
'all', 'warn', 'error' or 'none'.
Note that capitalization and extra letters in level are ignored, so
'Warn', 'warnings', etc. will all set the mode to 'warn'
'''
G2fil.G2SetPrintLevel(level)
global printLevel
for mode in 'all', 'warn', 'error', 'none':
if mode in level.lower():
printLevel = mode
return
[docs]
def installScriptingShortcut():
'''Creates a file named G2script in the current Python site-packages directory.
This is equivalent to the "Install GSASIIscriptable shortcut" command in the GUI's
File menu. Once this is done, a shortcut for calling GSASIIscriptable is created,
where the command:
>>> import G2script as G2sc
will provide access to GSASIIscriptable without changing the sys.path; also see
:ref:`ScriptingShortcut`.
Note that this only affects the current Python installation. If more than one
Python installation will be used with GSAS-II (for example because different
conda environments are used), this command should be called from within each
Python environment.
If more than one GSAS-II installation will be used with a Python installation,
this shortcut can only be used with one of them.
'''
f = GSASIIpath.makeScriptShortcut()
if f:
G2fil.G2Print(f'success creating {f}')
else:
raise G2ScriptException('error creating G2script')
[docs]
def ShowVersions():
'''Show the versions all of required Python packages, etc.
'''
out = ''
pkgList = [('Python',None), ('numpy',np)]
try:
import scipy
pkgList.append(('scipy',scipy))
except:
pass
try:
import IPython
pkgList.append(('IPython',IPython))
except:
pass
for s,m in pkgList:
msg = ''
if s == 'Python':
pkgver = platform.python_version()
#prefix = ''
msg = f"from {format(sys.executable)}"
else:
pkgver = m.__version__
#prefix = 'Package '
out += f" {s:12s}{pkgver}: {msg}\n"
out += GSASIIpath.getG2VersionInfo()
out += "\n\n"
try:
out += f"GSAS-II location: {GSASIIpath.path2GSAS2}\n"
except:
out += "GSAS-II location: not set\n"
try:
out += f"Binary location: {GSASIIpath.binaryPath}\n"
except:
out += "Binary location: not found\n"
return out
[docs]
def LoadG2fil():
'''Setup GSAS-II importers.
Delay importing this module when possible, it is slow.
Multiple calls are not. Only the first does anything.
'''
if len(Readers['Pwdr']) > 0: return
# initialize imports
Readers['Pwdr'] = G2fil.LoadImportRoutines("pwd", "Powder_Data")
Readers['Phase'] = G2fil.LoadImportRoutines("phase", "Phase")
Readers['Image'] = G2fil.LoadImportRoutines("img", "Image")
Readers['HKLF'] = G2fil.LoadImportRoutines('sfact','Struct_Factor')
# initialize exports
for obj in G2fil.LoadExportRoutines(None):
try:
obj.Writer
except AttributeError:
continue
for typ in obj.exporttype:
if typ not in exportersByExtension:
exportersByExtension[typ] = {obj.extension:obj}
elif obj.extension in exportersByExtension[typ]:
if type(exportersByExtension[typ][obj.extension]) is list:
exportersByExtension[typ][obj.extension].append(obj)
else:
exportersByExtension[typ][obj.extension] = [
exportersByExtension[typ][obj.extension],
obj]
else:
exportersByExtension[typ][obj.extension] = obj
[docs]
def LoadDictFromProjFile(ProjFile):
'''Read a GSAS-II project file and load items to dictionary
:param str ProjFile: GSAS-II project (name.gpx) full file name
:returns: Project,nameList, where
* Project (dict) is a representation of gpx file following the GSAS-II tree structure
for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
* nameList (list) has names of main tree entries & subentries used to reconstruct project file
Example for fap.gpx::
Project = { #NB:dict order is not tree order
'Phases':{'data':None,'fap':{phase dict}},
'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
'Rigid bodies':{'data': {rigid body dict}},
'Covariance':{'data':{covariance data dict}},
'Controls':{'data':{controls data dict}},
'Notebook':{'data':[notebook list]},
'Restraints':{'data':{restraint data dict}},
'Constraints':{'data':{constraint data dict}}]
}
nameList = [ #NB: reproduces tree order
['Notebook',],
['Controls',],
['Covariance',],
['Constraints',],
['Restraints',],
['Rigid bodies',],
['PWDR FAP.XRA Bank 1',
'Comments',
'Limits',
'Background',
'Instrument Parameters',
'Sample Parameters',
'Peak List',
'Index Peak List',
'Unit Cells List',
'Reflection Lists'],
['Phases', 'fap']
]
'''
# Let IOError be thrown if file does not exist
if not ospath.exists(ProjFile):
G2fil.G2Print ('\n*** Error attempt to open project file that does not exist: \n {}'.
format(ProjFile))
raise IOError('GPX file {} does not exist'.format(ProjFile))
try:
Project, nameList = G2stIO.GetFullGPX(ProjFile)
except Exception as msg:
raise IOError(msg)
return Project,nameList
[docs]
def SaveDictToProjFile(Project,nameList,ProjFile):
'''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
:param dict Project: representation of gpx file following the GSAS-II
tree structure as described for LoadDictFromProjFile
:param list nameList: names of main tree entries & subentries used to reconstruct project file
:param str ProjFile: full file name for output project.gpx file (including extension)
'''
file = open(ProjFile,'wb')
try:
for name in nameList:
data = []
item = Project[name[0]]
data.append([name[0],item['data']])
for item2 in name[1:]:
data.append([item2,item[item2]])
pickle.dump(data,file,1)
finally:
file.close()
G2fil.G2Print('gpx file saved as %s'%ProjFile)
[docs]
def PreSetup(data):
'''Create part of an initial (empty) phase dictionary
from GSASIIphsGUI.py, near end of UpdatePhaseData
Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
'''
if 'RBModels' not in data:
data['RBModels'] = {}
if 'MCSA' not in data:
data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
if 'dict' in str(type(data['MCSA']['Results'])):
data['MCSA']['Results'] = []
if 'Modulated' not in data['General']:
data['General']['Modulated'] = False
# if 'modulated' in data['General']['Type']:
# data['General']['Modulated'] = True
# data['General']['Type'] = 'nuclear'
[docs]
def SetupGeneral(data, dirname):
'''Initialize phase data.
'''
try:
G2elem.SetupGeneral(data,dirname)
except ValueError as msg:
raise G2ScriptException(msg)
[docs]
def make_empty_project(author=None, filename=None):
"""Creates an dictionary in the style of GSASIIscriptable, for an empty
project.
If no author name or filename is supplied, 'no name' and
<current dir>/test_output.gpx are used , respectively.
Returns: project dictionary, name list
Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
"""
if not filename:
filename = 'test_output.gpx'
filename = os.path.abspath(filename)
LoadG2fil()
controls_data = dict(G2obj.DefaultControls)
if author:
controls_data['Author'] = author
output = {'Constraints': {'data': {'Hist': [], 'HAP': [], 'Phase': [],
'Global': [],
'_seqmode':'auto-wildcard',
'_seqhist':0}},
'Controls': {'data': controls_data},
'Covariance': {'data': {}},
'Notebook': {'data': ['']},
'Restraints': {'data': {}},
'Rigid bodies': {'data':
{'Vector':{'AtInfo':{}},'Residue':{'AtInfo':{}},
"Spin": {},
'RBIds':{'Vector':[], 'Residue':[],'Spin':[]}}}
}
names = [['Notebook'], ['Controls'], ['Covariance'],
['Constraints'], ['Restraints'], ['Rigid bodies']]
return output, names
[docs]
def GenerateReflections(spcGrp,cell,Qmax=None,dmin=None,TTmax=None,wave=None):
"""Generates the crystallographically unique powder diffraction reflections
for a lattice and space group (see :func:`GSASIIlattice.GenHLaue`).
:param str spcGrp: A GSAS-II formatted space group (with spaces between
axial fields, e.g. 'P 21 21 21' or 'P 42/m m c'). Note that non-standard
space groups, such as 'P 21/n' or 'F -1' are allowed (see
:func:`GSASIIspc.SpcGroup`).
:param list cell: A list/tuple with six unit cell constants,
(a, b, c, alpha, beta, gamma) with values in Angstroms/degrees.
Note that the cell constants are not checked for consistency
with the space group.
:param float Qmax: Reflections up to this Q value are computed
(do not use with dmin or TTmax)
:param float dmin: Reflections with d-space above this value are computed
(do not use with Qmax or TTmax)
:param float TTmax: Reflections up to this 2-theta value are computed
(do not use with dmin or Qmax, use of wave is required.)
:param float wave: wavelength in Angstroms for use with TTmax (ignored
otherwise.)
:returns: a list of reflections, where each reflection contains four items:
h, k, l, d, where d is the d-space (Angstroms)
Example:
>>> import os,sys
>>> sys.path.insert(0,'/Users/toby/software/G2/GSASII')
>>> import GSASIIscriptable as G2sc
GSAS-II binary directory: /Users/toby/software/G2/GSASII/bin
17 values read from config file /Users/toby/software/G2/GSASII/config.py
>>> refs = G2sc.GenerateReflections('P 1',
... (5.,6.,7.,90.,90.,90),
... TTmax=20,wave=1)
>>> for r in refs: print(r)
...
[0, 0, 1, 7.0]
[0, 1, 0, 6.0]
[1, 0, 0, 5.0]
[0, 1, 1, 4.55553961419178]
[0, 1, -1, 4.55553961419178]
[1, 0, 1, 4.068667356033675]
[1, 0, -1, 4.068667356033674]
[1, 1, 0, 3.8411063979868794]
[1, -1, 0, 3.8411063979868794]
"""
if len(cell) != 6:
raise G2ScriptException("GenerateReflections: Invalid unit cell:" + str(cell))
opts = (Qmax is not None) + (dmin is not None) + (TTmax is not None)
if Qmax:
dmin = 2 * np.pi / Qmax
#print('Q,d',Qmax,dmin)
elif TTmax and wave is None:
raise G2ScriptException("GenerateReflections: specify a wavelength with TTmax")
elif TTmax:
dmin = wave / (2.0 * np.sin(np.pi*TTmax/360.))
#print('2theta,d',TTmax,dmin)
if opts != 1:
raise G2ScriptException("GenerateReflections: specify one Qmax, dmin or TTmax")
err,SGData = G2spc.SpcGroup(spcGrp)
if err != 0:
print('GenerateReflections space group error:',G2spc.SGErrors(err))
raise G2ScriptException("GenerateReflections: Invalid space group: " + str(spcGrp))
A = G2lat.cell2A(cell)
return G2lat.GenHLaue(dmin,SGData,A)
[docs]
class G2ImportException(Exception):
pass
[docs]
class G2ScriptException(Exception):
pass
[docs]
def import_generic(filename, readerlist, fmthint=None, bank=None):
"""Attempt to import a filename, using a list of reader objects.
Returns the first reader object which worked."""
# Translated from OnImportGeneric method in GSASII.py
primaryReaders, secondaryReaders = [], []
for reader in readerlist:
if fmthint is not None and fmthint not in reader.formatName: continue
flag = reader.ExtensionValidator(filename)
if flag is None:
secondaryReaders.append(reader)
elif flag:
primaryReaders.append(reader)
if not secondaryReaders and not primaryReaders:
raise G2ImportException("Could not read file: ", filename)
with open(filename, 'r'):
rd_list = []
for rd in primaryReaders + secondaryReaders:
# Initialize reader
rd.selections = []
if bank is None:
rd.selections = []
else:
try:
rd.selections = [i-1 for i in bank]
except TypeError:
rd.selections = [bank-1]
rd.dnames = []
rd.ReInitialize()
# Rewind file
rd.errors = ""
if not rd.ContentsValidator(filename):
# Report error
G2fil.G2Print("Warning: File {} has a validation error, continuing".format(filename))
#if len(rd.selections) > 1:
# raise G2ImportException("File {} has {} banks. Specify which bank to read with databank param."
# .format(filename,len(rd.selections)))
block = 0
rdbuffer = {}
repeat = True
while repeat:
repeat = False
block += 1
rd.objname = os.path.basename(filename)
try:
flag = rd.Reader(filename,buffer=rdbuffer, blocknum=block)
except Exception as msg:
if GSASIIpath.GetConfigValue('debug'):
print('Reader exception',msg)
flag = False
if flag:
# Omitting image loading special cases
rd.readfilename = filename
rd_list.append(copy.deepcopy(rd))
repeat = rd.repeat
else:
G2fil.G2Print("Warning: {} Reader failed to read {}".format(rd.formatName,filename))
if rd_list:
if rd.warnings:
G2fil.G2Print("Read warning by", rd.formatName, "reader:",
rd.warnings)
elif bank is None:
G2fil.G2Print("{} read by Reader {}"
.format(filename,rd.formatName))
else:
G2fil.G2Print("{} block # {} read by Reader {}"
.format(filename,bank,rd.formatName))
return rd_list
raise G2ImportException("No reader could read file: " + filename)
[docs]
def load_iprms(instfile, reader, bank=None):
"""Loads instrument parameters from a file, and edits the
given reader.
Returns a 2-tuple of (Iparm1, Iparm2) parameters
"""
LoadG2fil()
ext = os.path.splitext(instfile)[1]
if ext.lower() == '.instprm':
# New GSAS File, load appropriate bank
with open(instfile) as f:
lines = f.readlines()
if bank is None:
bank = reader.powderentry[2]
nbank,iparms = G2fil.ReadInstprm(lines, bank, reader.Sample)
reader.instfile = instfile
reader.instmsg = '{} (G2 fmt) bank {}'.format(instfile,bank)
return iparms
elif ext.lower() not in ('.prm', '.inst', '.ins'):
raise ValueError('Expected .prm file, found: ', instfile)
# It's an old GSAS file, load appropriately
Iparm = {}
with open(instfile, 'r') as fp:
for line in fp:
if '#' in line:
continue
Iparm[line[:12]] = line[12:-1]
ibanks = int(Iparm.get('INS BANK ', '1').strip())
if bank is not None:
# pull out requested bank # bank from the data, and change the bank to 1
Iparm,IparmC = {},Iparm
for key in IparmC:
if 'INS' not in key[:3]: continue #skip around rubbish lines in some old iparm
if key[4:6] == " ":
Iparm[key] = IparmC[key]
elif int(key[4:6].strip()) == bank:
Iparm[key[:4]+' 1'+key[6:]] = IparmC[key]
reader.instbank = bank
elif ibanks == 1:
reader.instbank = 1
else:
raise G2ImportException("Instrument parameter file has {} banks, select one with instbank param."
.format(ibanks))
reader.powderentry[2] = 1
reader.instfile = instfile
reader.instmsg = '{} bank {}'.format(instfile,reader.instbank)
return G2fil.SetPowderInstParms(Iparm, reader)
[docs]
def load_pwd_from_reader(reader, instprm, existingnames=[],bank=None):
"""Loads powder data from a reader object, and assembles it into a GSASII data tree.
:returns: (name, tree) - 2-tuple of the histogram name (str), and data
Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
"""
HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
HistName = G2obj.MakeUniqueLabel(HistName, existingnames)
try:
Iparm1, Iparm2 = instprm
except ValueError:
Iparm1, Iparm2 = load_iprms(instprm, reader, bank=bank)
G2fil.G2Print('Instrument parameters read:',reader.instmsg)
except TypeError: # instprm is None, get iparms from reader
Iparm1, Iparm2 = reader.pwdparms['Instrument Parameters']
if 'T' in Iparm1['Type'][0]:
if not reader.clockWd and reader.GSAS:
reader.powderdata[0] *= 100. #put back the CW centideg correction
cw = np.diff(reader.powderdata[0])
reader.powderdata[0] = reader.powderdata[0][:-1]+cw/2.
if reader.GSAS: #NB: old GSAS wanted intensities*CW even if normalized!
npts = min(len(reader.powderdata[0]),len(reader.powderdata[1]),len(cw))
reader.powderdata[1] = reader.powderdata[1][:npts]/cw[:npts]
reader.powderdata[2] = reader.powderdata[2][:npts]*cw[:npts]**2 #1/var=w at this point
else: #NB: from topas/fullprof type files
reader.powderdata[1] = reader.powderdata[1][:-1]
reader.powderdata[2] = reader.powderdata[2][:-1]
if 'Itype' in Iparm2:
Ibeg = np.searchsorted(reader.powderdata[0],Iparm2['Tminmax'][0])
Ifin = np.searchsorted(reader.powderdata[0],Iparm2['Tminmax'][1])
reader.powderdata[0] = reader.powderdata[0][Ibeg:Ifin]
YI,WYI = G2pwd.calcIncident(Iparm2,reader.powderdata[0])
reader.powderdata[1] = reader.powderdata[1][Ibeg:Ifin]/YI
var = 1./reader.powderdata[2][Ibeg:Ifin]
var += WYI*reader.powderdata[1]**2
var /= YI**2
reader.powderdata[2] = 1./var
reader.powderdata[1] = np.where(np.isinf(reader.powderdata[1]),0.,reader.powderdata[1])
reader.powderdata[3] = np.zeros_like(reader.powderdata[0])
reader.powderdata[4] = np.zeros_like(reader.powderdata[0])
reader.powderdata[5] = np.zeros_like(reader.powderdata[0])
Ymin = np.min(reader.powderdata[1])
Ymax = np.max(reader.powderdata[1])
valuesdict = {'wtFactor': 1.0,
'Dummy': False,
'ranId': ran.randint(0, sys.maxsize),
'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
'Yminmax': [Ymin, Ymax]}
# apply user-supplied corrections to powder data
if 'CorrectionCode' in Iparm1:
G2fil.G2Print('Applying corrections from instprm file')
corr = Iparm1['CorrectionCode'][0]
try:
exec(corr)
G2fil.G2Print('done')
except Exception as err:
print('error: {}'.format(err))
print('with commands -------------------')
print(corr)
print('---------------------------------')
finally:
del Iparm1['CorrectionCode']
reader.Sample['ranId'] = valuesdict['ranId']
# Ending keys:
# ['Reflection Lists',
# 'Limits',
# 'data',
# 'Index Peak List',
# 'Comments',
# 'Unit Cells List',
# 'Sample Parameters',
# 'Peak List',
# 'Background',
# 'Instrument Parameters']
Tmin = np.min(reader.powderdata[0])
Tmax = np.max(reader.powderdata[0])
Tmin1 = Tmin
if 'NT' in Iparm1['Type'][0] and G2lat.Pos2dsp(Iparm1,Tmin) < 0.4:
Tmin1 = G2lat.Dsp2pos(Iparm1,0.4)
default_background = [['chebyschev-1', False, 3, 1.0, 0.0, 0.0],
{'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0,
'peaksList': [],'background PWDR':['',1.0,False]}]
output_dict = {'Reflection Lists': {},
'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin1, Tmax]]),
'data': [valuesdict, reader.powderdata, HistName],
'Index Peak List': [[], []],
'Comments': reader.comments,
'Unit Cells List': [],
'Sample Parameters': reader.Sample,
'Peak List': {'peaks': [], 'sigDict': {}},
'Background': reader.pwdparms.get('Background', default_background),
'Instrument Parameters': [Iparm1, Iparm2],
}
names = ['Comments',
'Limits',
'Background',
'Instrument Parameters',
'Sample Parameters',
'Peak List',
'Index Peak List',
'Unit Cells List',
'Reflection Lists']
# TODO controls?? GSASII.py:1664-7
return HistName, [HistName] + names, output_dict
def _deep_copy_into(from_, into):
"""Helper function for reloading .gpx file. See G2Project.reload()
:author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
"""
if isinstance(from_, dict) and isinstance(into, dict):
combined_keys = set(from_.keys()).union(into.keys())
for key in combined_keys:
if key in from_ and key in into:
both_dicts = (isinstance(from_[key], dict)
and isinstance(into[key], dict))
both_lists = (isinstance(from_[key], list)
and isinstance(into[key], list))
if both_dicts or both_lists:
_deep_copy_into(from_[key], into[key])
else:
into[key] = from_[key]
elif key in from_:
into[key] = from_[key]
else: # key in into
del into[key]
elif isinstance(from_, list) and isinstance(into, list):
if len(from_) == len(into):
for i in range(len(from_)):
both_dicts = (isinstance(from_[i], dict)
and isinstance(into[i], dict))
both_lists = (isinstance(from_[i], list)
and isinstance(into[i], list))
if both_dicts or both_lists:
_deep_copy_into(from_[i], into[i])
else:
into[i] = from_[i]
else:
into[:] = from_
def _getCorrImage(ImageReaderlist,proj,imageRef):
'''Gets image & applies dark, background & flat background corrections.
based on :func:`GSASIIimgGUI.GetImageZ`. Expected to be for internal
use only.
:param list ImageReaderlist: list of Reader objects for images
:param object proj: references a :class:`G2Project` project
:param imageRef: A reference to the desired image in the project.
Either the Image tree name (str), the image's index (int) or
a image object (:class:`G2Image`)
:return: array sumImg: corrected image for background/dark/flat back
'''
ImgObj = proj.image(imageRef)
Controls = ImgObj.data['Image Controls']
formatName = Controls.get('formatName','')
imagefile = ImgObj.data['data'][1]
if isinstance(imagefile, tuple) or isinstance(imagefile, list):
imagefile, ImageTag = imagefile # fix for multiimage files
else:
ImageTag = None # single-image file
sumImg = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
if sumImg is None:
return []
sumImg = np.array(sumImg,dtype=np.int32)
darkImg = False
if 'dark image' in Controls:
darkImg,darkScale = Controls['dark image']
if darkImg:
dImgObj = proj.image(darkImg)
formatName = dImgObj.data['Image Controls'].get('formatName','')
imagefile = dImgObj.data['data'][1]
if type(imagefile) is tuple:
imagefile,ImageTag = imagefile
darkImage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
if darkImg is None:
raise Exception('Error reading dark image {}'.format(imagefile))
sumImg += np.array(darkImage*darkScale,dtype=np.int32)
if 'background image' in Controls:
backImg,backScale = Controls['background image']
if backImg: #ignores any transmission effect in the background image
bImgObj = proj.image(backImg)
formatName = bImgObj.data['Image Controls'].get('formatName','')
imagefile = bImgObj.data['data'][1]
ImageTag = None # fix this for multiimage files
backImage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
if backImage is None:
raise Exception('Error reading background image {}'.format(imagefile))
if darkImg:
backImage += np.array(darkImage*darkScale/backScale,dtype=np.int32)
else:
sumImg += np.array(backImage*backScale,dtype=np.int32)
if 'Gain map' in Controls:
gainMap = Controls['Gain map']
if gainMap:
gImgObj = proj.image(gainMap)
formatName = gImgObj.data['Image Controls'].get('formatName','')
imagefile = gImgObj.data['data'][1]
ImageTag = None # fix this for multiimage files
GMimage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
if GMimage is None:
raise Exception('Error reading Gain map image {}'.format(imagefile))
sumImg = sumImg*GMimage/1000
sumImg -= int(Controls.get('Flat Bkg',0))
Imax = np.max(sumImg)
Controls['range'] = [(0,Imax),[0,Imax]]
return np.asarray(sumImg,dtype=np.int32)
def _constr_type(var):
'''returns the constraint type based on phase/histogram use
in a variable
'''
if var.histogram and var.phase:
return 'HAP'
elif var.phase:
return 'Phase'
elif var.histogram:
return 'Hist'
else:
return 'Global'
[docs]
class G2ObjectWrapper(object):
"""Base class for all GSAS-II object wrappers.
The underlying GSAS-II format can be accessed as `wrapper.data`. A number
of overrides are implemented so that the wrapper behaves like a dictionary.
Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
"""
def __init__(self, datadict):
self.data = datadict
def __getitem__(self, key):
return self.data[key]
def __setitem__(self, key, value):
self.data[key] = value
def __contains__(self, key):
return key in self.data
def get(self, k, d=None):
return self.data.get(k, d)
def keys(self):
return self.data.keys()
def values(self):
return self.data.values()
def items(self):
return self.data.items()
[docs]
class G2Project(G2ObjectWrapper):
"""Represents an entire GSAS-II project. The object contains these
class variables:
* G2Project.filename: contains the .gpx filename
* G2Project.names: contains the contents of the project "tree" as a list
of lists. Each top-level entry in the tree is an item in the list. The
name of the top-level item is the first item in the inner list. Children
of that item, if any, are subsequent entries in that list.
* G2Project.data: contains the entire project as a dict. The keys
for the dict are the top-level names in the project tree (initial items
in the G2Project.names inner lists) and each top-level item is stored
as a dict.
* The contents of Top-level entries will be found in the item
named 'data', as an example, ``G2Project.data['Notebook']['data']``
* The contents of child entries will be found in the item
using the names of the parent and child, for example
``G2Project.data['Phases']['NaCl']``
:param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
creates an empty project.
:param str author: Author's name (not yet implemented)
:param str newgpx: The filename the project should be saved to in
the future. If both newgpx and gpxfile are present, the project is
loaded from the file named by gpxfile and then when saved will
be written to the file named by newgpx.
:param str filename: To be deprecated. Serves the same function as newgpx,
which has a somewhat more clear name.
(Do not specify both newgpx and filename).
There are two ways to initialize this object:
>>> # Load an existing project file
>>> proj = G2Project('filename.gpx')
>>> # Create a new project
>>> proj = G2Project(newgpx='new_file.gpx')
Histograms can be accessed easily.
>>> # By name
>>> hist = proj.histogram('PWDR my-histogram-name')
>>> # Or by index
>>> hist = proj.histogram(0)
>>> assert hist.id == 0
>>> # Or by random id
>>> assert hist == proj.histogram(hist.ranId)
Phases can be accessed the same way.
>>> phase = proj.phase('name of phase')
New data can also be loaded via :meth:`~G2Project.add_phase` and
:meth:`~G2Project.add_powder_histogram`.
>>> hist = proj.add_powder_histogram('some_data_file.chi',
'instrument_parameters.prm')
>>> phase = proj.add_phase('my_phase.cif', histograms=[hist])
Parameters for Rietveld refinement can be turned on and off at the project level
as well as described in
:meth:`~G2Project.set_refinement`, :meth:`~G2Project.iter_refinements` and
:meth:`~G2Project.do_refinements`.
"""
def __init__(self, gpxfile=None, author=None, filename=None, newgpx=None):
if filename is not None and newgpx is not None:
raise G2ScriptException('Do not use filename and newgpx together')
elif filename is not None:
G2fil.G2Print("Warning - recommending use of parameter newgpx rather than filename\n\twhen creating a G2Project")
elif newgpx:
filename = newgpx
if not gpxfile:
filename = os.path.abspath(os.path.expanduser(filename))
self.filename = filename
self.data, self.names = make_empty_project(author=author, filename=filename)
elif os.path.exists(os.path.expanduser(gpxfile)):
# TODO set author
self.data, self.names = LoadDictFromProjFile(gpxfile)
self.update_ids()
if filename:
filename = os.path.abspath(os.path.expanduser(filename))
dr = os.path.split(filename)[0]
if not os.path.exists(dr):
raise Exception("Directory {} for filename/newgpx does not exist".format(dr))
self.filename = filename
else:
self.filename = os.path.abspath(os.path.expanduser(gpxfile))
else:
raise ValueError("Not sure what to do with gpxfile {}. Does not exist?".format(gpxfile))
[docs]
@classmethod
def from_dict_and_names(cls, gpxdict, names, filename=None):
"""Creates a :class:`G2Project` directly from
a dictionary and a list of names. If in doubt, do not use this.
:returns: a :class:`G2Project`
"""
out = cls()
if filename:
filename = os.path.abspath(os.path.expanduser(filename))
out.filename = filename
gpxdict['Controls']['data']['LastSavedAs'] = filename
else:
try:
out.filename = gpxdict['Controls']['data']['LastSavedAs']
except KeyError:
out.filename = None
out.data = gpxdict
out.names = names
[docs]
def save(self, filename=None):
"""Saves the project, either to the current filename, or to a new file.
Updates self.filename if a new filename provided"""
controls_data = self.data['Controls']['data']
# place info in .gpx about GSAS-II
# report versions of modules; ignore what has not been used
modList = [np]
try:
modList += [mpl]
except NameError:
pass
try:
modList += [sp]
except NameError:
pass
controls_data['PythonVersions'] = G2fil.get_python_versions(modList)
# G2 version info
controls_data['LastSavedUsing'] = str(GSASIIpath.GetVersionNumber())
if GSASIIpath.HowIsG2Installed().startswith('git'):
try:
g2repo = GSASIIpath.openGitRepo(GSASIIpath.path2GSAS2)
commit = g2repo.head.commit
controls_data['LastSavedUsing'] += f" git {commit.hexsha[:6]} script"
except:
pass
# .gpx name
if filename:
filename = os.path.abspath(os.path.expanduser(filename))
controls_data['LastSavedAs'] = filename
self.filename = filename
elif not self.filename:
raise AttributeError("No file name to save to")
SaveDictToProjFile(self.data, self.names, self.filename)
[docs]
def add_powder_histogram(self, datafile, iparams=None, phases=[],
fmthint=None,
databank=None, instbank=None, multiple=False):
"""Loads a powder data histogram or multiple powder histograms
into the project.
Note that the data type (x-ray/CW neutron/TOF) for the histogram
will be set from the instrument parameter file. The instrument
geometry is assumed to be Debye-Scherrer except for
dual-wavelength x-ray, where Bragg-Brentano is assumed.
:param str datafile: A filename with the powder data file to read.
Note that in unix fashion, "~" can be used to indicate the
home directory (e.g. ~/G2data/data.fxye).
:param str iparams: A filenme for an instrument parameters file,
or a pair of instrument parameter dicts from :func:`load_iprms`.
This may be omitted for readers that provide the instrument
parameters in the file. (Only a few importers do this.)
:param list phases: A list of phases to link to the new histogram,
phases can be references by object, name, rId or number.
Alternately, use 'all' to link to all phases in the project.
:param str fmthint: If specified, only importers where the format name
(reader.formatName, as shown in Import menu) contains the
supplied string will be tried as importers. If not specified, all
importers consistent with the file extension will be tried
(equivalent to "guess format" in menu).
:param int databank: Specifies a dataset number to read, if file contains
more than set of data. This should be 1 to read the first bank in
the file (etc.) regardless of the number on the Bank line, etc.
Default is None which means the first dataset in the file is read.
When multiple is True, optionally a list of dataset numbers can
be supplied here.
:param int instbank: Specifies an instrument parameter set to read, if
the instrument parameter file contains more than set of parameters.
This will match the INS # in an GSAS type file so it will typically
be 1 to read the first parameter set in the file (etc.)
Default is None which means there should only be one parameter set
in the file.
:param bool multiple: If False (default) only one dataset is read, but if
specified as True, all selected banks of data (see databank)
are read in.
:returns: A :class:`G2PwdrData` object representing
the histogram, or if multiple is True, a list of :class:`G2PwdrData`
objects is returned.
"""
LoadG2fil()
datafile = os.path.abspath(os.path.expanduser(datafile))
try:
iparams = os.path.abspath(os.path.expanduser(iparams))
except:
pass
pwdrreaders = import_generic(datafile, Readers['Pwdr'],fmthint=fmthint,bank=databank)
if not multiple: pwdrreaders = pwdrreaders[0:1]
histlist = []
for r in pwdrreaders:
histname, new_names, pwdrdata = load_pwd_from_reader(r, iparams,
[h.name for h in self.histograms()],bank=instbank)
if histname in self.data:
G2fil.G2Print("Warning - redefining histogram", histname)
elif self.names[-1][0] == 'Phases':
self.names.insert(-1, new_names)
else:
self.names.append(new_names)
self.data[histname] = pwdrdata
self.update_ids()
if phases == 'all':
phases = self.phases()
for phase in phases:
phase = self.phase(phase)
self.link_histogram_phase(histname, phase)
histlist.append(self.histogram(histname))
if multiple:
return histlist
else:
return histlist[0]
[docs]
def clone_powder_histogram(self, histref, newname, Y, Yerr=None):
'''Creates a copy of a powder diffraction histogram with new Y values.
The X values are not changed. The number of Y values must match the
number of X values.
:param histref: The histogram object, the name of the histogram (str), or ranId
or histogram index.
:param str newname: The name to be assigned to the new histogram
:param list Y: A set of intensity values
:param list Yerr: A set of uncertainties for the intensity values (may be None,
sets all weights to unity)
:returns: the new histogram object (type G2PwdrData)
'''
hist = self.histogram(histref)
for i in self.names:
if i[0] == hist.name:
subkeys = i[1:]
break
else:
raise Exception("error in self.names, hist not found")
orighist = hist.name
newhist = 'PWDR '+newname
if len(Y) != len(self[orighist]['data'][1][0]):
raise Exception("clone error: length of Y does not match number of X values ({})"
.format(len(self[orighist]['data'][1][0])))
if Yerr is not None and len(Yerr) != len(self[orighist]['data'][1][0]):
raise Exception("clone error: length of Yerr does not match number of X values ({})"
.format(len(self[orighist]['data'][1][0])))
self[newhist] = copy.deepcopy(self[orighist])
# intensities
yo = self[newhist]['data'][1][1] = ma.MaskedArray(Y,mask=self[orighist]['data'][1][1].mask)
Ymin,Ymax = yo.min(),yo.max()
# set to zero: weights, calc, bkg, obs-calc
for i in [2,3,4,5]:
self[newhist]['data'][1][i] *= 0
# weights
if Yerr is not None:
self[newhist]['data'][1][2] += 1./np.array(Yerr)**2
else:
self[newhist]['data'][1][2] += 1 # set all weights to 1
self[newhist]['data'][0] = {'wtFactor': 1.0, 'Dummy': False,
'ranId': ran.randint(0, sys.maxsize),
'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
'Yminmax': [Ymin, Ymax]}
self[newhist]['Comments'].insert(0,'Cloned from '+orighist)
self[newhist]['Reflection Lists'] = {}
self[newhist]['Index Peak List'] = [[], []]
self[newhist]['Unit Cells List'] = []
self[newhist]['Peak List'] = {'peaks': [], 'sigDict': {}}
self.names.append([newhist]+subkeys)
self.update_ids()
return self.histogram(newhist)
[docs]
def add_simulated_powder_histogram(self, histname, iparams, Tmin, Tmax, Tstep=None,
wavelength=None, scale=None, phases=[], ibank=None,
Npoints=None):
"""Create a simulated powder data histogram for the project.
Requires an instrument parameter file.
Note that in unix fashion, "~" can be used to indicate the
home directory (e.g. ~/G2data/data.prm). The instrument parameter file
will determine if the histogram is x-ray, CW neutron, TOF, etc. as well
as the instrument type.
:param str histname: A name for the histogram to be created.
:param str iparams: The instrument parameters file, a filename.
:param float Tmin: Minimum 2theta or TOF (millisec) for dataset to be simulated
:param float Tmax: Maximum 2theta or TOF (millisec) for dataset to be simulated
:param float Tstep: Step size in 2theta or deltaT/T (TOF) for simulated dataset.
Default is to compute this from Npoints.
:param float wavelength: Wavelength for CW instruments, overriding the value
in the instrument parameters file if specified. For single-wavelength histograms,
this should be a single float value, for K alpha 1,2 histograms, this should
be a list or tuple with two values.
:param float scale: Histogram scale factor which multiplies the pattern. Note that
simulated noise is added to the pattern, so that if the maximum intensity is
small, the noise will mask the computed pattern. The scale needs to be a large
number for neutrons.
The default, None, provides a scale of 1 for x-rays, 10,000 for CW neutrons
and 100,000 for TOF.
:param list phases: Phases to link to the new histogram. Use proj.phases() to link to
all defined phases.
:param int ibank: provides a bank number for the instrument parameter file. The
default is None, corresponding to load the first bank.
:param int Νpoints: the number of data points to be used for computing the
diffraction pattern. Defaults as None, which sets this to 2500. Do not specify
both Npoints and Tstep. Due to roundoff the actual number of points used may differ
by +-1 from Npoints. Must be below 25,000.
:returns: A :class:`G2PwdrData` object representing the histogram
"""
LoadG2fil()
iparams = os.path.abspath(os.path.expanduser(iparams))
if not os.path.exists(iparams):
raise G2ScriptException("File does not exist:"+iparams)
rd = G2obj.ImportPowderData( # Initialize a base class reader
extensionlist=tuple(),
strictExtension=False,
formatName = 'Simulate dataset',
longFormatName = 'Compute a simulated pattern')
rd.powderentry[0] = '' # no filename
rd.powderentry[2] = 1 # only one bank
rd.comments.append('This is a dummy dataset for powder pattern simulation')
rd.idstring = histname
#Iparm1, Iparm2 = load_iprms(iparams, rd)
if Tmax < Tmin:
Tmin,Tmax = Tmax,Tmin
if Tstep is not None and Npoints is not None:
raise G2ScriptException("Error: Tstep and Npoints both specified")
elif Tstep is not None:
Tstep = abs(Tstep)
elif Npoints is None:
Npoints = 2500
Iparm1, Iparm2 = load_iprms(iparams, rd, bank=ibank)
#G2fil.G2Print('Instrument parameters read:',reader.instmsg)
if 'T' in Iparm1['Type'][0]:
# patch -- anticipate TOF values in microsec from buggy version
if Tmax > 200.:
print('Error: Tmax is too large. Note that input for TOF Tmin & Tmax has changed.')
print(' Tmin & Tmax are now in milliseconds not microsec. Step is now deltaT/T.')
raise G2ScriptException("Error: Tmax is too large")
if Npoints:
N = Npoints
Tstep = (np.log(Tmax)-np.log(Tmin))/N
else:
N = (np.log(Tmax)-np.log(Tmin))/Tstep
if N > 25000:
raise G2ScriptException("Error: Tstep is too small. Would need "+str(N)+" points.")
x = np.exp((np.arange(0,N))*Tstep+np.log(Tmin*1000.))
N = len(x)
unit = 'millisec'
else:
if Npoints:
N = Npoints
else:
N = int((Tmax-Tmin)/Tstep)+1
if N > 25000:
raise G2ScriptException("Error: Tstep is too small. Would need "+str(N)+" points.")
x = np.linspace(Tmin,Tmax,N,True)
N = len(x)
unit = 'degrees 2theta'
if N < 3:
raise G2ScriptException("Error: Range is too small or step is too large, <3 points")
G2fil.G2Print('Simulating {} points from {} to {} {}'.format(N,Tmin,Tmax,unit))
rd.powderdata = [
np.array(x), # x-axis values
np.zeros_like(x), # powder pattern intensities
np.ones_like(x), # 1/sig(intensity)^2 values (weights)
np.zeros_like(x), # calc. intensities (zero)
np.zeros_like(x), # calc. background (zero)
np.zeros_like(x), # obs-calc profiles
]
Tmin = rd.powderdata[0][0]
Tmax = rd.powderdata[0][-1]
histname, new_names, pwdrdata = load_pwd_from_reader(rd, iparams,
[h.name for h in self.histograms()],ibank)
if histname in self.data:
G2fil.G2Print("Warning - redefining histogram", histname)
elif self.names[-1][0] == 'Phases':
self.names.insert(-1, new_names)
else:
self.names.append(new_names)
if scale is not None:
pwdrdata['Sample Parameters']['Scale'][0] = scale
elif pwdrdata['Instrument Parameters'][0]['Type'][0].startswith('PNC'):
pwdrdata['Sample Parameters']['Scale'][0] = 10000.
elif pwdrdata['Instrument Parameters'][0]['Type'][0].startswith('PNT'):
pwdrdata['Sample Parameters']['Scale'][0] = 100000.
if wavelength is not None:
if 'Lam1' in pwdrdata['Instrument Parameters'][0]:
# have alpha 1,2 here
try:
pwdrdata['Instrument Parameters'][0]['Lam1'][0] = float(wavelength[0])
pwdrdata['Instrument Parameters'][0]['Lam1'][1] = float(wavelength[0])
pwdrdata['Instrument Parameters'][0]['Lam2'][0] = float(wavelength[1])
pwdrdata['Instrument Parameters'][0]['Lam2'][1] = float(wavelength[1])
except:
raise G2ScriptException("add_simulated_powder_histogram Error: only one wavelength with alpha 1+2 histogram?")
elif 'Lam' in pwdrdata['Instrument Parameters'][0]:
try:
pwdrdata['Instrument Parameters'][0]['Lam'][0] = float(wavelength)
pwdrdata['Instrument Parameters'][0]['Lam'][1] = float(wavelength)
except:
raise G2ScriptException("add_simulated_powder_histogram Error: invalid wavelength?")
else:
raise G2ScriptException("add_simulated_powder_histogram Error: can't set a wavelength for a non-CW dataset")
self.data[histname] = pwdrdata
self.update_ids()
for phase in phases:
phase = self.phase(phase)
self.link_histogram_phase(histname, phase)
self.set_Controls('cycles', 0)
self.data[histname]['data'][0]['Dummy'] = True
return self.histogram(histname)
[docs]
def add_phase(self, phasefile=None, phasename=None, histograms=[],
fmthint=None, mag=False,
spacegroup='P 1',cell=None):
"""Loads a phase into the project, usually from a .cif file
:param str phasefile: The CIF file (or other file type, see fmthint)
that the phase will be read from.
May be left as None (the default) if the phase will be constructed
a step at a time.
:param str phasename: The name of the new phase, or None for the
default. A phasename must be specified when a phasefile is not.
:param list histograms: The names of the histograms to associate with
this phase. Use proj.histograms() to add to all histograms.
:param str fmthint: If specified, only importers where the format name
(reader.formatName, as shown in Import menu) contains the
supplied string will be tried as importers. If not specified, all
importers consistent with the file extension will be tried
(equivalent to "guess format" in menu). Specifying this
is optional but is strongly encouraged.
:param bool mag: Set to True to read a magCIF
:param str spacegroup: The space group name as a string. The
space group must follow the naming rules used in
:func:`GSASIIspc.SpcGroup`. Defaults to 'P 1'. Note that
this is only used when phasefile is None.
:param list cell: a list with six unit cell constants
(a, b, c, alpha, beta and gamma in Angstrom/degrees).
:returns: A :class:`G2Phase` object representing the
new phase.
"""
LoadG2fil()
histograms = [self.histogram(h).name for h in histograms]
if phasefile is None: # create a phase, rather than reading it
if phasename is None:
raise Exception('add_phase: phasefile and phasename cannot both be None')
phaseNameList = [p.name for p in self.phases()]
phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
self.data['Phases'] = self.data.get('Phases', {'data': None})
err,SGData=G2spc.SpcGroup(spacegroup)
if err != 0:
print('Space group error:', G2spc.SGErrors(err))
raise Exception('Space group error')
self.data['Phases'][phasename] = G2obj.SetNewPhase(
Name=phasename,SGData=SGData)
self.data['Phases'][phasename]['General']['Name'] = phasename
for hist in histograms:
self.link_histogram_phase(hist, phasename)
for obj in self.names:
if obj[0] == 'Phases':
phasenames = obj
break
else:
phasenames = ['Phases']
self.names.append(phasenames)
phasenames.append(phasename)
data = self.data['Phases'][phasename]
if cell:
if len(cell) != 6:
raise Exception('Error: Unit cell must have 6 entries')
data['General']['Cell'][1:7] = cell
data['General']['Cell'][7] = G2lat.calc_V(G2lat.cell2A(cell))
SetupGeneral(data, None)
self.index_ids()
self.update_ids()
return self.phase(phasename)
phasefile = os.path.abspath(os.path.expanduser(phasefile))
if not os.path.exists(phasefile):
raise G2ImportException(f'File {phasefile} does not exist')
# TODO: handle multiple phases in a file
phasereaders = import_generic(phasefile, Readers['Phase'], fmthint=fmthint)
phasereader = phasereaders[0]
if phasereader.MPhase and mag:
phObj = phasereader.MPhase
else:
phObj = phasereader.Phase
phasename = phasename or phObj['General']['Name']
phaseNameList = [p.name for p in self.phases()]
phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
phObj['General']['Name'] = phasename
if 'Phases' not in self.data:
self.data['Phases'] = { 'data': None }
assert phasename not in self.data['Phases'], "phase names should be unique"
self.data['Phases'][phasename] = phObj
# process constraints, currently generated only from ISODISTORT CIFs
if phasereader.Constraints:
Constraints = self.data['Constraints']
for i in phasereader.Constraints:
if isinstance(i, dict):
if '_Explain' not in Constraints:
Constraints['_Explain'] = {}
Constraints['_Explain'].update(i)
else:
Constraints['Phase'].append(i)
data = self.data['Phases'][phasename]
# generalData = data['General']
# SGData = generalData['SGData']
# NShkl = len(G2spc.MustrainNames(SGData))
# NDij = len(G2spc.HStrainNames(SGData))
# Super = generalData.get('Super', 0)
# if Super:
# SuperVec = np.array(generalData['SuperVec'][0])
# else:
# SuperVec = []
# UseList = data['Histograms']
for hist in histograms:
self.link_histogram_phase(hist, phasename)
for obj in self.names:
if obj[0] == 'Phases':
phasenames = obj
break
else:
phasenames = ['Phases']
self.names.append(phasenames)
phasenames.append(phasename)
# TODO should it be self.filename, not phasefile?
SetupGeneral(data, os.path.dirname(phasefile))
self.index_ids()
self.update_ids()
return self.phase(phasename)
[docs]
def add_single_histogram(self, datafile, phase=None, fmthint=None):
"""Loads a powder data histogram or multiple powder histograms
into the project.
:param str datafile: A filename with the single crystal data file
to read. Note that in unix fashion, "~" can be used to indicate the
home directory (e.g. ~/G2data/data.hkl).
:param phases: A phase to link to the new histogram. A
phase can be referenced by object, name, rId or number.
If not specified, no phase will be linked.
:param str fmthint: If specified, only importers where the format name
(reader.formatName, as shown in Import menu) contains the
supplied string will be tried as importers. If not specified,
an error will be generated, as the file format will not distinguish
well between different data types.
:returns: A :class:`G2Single` object representing
the histogram
"""
LoadG2fil()
datafile = os.path.abspath(os.path.expanduser(datafile))
if fmthint is None:
choices = ', '.join(['"'+i.formatName+'"' for i in Readers['HKLF']])
raise ValueError("add_single_histogram: A value is required for the fmthint parameter"+
f'\n\nChoices are: {choices}')
readers = import_generic(datafile, Readers['HKLF'],fmthint=fmthint)
#histlist = []
for r in readers: # only expect 1
histname = 'HKLF ' + G2obj.StripUnicode(
os.path.split(r.readfilename)[1],'_')
#HistName = G2obj.MakeUniqueLabel(HistName, existingnames)
newnames = [histname,'Instrument Parameters','Reflection List']
if histname in self.data:
G2fil.G2Print("Warning - redefining histogram", histname)
elif self.names[-1][0] == 'Phases':
self.names.insert(-1, newnames) # add histogram into tree before phases
else:
self.names.append(newnames)
data = [{
'wtFactor': 1.0,
'Dummy': False,
'ranId': ran.randint(0, sys.maxsize),
'histTitle': ''},r.RefDict]
self.data[histname] = {
'data':data,
'Instrument Parameters':r.Parameters,
'Reflection List': {}
}
self.update_ids()
if phase is None: return
self.link_histogram_phase(histname, phase)
[docs]
def link_histogram_phase(self, histogram, phase):
"""Associates a given histogram and phase.
.. seealso::
:meth:`~G2Project.histogram`
:meth:`~G2Project.phase`"""
import GSASIImath as G2mth
hist = self.histogram(histogram)
phase = self.phase(phase)
generalData = phase['General']
if hist.name.startswith('HKLF '):
G2mth.UpdateHKLFvals(hist.name, phase, hist.data['data'][1])
elif hist.name.startswith('PWDR '):
hist['Reflection Lists'][generalData['Name']] = {}
UseList = phase['Histograms']
SGData = generalData['SGData']
NShkl = len(G2spc.MustrainNames(SGData))
NDij = len(G2spc.HStrainNames(SGData))
UseList[hist.name] = G2mth.SetDefaultDData(
'PWDR', hist.name, NShkl=NShkl, NDij=NDij)
UseList[hist.name]['hId'] = hist.id
for key, val in [('Use', True), ('LeBail', False),
('Babinet', {'BabA': [0.0, False],
'BabU': [0.0, False]})]:
if key not in UseList[hist.name]:
UseList[hist.name][key] = val
else:
raise RuntimeError("Unexpected histogram" + hist.name)
[docs]
def reload(self):
"""Reload self from self.filename"""
data, names = LoadDictFromProjFile(self.filename)
self.names = names
# Need to deep copy the new data file data into the current tree,
# so that any existing G2Phase, or G2PwdrData objects will still be
# valid
_deep_copy_into(from_=data, into=self.data)
[docs]
def refine(self, newfile=None, printFile=None, makeBack=False):
'''Invoke a refinement for the project. The project is written to
the currently selected gpx file and then either a single or sequential refinement
is performed depending on the setting of 'Seq Data' in Controls
(set in :meth:`get_Controls`).
'''
seqSetting = self.data['Controls']['data'].get('Seq Data',[])
if not seqSetting:
self.index_ids() # index_ids will automatically save the project
# TODO: migrate to RefineCore G2strMain does not properly use printFile
# G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
# calcControls,pawleyLookup,ifPrint,printFile,dlg)
# check that constraints are OK
errmsg, warnmsg = G2stIO.ReadCheckConstraints(self.filename)
if errmsg:
G2fil.G2Print('Constraint error',errmsg)
raise Exception('Constraint error')
if warnmsg:
G2fil.G2Print('\nNote these constraint warning(s):\n'+warnmsg)
G2fil.G2Print('Generated constraints\n'+G2mv.VarRemapShow([],True))
G2strMain.Refine(self.filename, makeBack=makeBack)
else:
self._seqrefine()
self.reload() # get file from GPX
def _seqrefine(self):
'''Perform a sequential refinement.
'''
self.data['Controls']['data']['ShowCell'] = True
# add to tree item to project, if not present
if 'Sequential results' not in self.data:
self.data['Sequential results'] = {'data':{}}
self.names.append(['Sequential results'])
self.index_ids() # index_ids will automatically save the project
#GSASIIpath.IPyBreak_base()
# check that constraints are OK
errmsg, warnmsg = G2stIO.ReadCheckConstraints(self.filename)
if errmsg:
G2fil.G2Print('Constraint error',errmsg)
raise Exception('Constraint error')
if warnmsg:
G2fil.G2Print('\nNote these constraint warning(s):\n'+warnmsg)
G2fil.G2Print('Generated constraints\n'+G2mv.VarRemapShow([],True))
OK,Msg = G2strMain.SeqRefine(self.filename,None)
[docs]
def histogram(self, histname):
"""Returns the histogram object associated with histname, or None
if it does not exist.
:param histname: The name of the histogram (str), or ranId or
(for powder) the histogram index.
:returns: A :class:`G2PwdrData` object, or :class:`G2Single` object, or
None if the histogram does not exist
.. seealso::
:meth:`~G2Project.histograms`
:meth:`~G2Project.phase`
:meth:`~G2Project.phases`
"""
if isinstance(histname, G2PwdrData) or isinstance(histname, G2Single):
if histname.proj == self:
return histname
else:
raise Exception(f'Histogram object {histname} (type {type(histname)}) is project {histname.proj}')
if histname in self.data:
if histname.startswith('HKLF '):
return G2Single(self.data[histname], self, histname)
elif histname.startswith('PWDR '):
return G2PwdrData(self.data[histname], self, histname)
try:
# see if histname is an id or ranId
histname = int(histname)
except ValueError:
return
for histogram in self.histograms():
if histogram.name.startswith('HKLF '):
if histogram.data['data'][0]['ranId'] == histname:
return histogram
elif histogram.name.startswith('PWDR '):
if histogram.id == histname or histogram.ranId == histname:
return histogram
[docs]
def histType(self, histname):
"""Returns the type for histogram object associated with histname, or
None if it does not exist.
:param histname: The name of the histogram (str), or ranId or
(for powder) the histogram index.
:returns: 'PWDR' for a Powder histogram,
'HKLF' for a single crystal histogram, or
None if the histogram does not exist
.. seealso::
:meth:`~G2Project.histogram`
"""
if isinstance(histname, G2PwdrData):
if histname.proj == self:
return 'PWDR'
else:
raise Exception(f'Histogram object {histname} (type {type(histname)}) is project {histname.proj}')
if isinstance(histname, G2Single):
if histname.proj == self:
return 'HKLF'
else:
raise Exception(f'Histogram object {histname} (type {type(histname)}) is project {histname.proj}')
if histname in self.data:
if histname.startswith('HKLF '):
return 'HKLF'
elif histname.startswith('PWDR '):
return 'PWDR'
try:
# see if histname is an id or ranId
histname = int(histname)
except ValueError:
return
for histogram in self.histograms():
if histogram.name.startswith('HKLF '):
if histogram.data['data'][0]['ranId'] == histname:
return 'HKLF'
elif histogram.name.startswith('PWDR '):
if histogram.id == histname or histogram.ranId == histname:
return 'PWDR'
[docs]
def histograms(self, typ=None):
"""Return a list of all histograms, as :class:`G2PwdrData` objects
For now this only finds Powder/Single Xtal histograms, since that is all that is
currently implemented in this module.
:param ste typ: The prefix (type) the histogram such as 'PWDR ' for
powder or 'HKLF ' for single crystal. If None
(the default) all known histograms types are found.
:returns: a list of objects
.. seealso::
:meth:`~G2Project.histogram`
:meth:`~G2Project.phase`
:meth:`~G2Project.phases`
"""
output = []
# loop through each tree entry. If it is more than one level (more than one item in the
# list of names). then it must be a histogram, unless labeled Phases or Restraints
if typ is None:
for obj in self.names:
if obj[0].startswith('PWDR ') or obj[0].startswith('HKLF '):
output.append(self.histogram(obj[0]))
else:
for obj in self.names:
if len(obj) > 1 and obj[0].startswith(typ):
output.append(self.histogram(obj[0]))
return output
[docs]
def phase(self, phasename):
"""
Gives an object representing the specified phase in this project.
:param str phasename: A reference to the desired phase. Either the phase
name (str), the phase's ranId, the phase's index (both int) or
a phase object (:class:`G2Phase`)
:returns: A :class:`G2Phase` object
:raises: KeyError
.. seealso::
:meth:`~G2Project.histograms`
:meth:`~G2Project.phase`
:meth:`~G2Project.phases`
"""
if isinstance(phasename, G2Phase):
if phasename.proj == self:
return phasename
phases = self.data['Phases']
if phasename in phases:
return G2Phase(phases[phasename], phasename, self)
try:
# phasename should be phase index or ranId
phasename = int(phasename)
except ValueError:
return
for phase in self.phases():
if phase.id == phasename or phase.ranId == phasename:
return phase
[docs]
def phases(self):
"""
Returns a list of all the phases in the project.
:returns: A list of :class:`G2Phase` objects
.. seealso::
:meth:`~G2Project.histogram`
:meth:`~G2Project.histograms`
:meth:`~G2Project.phase`
"""
for obj in self.names:
if obj[0] == 'Phases':
return [self.phase(p) for p in obj[1:]]
return []
def _images(self):
"""Returns a list of all the images in the project.
"""
return [i[0] for i in self.names if i[0].startswith('IMG ')]
[docs]
def image(self, imageRef):
"""
Gives an object representing the specified image in this project.
:param str imageRef: A reference to the desired image. Either the Image
tree name (str), the image's index (int) or
a image object (:class:`G2Image`)
:returns: A :class:`G2Image` object
:raises: KeyError
.. seealso::
:meth:`~G2Project.images`
"""
if isinstance(imageRef, G2Image):
if imageRef.proj == self:
return imageRef
else:
raise Exception("Image {} not in current selected project".format(imageRef.name))
if imageRef in self._images():
return G2Image(self.data[imageRef], imageRef, self)
try:
# imageRef should be an index
num = int(imageRef)
imageRef = self._images()[num]
return G2Image(self.data[imageRef], imageRef, self)
except ValueError:
raise Exception("imageRef {} not an object, name or image index in current selected project"
.format(imageRef))
except IndexError:
raise Exception("imageRef {} out of range (max={}) in current selected project"
.format(imageRef,len(self._images())-1))
[docs]
def images(self):
"""
Returns a list of all the images in the project.
:returns: A list of :class:`G2Image` objects
"""
return [G2Image(self.data[i],i,self) for i in self._images()]
def _pdfs(self):
"""Returns a list of all the PDF entries in the project.
"""
return [i[0] for i in self.names if i[0].startswith('PDF ')]
[docs]
def pdf(self, pdfRef):
"""
Gives an object representing the specified PDF entry in this project.
:param pdfRef: A reference to the desired image. Either the PDF
tree name (str), the pdf's index (int) or
a PDF object (:class:`G2PDF`)
:returns: A :class:`G2PDF` object
:raises: KeyError
.. seealso::
:meth:`~G2Project.pdfs`
:class:`~G2PDF`
"""
if isinstance(pdfRef, G2PDF):
if pdfRef.proj == self:
return pdfRef
else:
raise Exception(f"PDF {pdfRef.name} not in current selected project")
if pdfRef in self._pdfs():
return G2PDF(self.data[pdfRef], pdfRef, self)
try:
# pdfRef should be an index
num = int(pdfRef)
pdfRef = self._pdfs()[num]
return G2PDF(self.data[pdfRef], pdfRef, self)
except ValueError:
raise Exception(f"pdfRef {pdfRef} not an object, name or PDF index in current selected project")
except IndexError:
raise Exception(f"pdfRef {pdfRef} out of range (max={len(G2SmallAngle)-1}) in current selected project")
[docs]
def pdfs(self):
"""
Returns a list of all the PDFs in the project.
:returns: A list of :class:`G2PDF` objects
"""
return [G2PDF(self.data[i],i,self) for i in self._pdfs()]
[docs]
def copy_PDF(self, PDFobj, histogram):
'''Creates a PDF entry that can be used to compute a PDF
as a copy of settings in an existing PDF (:class:`G2PDF`)
object.
This places an entry in the project but :meth:`G2PDF.calculate`
must be used to actually perform the PDF computation.
:param PDFobj: A :class:`G2PDF` object which may be
in a separate project or the dict associated with the
PDF object (G2PDF.data).
:param histogram: A reference to a histogram,
which can be reference by object, name, or number.
:returns: A :class:`G2PDF` object for the PDF entry
'''
LoadG2fil()
PDFname = 'PDF ' + self.histogram(histogram).name[4:]
PDFdict = {'data':None}
for i in 'PDF Controls', 'PDF Peaks':
PDFdict[i] = copy.deepcopy(PDFobj[i])
self.names.append([PDFname]+['PDF Controls', 'PDF Peaks'])
self.data[PDFname] = PDFdict
for i in 'I(Q)','S(Q)','F(Q)','G(R)','g(r)':
self.data[PDFname]['PDF Controls'][i] = []
G2fil.G2Print('Adding "{}" to project'.format(PDFname))
return G2PDF(self.data[PDFname], PDFname, self)
[docs]
def add_PDF(self, prmfile, histogram):
'''Creates a PDF entry that can be used to compute a PDF.
Note that this command places an entry in the project,
but :meth:`G2PDF.calculate` must be used to actually perform
the computation.
:param str datafile: The powder data file to read, a filename.
:param histogram: A reference to a histogram,
which can be reference by object, name, or number.
:returns: A :class:`G2PDF` object for the PDF entry
'''
LoadG2fil()
PDFname = 'PDF ' + self.histogram(histogram).name[4:]
peaks = {'Limits':[1.,5.],'Background':[2,[0.,-0.2*np.pi],False],'Peaks':[]}
Controls = {
'Sample':{'Name':self.histogram(histogram).name,'Mult':1.0},
'Sample Bkg.':{'Name':'','Mult':-1.0,'Refine':False},
'Container':{'Name':'','Mult':-1.0,'Refine':False},
'Container Bkg.':{'Name':'','Mult':-1.0},
'ElList':{},
'Geometry':'Cylinder','Diam':1.0,'Pack':0.50,'Form Vol':0.0,'Flat Bkg':0,
'DetType':'Area detector','ObliqCoeff':0.2,'Ruland':0.025,'QScaleLim':[20,25],
'Lorch':False,'BackRatio':0.0,'Rmax':100.,'noRing':False,'IofQmin':1.0,'Rmin':1.0,
'I(Q)':[],'S(Q)':[],'F(Q)':[],'G(R)':[],'g(r)':[]}
fo = open(prmfile,'r')
S = fo.readline()
while S:
if '#' in S:
S = fo.readline()
continue
key,val = S.split(':',1)
try:
Controls[key] = eval(val)
except:
Controls[key] = val.strip()
S = fo.readline()
fo.close()
Controls['Sample']['Name'] = self.histogram(histogram).name
for i in 'Sample Bkg.','Container','Container Bkg.':
Controls[i]['Name'] = ''
PDFdict = {'data':None,'PDF Controls':Controls, 'PDF Peaks':peaks}
self.names.append([PDFname]+['PDF Controls', 'PDF Peaks'])
self.data[PDFname] = PDFdict
G2fil.G2Print('Adding "{}" to project'.format(PDFname))
return G2PDF(self.data[PDFname], PDFname, self)
[docs]
def seqref(self):
"""
Returns a sequential refinement results object, if present
:returns: A :class:`G2SeqRefRes` object or None if not present
"""
if 'Sequential results' not in self.data: return
return G2SeqRefRes(self.data['Sequential results']['data'], self)
[docs]
def update_ids(self):
"""Makes sure all phases and histograms have proper hId and pId"""
# Translated from GetUsedHistogramsAndPhasesfromTree,
# GSASIIdataGUI.py:4107
for i, h in enumerate(self.histograms()):
h.id = i
for i, p in enumerate(self.phases()):
p.id = i
[docs]
def do_refinements(self, refinements=[{}], histogram='all', phase='all',
outputnames=None, makeBack=False):
"""Conducts one or a series of refinements according to the
input provided in parameter refinements. This is a wrapper
around :meth:`iter_refinements`
:param list refinements: A list of dictionaries specifiying changes to be made to
parameters before refinements are conducted.
See the :ref:`Refinement_recipe` section for how this is defined.
If not specified, the default value is ``[{}]``, which performs a single
refinement step is performed with the current refinement settings.
:param str histogram: Name of histogram for refinements to be applied
to, or 'all'; note that this can be overridden for each refinement
step via a "histograms" entry in the dict.
:param str phase: Name of phase for refinements to be applied to, or
'all'; note that this can be overridden for each refinement
step via a "phases" entry in the dict.
:param list outputnames: Provides a list of project (.gpx) file names
to use for each refinement step (specifying None skips the save step).
See :meth:`save`.
Note that this can be overridden using an "output" entry in the dict.
:param bool makeBack: determines if a backup ).bckX.gpx) file is made
before a refinement is performed. The default is False.
To perform a single refinement without changing any parameters, use this
call:
.. code-block:: python
my_project.do_refinements([])
"""
for proj in self.iter_refinements(refinements, histogram, phase,
outputnames, makeBack):
pass
return self
[docs]
def iter_refinements(self, refinements, histogram='all', phase='all',
outputnames=None, makeBack=False):
"""Conducts a series of refinements, iteratively. Stops after every
refinement and yields this project, to allow error checking or
logging of intermediate results. Parameter use is the same as for
:meth:`do_refinements` (which calls this method).
>>> def checked_refinements(proj):
... for p in proj.iter_refinements(refs):
... # Track intermediate results
... log(p.histogram('0').residuals)
... log(p.phase('0').get_cell())
... # Check if parameter diverged, nonsense answer, or whatever
... if is_something_wrong(p):
... raise Exception("I need a human!")
"""
if outputnames:
if len(refinements) != len(outputnames):
raise ValueError("Should have same number of outputs to"
"refinements")
else:
outputnames = [None for r in refinements]
for output, refinedict in zip(outputnames, refinements):
if 'histograms' in refinedict:
hist = refinedict['histograms']
else:
hist = histogram
if 'phases' in refinedict:
ph = refinedict['phases']
else:
ph = phase
if 'output' in refinedict:
output = refinedict['output']
self.set_refinement(refinedict, hist, ph)
# Handle 'once' args - refinements that are disabled after this
# refinement
if 'once' in refinedict:
temp = {'set': refinedict['once']}
self.set_refinement(temp, hist, ph)
if output:
self.save(output)
if 'skip' not in refinedict:
self.refine(makeBack=makeBack)
yield self
# Handle 'once' args - refinements that are disabled after this
# refinement
if 'once' in refinedict:
temp = {'clear': refinedict['once']}
self.set_refinement(temp, hist, ph)
if 'call' in refinedict:
fxn = refinedict['call']
if callable(fxn):
fxn(*refinedict.get('callargs',[self]))
elif callable(eval(fxn)):
eval(fxn)(*refinedict.get('callargs',[self]))
else:
raise G2ScriptException("Error: call value {} is not callable".format(fxn))
[docs]
def set_refinement(self, refinement, histogram='all', phase='all'):
"""Set refinment flags at the project level to specified histogram(s)
or phase(s).
:param dict refinement: The refinements to be conducted
:param histogram: Specifies either 'all' (default), a single histogram or
a list of histograms. Histograms may be specified as histogram objects
(see :class:`G2PwdrData`), the histogram name (str) or the index number (int)
of the histogram in the project, numbered starting from 0.
Omitting the parameter or the string 'all' indicates that parameters in
all histograms should be set.
:param phase: Specifies either 'all' (default), a single phase or
a list of phases. Phases may be specified as phase objects
(see :class:`G2Phase`), the phase name (str) or the index number (int)
of the phase in the project, numbered starting from 0.
Omitting the parameter or the string 'all' indicates that parameters in
all phases should be set.
Note that refinement parameters are categorized as one of three types:
1. Histogram parameters
2. Phase parameters
3. Histogram-and-Phase (HAP) parameters
.. seealso::
:meth:`G2PwdrData.set_refinements`
:meth:`G2PwdrData.clear_refinements`
:meth:`G2Phase.set_refinements`
:meth:`G2Phase.clear_refinements`
:meth:`G2Phase.set_HAP_refinements`
:meth:`G2Phase.clear_HAP_refinements`
:meth:`G2Single.set_refinements`
"""
if histogram == 'all':
hists = self.histograms()
elif isinstance(histogram, list) or isinstance(histogram, tuple):
hists = []
for h in histogram:
if isinstance(h, str) or isinstance(h, int):
hists.append(self.histogram(h))
else:
hists.append(h)
elif isinstance(histogram, str) or isinstance(histogram, int):
hists = [self.histogram(histogram)]
else:
hists = [histogram]
if phase == 'all':
phases = self.phases()
elif isinstance(phase, list) or isinstance(phase, tuple):
phases = []
for ph in phase:
if isinstance(ph, str) or isinstance(ph, int):
phases.append(self.phase(ph))
else:
phases.append(ph)
elif isinstance(phase, str) or isinstance(phase, int):
phases = [self.phase(phase)]
else:
phases = [phase]
pwdr_set = {}
phase_set = {}
hap_set = {}
xtal_set = {}
# divide the refinement set commands across the various objects
for key, val in refinement.get('set', {}).items():
if G2PwdrData.is_valid_refinement_key(key):
pwdr_set[key] = val
elif G2Phase.is_valid_refinement_key(key):
phase_set[key] = val
elif G2Phase.is_valid_HAP_refinement_key(key):
if key == 'PhaseFraction': key = 'Scale'
hap_set[key] = val
elif G2Single.is_valid_refinement_key(key):
xtal_set[key] = val
else:
raise ValueError("Unknown refinement key", key)
for hist in hists:
if isinstance(hist, G2PwdrData):
hist.set_refinements(pwdr_set)
elif isinstance(hist, G2Single):
hist.set_refinements(xtal_set)
else:
raise KeyError(f"set_refinement error: unknown hist type {hist}")
for phase in phases:
phase.set_refinements(phase_set)
for phase in phases:
phase.set_HAP_refinements(hap_set, hists)
pwdr_clear = {}
phase_clear = {}
hap_clear = {}
xtal_clear = {}
for key, val in refinement.get('clear', {}).items():
# Clear refinement options
if G2PwdrData.is_valid_refinement_key(key):
pwdr_clear[key] = val
elif G2Phase.is_valid_refinement_key(key):
phase_clear[key] = val
elif G2Phase.is_valid_HAP_refinement_key(key):
if key == 'PhaseFraction': key = 'Scale'
hap_clear[key] = val # was _set, seems wrong
elif G2Single.is_valid_refinement_key(key):
xtal_clear[key] = val
else:
raise ValueError("Unknown refinement key", key)
for hist in hists:
if isinstance(hist, G2PwdrData):
hist.clear_refinements(pwdr_clear)
else:
hist.clear_refinements(xtal_clear)
for phase in phases:
phase.clear_refinements(phase_clear)
for phase in phases:
phase.clear_HAP_refinements(hap_clear, hists)
def index_ids(self):
self.save()
return G2stIO.GetUsedHistogramsAndPhases(self.filename)
def _sasd(self):
"""Returns a list of all the SASD entries in the project.
"""
return [i[0] for i in self.names if i[0].startswith('SASD ')]
[docs]
def SAS(self, sasRef):
"""
Gives an object representing the specified SAS entry in this project.
:param sasRef: A reference to the desired SASD entry. Either the SASD
tree name (str), the SASD's index (int) or
a SASD object (:class:`G2SmallAngle`)
:returns: A :class:`G2SmallAngle` object
:raises: KeyError
.. seealso::
:meth:`~G2Project.SASs`
:class:`~G2PDF`
"""
if isinstance(sasRef, G2SmallAngle):
if sasRef.proj == self:
return sasRef
else:
raise Exception(f"SASD {sasRef.name} not in current selected project")
if sasRef in self._sasd():
return G2SmallAngle(self.data[sasRef], sasRef, self)
try:
# sasRef should be an index
num = int(sasRef)
sasRef = self._sasd()[num]
return G2PDF(self.data[sasRef], sasRef, self)
except ValueError:
raise Exception(f"sasRef {sasRef} not an object, name or SAS index in current selected project")
except IndexError:
raise Exception(f"sasRef {sasRef} out of range (max={len(self._sasd())-1}) in current selected project")
[docs]
def SASs(self):
"""
Returns a list of all the Small Angle histograms in the project.
:returns: A list of :class:`G2SmallAngle` objects
"""
return [G2SmallAngle(self.data[i],i,self) for i in self._sasd()]
[docs]
def add_SmallAngle(self, datafile):
'''Placeholder for an eventual routine that will read a small angle dataset
from a file.
:param str datafile: The SASD data file to read, a filename.
:returns: A :class:`G2SmallAngle` object for the SASD entry
'''
raise G2ScriptException("Error: add_SmallAngle not yet implemented.")
#return G2SmallAngle(self.data[PDFname], PDFname, self)
[docs]
def get_Constraints(self,ctype):
'''Returns a list of constraints of the type selected.
:param str ctype: one of the following keywords: 'Hist', 'HAP', 'Phase', 'Global'
:returns: a list of constraints, see the
:ref:`constraint definition descriptions <Constraint_definitions_table>`. Note that
if this list is changed (for example by deleting elements or by changing them)
the constraints in the project are changed.
'''
if ctype in ('Hist', 'HAP', 'Phase', 'Global'):
return self.data['Constraints']['data'][ctype]
else:
raise Exception(("get_Constraints error: value of ctype ({})"
+" must be 'Hist', 'HAP', 'Phase', or 'Global'.")
.format(ctype))
[docs]
def add_HoldConstr(self,varlist,reloadIdx=True,override=False):
'''Set a hold constraint on a list of variables.
Note that this will cause the project to be saved if not
already done so. It will always save the .gpx file before
creating constraint(s) if reloadIdx is True.
:param list varlist: A list of variables to hold.
Each value in the list may be one of the following three items:
(A) a :class:`GSASIIobj.G2VarObj` object,
(B) a variable name (str), or
(C) a list/tuple of arguments for :meth:`make_var_obj`.
:param bool reloadIdx: If True (default) the .gpx file will be
saved and indexed prior to use. This is essential if atoms, phases
or histograms have been added to the project.
:param bool override: This routine looks up variables using
:func:`GSASIIobj.getDescr` (which is not comprehensive). If
not found, the routine will throw an exception, unless
override=True is specified.
Example::
gpx.add_HoldConstr(('0::A4','0:1:D12',':0:Lam'))
'''
if reloadIdx:
self.index_ids()
elif G2obj.TestIndexAll():
self.index_ids()
vardefwarn = False
for var in varlist:
# make var object
if isinstance(var, str):
var = self.make_var_obj(var,reloadIdx=False)
elif not isinstance(var, G2obj.G2VarObj):
var = self.make_var_obj(*var,reloadIdx=False)
if G2obj.getDescr(var.name) is None:
vardefwarn = True
print(f'Constraint var warning: No definition for variable {var.name}. Name correct?')
# make constraint
self.add_constraint_raw(_constr_type(var), [[1.0, var], None, None, 'h'])
if vardefwarn and not override:
raise Exception('Constraint var error: undefined variables.\n\tIf correct, use override=True in call.\n\tAlso, please let GSAS-II developers know so definition can be added.')
[docs]
def add_EquivConstr(self,varlist,multlist=[],reloadIdx=True,override=False):
'''Set a equivalence on a list of variables.
Note that this will cause the project to be saved if not
already done so. It will always save the .gpx file before
creating a constraint if reloadIdx is True.
:param list varlist: A list of variables to make equivalent to the
first item in the list.
Each value in the list may be one of the following three items:
(A) a :class:`GSASIIobj.G2VarObj` object,
(B) a variable name (str), or
(C) a list/tuple of arguments for :meth:`make_var_obj`.
:param list multlist: a list of multipliers for each variable in
varlist. If there are fewer values than supplied for varlist
then missing values will be set to 1. The default is [] which
means that all multipliers are 1.
:param bool reloadIdx: If True (default) the .gpx file will be
saved and indexed prior to use. This is essential if atoms, phases
or histograms have been added to the project.
:param bool override: This routine looks up variables using
:func:`GSASIIobj.getDescr` (which is not comprehensive). If
not found, the routine will throw an exception, unless
override=True is specified.
Examples::
gpx.add_EquivConstr(('0::AUiso:0','0::AUiso:1','0::AUiso:2'))
gpx.add_EquivConstr(('0::dAx:0','0::dAx:1'),[1,-1])
'''
if reloadIdx:
self.index_ids()
elif G2obj.TestIndexAll():
self.index_ids()
if len(varlist) < 2:
raise Exception('add_EquivConstr Error: varlist must have at least 2 variables')
constr = []
typ_prev = None
vardefwarn = False
for i,var in enumerate(varlist):
m = 1.
try:
m = float(multlist[i])
except IndexError:
pass
# make var object
if isinstance(var, str):
var = self.make_var_obj(var,reloadIdx=False)
elif not isinstance(var, G2obj.G2VarObj):
var = self.make_var_obj(*var,reloadIdx=False)
if G2obj.getDescr(var.name) is None:
vardefwarn = True
print(f'Constraint var warning: No definition for variable {var.name}. Name correct?')
# make constraint
constr.append([m,var])
typ = _constr_type(var)
if typ_prev is None:
typ_prev = typ
var_prev = var
if typ_prev != typ:
msg = 'Type ({}) for var {} is different from {} ({})'.format(typ,var,var_prev,typ_prev)
raise Exception('add_EquivConstr Error: '+msg)
typ_prev = typ
var_prev = var
if vardefwarn and not override:
raise Exception('Constraint var error: undefined variables.\n\tIf correct, use override=True in call.\n\tAlso, please let GSAS-II developers know so definition can be added.')
constr += [None, None, 'e']
self.add_constraint_raw(typ, constr)
[docs]
def add_EqnConstr(self,total,varlist,multlist=[],reloadIdx=True,override=False):
'''Set a constraint equation on a list of variables.
Note that this will cause the project to be saved if not
already done so. It will always save the .gpx file before
creating a constraint if reloadIdx is True.
:param float total: A value that the constraint must equal
:param list varlist: A list of variables to use in the equation.
Each value in the list may be one of the following three items:
(A) a :class:`GSASIIobj.G2VarObj` object,
(B) a variable name (str), or
(C) a list/tuple of arguments for :meth:`make_var_obj`.
:param list multlist: a list of multipliers for each variable in
varlist. If there are fewer values than supplied for varlist
then missing values will be set to 1. The default is [] which
means that all multipliers are 1.
:param bool reloadIdx: If True (default) the .gpx file will be
saved and indexed prior to use. This is essential if atoms, phases
or histograms have been added to the project.
:param bool override: This routine looks up variables using
:func:`GSASIIobj.getDescr` (which is not comprehensive). If
not found, the routine will throw an exception, unless
override=True is specified.
Example::
gpx.add_EqnConstr(1.0,('0::Ax:0','0::Ax:1'),[1,1])
'''
if reloadIdx:
self.index_ids()
elif G2obj.TestIndexAll():
self.index_ids()
if len(varlist) < 2:
raise Exception('Constraint var error: varlist must have at least 2 variables')
try:
float(total)
except:
raise Exception('Constraint var error: total be a valid float')
constr = []
typ_prev = None
vardefwarn = False
for i,var in enumerate(varlist):
m = 1.
try:
m = float(multlist[i])
except IndexError:
pass
# make var object
if isinstance(var, str):
var = self.make_var_obj(var,reloadIdx=False)
elif not isinstance(var, G2obj.G2VarObj):
var = self.make_var_obj(*var,reloadIdx=False)
if G2obj.getDescr(var.name) is None:
vardefwarn = True
print(f'Constraint var warning: No definition for variable {var.name}. Name correct?')
# make constraint
constr.append([m,var])
typ = _constr_type(var)
if typ_prev is None:
typ_prev = typ
var_prev = var
if typ_prev != typ:
msg = 'Type ({}) for var {} is different from {} ({})'.format(typ,var,var_prev,typ_prev)
raise Exception('add_EquivConstr Error: '+msg)
typ_prev = typ
var_prev = var
if vardefwarn and not override:
raise Exception('Constraint var error: undefined variables.\n\tIf correct, use override=True in call.\n\tAlso, please let GSAS-II developers know so definition can be added.')
constr += [float(total), None, 'c']
self.add_constraint_raw(typ, constr)
[docs]
def add_NewVarConstr(self,varlist,multlist=[],name=None,vary=False,
reloadIdx=True,override=False):
'''Set a new-variable constraint from a list of variables to
create a new parameter from two or more predefined parameters.
Note that this will cause the project to be saved, if not
already done so. It will always save the .gpx file before
creating a constraint if reloadIdx is True.
:param list varlist: A list of variables to use in the expression.
Each value in the list may be one of the following three items:
(A) a :class:`GSASIIobj.G2VarObj` object,
(B) a variable name (str), or
(C) a list/tuple of arguments for :meth:`make_var_obj`.
:param list multlist: a list of multipliers for each variable in
varlist. If there are fewer values than supplied for varlist
then missing values will be set to 1. The default is [] which
means that all multipliers are 1.
:param str name: An optional string to be supplied as a name for this
new parameter.
:param bool vary: Determines if the new variable should be flagged to
be refined.
:param bool reloadIdx: If True (default) the .gpx file will be
saved and indexed prior to use. This is essential if atoms, phases
or histograms have been added to the project.
:param bool override: This routine looks up variables using
:func:`GSASIIobj.getDescr` (which is not comprehensive). If
not found, the routine will throw an exception, unless
override=True is specified.
Examples::
gpx.add_NewVarConstr(('0::AFrac:0','0::AFrac:1'),[0.5,0.5],'avg',True)
gpx.add_NewVarConstr(('0::AFrac:0','0::AFrac:1'),[1,-1],'diff',False,False)
The example above is a way to treat two variables that are closely correlated.
The first variable, labeled as avg, allows the two variables to refine in tandem
while the second variable (diff) tracks their difference. In the initial stages of
refinement only avg would be refined, but in the final stages, it might be possible
to refine diff. The second False value in the second example prevents the
.gpx file from being saved.
'''
if reloadIdx:
self.index_ids()
elif G2obj.TestIndexAll():
self.index_ids()
if len(varlist) < 2:
raise Exception('add_EquivEquation Error: varlist must have at least 2 variables')
constr = []
if name is not None:
name = str(name)
typ_prev = None
vardefwarn = False
for i,var in enumerate(varlist):
m = 1.
try:
m = float(multlist[i])
except IndexError:
pass
# make var object
if isinstance(var, str):
var = self.make_var_obj(var,reloadIdx=False)
elif not isinstance(var, G2obj.G2VarObj):
var = self.make_var_obj(*var,reloadIdx=False)
if G2obj.getDescr(var.name) is None:
vardefwarn = True
print(f'Constraint var warning: No definition for variable {var.name}. Name correct?')
# make constraint
constr.append([m,var])
typ = _constr_type(var)
if typ_prev is None:
typ_prev = typ
var_prev = var
if typ_prev != typ:
msg = 'Type ({}) for var {} is different from {} ({})'.format(typ,var,var_prev,typ_prev)
raise Exception('add_EquivConstr Error: '+msg)
typ_prev = typ
var_prev = var
if vardefwarn and not override:
raise Exception('Constraint var error: undefined variables.\n\tIf correct, use override=True in call.\n\tAlso, please let GSAS-II developers know so definition can be added.')
constr += [name, bool(vary), 'f']
self.add_constraint_raw(typ, constr)
[docs]
def add_constraint_raw(self, cons_scope, constr):
"""Adds a constraint to the project.
:param str cons_scope: should be one of "Hist", "Phase", "HAP", or "Global".
:param list constr: a constraint coded with :class:`GSASIIobj.G2VarObj`
objects as described in the
:ref:`constraint definition descriptions <Constraint_definitions_table>`.
WARNING this function does not check the constraint is well-constructed.
Please use :meth:`G2Project.add_HoldConstr` or
:meth:`G2Project.add_EquivConstr` (etc.) instead, unless you are really
certain you know what you are doing.
"""
constrs = self.data['Constraints']['data']
if 'Global' not in constrs:
constrs['Global'] = []
constrs[cons_scope].append(constr)
[docs]
def hold_many(self, vars, ctype):
"""Apply holds for all the variables in vars, for constraint of a given type.
This routine has been superceeded by :meth:`add_Hold`
:param list vars: A list of variables to hold. Each may be a
:class:`GSASIIobj.G2VarObj` object, a variable name (str), or a
list/tuple of arguments for :meth:`make_var_obj`.
:param str ctype: A string constraint type specifier, passed directly to
:meth:`add_constraint_raw` as consType. Should be one of "Hist", "Phase",
or "HAP" ("Global" not implemented).
"""
G2fil.G2Print('G2Phase.hold_many Warning: replace calls to hold_many() with add_Hold()')
for var in vars:
if isinstance(var, str):
var = self.make_var_obj(var)
elif not isinstance(var, G2obj.G2VarObj):
var = self.make_var_obj(*var)
self.add_constraint_raw(ctype, [[1.0, var], None, None, 'h'])
[docs]
def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
reloadIdx=True):
"""Wrapper to create a G2VarObj. Takes either a string representation ("p:h:name:a")
or individual names of phase, histogram, varname, and atomId.
Automatically converts string phase, hist, or atom names into the ID required
by G2VarObj.
Note that this will cause the project to be saved if not
already done so.
"""
if reloadIdx:
self.index_ids()
elif G2obj.TestIndexAll():
self.index_ids()
# If string representation, short circuit
if hist is None and varname is None and atomId is None:
if isinstance(phase, str) and ':' in phase:
return G2obj.G2VarObj(phase)
# Get phase index
phaseObj = None
if isinstance(phase, G2Phase):
phaseObj = phase
phase = G2obj.PhaseRanIdLookup[phase.ranId]
elif phase in self.data['Phases']:
phaseObj = self.phase(phase)
phaseRanId = phaseObj.ranId
phase = G2obj.PhaseRanIdLookup[phaseRanId]
# Get histogram index
if isinstance(hist, G2PwdrData):
hist = G2obj.HistRanIdLookup[hist.ranId]
elif hist in self.data:
histRanId = self.histogram(hist).ranId
hist = G2obj.HistRanIdLookup[histRanId]
# Get atom index (if any)
if isinstance(atomId, G2AtomRecord):
atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
elif phaseObj:
atomObj = phaseObj.atom(atomId)
if atomObj:
atomRanId = atomObj.ranId
atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
return G2obj.G2VarObj(phase, hist, varname, atomId)
[docs]
def add_image(self, imagefile, fmthint=None, defaultImage=None,
indexList=None, cacheImage=False):
"""Load an image into a project
:param str imagefile: The image file to read, a filename.
:param str fmthint: If specified, only importers where the format name
(reader.formatName, as shown in Import menu) contains the
supplied string will be tried as importers. If not specified, all
importers consistent with the file extension will be tried
(equivalent to "guess format" in menu).
:param str defaultImage: The name of an image to use as a default for
setting parameters for the image file to read.
:param list indexList: specifies the image numbers (counting from zero)
to be used from the file when a file has multiple images. A value of
``[0,2,3]`` will cause the only first, third and fourth images in the file
to be included in the project.
:param bool cacheImage: When True, the image is cached to save
in rereading it later. Default is False (no caching).
:returns: a list of :class:`G2Image` object(s) for the added image(s)
"""
LoadG2fil()
imagefile = os.path.abspath(os.path.expanduser(imagefile))
readers = import_generic(imagefile, Readers['Image'], fmthint=fmthint)
objlist = []
for i,rd in enumerate(readers):
if indexList is not None and i not in indexList:
G2fil.G2Print("Image {} skipped".format(i))
continue
if rd.SciPy: #was default read by scipy; needs 1 time fixes
G2fil.G2Print('Warning: Image {} read by generic SciPy import. Image parameters likely wrong'.format(imagefile))
#see this: G2IO.EditImageParms(self,rd.Data,rd.Comments,rd.Image,imagefile)
rd.SciPy = False
rd.readfilename = imagefile
if rd.repeatcount == 1 and not rd.repeat: # skip image number if only one in set
rd.Data['ImageTag'] = None
else:
rd.Data['ImageTag'] = rd.repeatcount
rd.Data['formatName'] = rd.formatName
if rd.sumfile:
rd.readfilename = rd.sumfile
# Load generic metadata, as configured
G2fil.GetColumnMetadata(rd)
# Code from G2IO.LoadImage2Tree(rd.readfilename,self,rd.Comments,rd.Data,rd.Npix,rd.Image)
Imax = np.amax(rd.Image)
ImgNames = [i[0] for i in self.names if i[0].startswith('IMG ')]
TreeLbl = 'IMG '+os.path.basename(imagefile)
ImageTag = rd.Data.get('ImageTag')
if ImageTag:
TreeLbl += ' #'+'%04d'%(ImageTag)
imageInfo = (imagefile,ImageTag)
else:
imageInfo = imagefile
TreeName = G2obj.MakeUniqueLabel(TreeLbl,ImgNames)
# MT dict to contain image info
ImgDict = {}
ImgDict['data'] = [rd.Npix,imageInfo]
ImgDict['Comments'] = rd.Comments
if defaultImage:
if isinstance(defaultImage, G2Image):
if defaultImage.proj == self:
defaultImage = self.data[defaultImage.name]['data']
else:
raise Exception("Image {} not in current selected project".format(defaultImage.name))
elif defaultImage in self._images():
defaultImage = self.data[defaultImage]['data']
else:
defaultImage = None
Data = rd.Data
if defaultImage:
Data = copy.copy(defaultImage)
Data['showLines'] = True
Data['ring'] = []
Data['rings'] = []
Data['cutoff'] = 10.
Data['pixLimit'] = 20
Data['edgemin'] = 100000000
Data['calibdmin'] = 0.5
Data['calibskip'] = 0
Data['ellipses'] = []
Data['calibrant'] = ''
Data['GonioAngles'] = [0.,0.,0.]
Data['DetDepthRef'] = False
else:
Data['type'] = 'PWDR'
Data['color'] = GSASIIpath.GetConfigValue('Contour_color','Paired')
if 'tilt' not in Data: #defaults if not preset in e.g. Bruker importer
Data['tilt'] = 0.0
Data['rotation'] = 0.0
Data['pixLimit'] = 20
Data['calibdmin'] = 0.5
Data['cutoff'] = 10.
Data['showLines'] = False
Data['calibskip'] = 0
Data['ring'] = []
Data['rings'] = []
Data['edgemin'] = 100000000
Data['ellipses'] = []
Data['GonioAngles'] = [0.,0.,0.]
Data['DetDepth'] = 0.
Data['DetDepthRef'] = False
Data['calibrant'] = ''
Data['IOtth'] = [5.0,50.0]
Data['LRazimuth'] = [0.,180.]
Data['azmthOff'] = 0.0
Data['outChannels'] = 2500
Data['outAzimuths'] = 1
Data['centerAzm'] = False
Data['fullIntegrate'] = GSASIIpath.GetConfigValue('fullIntegrate',True)
Data['setRings'] = False
Data['background image'] = ['',-1.0]
Data['dark image'] = ['',-1.0]
Data['Flat Bkg'] = 0.0
Data['Oblique'] = [0.5,False]
Data['varyList'] = {'dist':True,'det-X':True,'det-Y':True,'tilt':True,'phi':True,'dep':False,'wave':False}
Data['setDefault'] = False
Data['range'] = [(0,Imax),[0,Imax]]
ImgDict['Image Controls'] = Data
ImgDict['Masks'] = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],
'Frames':[],'Thresholds':[(0,Imax),[0,Imax]],'SpotMask':{'esdMul':3.,'spotMask':None},}
ImgDict['Stress/Strain'] = {'Type':'True','d-zero':[],'Sample phi':0.0,
'Sample z':0.0,'Sample load':0.0}
self.names.append([TreeName]+['Comments','Image Controls','Masks','Stress/Strain'])
self.data[TreeName] = ImgDict
if cacheImage:
objlist.append(G2Image(self.data[TreeName], TreeName, self, image=rd.Image))
else:
objlist.append(G2Image(self.data[TreeName], TreeName, self))
del rd.Image
return objlist
[docs]
def imageMultiDistCalib(self,imageList=None,verbose=False):
'''Invokes a global calibration fit (same as Image Controls/Calibration/Multi-distance Recalibrate
menu command) with images as multiple distance settings.
Note that for this to work properly, the initial calibration parameters
(center, wavelength, distance & tilts) must be close enough to converge.
This may produce a better result if run more than once.
See :ref:`MultiDist_Example` for example code.
:param str imageList: the images to include in the fit, if not specified
all images in the project will be included.
:returns: parmDict,covData where parmDict has the refined parameters
and their values and covData is a dict containing the covariance matrix ('covMatrix'),
the number of ring picks ('obs') the reduced Chi-squared ('chisq'),
the names of the variables ('varyList') and their values ('variables')
'''
if imageList is None:
imageList = self.images()
# code based on GSASIIimgGUI..OnDistRecalib
obsArr = np.array([]).reshape(0,4)
parmDict = {}
varList = []
HKL = {}
for img in imageList:
name = img.name
G2fil.G2Print ('getting rings for',name)
Data = img.data['Image Controls']
key = str(int(Data['setdist']))
# create a parameter dict for combined fit
if 'wavelength' not in parmDict:
parmDict['wavelength'] = Data['wavelength']
if Data['varyList']['wave']:
varList += ['wavelength']
if Data['varyList']['dist']:
raise Exception(
'You cannot vary individual detector positions and the global wavelength.\n\nChange flags for 1st image.',
'Conflicting vars')
parmDict['dep'] = Data['DetDepth']
if Data['varyList']['dep']:
varList += ['dep']
# distance flag determines if individual values are refined
if not Data['varyList']['dist']:
# starts as zero, single variable, always refined
parmDict['deltaDist'] = 0.
varList += ['deltaDist']
parmDict['phi'] = Data['rotation']
if Data['varyList']['phi']:
varList += ['phi']
parmDict['tilt'] = Data['tilt']
if Data['varyList']['tilt']:
varList += ['tilt']
ImageZ = _getCorrImage(Readers['Image'],self,img)
Data['setRings'] = True
Masks = img.data['Masks']
result = G2img.ImageRecalibrate(None,ImageZ,Data,Masks,getRingsOnly=True)
if not len(result):
raise Exception('calibrant missing from local image calibrants files')
rings,HKL[key] = result
# add detector set dist into data array, create a single really large array
distarr = np.zeros_like(rings[:,2:3])
if 'setdist' not in Data:
raise Exception('Distance (setdist) not in image metadata')
distarr += Data['setdist']
obsArr = np.concatenate((
obsArr,
np.concatenate((rings[:,0:2],distarr,rings[:,2:3]),axis=1)),axis=0)
if 'deltaDist' not in parmDict:
# starts as zero, variable refined for each image
parmDict['delta'+key] = 0
varList += ['delta'+key]
for i,z in enumerate(['X','Y']):
v = 'det-'+z
if v+key in parmDict:
raise Exception('Error: two images with setdist ~=',key)
parmDict[v+key] = Data['center'][i]
if Data['varyList'][v]:
varList += [v+key]
#GSASIIpath.IPyBreak()
G2fil.G2Print('\nFitting',obsArr.shape[0],'ring picks and',len(varList),'variables...')
result = G2img.FitMultiDist(obsArr,varList,parmDict,covar=True,Print=verbose)
for img in imageList: # update GPX info with fit results
name = img.name
#print ('updating',name)
Data = img.data['Image Controls']
Data['wavelength'] = parmDict['wavelength']
key = str(int(Data['setdist']))
Data['center'] = [parmDict['det-X'+key],parmDict['det-Y'+key]]
if 'deltaDist' in parmDict:
Data['distance'] = Data['setdist'] - parmDict['deltaDist']
else:
Data['distance'] = Data['setdist'] - parmDict['delta'+key]
Data['rotation'] = np.mod(parmDict['phi'],360.0)
Data['tilt'] = parmDict['tilt']
Data['DetDepth'] = parmDict['dep']
N = len(Data['ellipses'])
Data['ellipses'] = [] #clear away individual ellipse fits
for H in HKL[key][:N]:
ellipse = G2img.GetEllipse(H[3],Data)
Data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
covData = {'title':'Multi-distance recalibrate','covMatrix':result[3],
'obs':obsArr.shape[0],'chisq':result[0],
'varyList':varList,'variables':result[1]}
return parmDict,covData
[docs]
def set_Frozen(self, variable=None, histogram=None, mode='remove'):
'''Removes one or more Frozen variables (or adds one)
(See :ref:`Parameter Limits<ParameterLimits>` description.)
Note that use of this
will cause the project to be saved if not already done so.
:param str variable: a variable name as a str or
(as a :class:`GSASIIobj.G2VarObj` object). Should
not contain wildcards.
If None (default), all frozen variables are deleted
from the project, unless a sequential fit and
a histogram is specified.
:param histogram: A reference to a histogram,
which can be reference by object, name, or number.
Used for sequential fits only.
:param str mode: The default mode is to remove variables
from the appropriate Frozen list, but if the mode
is specified as 'add', the variable is added to the
list.
:returns: True if the variable was added or removed, False
otherwise. Exceptions are generated with invalid requests.
'''
Controls = self.data['Controls']['data']
for key in ('parmMinDict','parmMaxDict','parmFrozen'):
if key not in Controls: Controls[key] = {}
if G2obj.TestIndexAll(): self.index_ids()
G2obj.patchControls(Controls)
if mode == 'remove':
if variable is None and histogram is None:
Controls['parmFrozen'] = {}
return True
elif mode == 'add':
if variable is None:
raise Exception('set_Frozen error: variable must be specified')
else:
raise Exception('Undefined mode ({}) in set_Frozen'.format(mode))
if histogram is None:
h = 'FrozenList'
else:
hist = self.histogram(histogram)
if hist:
h = hist.name
if h not in Controls['parmFrozen']: # no such list, not found
return False
else:
raise Exception('set_Frozen error: histogram {} not found'.format(histogram))
if mode == 'remove':
if variable is None:
Controls['parmFrozen'][h] = []
return True
if h not in Controls['parmFrozen']:
return True
delList = []
for i,v in enumerate(Controls['parmFrozen'][h]):
if v == variable: delList.append(i)
if delList:
for i in reversed(delList):
del Controls['parmFrozen'][h][i]
return True
return False
elif mode == 'add':
if type(variable) is str:
variable = G2obj.G2VarObj(variable)
elif type(v) is not G2obj.G2VarObj:
raise Exception(
'set_Frozen error: variable {} wrong type ({})'
.format(variable,type(variable)))
if h not in Controls['parmFrozen']: Controls['parmFrozen'][h] = []
Controls['parmFrozen'][h].append(variable)
return True
[docs]
def get_Frozen(self, histogram=None):
'''Gets a list of Frozen variables.
(See :ref:`Parameter Limits<ParameterLimits>` description.)
Note that use of this
will cause the project to be saved if not already done so.
:param histogram: A reference to a histogram,
which can be reference by object, name, or number. Used
for sequential fits only. If left as the default (None)
for a sequential fit, all Frozen variables in all
histograms are returned.
:returns: a list containing variable names, as str values
'''
Controls = self.data['Controls']['data']
for key in ('parmMinDict','parmMaxDict','parmFrozen'):
if key not in Controls: Controls[key] = {}
if G2obj.TestIndexAll(): self.index_ids()
G2obj.patchControls(Controls)
if histogram:
hist = self.histogram(histogram)
if hist:
h = hist.name
return [str(i) for i in Controls['parmFrozen'][h]]
elif histogram.lower() == 'all':
names = set()
for h in self.histograms():
if h.name not in Controls['parmFrozen']: continue
names = names.union([str(i) for i in Controls['parmFrozen'][h.name]])
return list(names)
else:
raise Exception('histogram {} not recognized'.format(histogram))
elif 'FrozenList' in Controls['parmFrozen']:
return [str(i) for i in Controls['parmFrozen']['FrozenList']]
else:
return []
[docs]
def get_Controls(self, control, variable=None):
'''Return project controls settings
:param str control: the item to be returned. See below for allowed values.
:param str variable: a variable name as a str or
(as a :class:`GSASIIobj.G2VarObj` object).
Used only with control set to "parmMin" or "parmMax".
:returns: The value for the control.
Allowed values for parameter control:
* cycles: the maximum number of cycles (returns int)
* sequential: the histograms used for a sequential refinement as a list
of histogram names or an empty list when in non-sequential mode.
* Reverse Seq: returns True or False. True indicates that fitting of the
sequence of histograms proceeds in reversed order.
* seqCopy: returns True or False. True indicates that results from
each sequential fit are used as the starting point for the next
histogram.
* parmMin & parmMax: retrieves a maximum or minimum value for
a refined parameter. Note that variable will be a GSAS-II
variable name, optionally with * specified for a histogram
or atom number. Return value will be a float.
(See :ref:`Parameter Limits<ParameterLimits>` description.)
* Anything else returns the value in the Controls dict, if present. An
exception is raised if the control value is not present.
.. seealso::
:meth:`set_Controls`
'''
if control.startswith('parmM'):
if not variable:
raise Exception('set_Controls requires a value for variable for control=parmMin or parmMax')
for key in ('parmMinDict','parmMaxDict','parmFrozen'):
if key not in self.data['Controls']['data']: self.data['Controls']['data'][key] = {}
if G2obj.TestIndexAll(): self.index_ids()
G2obj.patchControls(self.data['Controls']['data'])
if control == 'cycles':
return self.data['Controls']['data']['max cyc']
elif control == 'sequential':
return self.data['Controls']['data']['Seq Data']
elif control == 'Reverse Seq':
return self.data['Controls']['data']['Reverse Seq']
elif control == 'seqCopy':
return self.data['Controls']['data']['Copy2Next']
elif control == 'parmMin' or control == 'parmMax':
key = G2obj.G2VarObj(variable)
return G2obj.prmLookup(
variable,
self.data['Controls']['data'][control+'Dict'])[1]
elif control in self.data['Controls']['data']:
return self.data['Controls']['data'][control]
else:
G2fil.G2Print('Defined Controls:',self.data['Controls']['data'].keys())
raise Exception('{} is an invalid control value'.format(control))
[docs]
def set_Controls(self, control, value, variable=None):
'''Set project controls.
Note that use of this with control set to parmMin or parmMax
will cause the project to be saved if not already done so.
:param str control: the item to be set. See below for allowed values.
:param value: the value to be set.
:param str variable: used only with control set to "parmMin" or "parmMax"
Allowed values for *control* parameter:
* ``'cycles'``: sets the maximum number of cycles (value must be int)
* ``'sequential'``: sets the histograms to be used for a sequential refinement.
Use an empty list to turn off sequential fitting.
The values in the list may be the name of the histogram (a str), or
a ranId or index (int values), see :meth:`histogram`.
* ``'seqCopy'``: when True, the results from each sequential fit are used as
the starting point for the next. After each fit is is set to False.
Ignored for non-sequential fits.
* ``'Reverse Seq'``: when True, sequential refinement is performed on the
reversed list of histograms.
* ``'parmMin'`` & ``'parmMax'``: set a maximum or minimum value for a refined
parameter. Note that variable will be a GSAS-II variable name,
optionally with * specified for a histogram or atom number and
value must be a float.
(See :ref:`Parameter Limits<ParameterLimits>` description.)
.. seealso::
:meth:`get_Controls`
'''
if control.startswith('parmM'):
if not variable:
raise Exception('set_Controls requires a value for variable for control=parmMin or parmMax')
for key in ('parmMinDict','parmMaxDict','parmFrozen'):
if key not in self.data['Controls']['data']: self.data['Controls']['data'][key] = {}
if G2obj.TestIndexAll(): self.index_ids()
G2obj.patchControls(self.data['Controls']['data'])
if control == 'cycles':
self.data['Controls']['data']['max cyc'] = int(value)
elif control == 'seqCopy':
self.data['Controls']['data']['Copy2Next'] = bool(value)
elif control == 'Reverse Seq':
self.data['Controls']['data']['Reverse Seq'] = bool(value)
elif control == 'sequential':
histlist = []
for i,j in enumerate(value):
h = self.histogram(j)
if h:
histlist.append(h.name)
else:
raise Exception('item #{} ({}) is an invalid histogram value'
.format(i,j))
self.data['Controls']['data']['Seq Data'] = histlist
elif control == 'parmMin' or control == 'parmMax':
key = G2obj.G2VarObj(variable)
self.data['Controls']['data'][control+'Dict'][key] = float(value)
else:
raise Exception('{} is an invalid control value'.format(control))
[docs]
def copyHistParms(self,sourcehist,targethistlist='all',modelist='all'):
'''Copy histogram information from one histogram to others
:param sourcehist: is a histogram object (:class:`G2PwdrData`) or
a histogram name or the index number of the histogram
:param list targethistlist: a list of histograms where each item in the
list can be a histogram object (:class:`G2PwdrData`),
a histogram name or the index number of the histogram.
if the string 'all' (default value), then all histograms in
the project are used.
:param list modelist: May be a list of sections to copy, which may
include 'Background', 'Instrument Parameters', 'Limits' and
'Sample Parameters' (items may be shortened to uniqueness and
capitalization is ignored, so ['b','i','L','s'] will work.)
The default value, 'all' causes the listed sections to
'''
sections = ('Background','Instrument Parameters','Limits',
'Sample Parameters')
hist_in = self.histogram(sourcehist)
if not hist_in:
raise Exception('{} is not a valid histogram'.format(sourcehist))
if targethistlist == "all":
targethistlist = self.histograms()
if 'all' in modelist:
copysections = sections
else:
copysections = set()
for s in sections:
for m in modelist:
if s.lower().startswith(m.lower()):
copysections.add(s)
for h in targethistlist:
hist_out = self.histogram(h)
if not hist_out:
raise Exception('{} is not a valid histogram'.format(h))
for key in copysections:
hist_out[key] = copy.deepcopy(hist_in[key])
[docs]
def get_VaryList(self):
'''Returns a list of the refined variables in the
last refinement cycle
:returns: a list of variables or None if no refinement has been
performed.
'''
try:
return self['Covariance']['data']['varyList']
except:
return
[docs]
def get_ParmList(self):
'''Returns a list of all the parameters defined in the
last refinement cycle
:returns: a list of parameters or None if no refinement has been
performed.
'''
if 'parmDict' not in self['Covariance']['data']:
raise G2ScriptException('No parameters found in project, has a refinement been run?')
try:
return list(self['Covariance']['data']['parmDict'].keys())
except:
return
[docs]
def get_Variable(self,var):
'''Returns the value and standard uncertainty (esd) for a variable
parameters, as defined in the last refinement cycle
:param str var: a variable name of form '<p>:<h>:<name>', such as
':0:Scale'
:returns: (value,esd) if the parameter is refined or
(value, None) if the variable is in a constraint or is not
refined or None if the parameter is not found.
'''
if 'parmDict' not in self['Covariance']['data']:
raise G2ScriptException('No parameters found in project, has a refinement been run?')
if var not in self['Covariance']['data']['parmDict']:
return None
val = self['Covariance']['data']['parmDict'][var]
try:
pos = self['Covariance']['data']['varyList'].index(var)
esd = np.sqrt(self['Covariance']['data']['covMatrix'][pos,pos])
return (val,esd)
except ValueError:
return (val,None)
[docs]
def get_Covariance(self,varList):
'''Returns the values and covariance matrix for a series of variable
parameters. as defined in the last refinement cycle
:param tuple varList: a list of variable names of form '<p>:<h>:<name>'
:returns: (valueList,CovMatrix) where valueList contains the (n) values
in the same order as varList (also length n) and CovMatrix is a
(n x n) matrix. If any variable name is not found in the varyList
then None is returned.
Use this code, where sig provides standard uncertainties for
parameters and where covArray provides the correlation between
off-diagonal terms::
sig = np.sqrt(np.diag(covMatrix))
xvar = np.outer(sig,np.ones_like(sig))
covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
'''
missing = [i for i in varList if i not in self['Covariance']['data']['varyList']]
if missing:
G2fil.G2Print('Warning: Variable(s) {} were not found in the varyList'.format(missing))
return None
if 'parmDict' not in self['Covariance']['data']:
raise G2ScriptException('No parameters found in project, has a refinement been run?')
vals = [self['Covariance']['data']['parmDict'][i] for i in varList]
import GSASIImath as G2mth
cov = G2mth.getVCov(varList,
self['Covariance']['data']['varyList'],
self['Covariance']['data']['covMatrix'])
return (vals,cov)
[docs]
def ComputeWorstFit(self):
'''Computes the worst-fit parameters in a model.
:returns: (keys, derivCalcs, varyList) where:
* keys is a list of parameter names
where the names are ordered such that first entry in the list
will produce the largest change in the fit if refined and the last
entry will have the smallest change;
* derivCalcs is a dict where the key is a variable name and the
value is a list with three partial derivative values for
d(Chi**2)/d(var) where the derivatives are computed
for values v-d to v; v-d to v+d; v to v+d where v is
the current value for the variable and d is a small delta
value chosen for that variable type;
* varyList is a list of the parameters that are currently set to
be varied.
'''
self.save()
derivCalcs,varyList = G2strMain.Refine(self.filename,None,allDerivs=True)
keys = sorted(derivCalcs,key=lambda x:abs(derivCalcs[x][1]),reverse=True)
return (keys,derivCalcs,varyList)
[docs]
class G2AtomRecord(G2ObjectWrapper):
"""Wrapper for an atom record. Allows many atom properties to be access
and changed. See the :ref:`Atom Records description <Atoms_table>`
for the details on what information is contained in an atom record.
Scripts should not try to create a :class:`G2AtomRecord` object directly as
these objects are created via access from a :class:`G2Phase` object.
Example showing some uses of :class:`G2AtomRecord` methods:
>>> atom = some_phase.atom("O3")
>>> # We can access the underlying data structure (a list):
>>> atom.data
['O3', 'O-2', '', ... ]
>>> # We can also use wrapper accessors to get or change atom info:
>>> atom.coordinates
(0.33, 0.15, 0.5)
>>> atom.coordinates = [1/3, .1, 1/2]
>>> atom.coordinates
(0.3333333333333333, 0.1, 0.5)
>>> atom.refinement_flags
'FX'
>>> atom.ranId
4615973324315876477
>>> atom.occupancy
1.0
"""
def __init__(self, data, indices, proj):
self.data = data
self.cx, self.ct, self.cs, self.cia = indices
self.proj = proj
# @property
# def X(self):
# '''Get or set the associated atom's X.
# Use as ``x = atom.X`` to obtain the value and
# ``atom.X = x`` to set the value.
# '''
# pass
# @X.setter
# def X(self, val):
# pass
@property
def label(self):
'''Get the associated atom's label.
Use as ``x = atom.label`` to obtain the value and
``atom.label = x`` to set the value.
'''
return self.data[self.ct-1]
@label.setter
def label(self, val):
self.data[self.ct-1] = str(val)
@property
def type(self):
'''Get or set the associated atom's type. Call as a variable
(``x = atom.type``) to obtain the value or use
``atom.type = x`` to change the type. It is the user's
responsibility to make sure that the atom type is valid;
no checking is done here.
.. seealso::
:meth:`element`
'''
return self.data[self.ct]
@type.setter
def type(self, val):
# TODO: should check if atom type is defined
self.data[self.ct] = str(val)
@property
def element(self):
'''Parses element symbol from the atom type symbol for the atom
associated with the current object.
.. seealso::
:meth:`type`
'''
import re
try:
return re.match('^([A-Z][a-z]?)',self.data[self.ct]).group(1)
except:
raise Exception("element parse error with type {}".
format(self.data[self.ct]))
@property
def refinement_flags(self):
'''Get or set refinement flags for the associated atom.
Use as ``x = atom.refinement_flags`` to obtain the flags and
``atom.refinement_flags = "XU"`` (etc) to set the value.
'''
return self.data[self.ct+1]
@refinement_flags.setter
def refinement_flags(self, other):
# Automatically check it is a valid refinement
for c in other:
if c not in ' FXU':
raise ValueError("Invalid atom refinement: ", other)
self.data[self.ct+1] = other
@property
def coordinates(self):
'''Get or set the associated atom's coordinates.
Use as ``x = atom.coordinates`` to obtain a tuple with
the three (x,y,z) values and ``atom.coordinates = (x,y,z)``
to set the values.
Changes needed to adapt for changes in site symmetry have not yet been
implemented:
'''
return tuple(self.data[self.cx:self.cx+3])
@coordinates.setter
def coordinates(self, val):
if len(val) != 3:
raise ValueError(f"coordinates are of wrong length {val}")
try:
self.data[self.cx:self.cx+3] = [float(i) for i in val]
except:
raise ValueError(f"conversion error with coordinates {val}")
# TODO: should recompute the atom site symmetries here
@property
def occupancy(self):
'''Get or set the associated atom's site fraction.
Use as ``x = atom.occupancy`` to obtain the value and
``atom.occupancy = x`` to set the value.
'''
return self.data[self.cx+3]
@occupancy.setter
def occupancy(self, val):
self.data[self.cx+3] = float(val)
@property
def mult(self):
'''Get the associated atom's multiplicity value. Should not be
changed by user.
'''
return self.data[self.cs+1]
@property
def ranId(self):
'''Get the associated atom's Random Id number. Don't change this.
'''
return self.data[self.cia+8]
@property
def adp_flag(self):
'''Get the associated atom's iso/aniso setting. The value
will be 'I' or 'A'. No API provision is offered to change
this.
'''
# Either 'I' or 'A'
return self.data[self.cia]
@property
def ADP(self):
'''Get or set the associated atom's Uiso or Uaniso value(s).
Use as ``x = atom.ADP`` to obtain the value(s) and
``atom.ADP = x`` to set the value(s). For isotropic atoms
a single float value is returned (or used to set). For
anisotropic atoms a list of six values is used.
.. seealso::
:meth:`adp_flag`
:meth:`uiso`
'''
if self.adp_flag == 'I':
return self.data[self.cia+1]
else:
return self.data[self.cia+2:self.cia+8]
@ADP.setter
def ADP(self, value):
if self.adp_flag == 'I':
self.data[self.cia+1] = float(value)
else:
assert len(value) == 6
self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
@property
def uiso(self):
'''A synonym for :meth:`ADP` to be used for Isotropic atoms.
Get or set the associated atom's Uiso value.
Use as ``x = atom.uiso`` to obtain the value and
``atom.uiso = x`` to set the value. A
single float value is returned or used to set.
.. seealso::
:meth:`adp_flag`
:meth:`ADP`
'''
if self.adp_flag != 'I':
raise G2ScriptException(f"Atom {self.label} is not isotropic")
return self.ADP
@uiso.setter
def uiso(self, value):
if self.adp_flag != 'I':
raise G2ScriptException(f"Atom {self.label} is not isotropic")
self.ADP = value
[docs]
class G2PwdrData(G2ObjectWrapper):
"""Wraps a Powder Data Histogram.
The object contains these class variables:
* G2PwdrData.proj: contains a reference to the :class:`G2Project`
object that contains this histogram
* G2PwdrData.name: contains the name of the histogram
* G2PwdrData.data: contains the histogram's associated data in a dict,
as documented for the :ref:`Powder Diffraction Tree<Powder_table>`.
The actual histogram values are contained in the 'data' dict item,
as documented for Data.
Scripts should not try to create a :class:`G2PwdrData` object directly as
:meth:`G2PwdrData.__init__` should be invoked from inside :class:`G2Project`.
"""
def __init__(self, data, proj, name):
self.data = data
self.proj = proj
self.name = name
@staticmethod
def is_valid_refinement_key(key):
valid_keys = ['Limits', 'Sample Parameters', 'Background',
'Instrument Parameters']
return key in valid_keys
#@property
#def name(self):
# return self.data['data'][-1]
@property
def ranId(self):
return self.data['data'][0]['ranId']
@property
def residuals(self):
'''Provides a dictionary with with the R-factors for this histogram.
Includes the weighted and unweighted profile terms (R, Rb, wR, wRb, wRmin)
as well as the Bragg R-values for each phase (ph:H:Rf and ph:H:Rf^2).
'''
data = self.data['data'][0]
return {key: data[key] for key in data
if key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']
or ':' in key}
@property
def InstrumentParameters(self):
'''Provides a dictionary with with the Instrument Parameters
for this histogram.
'''
return self.data['Instrument Parameters'][0]
@property
def SampleParameters(self):
'''Provides a dictionary with with the Sample Parameters
for this histogram.
'''
return self.data['Sample Parameters']
@property
def Background(self):
'''Provides a list with with the Background parameters
for this histogram.
:returns: list containing a list and dict with background values
'''
return self.data['Background']
[docs]
def add_back_peak(self,pos,int,sig,gam,refflags=[]):
'''Adds a background peak to the Background parameters
:param float pos: position of peak, a 2theta or TOF value
:param float int: integrated intensity of background peak, usually large
:param float sig: Gaussian width of background peak, usually large
:param float gam: Lorentzian width of background peak, usually unused (small)
:param list refflags: a list of 1 to 4 boolean refinement flags for
pos,int,sig & gam, respectively (use [0,1] to refine int only).
Defaults to [] which means nothing is refined.
'''
if 'peaksList' not in self.Background[1]:
self.Background[1]['peaksList'] = []
flags = 4*[False]
for i,f in enumerate(refflags):
if i>3: break
flags[i] = bool(f)
bpk = []
for i,j in zip((pos,int,sig,gam),flags):
bpk += [float(i),j]
self.Background[1]['peaksList'].append(bpk)
self.Background[1]['nPeaks'] = len(self.Background[1]['peaksList'])
[docs]
def del_back_peak(self,peaknum):
'''Removes a background peak from the Background parameters
:param int peaknum: the number of the peak (starting from 0)
'''
npks = self.Background[1].get('nPeaks',0)
if peaknum >= npks:
raise Exception('peak {} not found in histogram {}'.format(peaknum,self.name))
del self.Background[1]['peaksList'][peaknum]
self.Background[1]['nPeaks'] = len(self.Background[1]['peaksList'])
[docs]
def ref_back_peak(self,peaknum,refflags=[]):
'''Sets refinement flag for a background peak
:param int peaknum: the number of the peak (starting from 0)
:param list refflags: a list of 1 to 4 boolean refinement flags for
pos,int,sig & gam, respectively. If a flag is not specified
it defaults to False (use [0,1] to refine int only).
Defaults to [] which means nothing is refined.
'''
npks = self.Background[1].get('nPeaks',0)
if peaknum >= npks:
raise Exception('peak {} not found in histogram {}'.format(peaknum,self.name))
flags = 4*[False]
for i,f in enumerate(refflags):
if i>3: break
flags[i] = bool(f)
for i,f in enumerate(flags):
self.Background[1]['peaksList'][peaknum][2*i+1] = f
@property
def id(self):
self.proj.update_ids()
return self.data['data'][0]['hId']
@id.setter
def id(self, val):
self.data['data'][0]['hId'] = val
[docs]
def Limits(self,typ,value=None):
'''Used to obtain or set the histogram limits. When a value is
specified, the appropriate limit is set. Otherwise, the value is
returned. Note that this provides an alternative to setting
histogram limits with the :meth:`G2Project:do_refinements` or
:meth:`G2PwdrData.set_refinements` methods.
:param str typ: a string which must be either 'lower'
(for 2-theta min or TOF min) or 'upper' (for 2theta max or TOF max).
Anything else produces an error.
:param float value: the number to set the limit (in units of degrees
or TOF (microsec.). If not specified, the command returns the
selected limit value rather than setting it.
:returns: The current value of the requested limit
(when ``value=None``). Units are
2-theta (degrees) or TOF (microsec).
Examples::
h = gpx.histogram(0)
val = h.Limits('lower')
h.Limits('upper',75)
'''
if value == None:
if typ == 'lower':
return self.data['Limits'][1][0]
elif typ == 'upper':
return self.data['Limits'][1][1]
else:
raise G2ScriptException('Limit errr: must be "lower" or "upper"')
else:
if typ == 'lower':
self.data['Limits'][1][0] = float(value)
elif typ == 'upper':
self.data['Limits'][1][1] = float(value)
else:
raise G2ScriptException('G2PwdData.Limit error: must be "lower" or "upper"')
[docs]
def Excluded(self,value=None):
'''Used to obtain or set the excluded regions for a histogram.
When a value is specified, the excluded regions are set.
Otherwise, the list of excluded region pairs is returned.
Note that excluded regions may be an empty list or a list
of regions to be excluded, where each region is provided
as pair of numbers, where the lower limit comes first.
Some sample excluded region lists are::
[[4.5, 5.5], [8.0, 9.0]]
[[130000.0, 140000.0], [160000.0, 170000.0]]
[]
The first above describes two excluded regions from 4.5-5.5 and 8-9
degrees 2-theta. The second is for a TOF pattern and also
describes two excluded regions, for 130-140 and 160-170 milliseconds.
The third line would be the case where there are no excluded
regions.
:param list value: A list of pairs of excluded region numbers
(as two-element lists). Some error checking/reformatting is done,
but users are expected to get this right. Use the GUI to
create examples or check input. Numbers in the list
are in units of degrees or TOF (microsec.).
If a value is not specified, the command returns the list of excluded regions.
:returns: The list of excluded regions (when ``value=None``). Units are
2-theta (degrees) or TOF (microsec).
Example 1::
h = gpx.histogram(0) # adds an excluded region (11-13 degrees)
h.Excluded(h.Excluded() + [[11,13]])
Example 2::
h = gpx.histogram(0) # changes the range of the first excluded region
excl = h.Excluded()
excl[0] = [120000.0, 160000.0] # microsec
h.Excluded(excl)
Example 3::
h = gpx.histogram(0) # deletes all excluded regions
h.Excluded([])
'''
if value is None:
return self.data['Limits'][2:]
# set up the excluded regions. Make sure we have a list of list pairs
s = ''
listValues = []
try:
for i,(l,h) in enumerate(value): # quick bit of validation
listValues.append([
min(float(l),float(h)),
max(float(l),float(h)),
])
if float(l) < self.data['Limits'][0][0]:
s += f'Lower limit {l} too low. '
if float(h) > self.data['Limits'][0][1]:
s += f'Upper limit {h} too high. '
except:
raise G2ScriptException(f'G2PwdData.Excluded error: incorrectly formatted list or\n\tinvalid value. In: {value}')
if s:
raise G2ScriptException(f'G2PwdData.Excluded error(s): {s}')
self.data['Limits'][2:] = listValues
[docs]
def fit_fixed_points(self):
"""Attempts to apply a background fit to the fixed points currently specified."""
def SetInstParms(Inst):
dataType = Inst['Type'][0]
insVary = []
insNames = []
insVals = []
for parm in Inst:
insNames.append(parm)
insVals.append(Inst[parm][1])
if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
Inst[parm][2] = False
instDict = dict(zip(insNames, insVals))
instDict['X'] = max(instDict['X'], 0.01)
instDict['Y'] = max(instDict['Y'], 0.01)
if 'SH/L' in instDict:
instDict['SH/L'] = max(instDict['SH/L'], 0.002)
return dataType, instDict, insVary
import scipy.interpolate as si
bgrnd = self.data['Background']
# Need our fixed points in order
bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
X = [x for x, y in bgrnd[1]['FixedPoints']]
Y = [y for x, y in bgrnd[1]['FixedPoints']]
limits = self.data['Limits'][1]
if X[0] > limits[0]:
X = [limits[0]] + X
Y = [Y[0]] + Y
if X[-1] < limits[1]:
X += [limits[1]]
Y += [Y[-1]]
# Some simple lookups
controls = self.proj['Controls']['data']
inst, inst2 = self.data['Instrument Parameters']
pwddata = self.data['data'][1]
# Construct the data for background fitting
xBeg = np.searchsorted(pwddata[0], limits[0])
xFin = np.searchsorted(pwddata[0], limits[1])
xdata = pwddata[0][xBeg:xFin]
ydata = si.interp1d(X,Y)(ma.getdata(xdata))
W = [1]*len(xdata)
Z = [0]*len(xdata)
dataType, insDict, insVary = SetInstParms(inst)
bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
# Do the fit
data = np.array([xdata, ydata, W, Z, Z, Z])
G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
prevVaryList=bakVary, controls=controls)
# Post-fit
parmDict = {}
bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
parmDict.update(bakDict)
parmDict.update(insDict)
pwddata[3][xBeg:xFin] *= 0
pwddata[5][xBeg:xFin] *= 0
pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
# TODO adjust pwddata? GSASIIpwdGUI.py:1041
# TODO update background
self.proj.save()
[docs]
def getdata(self,datatype):
'''Provides access to the histogram data of the selected data type
:param str datatype: must be one of the following values
(case is ignored):
* 'X': the 2theta or TOF values for the pattern
* 'Q': the 2theta or TOF values for the pattern transformed to Q
* 'd': the 2theta or TOF values for the pattern transformed to d-space
* 'Yobs': the observed intensity values
* 'Yweight': the weights for each data point (1/sigma**2)
* 'Ycalc': the computed intensity values
* 'Background': the computed background values
* 'Residual': the difference between Yobs and Ycalc (obs-calc)
:returns: an numpy MaskedArray with data values of the requested type
'''
enums = ['x', 'yobs', 'yweight', 'ycalc', 'background', 'residual', 'q', 'd']
if datatype.lower() not in enums:
raise G2ScriptException("Invalid datatype = "+datatype+" must be one of "+str(enums))
if datatype.lower() == 'q':
Inst,Inst2 = self.data['Instrument Parameters']
return 2 * np.pi / G2lat.Pos2dsp(Inst,self.data['data'][1][0])
elif datatype.lower() == 'd':
Inst,Inst2 = self.data['Instrument Parameters']
return G2lat.Pos2dsp(Inst,self.data['data'][1][0])
else:
return self.data['data'][1][enums.index(datatype.lower())]
[docs]
def y_calc(self):
'''Returns the calculated intensity values; better to
use :meth:`getdata`
'''
return self.data['data'][1][3]
[docs]
def reflections(self):
'''Returns a dict with an entry for every phase in the
current histogram. Within each entry is a dict with keys
'RefList' (reflection list, see
:ref:`Powder Reflections <PowderRefl_table>`),
'Type' (histogram type), 'FF'
(form factor information), 'Super' (True if this is superspace
group).
'''
return self.data['Reflection Lists']
[docs]
def Export(self,fileroot,extension,fmthint=None):
'''Write the histogram into a file. The path is specified by fileroot and
extension.
:param str fileroot: name of the file, optionally with a path (extension is
ignored)
:param str extension: includes '.', must match an extension in global
exportersByExtension['powder'] or a Exception is raised.
:param str fmthint: If specified, the first exporter where the format
name (obj.formatName, as shown in Export menu) contains the
supplied string will be used. If not specified, an error
will be generated showing the possible choices.
:returns: name of file that was written
'''
LoadG2fil()
if extension not in exportersByExtension.get('powder',[]):
print('Defined exporters are:')
print(' ',list(exportersByExtension.get('powder',[])))
raise G2ScriptException('No Writer for file type = "'+extension+'"')
fil = os.path.abspath(os.path.splitext(fileroot)[0]+extension)
obj = exportersByExtension['powder'][extension]
if type(obj) is list:
if fmthint is None:
print('Defined ',extension,'exporters are:')
for o in obj:
print('\t',o.formatName)
raise G2ScriptException('No format hint for file type = "'+extension+'"')
for o in obj:
if fmthint.lower() in o.formatName.lower():
obj = o
break
else:
print('Hint ',fmthint,'not found. Defined ',extension,'exporters are:')
for o in obj:
print('\t',o.formatName)
raise G2ScriptException('Bad format hint for file type = "'+extension+'"')
self._SetFromArray(obj)
obj.Writer(self.name,fil)
return fil
def _SetFromArray(self,expObj):
'''Load a histogram into the exporter in preparation for use of
the .Writer method in the object.
:param Exporter expObj: Exporter object
'''
expObj.Histograms[self.name] = {}
expObj.Histograms[self.name]['Data'] = self.data['data'][1]
for key in 'Instrument Parameters','Sample Parameters','Reflection Lists':
expObj.Histograms[self.name][key] = self.data[key]
def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
try:
import matplotlib.pyplot as plt
except ImportError:
G2fil.G2Print('Warning: Unable to import matplotlib, skipping plot')
return
data = self.data['data'][1]
if Yobs:
plt.plot(data[0], data[1], label='Yobs')
if Ycalc:
plt.plot(data[0], data[3], label='Ycalc')
if Background:
plt.plot(data[0], data[4], label='Background')
if Residual:
plt.plot(data[0], data[5], label="Residual")
[docs]
def get_wR(self):
"""returns the overall weighted profile R factor for a histogram
:returns: a wR value as a percentage or None if not defined
"""
return self['data'][0].get('wR')
def _decodeHist(self,hist):
'''Convert a histogram reference to a histogram name string
'''
if isinstance(hist, G2PwdrData):
return hist.name
elif isinstance(hist, G2Single):
return hist.name
elif hist in [h.name for h in self.proj.histograms()]:
return hist
elif type(hist) is int:
return self.proj.histograms()[hist].name
else:
raise G2ScriptException("Invalid histogram reference: "+str(hist))
[docs]
def set_background(self, key, value):
'''Set background parameters (this serves a similar function as in
:meth:`set_refinements`, but with a simplified interface).
:param str key: a string that defines the background parameter that will
be changed. Must appear in the table below.
================= ============== ===========================================
key name type of value meaning of value
================= ============== ===========================================
fixedHist int, str, reference to a histogram in the current
None or project or None to remove the reference.
G2PwdrData
fixedFileMult float multiplier applied to intensities in
the background histogram where a value
of -1.0 means full subtraction of
the background histogram.
================= ============== ===========================================
:param value: a value to set the selected background parameter. The meaning
and type for this parameter is listed in the table above.
'''
bkgPrms, bkgDict = self.data['Background']
if key == 'fixedHist':
if value is None:
bkgDict['background PWDR'][0] = ''
return
bkgDict['background PWDR'][0] = self._decodeHist(value)
elif key == 'fixedFileMult':
bkgDict['background PWDR'][1] = float(value)
else:
raise ValueError("Invalid key in set_background:", key)
[docs]
def set_refinements(self, refs):
"""Sets the PWDR histogram refinement parameter 'key' to the specification 'value'.
:param dict refs: A dictionary of the parameters to be set. See the
:ref:`Histogram_parameters_table` table for a description of
what these dictionaries should be.
:returns: None
"""
do_fit_fixed_points = False
for key, value in refs.items():
if key == 'Limits':
old_limits = self.data['Limits'][1]
new_limits = value
if isinstance(new_limits, dict):
if 'low' in new_limits:
old_limits[0] = new_limits['low']
if 'high' in new_limits:
old_limits[1] = new_limits['high']
else:
old_limits[0], old_limits[1] = new_limits
elif key == 'Sample Parameters':
sample = self.data['Sample Parameters']
for sparam in value:
if sparam not in sample:
raise ValueError("Unknown refinement parameter, "
+ str(sparam))
sample[sparam][1] = True
elif key == 'Background':
bkg, peaks = self.data['Background']
# If True or False, just set the refine parameter
if value in (True, False):
bkg[1] = value
return
if 'type' in value:
bkg[0] = value['type']
if 'refine' in value:
bkg[1] = value['refine']
if 'no. coeffs' in value:
cur_coeffs = bkg[2]
n_coeffs = value['no. coeffs']
if n_coeffs > cur_coeffs:
for x in range(n_coeffs - cur_coeffs):
bkg.append(0.0)
else:
for _ in range(cur_coeffs - n_coeffs):
bkg.pop()
bkg[2] = n_coeffs
if 'coeffs' in value:
bkg[3:] = value['coeffs']
if 'FixedPoints' in value:
peaks['FixedPoints'] = [(float(a), float(b))
for a, b in value['FixedPoints']]
if value.get('fit fixed points', False):
do_fit_fixed_points = True
if 'peaks' in value:
for i,flags in enumerate(value['peaks']):
self.ref_back_peak(i,flags)
elif key == 'Instrument Parameters':
instrument, secondary = self.data['Instrument Parameters']
for iparam in value:
try:
instrument[iparam][2] = True
except IndexError:
raise ValueError("Invalid key:", iparam)
else:
raise ValueError("Unknown key:", key)
# Fit fixed points after the fact - ensure they are after fixed points
# are added
if do_fit_fixed_points:
# Background won't be fit if refinement flag not set
orig = self.data['Background'][0][1]
self.data['Background'][0][1] = True
self.fit_fixed_points()
# Restore the previous value
self.data['Background'][0][1] = orig
[docs]
def clear_refinements(self, refs):
"""Clears the PWDR refinement parameter 'key' and its associated value.
:param dict refs: A dictionary of parameters to clear.
See the :ref:`Histogram_parameters_table` table for what can be specified.
"""
for key, value in refs.items():
if key == 'Limits':
old_limits, cur_limits = self.data['Limits']
cur_limits[0], cur_limits[1] = old_limits
elif key == 'Sample Parameters':
sample = self.data['Sample Parameters']
for sparam in value:
sample[sparam][1] = False
elif key == 'Background':
bkg, peaks = self.data['Background']
# If True or False, just set the refine parameter
if value in (True, False):
bkg[1] = False
return
bkg[1] = False
if 'FixedPoints' in value:
if 'FixedPoints' in peaks:
del peaks['FixedPoints']
if 'peaks' in value:
for i in range(len(self.Background[1]['peaksList'])):
self.ref_back_peak(i,[])
elif key == 'Instrument Parameters':
instrument, secondary = self.data['Instrument Parameters']
for iparam in value:
instrument[iparam][2] = False
else:
raise ValueError("Unknown key:", key)
[docs]
def add_peak(self,area,dspace=None,Q=None,ttheta=None):
'''Adds a single peak to the peak list
:param float area: peak area
:param float dspace: peak position as d-space (A)
:param float Q: peak position as Q (A-1)
:param float ttheta: peak position as 2Theta (deg)
Note: only one of the parameters: dspace, Q or ttheta may be specified.
See :ref:`PeakRefine` for an example.
'''
import GSASIImath as G2mth
if (not dspace) + (not Q) + (not ttheta) != 2:
G2fil.G2Print('add_peak error: too many or no peak position(s) specified')
return
pos = ttheta
Parms,Parms2 = self.data['Instrument Parameters']
if Q:
pos = G2lat.Dsp2pos(<