Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88685176
test_surface.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Oct 20, 04:07
Size
4 KB
Mime Type
text/x-python
Expires
Tue, Oct 22, 04:07 (2 d)
Engine
blob
Format
Raw Data
Handle
21806658
Attached To
rTAMAAS tamaas
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
numpy
as
np
import
tamaas
as
tm
def
compute_mse
(
values
,
mean
):
return
np
.
mean
((
values
-
mean
)
**
2
)
/
values
.
size
def
test_surface
():
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
.
q0
=
k_l
iso
.
q1
=
k_r
iso
.
q2
=
k_s
iso
.
hurst
=
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
.
random_seed
=
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
):
assert
mse
<
tol
tm
.
finalize
()
if
__name__
==
"__main__"
:
test_surface
()
Event Timeline
Log In to Comment