diff --git a/CHANGELOG.md b/CHANGELOG.md
index 03175ce..cab0b98 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,327 +1,328 @@
# Changelog
All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html) for final versions and [PEP440](https://www.python.org/dev/peps/pep-0440/) in case intermediate versions need to be released (e.g. development version `2.2.3.dev1` or release candidates `2.2.3rc1`), or individual commits are packaged.
## Unrelease
### Added
- Added `tamaas.utils.publications` to print a list of publications relevant to
the parts of Tamaas used in a script. Call at the end of a script for an
exhaustive list.
- Added `tamaas.mpi.gather` and `tamaas.mpi.scatter` to help handling of 2D data
in MPI contexts.
- Added `getBoundaryFields()/boundary_fields` function/property to `Model`.
- Added Python constructor to `Model`, which is an alias to
`ModelFactory.createModel`.
- Added convenience `tamaas.dumpers._helper.netCDFtoParaview` function to
convert model sequence to Paraview data file.
- Added dump of surface fields as `FieldData` in VTK files
- Added `scipy.sparse.linalg.LinearOperator` interface to `Westergaard`. A
`LinearOperator` can be created by calling
`scipy.sparse.linalg.aslinearoperator` on the operator object.
- Allowed sub-classing of `ContactSolver` in Python.
- Exposed elastic contact functionals in Python to help with custom solver
implementation.
- Added an example of custom solver with penalty method using Scipy.
- Added `graphArea` function in statistics to compute area of function graph
+- Added `tamaas.utils.radial_average` to radialy average a 2D field
## Changed
- Filter for boundary fields in dumper should be smarter.
- `read` functions in dumpers are now `@classmethod`
- Changed sign of plastic correction in KatoSaturated solver to match the
EPICSolver sign.
- Plastic correction of KatoSaturated is registered in the model as
`KatoSaturated::residual_displacement`.
- Changed to modern Python packaging (PEP517) with `setup.cfg` and
`pyproject.toml`. A `setup.py` script is still provided for editable
installations, see [gh#7953](https://github.com/pypa/pip/issues/7953). Once
setuptools has stable support for
[PEP621](https://peps.python.org/pep-0621/), `setup.cfg` will be merged into
`pyproject.toml`. Note that source distributions that can be built with
[`build`](https://pypi.org/project/build/) are not true source distributions,
since they contain the compiled Python module of Tamaas.
## Fixed
- Fixed `tamaas.dumpers.NetCDFDumper` to dump in MPI context.
## Deprecated
- SCons versions < 3 and Python < 3.6 are explicitely no longer supported.
## v2.4.0 -- 2022-03-22
### Added
- Added a `tamaas.utils` module.
- Added `tamaas.utils.load_path` generator, which yields model objects for a
sequence of applied loads.
- Added `tamaas.utils.seeded_surfaces` generator, which yields surfaces for a
sequence of random seeds.
- Added `tamaas.utils.hertz_surface` which generates a parabolic (Hertzian)
surface.
- Added `tamaas.utils.log_context` context manager which locally sets the log
level for Tamaas' logger.
- Added `tamaas.compute.from_voigt` to convert Voigt/Mendel fields to dense
tensor fields.
- Added a deep-copy function to `ModelFactory` to copy `Model` objects. Use
`model_copy = copy.deepcopy(model)` in Python to make a copy. Currently only
copies registered fields, same as dumpers/readers.
- Added a `read_sequence` method for dumpers to read model frames.
- Automatic draft release on Zenodo.
### Changed
- `*_ROOT` variables for vendored dependencies GoogleTest and pybind11 now
default to empty strings, so that the vendored trees in `third-party/` are not
selected by default. This is so system packages are given priority. Vendored
dependency submodules will eventually be deprecated.
### Fixed
- Fixed an issue with `scons -c` when GTest was missing.
## v2.3.1 -- 2021-11-08
### Added
- Now using `clang-format`, `clang-tidy` and `flake8` for linting
- Pre-built Tamaas images can be pulled with `docker pull
registry.gitlab.com/tamaas/tamaas`
- Added a `--version` flag to the `tamaas` command
### Changed
- The root `Dockerfile` now compiles Tamaas, so using Tamaas in Docker is easier.
- The model attributes dumped to Numpy files are now written in a JSON-formatted
string to avoid unsafe loading/unpickling of objects.
- Removed the `build_doc` build option: now the doc targets are automatically
added if the dependencies are met, and built if `scons doc` is called.
- Removed the `use_googletest` build option: if tests are built and gtest is
present, the corresponding tests will be built
### Fixed
- The command `tamaas plot` gracefully fails if Matplotlib is missing
- Better heuristic to guess which fields are defined on boundary in dumpers
(still not perfect and may give false negatives)
## v2.3.0 -- 2021-06-15
### Added
- Added `read()` method to dumpers to create a model from a dump file
- `getClusters()` can be called in MPI contact with partial contact maps
- Added a JSON encoder class for models and a JSON dumper
- CUDA compatibility is re-established, but has not been tested
- Docstrings in the Python bindings for many classes/methods
### Changed
- Tamaas version numbers are now managed by
[versioneer](https://github.com/python-versioneer/python-versioneer). This
means that Git tags prefixed with `v` (e.g. `v2.2.3`) carry meaning and
determine the version. When no tag is set, versioneer uses the last tag,
specifies the commit short hash and the distance to the last tag (e.g.
`2.2.2+33.ge314b0e`). This version string is used in the compiled library, the
`setup.py` script and the `__version__` variable in the python module.
- Tamaas migrated to [GitLab](https://gitlab.com/tamaas/tamaas)
- Continuous delivery has been implemented:
- the `master` branch will now automatically build and publish Python wheels
to `https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple`. These
"nightly" builds can be installed with:
pip install \
--extra-index-url https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple \
tamaas
- version tags pushed to `master` will automatically publish the wheels to
[PyPI](https://pypi.org/project/tamaas/)
### Deprecated
- The `finalize()` function is now deprecated, since it is automatically called
when the process terminates
- Python versions 3.5 and below are not supported anymore
### Fixed
- Fixed a host of dump read/write issues when model type was not `volume_*d`.
Dumper tests are now streamlined and systematic.
- Fixed a bug where `Model::solveDirichlet` would not compute correctly
- Fixed a bug where `Statistics::contact` would not normalize by the global
number of surface points
## v2.2.2 -- 2021-04-02
### Added
- Entry-point `tamaas` defines a grouped CLI for `examples/pipe_tools`. Try
executing `tamaas surface -h` from the command-line!
### Changed
- `CXXFLAGS` are now passed to the linker
- Added this changelog
- Using absolute paths for environmental variables when running `scons test`
- Reorganized documentation layout
- Gave the build system a facelift (docs are now generated directly with SCons
instead of a Makefile)
### Deprecated
- Python 2 support is discontinued. Version `v2.2.1` is the last PyPi build with
a Python 2 wheel.
- The scripts in `examples/pipe_tools` have been replaced by the `tamaas` command
### Fixed
- `UVWDumper` no longer imports `mpi4py` in sequential
- Compiling with different Thrust/FFTW backends
## v2.2.1 -- 2021-03-02
### Added
- Output registered fields and dumpers in `print(model)`
- Added `operator[]` to the C++ model class (for fields)
- Added `traction` and `displacement` properties to Python model bindings
- Added `operators` property to Python model bindings, which provides a
dict-like access to registered operators
- Added `shape` and `spectrum` to properties to Python surface generator
bindings
- Surface generator constructor accepts surface global shape as argument
- Choice of FFTW thread model
### Changed
- Tests use `/tmp` for temporary files
- Updated dependency versions (Thrust, Pybind11)
### Deprecated
- Most `get___()` and `set___()` in Python bindings have been deprecated. They
will generate a `DeprecationWarning`.
### Removed
- All legacy code
## v2.2.0 -- 2020-12-31
### Added
- More accurate function for computation of contact area
- Function to compute deviatoric of tensor fields
- MPI implementation
- Convenience `hdf5toVTK` function
- Readonly properties `shape`, `global_shape`, `boundary_shape` on model to give
shape information
### Changed
- Preprocessor defined macros are prefixed with `TAMAAS_`
- Moved `tamaas.to_voigt` to `tamaas.compute.to_voigt`
### Fixed
- Warning about deprecated constructors with recent GCC versions
- Wrong computation of grid strides
- Wrong computation of grid sizes in views
## v2.1.4 -- 2020-08-07
### Added
- Possibility to generate a static `libTamaas`
- C++ implementation of DFSANE solver
- Allowing compilation without OpenMP
### Changed
- NetCDF dumper writes frames to a single file
### Fixed
- Compatibility with SCons+Python 3
## v2.1.3 -- 2020-07-27
### Added
- Version number to `TamaasInfo`
### Changed
- Prepending root directory when generating archive
## v2.1.2 -- 2020-07-24
This release changes some core internals related to discrete Fourier transforms
for future MPI support.
### Added
- Caching `CXXFLAGS` in SCons build
- SCons shortcut to create code archive
- Test of the elastic-plastic contact solver
- Paraview data dumper (`.pvd` files)
- Compression for UVW dumper
- `__contains__` and `__iter__` Python bindings of model
- Warning message of possible overflow in Kelvin
### Changed
- Simplified `tamaas_info.cpp`, particularly the diff part
- Using a new class `FFTEngine` to manage discrete Fourier transforms. Plans are
re-used as much as possible with different data with the same shape. This is
in view of future MPI developments
- Redirecting I/O streams in solve functions so they can be used from Python
(e.g. in Jupyter notebooks)
- Calling `initialize()` and `finalize()` is no longer necessary
### Fixed
- Convergence issue with non-linear solvers
- Memory error in volume potentials
## v2.1.1 -- 2020-04-22
### Added
- SCons shortcut to run tests
### Fixed
- Correct `RPATH` for shared libraries
- Issues with SCons commands introduced in v2.1.0
- Tests with Python 2.7
## v2.1.0 -- 2020-04-17
### Added
- SCons shortcuts to build/install Tamaas and its components
- Selection of integration method for Kelvin operator
- Compilation option to remove the legacy part of Tamaas
- NetCDF dumper
### Fixed
- Link bug with clang
- NaNs in Kato saturated solver
## v2.0.0 -- 2019-11-11
First public release. Contains relatively mature elastic-plastic contact code.
diff --git a/python/tamaas/utils.py b/python/tamaas/utils.py
index b2d6e01..ccbee20 100644
--- a/python/tamaas/utils.py
+++ b/python/tamaas/utils.py
@@ -1,295 +1,310 @@
# -*- mode: python; coding: utf-8 -*-
# vim: set ft=python:
#
# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
"""Convenience utilities."""
import inspect
from collections import namedtuple
from contextlib import contextmanager
from typing import Iterable, Union
from functools import partial
from operator import contains
import numpy as np
from . import (
ContactSolver,
Model,
SurfaceGenerator1D,
SurfaceGenerator2D,
dtype,
set_log_level,
get_log_level,
LogLevel,
Logger,
)
__all__ = [
"log_context",
"publications",
"load_path",
"seeded_surfaces",
"hertz_surface",
]
class NoConvergenceError(RuntimeError):
"""Convergence not reached exception."""
@contextmanager
def log_context(log_level: LogLevel):
"""Context manager to easily control Tamaas' logging level."""
current = get_log_level()
set_log_level(log_level)
try:
yield
finally:
set_log_level(current)
def publications(format_str="{pub.citation}\n\t{pub.doi}"):
"""Print publications associated with objects in use."""
Publication = namedtuple("Publication", ["citation", "doi"])
joss = Publication(
citation=(
"Frérot, L., Anciaux, G., Rey, V., Pham-Ba, S. & Molinari, J.-F."
" Tamaas: a library for elastic-plastic contact of periodic rough"
" surfaces. Journal of Open Source Software 5, 2121 (2020)."),
doi="10.21105/joss.02121",
)
zenodo = Publication(
citation=(
"Frérot, Lucas, Anciaux, Guillaume, Rey, Valentine, Pham-Ba, Son, "
"& Molinari, Jean-François. (2022). Tamaas, a high-performance "
"library for periodic rough surface contact (2.4.0). Zenodo."),
doi="10.5281/zenodo.3479236",
)
_publications = {
k: v
for keys, v in [
(
(
"SurfaceGeneratorRandomPhase1D",
"SurfaceGeneratorRandomPhase2D",
),
[
Publication(
citation=(
"Wu, J.-J. Simulation of rough surfaces with FFT."
" Tribology International 33, 47–58 (2000)."),
doi="10.1016/S0301-679X(00)00016-5",
),
],
),
(
("SurfaceGeneratorFilter1D", "SurfaceGeneratorFilter2D"),
[
Publication(
citation=(
"Hu, Y. Z. & Tonder, K. Simulation of 3-D random"
" rough surface by 2-D digital filter and fourier"
" analysis. International Journal of Machine Tools"
" and Manufacture 32, 83–90 (1992)."),
doi="10.1016/0890-6955(92)90064-N",
),
],
),
(
("PolonskyKeerRey", ),
[
Publication(
citation=(
"Polonsky, I. A. & Keer, L. M. A numerical method"
" for solving rough contact problems based on the"
" multi-level multi-summation and conjugate"
" gradient techniques. Wear 231, 206–219 (1999)."),
doi="10.1016/S0043-1648(99)00113-1",
),
Publication(
citation=(
"Rey, V., Anciaux, G. & Molinari, J.-F. Normal"
" adhesive contact on rough surfaces: efficient"
" algorithm for FFT-based BEM resolution. Comput"
" Mech 1–13 (2017)."),
doi="10.1007/s00466-017-1392-5",
),
],
),
(
("DFSANESolver", "DFSANECXXSolver"),
[
Publication(
citation=(
"La Cruz, W., Martínez, J. & Raydan, M. Spectral"
" residual method without gradient information for"
" solving large-scale nonlinear systems of"
" equations. Math. Comp. 75, 1429–1448 (2006)."),
doi="10.1090/S0025-5718-06-01840-0",
),
],
),
(
("Residual", ),
[
Publication(
citation=("Frérot, L., Bonnet, M., Molinari, J.-F. &"
" Anciaux, G. A Fourier-accelerated volume"
" integral method for elastoplastic contact."
" Computer Methods in Applied Mechanics and"
" Engineering 351, 951–976 (2019)."),
doi="10.1016/j.cma.2019.04.006",
),
],
),
(
("EPICSolver", ),
[
Publication(
citation=("Frérot, L., Bonnet, M., Molinari, J.-F. &"
" Anciaux, G. A Fourier-accelerated volume"
" integral method for elastoplastic contact."
" Computer Methods in Applied Mechanics and"
" Engineering 351, 951–976 (2019)."),
doi="10.1016/j.cma.2019.04.006",
),
Publication(
citation=(
"Jacq, C., Nélias, D., Lormand, G. & Girodin, D."
" Development of a Three-Dimensional"
" Semi-Analytical Elastic-Plastic Contact Code."
" Journal of Tribology 124, 653 (2002)."),
doi="10.1115/1.1467920",
),
],
),
(
("Condat", ),
[
Publication(
citation=(
"Condat, L. A Primal–Dual Splitting Method for"
" Convex Optimization Involving Lipschitzian,"
" Proximable and Linear Composite Terms. J Optim"
" Theory Appl 158, 460–479 (2012)."),
doi="10.1007/s10957-012-0245-9",
),
],
),
(
("KatoSaturated", ),
[
Publication(
citation=(
"Stanley, H. M. & Kato, T. An FFT-Based Method for"
" Rough Surface Contact. J. Tribol 119, 481–485"
" (1997)."),
doi="10.1115/1.2833523",
),
Publication(
citation=(
"Almqvist, A., Sahlin, F., Larsson, R. &"
" Glavatskih, S. On the dry elasto-plastic contact"
" of nominally flat surfaces. Tribology"
" International 40, 574–579 (2007)."),
doi="10.1016/j.triboint.2005.11.008",
),
],
),
] for k in keys
}
frame = inspect.stack()[1][0]
caller_vars = (frame.f_globals | frame.f_locals)
citable = filter(
partial(contains, _publications.keys()),
map(lambda x: type(x).__name__, caller_vars.values()),
)
citable = [joss, zenodo] \
+ list({pub for k in citable for pub in _publications[k]})
msg = "Please cite the following publications:\n\n"
msg += "\n".join(format_str.format(pub=pub) for pub in citable)
Logger().get(LogLevel.info) << msg
return citable
def load_path(
solver: ContactSolver,
loads: Iterable[Union[float, np.ndarray]],
verbose: bool = False,
callback=None,
) -> Iterable[Model]:
"""
Generate model objects solutions for a sequence of applied loads.
:param solver: a contact solver object
:param loads: an iterable sequence of loads
:param verbose: print info output of solver
:param callback: a callback executed after the yield
"""
log_level = LogLevel.info if verbose else LogLevel.warning
with log_context(log_level):
for load in loads:
if solver.solve(load) > solver.tolerance:
raise NoConvergenceError("Solver error exceeded tolerance")
yield solver.model
if callback is not None:
callback()
def seeded_surfaces(
generator: Union[SurfaceGenerator1D, SurfaceGenerator2D],
seeds: Iterable[int],
) -> Iterable[np.ndarray]:
"""
Generate rough surfaces with a prescribed seed sequence.
:param generator: surface generator object
:param seeds: random seed sequence
"""
for seed in seeds:
generator.random_seed = seed
yield generator.buildSurface()
def hertz_surface(system_size: Iterable[float], shape: Iterable[int],
radius: float) -> np.ndarray:
"""
Construct a parabolic surface.
:param system_size: size of the domain in each direction
:param shape: number of points in each direction
:param radius: radius of surface
"""
coords = map(
lambda L, N: np.linspace(0, L, N, endpoint=False, dtype=dtype),
system_size,
shape,
)
coords = np.meshgrid(*coords, "ij", sparse=True)
surface = (-1 / (2 * radius)) * sum(
(x - L / 2)**2 for x, L in zip(coords, system_size))
return surface
+
+
+def radial_average(x, y, values, r, theta, method='linear', endpoint=False):
+ """Compute the radial average of a 2D field."""
+ try:
+ from scipy.interpolate import RegularGridInterpolator
+ except ImportError:
+ raise ImportError("Install scipy to use tamaas.utils.radial_average")
+
+ interpolator = RegularGridInterpolator((x, y), values, method=method)
+ rr, tt = np.meshgrid(r, theta, indexing='ij', sparse=True)
+ x, y = rr * np.cos(tt), rr * np.sin(tt)
+ X = np.vstack((x.flatten(), y.flatten())).T
+ return interpolator(X).reshape(x.shape).sum(axis=1) \
+ / (theta.size - int(not endpoint))
diff --git a/site_scons/version.py b/site_scons/version.py
index aaaf989..94bbd88 100644
--- a/site_scons/version.py
+++ b/site_scons/version.py
@@ -1,51 +1,52 @@
# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
import subprocess
import base64
from zlib import compress
def git(*args):
"""Run a git command"""
return bytes(subprocess.check_output(["git"] + list(args),
stderr=subprocess.STDOUT))
+
def get_git_subst():
"""Get info about state of git repository"""
try:
branch = git("rev-parse", "--abbrev-ref", "HEAD")[:-1]
commit = git("rev-parse", branch)[:-1]
diff = git("diff", "HEAD")
remotes = git('remote', '-v')
except subprocess.CalledProcessError as err:
if b'not a git repository' in err.output:
print("Not in a git repository: skipping info gathering")
branch, commit, diff, remotes = [bytes()] * 4
if remotes != b"":
remotes_string = remotes[:-1].decode().replace('\n', '\\\\n')
else:
remotes_string = '""'
return {
'@commit@': commit.decode(),
'@branch@': branch.decode(),
'@diff@': base64.b64encode(compress(diff, 9)).decode(),
'@remotes@': remotes_string,
}
diff --git a/tests/test_surface.py b/tests/test_surface.py
index 368ad91..10feaf3 100644
--- a/tests/test_surface.py
+++ b/tests/test_surface.py
@@ -1,146 +1,140 @@
# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
from __future__ import print_function, division
+import pytest
import numpy as np
import tamaas as tm
-import pytest
-import scipy.integrate as integrate
-from scipy.special import jv, ellipe as E
+
from numpy import pi
-from scipy.interpolate import interp2d
+
+integrate = pytest.importorskip('scipy.integrate')
+special = pytest.importorskip('scipy.special')
+
+jv = special.jv
+E = special.ellipe
def spectrum(r, isopowerlaw):
kl, kr, ks = isopowerlaw.q0, isopowerlaw.q1, isopowerlaw.q2
H = isopowerlaw.hurst
C = 1.
if r > ks:
return 0
elif r > kr:
return C * (r / float(kr))**(-2. * (H + 1))
elif r > kl:
return C
else:
return 0
def integrate_bessel(x, spectrum, n):
return 2 * pi * integrate.quad(
lambda y: y * jv(0, y * x) * spectrum(n * y), 0, np.inf, limit=200)[0]
-def radial_average(field):
- x = np.linspace(-1, 1, field.shape[0])
- y = np.linspace(-1, 1, field.shape[1])
- interpolation = interp2d(x, y, field)
-
- n = field.shape[0] // 2
- radii = np.linspace(0, 1, n)
- thetas = np.linspace(0, 2 * pi, endpoint=False)
-
- sum = np.zeros_like(radii)
-
- for Θ in thetas:
- for i, r in enumerate(radii):
- sum[i] += interpolation(r * np.cos(Θ), r * np.sin(Θ))
- return sum / thetas.size
-
-
integrate_bessel = np.vectorize(integrate_bessel, otypes=[tm.dtype])
@pytest.fixture(
scope="module",
params=[tm.SurfaceGeneratorFilter2D, tm.SurfaceGeneratorRandomPhase2D])
def stats(request):
iso = tm.Isopowerlaw2D()
iso.q0 = 16
iso.q1 = 32
iso.q2 = 64
iso.hurst = 0.8
N = 243
sg = request.param()
sg.spectrum = iso
sg.shape = [N, N]
analytical_stats = {
"rms_heights": iso.rmsHeights(),
"rms_slopes_spectral": iso.rmsSlopes(),
"moments": iso.moments(),
"alpha": iso.alpha()
}
stats = {k: [] for k in analytical_stats}
acf = np.zeros(sg.shape, dtype=tm.dtype)
realizations = 1000
for i in range(realizations):
sg.random_seed = i
s = sg.buildSurface()
stats['rms_heights'].append(tm.Statistics2D.computeRMSHeights(s))
stats['rms_slopes_spectral'].append(
tm.Statistics2D.computeSpectralRMSSlope(s))
stats['moments'].append(tm.Statistics2D.computeMoments(s))
acf += tm.Statistics2D.computeAutocorrelation(s) / realizations
L = 1 * sg.spectrum.q0
n = sg.shape[0] // 2
x = np.linspace(0, L / 2, n)
true_acf = (L / (2 * pi))**2 \
* integrate_bessel(x, lambda x: spectrum(x, sg.spectrum), L / (2 * pi))
return analytical_stats, stats, acf, true_acf
@pytest.mark.parametrize('key', ['rms_heights', 'rms_slopes_spectral'])
def test_statistics(stats, key):
analytical_stats, stats, _, _ = stats
mu = np.mean(stats[key])
error = np.abs(mu - analytical_stats[key]) / analytical_stats[key]
assert error < 3e-3
def test_autocorrelation(stats):
+ from tamaas.utils import radial_average
+
_, _, acf, true_acf = stats
- acf = radial_average(np.fft.fftshift(acf))
- acf_error = np.sum((acf - true_acf)**2) / np.sum(true_acf**2)
+ x, y = (np.linspace(-1, 1, acf.shape[i]) for i in range(2))
+ r = np.linspace(0, 1, true_acf.shape[0])
+ theta = np.linspace(0, 2 * np.pi, endpoint=False)
+
+ r_acf = radial_average(x, y, np.fft.fftshift(acf),
+ r, theta, endpoint=False)
+ acf_error = np.sum((r_acf - true_acf)**2) / np.sum(true_acf**2)
assert acf_error < 2e-3
def test_graph_area():
N = 100
x = np.linspace(0, 1, N, endpoint=False)
f = np.sin(2 * np.pi * x)
df = 2 * np.pi * np.cos(2 * np.pi * x)
num_area = np.sqrt(1 + df**2).mean()
area = 2 * np.sqrt(1 + 4 * np.pi**2) / np.pi * E(4 * np.pi**2 /
(1 + 4 * np.pi**2))
fourier_area = tm.Statistics1D.graphArea(f)
assert np.abs(num_area - fourier_area) / num_area < 1e-10
assert np.abs(area - fourier_area) / area < 2e-10