Commit ed42c6f7 authored by Paul Erhart's avatar Paul Erhart

BLD: prepared dynasor for distribution

* upgraded from Python2 to Python3
* prepared setup.py for pip installation
* fixed docstrings in filon and trajectory
* improved Python interface documentation in user guide
parent ed94a3d2
......@@ -5,7 +5,10 @@ build_linux:
tags:
- linux
script:
- python ./setup.py install --prefix=./opt
- INSTDIR=$PWD/opt/lib/python
- mkdir -p $INSTDIR
- export PYTHONPATH=$INSTDIR:$PYTHONPATH
- python3 ./setup.py install --home=$PWD/opt
artifacts:
paths:
- opt
......@@ -18,6 +21,11 @@ run_tests:
dependencies:
- build_linux
script:
- export PATH=${PATH}:${HOME}/.local/bin
# Python path used throughout
- INSTDIR=$PWD/opt/lib/python
- export PYTHONPATH=$INSTDIR:$PYTHONPATH
# run tests
- cd tests
- echo Do stuff, e.g.,
- echo PYTHONPATH=$PWD/../opt python run-tests.py
......@@ -30,15 +38,20 @@ pages:
dependencies:
- build_linux
script:
- export PYTHONPATH=$PWD/opt/lib/python2.7/site-packages
- sphinx-apidoc -e -o doc/ dsf/
# build documentation
- INSTDIR=$PWD/opt/lib/python
- export PYTHONPATH=$INSTDIR:$PYTHONPATH
- sphinx-build doc/ public/
# site verification
- mkdir -p public/.well-known/acme-challenge/
- fname=`awk '{split($0, s, "."); print s[1]}' doc/letsencrypt-setup.html`
- mv doc/letsencrypt-setup.html public/.well-known/acme-challenge/$fname
# output for error tracking
- ls -l public/.well-known/acme-challenge/
- ls -l public/
# check settings on public directory
- chmod go-rwX -R public/
- chmod u-w -R public/
artifacts:
paths:
- public
......
# Base image
FROM ubuntu:17.10
FROM ubuntu:18.04
# Install required packages
RUN apt-get update -qy
RUN apt-get upgrade -qy
RUN apt-get install -qy python3-dev \
python3-pip \
python3-numpy \
python3-scipy \
python3-h5py \
python3-sphinx \
pyflakes \
pep8 \
graphviz
RUN apt update -qy
RUN apt upgrade -qy
RUN apt install -qy tzdata
RUN apt install -qy python3-dev \
python3-pip \
python3-numpy \
python3-scipy \
python3-h5py \
python3-sphinx \
pyflakes \
pep8 \
graphviz
# Set up some Python3 packages via pip
RUN pip3 install --upgrade pip
......
include build_config.py
include dynasor.1
include MANIFEST.in
graft examples
\ No newline at end of file
......@@ -11,7 +11,8 @@
# Example: Explicitly use gcc
#
local_compiler = 'gcc'
extra_compile_args = ['-fPIC', '-fopenmp', '-Ofast', '-march=native', '-std=c99']
extra_compile_args = ['-fPIC', '-fopenmp', '-Ofast', '-march=native',
'-std=c99']
local_linker = local_compiler
local_link_shared = ['-shared']
......@@ -21,7 +22,8 @@ extra_link_args = ['-fopenmp']
# # Example: Use icc instead of the default compiler
# #
# local_compiler = 'icc'
# extra_compile_args = ['-openmp', '-xHOST', '-O3', '-fno-alias', '-fPIC', '-std=c99']
# extra_compile_args = ['-openmp', '-xHOST', '-O3',
# '-fno-alias', '-fPIC', '-std=c99']
# local_linker = local_compiler
# local_link_shared = ['-shared']
......@@ -41,7 +43,8 @@ extra_link_args = ['-fopenmp']
# # Example: Use pgcc to generate GPU-code using OpenACC directives
# #
# local_compiler = 'pgcc'
# extra_compile_args = ['-acc', '-ta=nvidia', '-O4', '-g', '-Minfo=acc', '-fPIC']
# extra_compile_args = ['-acc', '-ta=nvidia', '-O4',
# '-g', '-Minfo=acc', '-fPIC']
# local_linker = local_compiler
# local_link_shared = ['-shared']
......
......@@ -9,10 +9,13 @@ sys.path.insert(0, os.path.abspath('../dsf/'))
extensions = [
'sphinx.ext.autodoc',
'sphinx.ext.todo',
'sphinx_autodoc_typehints',
'sphinx.ext.mathjax',
'sphinx.ext.todo',
'sphinx.ext.napoleon',
'sphinxcontrib.bibtex',
'sphinx_sitemap',
'breathe',
]
site_url = 'https://dynasor.materialsmodeling.org/'
......
......@@ -3,7 +3,31 @@
Python interface
****************
.. toctree::
:maxdepth: 4
This section includes documentation of some of the internal methods and classes
can be useful for further analysis.
dsf
.. index::
single: Function reference; Fourier transformations
single: Function reference; Filon's method
Fourier transformations
-----------------------
.. automodule:: dsf.filon
:members:
:undoc-members:
:inherited-members:
.. index::
single: Function reference; Trajectory handling
Trajectory handling
===================
.. automodule:: dsf.trajectory
:members:
:undoc-members:
:inherited-members:
# -*- coding: utf-8 -*-
'''
"""
dynasor module.
'''
"""
__project__ = 'dynasor'
__description__ = 'A tool for calculating dynamical structure factors'
......@@ -11,7 +11,7 @@ __authors__ = ['Mattias Slabanja',
__copyright__ = '2018'
__license__ = 'GPL2+'
__credits__ = ['The dynasor developers team']
__version__ = '0.2'
__version__ = '0.1'
__maintainer__ = 'The dynasor developers team'
__maintainer_email__ = 'dynasor@materialsmodeling.org'
__status__ = 'beta-version'
......
......@@ -20,38 +20,24 @@ __all__ = ['fixed_bin_averager']
import numpy as np
import logging
logger = logging.getLogger('dynasor')
class fixed_bin_averager:
"""Class for averaging data sets of points (x,y) with
a priori decided positions (x), over a pre-defined number
of equally sized, linearely distriuted bins.
This is used to get an average of y for a set of
lineary spaced x values when a possibly large set
of {y(x)} is given.
y ^ *
| *
| * * * *
| * * **
| * * *
| * * *
| * * * *
+--------------------------------------->
x
| | | | | | | | | |
+---+---+---+---+---+---+---+---+---+
x1 x2 x3 x4 x5 x6 x7 x8 x9
9 linearely distributed "bins" / ranges.
The result will be an array containing,
for each bin, the average y-value over the data points
having its x-value within the respective bin x-range.
"""
Class for averaging data sets of points (x,y) with a priori decided
positions (x), over a pre-defined number of equally sized, linearely
distriuted bins.
This is used to get an average of y for a set of lineary spaced x values
when a possibly large set of {y(x)} is given.
The result will be an array containing, for each bin, the average y-value
over the data points having its x-value within the respective bin x-range.
"""
def __init__(self, x_max, x_bins, x_distances, x_min=0.0):
assert x_max > x_min
assert x_bins > 1
......@@ -63,14 +49,14 @@ class fixed_bin_averager:
range=x_range)
self.x_linspace = 0.5 * (edges[1:]+edges[:-1])
I = np.nonzero(bin_count)
self.bin_count = bin_count[I]
self.x = self.x_linspace[I]
m = np.nonzero(bin_count)
self.bin_count = bin_count[m]
self.x = self.x_linspace[m]
self.input_length = len(x_distances)
self.bins = len(self.x)
if self.bins != x_bins:
logger.info('Ignoring %d bins without coverage' % (
x_bins-self.bins))
logger.info('Ignoring {} bins without coverage'
.format(x_bins-self.bins))
def bin(self, y, axis=0):
y = np.require(y)
......
......@@ -16,53 +16,77 @@
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA.
"""Python implementation of Filon's integration formula.
"""
This module provides an implementation of Filon's integration formula.
For information about Filon's formula, see e.g.
Abramowitz, Stegun, Handbook of Mathematical Functions, section 25,
http://mathworld.wolfram.com/FilonsIntegrationFormula.html,
Allen, Tildesley, Computer Simulation of Liquids, Appendix D.
`Abramowitz and Stegun, Handbook of Mathematical Functions,
section 25 <http://mathworld.wolfram.com/FilonsIntegrationFormula.html>`_ or
Allen and Tildesley, *Computer Simulation of Liquids*, Appendix D.
Integration is performed along one dimension (default ``axis=0``), e.g.,
Integration is performed along one dimension (default axis=0), e.g.
.. code::
[F0[0] F1[0] .. FN[0] ] [f0[0] f1[0] .. fN[0] ]
[ . . . ] [ . . . ]
[F0[.] F1[.] .. FN[.] ] = I([f0[.] f1[.] .. fN[.] ], dx, [k[0] .. k[Nk]])
[ . . . ] [ . . . ]
[F0[Nk] F1[Nk] .. FN[Nk]] [f0[Nx] f1[Nx] .. fN[Nx]]
[F0[0] F1[0] .. FN[0] ] [f0[0] f1[0] .. fN[0] ]
[ . . . ] [ . . . ]
[F0[.] F1[.] .. FN[.] ] = I([f0[.] f1[.] .. fN[.] ], dx, [k[0] .. k[Nk]])
[ . . . ] [ . . . ]
[F0[Nk] F1[Nk] .. FN[Nk]] [f0[Nx] f1[Nx] .. fN[Nx]]
where k and Fj have end index Nk, and fj have end index Nx.
Nk is arbitray (implicitly derived from the length of k).
Due to the algorithm, fj[Nx] must be of odd length (Nx must be an even number),
and should correspond to a linearly spaced set of data points (separated by
dx along the integration axis).
where ``k`` and ``Fj`` have end index ``Nk``, and ``fj`` has end index ``Nx``.
``Nk`` is automatically set by the length of ``k``. Due to the algorithm,
``fj[Nx]`` must be of odd length (``Nx`` must be an even number), and should
correspond to a linearly spaced set of data points (separated by `dx` along the
integration axis).
sin_integral and cos_integral allows for shifted
integration intervals by the optional argument x0.
:func:`sin_integral` and :func:`cos_integral` allow for shifted integration
intervals by the optional argument ``x0``.
"""
__all__ = ['fourier_cos', 'sin_integral', 'cos_integral']
from numpy import sin, cos, pi, mod, \
zeros, arange, linspace, require, \
where, sum
from numpy import (sin, cos, pi, mod, zeros, arange,
linspace, require, where, sum)
def fourier_cos(f, dx, k=None, axis=0):
"""Calculate a direct fourier cosine transform of function f(x) using
Filon's integration method
k, F = fourier_cos(f, dx)
Array values f[0]..f[2n] is expected to correspond to f(0.0)..f(2n*dx),
hence, it should contain an odd number of elements.
The transform is approximated with the integral
2*\int_{0}^{xmax}
where xmax = 2n*dx
If k is not provided, linspace(0.0, 2*pi/dx, f.shape[axis]),
will be used.
"""Calculates the direct Fourier cosine transform :math:`F(k)` of a
function :math:`f(x)` using Filon's integration method.
The array values ``f[0]..f[2n]`` are expected to correspond to
:math:`f(0.0)\ldots f(2n\Delta x)`. Hence, ``f`` should contain an odd
number of elements.
The transform is approximated by the integral
:math:`F(k) = 2\int_{0}^{x_{max}} f(x) \cos(k x) dx`,
where :math:`x_{max} = 2n \Delta x`.
Parameters
----------
f : array
function values; must contain an odd number of elements
dx : float
spacing of x-axis (:math:`\Delta x`)
k : array
values of reciprocal axis, at which to evaluate transform;
if ``k`` is not provided, ``linspace(0.0, 2*pi/dx, f.shape[axis])``,
will be used.
axis : int
axis along which to carry out integration
Returns
-------
(array, array)
tuple of ``k`` and ``F`` values
Example
-------
A common use case is
.. code-block:: python
k, F = fourier_cos(f, dx)
"""
if k is None:
......@@ -72,19 +96,59 @@ def fourier_cos(f, dx, k=None, axis=0):
def cos_integral(f, dx, k, x0=0.0, axis=0):
"""\int_{x0}^{2n*dx} f(x)*cos(k x) dx
f must have length 2n+1 along the integration axis.
"""Calculates the integral
:math:`\int_{x_0}^{2n\Delta x} f(x) \cos(k x) dx`.
Parameters
----------
f : array
function values; must contain an odd number of elements
dx : float
spacing of x-axis (:math:`\Delta x`)
k : array
values of reciprocal axis, at which to evaluate transform;
if ``k`` is not provided, ``linspace(0.0, 2*pi/dx, f.shape[axis])``,
will be used.
x0 : float
offset for integration interval
axis : int
axis along which to carry out integration
Returns
-------
float
tuple of ``k`` and ``F`` values
"""
return _gen_sc_int(f, dx, k, x0, axis, cos)
def sin_integral(f, dx, k, x0=0.0, axis=0):
"""\int_{x0}^{2n*dx} f(x)*sin(k x) dx
f must have length 2n+1 along the integration axis.
def sin_integral(f, dx, k, x0=0.0, axis=0):
"""Calculates the integral
:math:`\int_{x_0}^{2n\Delta x} f(x) \sin(k x) dx`.
Parameters
----------
f : array
function values; must contain an odd number of elements
dx : float
spacing of x-axis (:math:`\Delta x`)
k : array
values of reciprocal axis, at which to evaluate transform;
if ``k`` is not provided, ``linspace(0.0, 2*pi/dx, f.shape[axis])``,
will be used.
x0 : float
offset for integration interval
axis : int
axis along which to carry out integration
Returns
-------
float
tuple of ``k`` and ``F`` values
"""
return _gen_sc_int(f, dx, k, x0, axis, sin)
def _gen_sc_int(f, dx, k, x0, axis, sc):
f = require(f)
......@@ -93,11 +157,10 @@ def _gen_sc_int(f, dx, k, x0, axis, sc):
try:
axis = range(f.ndim)[axis]
except (IndexError, TypeError):
print('Error: axis(=%s) is invalid' % str(axis))
raise
raise Exception('axis(={}) is invalid'.format(str(axis)))
if k.ndim != 1:
raise ValueError('k is not one dimensional')
raise ValueError('k must be one-dimensional')
Nk = len(k)
Nx = f.shape[axis]
......@@ -106,7 +169,7 @@ def _gen_sc_int(f, dx, k, x0, axis, sc):
k_shape[axis] = Nk #
if mod((Nx - 1), 2) != 0 or Nx < 3:
raise ValueError('f must have an odd length, >=3, along its integration axis')
raise ValueError('len f must be odd and >2 along its integration axis')
s = (slice(None),) * f.ndim
odd_index = s + (slice(1, None, 2),)
......@@ -114,11 +177,13 @@ def _gen_sc_int(f, dx, k, x0, axis, sc):
first_index = s + ((0,),)
last_index = s + ((-1,),)
alpha, beta, gamma = [x.reshape(k_shape[:-1]) for x in _alpha_beta_gamma(dx * k)]
alpha, beta, gamma = [x.reshape(k_shape[:-1])
for x in _alpha_beta_gamma(dx * k)]
x = (x0 + dx * arange(0.0, Nx)).reshape(x_shape)
k = k.reshape(k_shape)
# Add an extra dimension to f, and transpose it to put the x-dimension at axis=-1
# Add an extra dimension to f, and transpose it to put the x-dimension at
# axis=-1
t = list(range(f.ndim + 1))
t[axis], t[-1] = t[-1], t[axis]
f = f.reshape(f.shape + (1,)).transpose(t)
......@@ -129,20 +194,19 @@ def _gen_sc_int(f, dx, k, x0, axis, sc):
if sc == sin:
return dx * (alpha * sum((f[first_index] * cos(k * x0) -
f[last_index] * cos(k * x[last_index])), axis=-1) +
beta * sum(f[even_index] * sc_k_x[even_index], axis=-1) +
gamma * sum(f[odd_index] * sc_k_x[odd_index], axis=-1))
f[last_index] * cos(k * x[last_index])),
axis=-1) +
beta * sum(f[even_index] * sc_k_x[even_index], axis=-1) +
gamma * sum(f[odd_index] * sc_k_x[odd_index], axis=-1))
elif sc == cos:
return dx * (alpha * sum((f[last_index] * sin(k * x[last_index]) -
f[first_index] * sin(k * x0)), axis=-1) +
beta * sum(f[even_index] * sc_k_x[even_index], axis=-1) +
gamma * sum(f[odd_index] * sc_k_x[odd_index], axis=-1))
f[first_index] * sin(k * x0)), axis=-1) +
beta * sum(f[even_index] * sc_k_x[even_index], axis=-1) +
gamma * sum(f[odd_index] * sc_k_x[odd_index], axis=-1))
raise RuntimeError('Internal error, this should not happen')
def _alpha_beta_gamma(theta):
# From theta, calculate alpha, beta, and gamma
......@@ -170,4 +234,4 @@ def _alpha_beta_gamma(theta):
beta[I_nz] = 2 * itheta3 * (theta * (1 + cos2_t) - 2 * sin_t * cos_t)
gamma[I_nz] = 4 * itheta3 * (sin_t - theta * cos_t)
return (alpha, beta, gamma)
return alpha, beta, gamma
#
# Taken from the SciPy Cookbook, http://www.scipy.org/Cookbook/Multithreading
#
import sys
import time
import threading
from itertools import izip, count
from itertools import count
def foreach(f,l,threads=3,return_=False):
def foreach(f, l, threads=3, return_=False):
"""
Apply f to each element of l, in parallel
Apply f to each element of l, in parallel.
"""
if threads>1:
if threads > 1:
iteratorlock = threading.Lock()
exceptions = []
if return_:
n = 0
d = {}
i = izip(count(),l.__iter__())
i = zip(count(), l.__iter__())
else:
i = l.__iter__()
def runall():
while True:
iteratorlock.acquire()
......@@ -39,30 +36,30 @@ def foreach(f,l,threads=3,return_=False):
return
try:
if return_:
n,x = v
n, x = v
d[n] = f(x)
else:
f(v)
except:
except Exception:
e = sys.exc_info()
iteratorlock.acquire()
try:
exceptions.append(e)
finally:
iteratorlock.release()
threadlist = [threading.Thread(target=runall) for j in xrange(threads)]
threadlist = [threading.Thread(target=runall) for j in range(threads)]
for t in threadlist:
t.start()
for t in threadlist:
t.join()
if exceptions:
a, b, c = exceptions[0]
raise a, b, c
a = exceptions[0]
raise a
if return_:
r = d.items()
r.sort()
return [v for (n,v) in r]
return [v for (n, v) in r]
else:
if return_:
return [f(v) for v in l]
......@@ -71,6 +68,6 @@ def foreach(f,l,threads=3,return_=False):
f(v)
return
def parallel_map(f,l,threads=3):
return foreach(f,l,threads=threads,return_=True)
def parallel_map(f, l, threads=3):
return foreach(f, l, threads=threads, return_=True)
......@@ -16,12 +16,14 @@
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA.
import re
import numpy as np
import re
from os.path import isfile
section_re = re.compile(r'^ *\[ *([a-zA-Z0-9_.-]+) *\] *$')
class section_index:
"""Read an ini-style gromacs index file
......@@ -40,32 +42,36 @@ class section_index:
m = section_re.match(L)
if m:
if members and name:
sections.append((name, np.unique(np.concatenate(members))-1))
sections.append(
(name, np.unique(np.concatenate(members))-1))
name = m.group(1)
members = []
elif not L.isspace():
members.append(np.fromstring(L, dtype=int, sep=' '))
if members and name:
sections.append((name, np.unique(np.concatenate(members))-1))
sections.append(
(name, np.unique(np.concatenate(members))-1))
else:
sections = [("all", np.arange(max_index, dtype=int))]
self.sections = sections
if not self.valid_index_limits(max_index):
raise ValueError("section_index: Index file seems to contain one or more invalid indices. " + \
("For the provided trajectory file, indices must be in range [1, %d]" % max_index))
raise ValueError('section_index: Index file seems to contain one'
' or more invalid indices. For the provided'
' trajectory file, indices must be in range'
' [1, {}]'.format(max_index))
def valid_index_limits(self, N):
for _,I in self.sections:
for _, I in self.sections:
if I[0] < 0 or I[-1] >= N:
return False
return True
def get_section_names(self):
return [n for n,_ in self.sections]
return [n for n, _ in self.sections]
def get_section_indices(self):
return [i for _,i in self.sections]
return [i for _, i in self.sections]
def N_sections(self):
return len(self.sections)
......@@ -75,11 +81,12 @@ class section_index:
Split x/v into list of xs/vs in accordance with specified sections.
"""
indices = [I for _,I in self.sections]
indices = [I for _, I in self.sections]
def fun(frame):
frame = frame.copy()
frame['xs'] = [frame['x'][:,I] for I in indices]
frame['xs'] = [frame['x'][:, I] for I in indices]
if 'v' in frame:
frame['vs'] = [frame['v'][:,I] for I in indices]
frame['vs'] = [frame['v'][:, I] for I in indices]
return frame
return fun
......@@ -20,17 +20,15 @@ __all__ = ['create_mfile', 'create_pfile']
import numpy as np
import logging
from itertools import repeat
logger = logging.getLogger('dynasor')
def create_mfile(filename, output, comment=None):
"""Create matlab style m-file.
"""
""" Creates matlab style m-file. """
with open(filename, 'w') as fh:
popts = np.get_printoptions()
np.set_printoptions(threshold='inf',
linewidth='inf')
np.set_printoptions(threshold='inf', linewidth='inf')
if comment is not None:
fh.write('%%%\n')
......@@ -41,12 +39,12 @@ def create_mfile(filename, output, comment=None):
fh.write("\n%% %s\n%s = ...\n%s;\n" % (desc, n, str(v)))
np.set_printoptions(**popts)
logger.info('Wrote Matlab-style output to %s' % fh.name)
logger.info('Wrote Matlab-style output to {}'.format(fh.name))
def create_pfile(filename, output, comment=None):
"""Create python pickle file
"""
""" Creates python pickle file. """
import cPickle
with open(filename, 'w') as fh:
cPickle.dump(output, fh)
logger.info('Wrote pickled output to %s' % fh.name)
logger.info('Wrote pickled output to {}'.format(fh.name))
......@@ -19,11 +19,13 @@
__all__ = ['reciprocal_isotropic', 'reciprocal_line']
from os.path import dirname, join
import numpy as np
import logging
from numpy import linalg, array, arange, require, nonzero, pi, sqrt, prod
import numpy as np
from ctypes import cdll, c_int
from distutils.sysconfig import get_config_var
from numpy import linalg, array, arange, require, nonzero, pi, sqrt, prod
from os.path import dirname, join
logger = logging.getLogger('dynasor')
......@@ -38,14 +40,16 @@ _rho_j_k = {}
ndp_f_2d_r = {}
ndp_c_1d_rw = {}
ndp_c_2d_rw = {}
for t in "ds":
ndp_f_2d_r[t] = np_ndp(dtype=np_f[t], ndim=2, flags='f_contiguous, aligned')
for t in 'ds':
ndp_f_2d_r[t] = np_ndp(dtype=np_f[t], ndim=2,
flags='f_contiguous, aligned')
ndp_c_1d_rw[t] = np_ndp(dtype=np_c[t], ndim=1,
flags='f_contiguous, aligned, writeable')
ndp_c_2d_rw[t] = np_ndp(dtype=np_c[t], ndim=2,
flags='f_contiguous, aligned, writeable')