Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_surface.py
2805b1e8fc1ffeatures/netcdf-mpi

Last Change
Raw File

test_surface.py
View Options

#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @author Lucas Frérot <lucas.frerot@epfl.ch>
#
# @section LICENSE
#
# Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
# Solides)
#
# Tamaas is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Tamaas 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
from __future__ import print_function, division
import sys
import numpy as np
import tamaas as tm
def compute_mse(values, mean):
return np.mean((values - mean)**2) / values.size
def main():
tm.initialize(1)
# Spectrum parameters in wavelengths
disctretization_size = 1.
lambda_s = 8.
lambda_r = 64.
lambda_l = 64.
box_size = 512
# Independent parameters
Hurst = 0.8
hrms = 1.
# Spectrum parameters in wavenumbers
k_l = int(box_size // lambda_l)
k_r = int(box_size // lambda_r)
k_s = int(box_size // lambda_s)
N = int(box_size // disctretization_size)
# Output setup
realizations = 100
hrms_out = np.zeros(realizations)
alpha_out = np.zeros(realizations)
rms_slopes_out = {"spectral": np.zeros(realizations),
"geometric": np.zeros(realizations)}
iso = tm.Isopowerlaw2D()
iso.getQ0().assign(k_l)
iso.getQ1().assign(k_r)
iso.getQ2().assign(k_s)
iso.getHurst().assign(Hurst)
sg = tm.SurfaceGeneratorFilter2D()
sg.setFilter(iso)
sg.setSizes([N, N])
# Analytical moments
m = iso.moments()
alpha = iso.alpha()
rms_slopes = iso.rmsSlopes()
for i in range(realizations):
# Surface generator setup
sg.getRandomSeed().assign(i)
# Generating surface
surface = sg.buildSurface() / iso.rmsHeights()
# Computing RMS heights
# We know surface.mean() = 0, so we get an unbiased estimator here
hrms_out[i] = np.sqrt(1 / (N**2-1) * np.sum(surface**2))
# Computing spectrum bandwidth
surface /= box_size**2 # Normalize so that the computed PSD does not contain the box size
moments = tm.SurfaceStatistics.computeMoments(surface)
m0 = hrms_out[i]**2 / box_size**4
m2 = moments[0]
m4 = moments[1]
alpha_out[i] = m0*m4/m2**2
# Computing RMS slopes
surface *= N * iso.rmsHeights() * 2 # TODO: Must find out why it is not exact (factor 2?)
rms_slopes_out["spectral"][i] = tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
rms_slopes_out["geometric"][i] = tm.SurfaceStatistics.computeRMSSlope(surface)
print(rms_slopes_out['spectral'][i])
# print(rms_slopes_out)
# print(rms_slopes)
# Computing Mean Square Error
print("""========= Analytical values =========
hrms = {}
alpha = {}
m0, m2, m4 = {}
rms_slopes = {}""".format(hrms, alpha, m, rms_slopes))
print("========= MSE Calculations =========")
data_association = [
(hrms_out, hrms, "hrms"),
(alpha_out, alpha, "alpha"),
(rms_slopes_out["spectral"], rms_slopes, "RMS Slopes Spectral"),
(rms_slopes_out["geometric"], rms_slopes, "RMS Slopes Geometric")
]
# Computing and printing MSE values
mse_vals = []
for out, analytical, string in data_association:
mse = np.sqrt(compute_mse(out, analytical)) / analytical
print("Normalized error (MSE) for {}: {}".format(string, mse))
mse_vals.append(mse)
# Checking tolerances
tols = [5e-3, 5e-3, 1e-1, 1e-1]
for mse, tol in zip(mse_vals, tols):
if mse >= tol:
print("{} >= {}".format(mse, tol))
return 1
tm.finalize()
return 0
if __name__ == "__main__":
sys.exit(main())
c4science · Help