Homec4science
Diffusion tamaas (ci-fix-static-checks) Browse

tamaas/tests/test_autocorrelation.py
6de98daa963eci-fix-static-checks

Last Change
Raw File

test_autocorrelation.py
View Options

# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-19 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 <https://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()
c4science · Help