Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78564429
test_autocorrelation.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
Wed, Aug 21, 15:53
Size
3 KB
Mime Type
text/x-python
Expires
Fri, Aug 23, 15:53 (2 d)
Engine
blob
Format
Raw Data
Handle
19953638
Attached To
rTAMAAS tamaas
test_autocorrelation.py
View Options
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# @author Valentine Rey <valentine.rey@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
import
tamaas
as
tm
import
numpy
as
np
from
scipy.special
import
jv
import
scipy.integrate
as
integrate
def
spectrum
(
r
,
kl
,
kr
,
ks
,
C
,
H
):
if
r
>
ks
:
return
0
elif
r
>
kr
:
return
C
*
(
r
/
float
(
kr
))
**
(
-
2.
*
(
H
+
1
))
elif
r
>
kl
:
return
C
else
:
return
0
# n is normalization parameter
def
integrate_bessel
(
x
,
spectrum
,
n
):
result
=
2
*
np
.
pi
*
integrate
.
quad
(
lambda
y
:
y
*
jv
(
0
,
y
*
x
)
*
spectrum
(
n
*
y
),
0
,
np
.
inf
)[
0
]
return
result
def
test_autocorrelation
():
tm
.
initialize
()
grid_size
=
256
lambda_l
=
8.
lambda_r
=
8.
lambda_s
=
1.
L
=
4
*
lambda_l
kl
=
int
(
L
/
lambda_l
)
kr
=
int
(
L
/
lambda_r
)
ks
=
int
(
L
/
lambda_s
)
# Computing ACF from integral with Bessel function
l
=
grid_size
//
2
# points on curve
x
=
np
.
linspace
(
0
,
L
/
2
,
l
)
bessel_acf
=
np
.
zeros
(
l
)
integrate_bessel_v
=
np
.
vectorize
(
integrate_bessel
,
otypes
=
[
np
.
float
])
bessel_acf
=
(
L
/
(
2
*
np
.
pi
)
)
**
2
*
integrate_bessel_v
(
x
,
lambda
x
:
spectrum
(
x
,
kl
,
kr
,
ks
,
1.
,
0.8
),
L
/
(
2
*
np
.
pi
))
nsamples
=
1000
etages
=
np
.
zeros
((
grid_size
,
grid_size
,
nsamples
))
# Computing ACF with Hu & Tonder algo reimplemented
for
i
in
range
(
nsamples
):
SG
=
tm
.
SurfaceGeneratorFilterFFT
()
# generate surface
SG
.
getGridSize
()
.
assign
(
grid_size
)
SG
.
getHurst
()
.
assign
(
0.8
)
SG
.
getQ0
()
.
assign
(
kl
);
SG
.
getQ1
()
.
assign
(
kr
);
SG
.
getQ2
()
.
assign
(
ks
);
SG
.
getRandomSeed
()
.
assign
(
i
);
SG
.
Init
()
a
=
SG
.
buildSurface
()
etages
[:,:,
i
]
=
a
[:,
:]
cov
=
np
.
zeros
(
l
)
for
i
in
np
.
arange
(
l
,
grid_size
,
1
):
cov
[
i
-
l
]
=
np
.
sum
((
etages
[
l
,
i
,:]
-
np
.
mean
(
etages
[
l
,
i
,:]))
*
(
etages
[
l
,
l
,:]
-
np
.
mean
(
etages
[
l
,
l
,:])))
/
nsamples
error
=
np
.
sum
((
cov
-
bessel_acf
)
*
(
cov
-
bessel_acf
))
/
np
.
sum
(
bessel_acf
*
bessel_acf
)
print
(
error
)
assert
error
<
1e-1
if
__name__
==
"__main__"
:
test_autocorrelation
()
Event Timeline
Log In to Comment