Skip to content

Commit 155807b

Browse files
pavoljuhaspavoljuhas
committed
Merge branch add-sfaverage.
2 parents 8f7408d + a135e24 commit 155807b

7 files changed

Lines changed: 376 additions & 13 deletions

File tree

‎diffpy/srreal/scatteringfactortable.py‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,21 +14,23 @@
1414
##############################################################################
1515

1616

17-
"""class ScatteringFactorTable -- scattering factors for atoms, ions and
18-
isotopes.
17+
"""\
18+
class ScatteringFactorTable -- scattering factors for atoms, ions and isotopes.
1919
"""
2020

2121

2222
# exported items, these also makes them show in pydoc.
2323
__all__ = ['ScatteringFactorTable',
24-
'SFTXray', 'SFTElectron', 'SFTNeutron', 'SFTElectronNumber']
24+
'SFTXray', 'SFTElectron', 'SFTNeutron', 'SFTElectronNumber',
25+
'SFAverage']
2526

2627
from diffpy.srreal.srreal_ext import ScatteringFactorTable
2728
from diffpy.srreal.srreal_ext import SFTXray
2829
from diffpy.srreal.srreal_ext import SFTElectron
2930
from diffpy.srreal.srreal_ext import SFTNeutron
3031
from diffpy.srreal.srreal_ext import SFTElectronNumber
3132
from diffpy.srreal.wraputils import _pickle_getstate, _pickle_setstate
33+
from diffpy.srreal.sfaverage import SFAverage
3234

3335
# Pickling Support -----------------------------------------------------------
3436

‎diffpy/srreal/sfaverage.py‎

Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
#!/usr/bin/env python
2+
##############################################################################
3+
#
4+
# diffpy.srreal Complex Modeling Initiative
5+
# (c) 2016 Brookhaven Science Associates,
6+
# Brookhaven National Laboratory.
7+
# All rights reserved.
8+
#
9+
# File coded by: Pavol Juhas
10+
#
11+
# See AUTHORS.txt for a list of people who contributed.
12+
# See LICENSE.txt for license information.
13+
#
14+
##############################################################################
15+
16+
17+
"""\
18+
Compositional averaging of atom scattering factors.
19+
20+
21+
Examples
22+
--------
23+
::
24+
25+
import numpy as np
26+
from diffpy.srreal.scatteringfactortable import SFTXray
27+
xtb = SFTXray()
28+
qa = np.linspace(0, 25)
29+
sfavg1 = SFAverage.fromComposition(xtb, {'Na' : 1, 'Cl' : 1}, qa)
30+
sfavg2 = SFAverage.fromComposition(xtb, [('Mn', 1), ('O', 2)], qa)
31+
32+
from diffpy.Structure import loadStructure
33+
from diffpy.srreal.pdfcalculator import DebyePDFCalculator
34+
dpdfc = DebyePDFCalculator()
35+
dpdfc(loadStructure('rutile.cif'))
36+
sfavg3 = SFAverage.fromStructure(dpdfc.scatteringfactortable,
37+
dpdfc.getStructure(), dpdfc.qgrid)
38+
"""
39+
40+
41+
# class SFAverage ------------------------------------------------------------
42+
43+
class SFAverage(object):
44+
"""\
45+
Calculate compositional statistics of atom scattering factors.
46+
47+
Compositional averages can be calculated for an array of Q-values.
48+
Results are stored in the class attributes.
49+
50+
Attributes
51+
----------
52+
f1sum :
53+
Sum of scattering factors from all atoms.
54+
Float or NumPy array.
55+
f2sum :
56+
Sum of squared scattering factors from all atoms.
57+
Float or NumPy array.
58+
count :
59+
Total number of atoms. Can be non-integer in case of
60+
fractional occupancies.
61+
f1avg :
62+
Compositional average of scattering factors.
63+
Float or NumPy array.
64+
f2avg :
65+
Compositional average of squared scattering factors.
66+
Float or NumPy array.
67+
composition :
68+
Dictionary of atom symbols and their total abundancies.
69+
"""
70+
71+
f1sum = 0
72+
f2sum = 0
73+
count = 0
74+
f1avg = 0
75+
f2avg = 0
76+
composition = None
77+
78+
@classmethod
79+
def fromStructure(cls, sftb, stru, q=0):
80+
"""\
81+
Calculate average scattering factors from a structure object.
82+
83+
Parameters
84+
----------
85+
sftb : ScatteringFactorTable
86+
The ScatteringFactorTable object for looking up the values.
87+
stru : diffpy Structure or pyobjcryst Crystal or StructureAdapter
88+
The structure object that stores the atom species and their
89+
occupancies. Can be any type with a registered conversion
90+
to the StructureAdapter class.
91+
q : float or NumPy array (optional)
92+
The Q value in inverse Angstroms for which to lookup
93+
the scattering factor values.
94+
95+
See also
96+
--------
97+
RegisterStructureAdapter : to add support for more structure types.
98+
99+
Returns
100+
-------
101+
SFAverage
102+
The calculated scattering factor averages.
103+
"""
104+
# a bit of duck-typing for faster handling of diffpy.Structure
105+
if hasattr(type(stru), 'composition'):
106+
composition = stru.composition
107+
if isinstance(composition, dict):
108+
return cls.fromComposition(sftb, composition, q)
109+
# otherwise let's convert to a known structure type
110+
from diffpy.srreal.structureadapter import createStructureAdapter
111+
adpt = createStructureAdapter(stru)
112+
composition = {}
113+
for i in range(adpt.countSites()):
114+
smbl = adpt.siteAtomType(i)
115+
cnt = adpt.siteOccupancy(i) * adpt.siteMultiplicity(i)
116+
composition[smbl] = composition.get(smbl, 0) + cnt
117+
return cls.fromComposition(sftb, composition, q)
118+
119+
120+
@classmethod
121+
def fromComposition(cls, sftb, composition, q=0):
122+
"""\
123+
Calculate average scattering factors from atom concentrations.
124+
125+
Parameters
126+
----------
127+
sftb : ScatteringFactorTable
128+
The ScatteringFactorTable object for looking up the values.
129+
composition : dictionary or a list of (symbol, amount) pairs.
130+
The chemical composition for evaluating the average. Atom
131+
symbols may repeat when it is a list of (symbol, amount) pairs.
132+
q : float or NumPy array (optional)
133+
The Q value in inverse Angstroms for which to lookup
134+
the scattering factor values.
135+
136+
Returns
137+
-------
138+
SFAverage
139+
The calculated scattering factor averages.
140+
"""
141+
sfa = cls()
142+
sfa.composition = {}
143+
if isinstance(composition, dict):
144+
sfa.composition.update(composition)
145+
else:
146+
for smbl, cnt in composition:
147+
if not smbl in sfa.composition:
148+
sfa.composition[smbl] = 0
149+
sfa.composition[smbl] += cnt
150+
sfa.f1sum = 0.0 * q
151+
sfa.f2sum = 0.0 * q
152+
for smbl, cnt in sfa.composition.items():
153+
sfq = sftb.lookup(smbl, q)
154+
sfa.f1sum += cnt * sfq
155+
sfa.f2sum += cnt * sfq**2
156+
sfa.count += cnt
157+
denom = sfa.count if sfa.count > 0 else 1
158+
sfa.f1avg = sfa.f1sum / denom
159+
sfa.f2avg = sfa.f2sum / denom
160+
return sfa
161+
162+
# End of class SFAverage

‎diffpy/srreal/tests/testscatteringfactortable.py‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66

77
import unittest
88
import cPickle
9+
import numpy
910

1011
from diffpy.srreal.scatteringfactortable import ScatteringFactorTable
1112

@@ -130,6 +131,23 @@ def test_derived_clone(self):
130131
self.assertEqual(set(['Xy']), lsft2.getCustomSymbols())
131132
return
132133

134+
def test_lookup(self):
135+
"""Check ScatteringFactorTable.lookup handling of array arguments.
136+
"""
137+
qa = numpy.linspace(0, 50)
138+
sftx = self.sftx
139+
fmn0 = numpy.array([sftx.lookup('Mn', x) for x in qa])
140+
fmn1 = sftx.lookup('Mn', qa)
141+
self.assertTrue(numpy.array_equal(fmn0, fmn1))
142+
self.assertTrue(numpy.array_equal(
143+
fmn0.reshape(5, 10), sftx.lookup('Mn', qa.reshape(5, 10))))
144+
self.assertTrue(numpy.array_equal(
145+
fmn0.reshape(5, 2, 5), sftx.lookup('Mn', qa.reshape(5, 2, 5))))
146+
self.assertTrue(numpy.array_equal(fmn0, sftx.lookup('Mn', list(qa))))
147+
self.assertRaises(TypeError, sftx.lookup, 'Na', 'asdf')
148+
self.assertRaises(TypeError, sftx.lookup, 'Na', {})
149+
return
150+
133151
# End of class TestScatteringFactorTable
134152

135153
if __name__ == '__main__':

‎diffpy/srreal/tests/testsfaverage.py‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#!/usr/bin/env python
2+
3+
"""Unit tests for diffpy.srreal.sfaverage
4+
"""
5+
6+
7+
import unittest
8+
import numpy
9+
10+
from diffpy.srreal.sfaverage import SFAverage
11+
from diffpy.srreal.scatteringfactortable import ScatteringFactorTable
12+
from diffpy.srreal.tests.testutils import loadDiffPyStructure
13+
from diffpy.srreal.tests.testutils import TestCaseObjCrystOptional
14+
from diffpy.srreal.tests.testutils import loadObjCrystCrystal
15+
16+
17+
##############################################################################
18+
class TestSFAverage(unittest.TestCase):
19+
20+
def setUp(self):
21+
self.sftx = ScatteringFactorTable.createByType('X')
22+
self.sftn = ScatteringFactorTable.createByType('N')
23+
return
24+
25+
def tearDown(self):
26+
return
27+
28+
29+
def test_fromStructure_CdSe(self):
30+
"""check SFAverage.fromStructure() for CdSe
31+
"""
32+
cdse = loadDiffPyStructure('CdSe_cadmoselite.cif')
33+
sfavg = SFAverage.fromStructure(self.sftx, cdse)
34+
fcd = self.sftx.lookup('Cd')
35+
fse = self.sftx.lookup('Se')
36+
self.assertTrue(isinstance(sfavg.f1sum, float))
37+
self.assertAlmostEqual(0.5 * (fcd + fse), sfavg.f1avg)
38+
self.assertAlmostEqual(0.5 * (fcd**2 + fse**2), sfavg.f2avg)
39+
self.assertEqual(4, sfavg.count)
40+
self.assertEqual(cdse.composition, sfavg.composition)
41+
qa = numpy.arange(0, 25, 0.1)
42+
sfavg2 = SFAverage.fromStructure(self.sftx, cdse, qa)
43+
self.assertTrue(isinstance(sfavg2.f1sum, numpy.ndarray))
44+
self.assertNotEqual(sfavg2.f1sum[0], sfavg2.f1sum[-1])
45+
self.assertEqual(sfavg.f1sum, sfavg2.f1sum[0])
46+
self.assertEqual(sfavg.f2sum, sfavg2.f2sum[0])
47+
sfavg3 = SFAverage.fromStructure(self.sftn, cdse, qa)
48+
self.assertEqual(sfavg3.f1sum[0], sfavg3.f1sum[-1])
49+
self.assertRaises(TypeError, SFAverage.fromStructure,
50+
self.sftx, 'notastructure')
51+
return
52+
53+
54+
def test_fromComposition(self):
55+
"""check SFAverage.fromComposition()
56+
"""
57+
sfavg1 = SFAverage.fromComposition(self.sftx, {'Na' : 1, 'Cl' : 1})
58+
fm = ['Na', 0.25, 'Cl', 0.75, 'Cl', 0.25, 'Na', 0.5, 'Na', 0.25]
59+
smblcnts = zip(fm[0::2], fm[1::2])
60+
sfavg2 = SFAverage.fromComposition(self.sftx, smblcnts)
61+
self.assertEqual(sfavg1.f1sum, sfavg2.f1sum)
62+
self.assertEqual(sfavg1.f2sum, sfavg2.f2sum)
63+
self.assertEqual(sfavg1.count, sfavg2.count)
64+
sfempty = SFAverage.fromComposition(self.sftx, [])
65+
self.assertEqual(0, sfempty.f1avg)
66+
self.assertEqual(0, sfempty.f2avg)
67+
self.assertEqual(0, sfempty.count)
68+
return
69+
70+
# End of class TestSFAverage
71+
72+
##############################################################################
73+
class TestSFAverageObjCryst(TestCaseObjCrystOptional):
74+
75+
def setUp(self):
76+
self.sftx = ScatteringFactorTable.createByType('X')
77+
return
78+
79+
def tearDown(self):
80+
return
81+
82+
83+
def test_from_rutile(self):
84+
"""check SFAverage.fromStructure for pyobjcryst Crystal of rutile.
85+
"""
86+
rutile = loadObjCrystCrystal('rutile.cif')
87+
qa = numpy.arange(0, 25, 0.1)
88+
sfavg = SFAverage.fromStructure(self.sftx, rutile, qa)
89+
fti = self.sftx.lookup('Ti', qa)
90+
fo = self.sftx.lookup('O', qa)
91+
self.assertTrue(numpy.allclose((fti + 2 * fo) / 3.0, sfavg.f1avg))
92+
fti2, fo2 = fti**2, fo**2
93+
self.assertTrue(numpy.allclose((fti2 + 2 * fo2) / 3.0, sfavg.f2avg))
94+
self.assertEqual(6, sfavg.count)
95+
96+
# End of class TestSFAverageObjCryst
97+
98+
if __name__ == '__main__':
99+
unittest.main()
100+
101+
# End of file

‎srrealmodule/srreal_converters.cpp‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,22 @@ NumPyArray_DoublePtr createNumPyDoubleArray(int dim, const int* sz)
106106
}
107107

108108

109+
/// helper for creating numpy array of the same shape as the argument
110+
NumPyArray_DoublePtr createNumPyDoubleArrayLike(boost::python::object& obj)
111+
{
112+
PyObject* arr = obj.ptr();
113+
int dim = PyArray_NDIM(arr);
114+
npy_intp* shape = PyArray_DIMS(arr);
115+
// create numpy array
116+
boost::python::object rvobj(
117+
boost::python::handle<>(
118+
PyArray_SimpleNew(dim, shape, PyArray_DOUBLE)));
119+
double* rvdata = static_cast<double*>(PyArray_DATA(rvobj.ptr()));
120+
NumPyArray_DoublePtr rv(rvobj, rvdata);
121+
return rv;
122+
}
123+
124+
109125
/// helper for creating a numpy array view on a double array
110126
boost::python::object
111127
createNumPyDoubleView(double* data, int dim, const int* sz)
@@ -248,6 +264,24 @@ extractQuantityType(
248264
}
249265

250266

267+
/// efficient conversion of Python object to a numpy array of doubles
268+
NumPyArray_DoublePtr extractNumPyDoubleArray(::boost::python::object& obj)
269+
{
270+
PyObject* arr = PyArray_ContiguousFromAny(obj.ptr(), PyArray_DOUBLE, 0, 0);
271+
if (!arr)
272+
{
273+
const char* emsg = "Cannot convert this object to numpy array.";
274+
PyErr_SetString(PyExc_TypeError, emsg);
275+
boost::python::throw_error_already_set();
276+
assert(false);
277+
}
278+
boost::python::object rvobj((boost::python::handle<>(arr)));
279+
double* rvdata = static_cast<double*>(PyArray_DATA(arr));
280+
NumPyArray_DoublePtr rv(rvobj, rvdata);
281+
return rv;
282+
}
283+
284+
251285
/// extract integer with a support for numpy.int types
252286
int extractint(boost::python::object obj)
253287
{

0 commit comments

Comments
 (0)