Commit f19273bb authored by Pierre Navaro's avatar Pierre Navaro
Browse files

mv helper functions from src to test

parent 0c9dbd4c
......@@ -7,6 +7,7 @@ version = "0.1.0"
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
......
%% Cell type:code id: tags:
``` julia
using Plots
pyplot()
```
%% Output
Plots.PyPlotBackend()
%% Cell type:markdown id: tags:
HermiteGf is a julia package, install it with
```julia
using Pkg
Pkg.add("https://gitlab.mpcdf.mpg.de/clapp/HermiteGF.jl")
```
%% Cell type:code id: tags:
``` julia
using HermiteGF
```
%% Output
┌ Info: Recompiling stale cache file /Users/navaro/.julia/compiled/v1.0/HermiteGF/V9Jve.ji for HermiteGF [6d312152-d853-11e8-033d-978d5c49d094]
└ @ Base loading.jl:1187
%% Cell type:code id: tags:
``` julia
import HermiteGF:trapz
```
%% Cell type:code id: tags:
``` julia
xmin, xmax, nx = -1, 1, 32
ϵ, γ = 0.1, 3.0
interp = Hermite(Chebyshev(xmin, xmax, nx), 0.1, 3.0)
xk = interp.xk
f = cos.(xk.^2)
xe = collect(range(xmin, stop=xmax, length=nx))
dx = xe[2]-xe[1]
fe = cos.(xe.^2)
s = interpolate( interp, f, xe )
sqrt(trapz((s - fe).^2, dx)), maximum(abs.(s .- fe))
```
%% Output
MethodError: no method matching Hermite(::Chebyshev, ::Float64, ::Int64)
MethodError: no method matching Hermite(::Chebyshev, ::Float64, ::Float64)
Closest candidates are:
Hermite(!Matched::Array{HermiteGF.NodesType,1}, ::Real, ::Real) at /Users/navaro/.julia/packages/HermiteGF/6sCXJ/src/hermite.jl:17
Stacktrace:
[1] top-level scope at In[4]:2
[1] top-level scope at In[4]:3
%% Cell type:code id: tags:
``` julia
xmin, xmax, nx = -1, 1, 8
xe = collect(range(xmin, stop=xmax, length=nx))
fe = cos.(xe.^2)
```
%% Output
8-element Array{Float64,1}:
0.5403023058681398
0.8726448614622705
0.9831793964382022
0.9997917606637399
0.9997917606637399
0.9831793964382022
0.8726448614622705
0.5403023058681398
%% Cell type:code id: tags:
``` julia
hermite_chebyshev = Hermite(Chebyshev(xmin, xmax, nx), 0.1, 3)
xk = hermite_chebyshev.xk
f = cos.(xk.^2)
scatter(xe, interpolate( hermite_chebyshev, f, xe );
label = "chebyshev")
plot!(xe, fe; label="f", title="Hermite-Chebyshev")
```
%% Output
MethodError: no method matching Hermite(::Chebyshev, ::Float64, ::Int64)
Stacktrace:
[1] top-level scope at In[6]:1
%% Cell type:code id: tags:
``` julia
hermite_uniform = Hermite(Uniform( xmin, xmax, nx), 0.1, 3)
xk = hermite_uniform.xk
f = cos.(xk.^2)
scatter!(xe, interpolate( hermite_uniform, f, xe );
label="uniform")
plot!(xe, fe; label="f", title="Hermite-Uniform")
```
%% Output
MethodError: no method matching Hermite(::Uniform, ::Float64, ::Int64)
Stacktrace:
[1] top-level scope at In[7]:1
%% Cell type:code id: tags:
``` julia
```
......
......@@ -4,8 +4,6 @@ module HermiteGF
export Chebyshev, Uniform
export interpolate
include("helper_functions.jl")
include("ndgrid.jl")
include("evaluate_tensor_product.jl")
include("evaluate_hermite.jl")
......@@ -19,11 +17,6 @@ module HermiteGF
include("uniform.jl")
include("hermite.jl")
include("radial.jl")
include("interpolation1d.jl")
include("interpolation2d.jl")
include("interpolation3d.jl")
include("interpolation4d.jl")
include("interpolation5d.jl")
include("interpolation.jl")
end
......@@ -6,18 +6,18 @@ Hermite inteprolation
"""
mutable struct Hermite <: InterpolationType
nx :: Int64
nnodes :: Int64
epsilon :: Float64
gamma :: Float64
xk :: Vector{Array{Float64,1}}
colloc_mat :: AbstractVecOrMat{Float64}
nodes :: Array{Float64,1}
colloc_mat :: Array{Float64,2}
function Hermite( nodes::Vector{NodesType}, epsilon::Real, gamma::Real )
function Hermite( nodes_t::NodesType, epsilon::Float64, gamma::Float64 )
xk = nodes.xk
nx = nodes.nx
colloc_mat = evaluate_hermite(xk, nx, epsilon, gamma)
new( nx, epsilon, gamma, xk, colloc_mat )
nodes = nodes_t.xk
nnodes = nodes_t.nx
colloc_mat = evaluate_hermite(nodes, nnodes, epsilon, gamma)
new( nnodes, epsilon, gamma, nodes, colloc_mat )
end
......@@ -26,14 +26,43 @@ end
"""
Hermite( xe )
returns evaluation matrix
This function returns the matrix He of the values of the HermiteGF basis
functions. The computation is done via three term recurrence for Hermite
functions with an argument gamma*x and then appropriate exponential scaling.
More details can be found in Section 5.1 of the paper
*STABLE EVALUATION OF GAUSSIAN RADIAL BASIS
FUNCTIONS USING HERMITE POLYNOMIALS*
by Anna Yurova and Katharina Kormann.
"""
function (interp::Hermite)( xe::Array{Float64,1} )
function (interp::Hermite)( x::Array{Float64,1} )
nx = interp.nx
nx = interp.nnodes
epsilon = interp.epsilon
gamma = interp.gamma
evaluate_hermite( xe, nx, epsilon, gamma)
eval_mat = zeros(Int(length(x)), nx)
# Write the values of the first two Hermite functions to initialize the
# three-term recurrence
eval_mat[:,1] = exp.(-0.5.*((gamma*x).^2))
eval_mat[:,2] = gamma*x.*sqrt(2).*exp.(-0.5.*((gamma*x).^2))
# Three term recurrence for the Hermite functions with argument gamma*x
for i = 3:nx
eval_mat[:,i] = sqrt(2.0/(i-1)) .* (gamma*x) .* eval_mat[:,i-1] .- sqrt((i-2)/(i-1)) .* eval_mat[:,i-2]
end
# Scaling the Hermite functions with the exponential factor
# exp(-epsilon^2 x^2 + (gamma x)^2/2)
for i=1:nx
eval_mat[:,i] = (pi^(1/4)) * eval_mat[:,i] .* exp.((x.^2) .* (gamma*gamma*0.5 - epsilon^2))
end
eval_mat / interp.colloc_mat
end
function interpolate(interp::InterpolationType,
f::Array{Float64,1}, xe::Array{Float64,1} )
"""
ϵ = interp.epsilon
nx = interp.nx
xk = interp.xk
interpolate( interp :: InterpolationType,
f :: Array{Float64,1},
xe :: Array{Float64,1})
# Precompute 1D interpolation matrix
x_all = interp(xe) / interp.colloc_mat
Computes interpolation of a 1D function via Hermite-tensor method.
# Compute the interpolated values
ne = Int(length(xe))
Arguments:
- interp : inteprolation type (Hermite or Radial)
- f : 2d array containing values at interpolation nodes (Chebyshev or Uniform)
- xe : vector of evaluation points
# Initialize the result tensor
s = similar(xe)
© Anna Yurova, 2017
"""
function interpolate(interp :: InterpolationType,
f :: Array{Float64,1},
xe :: Array{Float64,1} )
x = interp(xe)
@tensor begin
s[e] = s[e] + x_all[e, c]*f[c]
s[e] := x[e, c]*f[c]
end
s
end
function interpolate(interp_x :: InterpolationType,
interp_y :: InterpolationType,
f :: Array{Float64,2},
xe :: Array{Float64,1},
ye :: Array{Float64,1})
x = interp_x(xe)
y = interp_y(ye)
@tensor begin
s[e1,e2] := x[e1,c1]*y[e2,c2]*f[c1,c2]
end
s
......
"""
interpolate_1D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
Computes an interpolant of a 1D function via Hermite-tensor method.
Arguments:
- function_name : the name of the tested function. Options_available: "f_3"
- nodes_types : an array of two strings with the spacing of nodes in corresponding directions.
Options available: :Uniform, :Chebyshev
- epsilon : a vector of the values of shape parameters in corresponding directions.
The values cannot be larger than 1!
- N : a vector with the numbers of basis functions in corresponding directions.
- Ne : a vector with the numbers of evaluation points in corresponding directions.
- interpolation_type - a string with the name of preferred basis.
Options available: :Hermite, :Radial
gamma - a factor allowing to scale the domain and the function accordingly.
© Anna Yurova, 2017
"""
function interpolate_1D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
# Set the nodes types in all directions.
nodes_type_x = nodes_types[1]
## ----------------------- TEST CASES -------------------------------------
if function_name == :f_3
# Setup the interpolated function.
f(x) = cos.(x.^2)
# Setup the bounds of the interpolation domain.
xmin = -1
xmax = 1
else
@error ("Wrong name of the test case! Possible options: :f_3")
end
## -------------- COLLOCATION AND EVALUATION MATRICES ---------------------
# Setup the basis used for the interpolation.
# Declare the functions for the generating of the collocation and evaluation
# matrices according to the chosen basis.
if interpolation_type == :Hermite
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
elseif interpolation_type == :Radial
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
else
@error ("Wrong basis type! Possible options are: Radial, Hermite.")
end
## -------------------- NODES TYPES ---------------------------------------
# Declare the nodes generating functions in all directions according to the specified by the user.
# X - direction
if nodes_type_x == :Chebyshev
generate_nodes_x = (Nx) -> generate_chebyshev_nodes(Nx, xmin, xmax)
elseif nodes_type_x == :Uniform
generate_nodes_x = (Nx) -> generate_uniform_nodes(Nx, xmin, xmax)
else
@error ("Wrong node type! Possible options: Chebyshev, Uniform")
end
# ---------------- INTERPOLATION PARAMETERS ------------------------
# Setup the value of epsilon
ep1 = epsilon[1]
# Setup the number of collocation points
Nx = N[1]
#Setup the number of evaluation points
Ne1 = Ne[1]
# Generate node points
xk = generate_nodes_x(Nx)
min_dist = minimum(abs.(xk[2:end] .- xk[1:end-1]))
xx_k = xk
# Generate evaluation points
xe = range(xmin, stop=xmax, length=Ne1)
xx_e = xe
## ------------------- ACTUAL INTERPOLATION -------------------------------
# Calculate collocation and evaluation matrices
X_collocation = generate_collocation_matrix_x(xk, Nx, ep1)
X_evaluation = generate_evaluation_matrix_x(xk, xe, Nx, ep1)
# Precompute 1D interpolation matrix
X_all = X_evaluation/X_collocation
# Evaluate the function f at collocation points and evaluation points
F = f(xx_k)
f_vals = f(xx_e)
# Compute the interpolated values
s = evaluate_s_1D(X_all, F)
## ----------------- OBSERVABLES COMPUTATION ------------------------------
# Maximum error
Hermite_max_error = maximum(abs.(s .- f_vals))
# L2 error. Using the fact that the evaluation grid is uniform
dx = abs(xe[2] - xe[1])
Hermite_L2_error = sqrt(trapz_1D((s - f_vals).^2, dx))
[Hermite_max_error, Hermite_L2_error]
end
"""
interpolate_2D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
Computes an interpolant of a 2D function via Hermite-tensor method.
Arguments:
- function_name : the name of the tested function.
Options_available: "f_3"
- nodes_types : an array of two strings with the spacing of nodes in corresponding directions.
Options available: "Uniform", "Chebyshev"
- epsilon : a vector of the values of shape parameters in corresponding directions.
The values cannot be larger than 1!
- N : a vector with the numbers of basis functions in corresponding directions.
- Ne : a vector with the numbers of evaluation points in corresponding directions.
- interpolation_type : a string with the name of preferred basis.
Options available: "Hermite", "Radial"
gamma - a factor allowing to scale the domain and the function accordingly.
© Anna Yurova, 2017
"""
function interpolate_2D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
# Set the nodes types in all directions.
nodes_type_x = nodes_types[1]
nodes_type_y = nodes_types[2]
## ----------------------- TEST CASES -------------------------------------
if function_name == "f_3"
# Setup the interpolated function.
f(x, y) = cos.(x.^2 + y.^2)
# Setup the bounds of the interpolation domain.
xmin = -1
xmax = 1
ymin = -1
ymax = 1
else
error("Wrong name of the test case! Possible options: f_3")
end
## -------------- COLLOCATION AND EVALUATION MATRICES ---------------------
# Setup the basis used for the interpolation.
# Declare the functions for the generating of the collocation and evaluation
# matrices according to the chosed basis.
if interpolation_type == "Hermite"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
elseif interpolation_type == "Radial"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
else
error("Wrong basis type! Possible options are: Radial, Hermite.")
end
## -------------------- NODES TYPES ---------------------------------------
# Declare the nodes generating functions in all directions according to the specified by the user.
# X - direction
if nodes_type_x == "Chebyshev"
generate_nodes_x = (Nx) -> generate_chebyshev_nodes(Nx, xmin, xmax)
elseif nodes_type_x == "Uniform"
generate_nodes_x = (Nx) -> generate_uniform_nodes(Nx, xmin, xmax)
else
@error("Wrong node type! Possible options: :Chebyshev, :Uniform")
end
# Y - direction
if nodes_type_y == "Chebyshev"
generate_nodes_y = (Ny) -> generate_chebyshev_nodes(Ny, ymin, ymax)
elseif nodes_type_y == "Uniform"
generate_nodes_y = (Ny) -> generate_uniform_nodes(Ny, ymin, ymax)
else
@error("Wrong node type! Possible options: :Chebyshev, :Uniform")
end
## ---------------- INTERPOLATION PARAMETERS ------------------------
# Setup the value of epsilon
ep1 = epsilon[1]
ep2 = epsilon[2]
# Setup the number of collocation points.
Nx = N[1]
Ny = N[2]
#Setup the number of evaluation points.
Ne1 = Ne[1]
Ne2 = Ne[2]
# Generate node points.
xk = generate_nodes_x(Nx)
yk = generate_nodes_y(Ny)
grid_col = ndgrid(xk, yk)
xx_k = grid_col[1]
yy_k = grid_col[2]
# Generate evaluation points
xe = linspace(xmin, xmax, Ne1)
ye = linspace(ymin, ymax, Ne2)
grid_eval = ndgrid(xe, ye)
xx_e = grid_eval[1]
yy_e = grid_eval[2]
## ------------------- ACTUAL INTERPOLATION -------------------------------
# Calculate collocation and evaluation matrices
X_collocation = generate_collocation_matrix_x(xk, Nx, ep1)
X_evaluation = generate_evaluation_matrix_x(xk, xe, Nx, ep1)
Y_collocation = generate_collocation_matrix_y(yk, Ny, ep2)
Y_evaluation = generate_evaluation_matrix_y(yk, ye, Ny, ep2)
# Precompute 1D interpolation matrices
X_all = X_evaluation/X_collocation
Y_all = Y_evaluation/Y_collocation
# Evaluate the function f at collocation points and evaluation points
F = f(xx_k, yy_k)
f_vals = f(xx_e, yy_e)
# Compute the interpolated values
s = evaluate_s_2D(X_all, Y_all, F)
## ----------------- OBSERVABLES COMPUTATION ------------------------------
# Maximum error
Hermite_max_error = maximum(abs.(s .- f_vals))
dx = abs(xe[2] - xe[1])
dy = abs(ye[2] - ye[1])
# L2 error. Using the fact that the evaluation grid is uniform
Hermite_L2_error = sqrt(trapz_2D((s .- f_vals).^2, dx, dy))
[Hermite_max_error, Hermite_L2_error]
end
"""
interpolate_3D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
Computes an interpolant of a 3D function via Hermite-tensor method.
Arguments:
- function_name : the name of the tested function.
Options_available: "f_3"
- nodes_types : an array of three strings with the spacing of nodes in corresponding directions.
Options available: "Uniform", "Chebyshev"
- epsilon : a vector of the values of shape parameters in corresponding directions.
The values cannot be larger than 1!
- N : a vector with the numbers of basis functions in corresponding directions.
- Ne : a vector with the numbers of evaluation points in corresponding directions.
- interpolation_type : a string with the name of preferred basis.
Options available: :Hermite, :Radial
- gamma : a factor allowing to scale the domain and the function accordingly.
© Anna Yurova, 2017
"""
function interpolate_3D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
# Set the nodes types in all directions.
nodes_type_x = nodes_types[1]
nodes_type_y = nodes_types[2]
nodes_type_z = nodes_types[3]
## ----------------------- TEST CASES -------------------------------------
if function_name == "f_3"
# Setup the interpolated function.
f(x, y, z) = cos(x.^2 + y.^2 + z.^2)
# Setup the bounds of the interpolation domain.
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
else
error("Wrong name of the test case! Possible options: f_3")
end
## -------------- COLLOCATION AND EVALUATION MATRICES ---------------------
# Setup the basis used for the interpolation.
# Declare the functions for the generating of the collocation and evaluation
# matrices according to the chosed basis.
if interpolation_type == "Hermite"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
generate_collocation_matrix_z = (zk, N, epsilon) -> evaluate_hermite(zk, N, epsilon, gamma)
generate_evaluation_matrix_z = (zk, ze, N, epsilon) -> evaluate_hermite(ze, N, epsilon, gamma)
elseif interpolation_type == "Radial"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
generate_collocation_matrix_z = (zk, N, epsilon) -> evaluate_radial(zk, zk, N, epsilon)
generate_evaluation_matrix_z = (zk, ze, N, epsilon) -> evaluate_radial(zk, ze, N, epsilon)
else
@error ("Wrong basis type! Possible options are: Radial, Hermite.")
end
## -------------------- NODES TYPES ---------------------------------------
# Declare the nodes generating functions in all directions according to the specified by the user.
# X - direction
if nodes_type_x == "Chebyshev"
generate_nodes_x = (Nx) -> generate_chebyshev_nodes(Nx, xmin, xmax)
elseif nodes_type_x == "Uniform"
generate_nodes_x = (Nx) -> generate_uniform_nodes(Nx, xmin, xmax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")
end
# Y - direction
if nodes_type_y == "Chebyshev"
generate_nodes_y = (Ny) -> generate_chebyshev_nodes(Ny, ymin, ymax)
elseif nodes_type_y == "Uniform"
generate_nodes_y = (Ny) -> generate_uniform_nodes(Ny, ymin, ymax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")
end
# Z - direction
if nodes_type_z == :Chebyshev
generate_nodes_z = (Nz) -> generate_chebyshev_nodes(Nz, zmin, zmax)
elseif nodes_type_z == :Uniform
generate_nodes_z = (Nz) -> generate_uniform_nodes(Nz, zmin, zmax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")
end
## ---------------- INTERPOLATION PARAMETERS ------------------------
# Setup the value of epsilon
ep1 = epsilon[1]
ep2 = epsilon[2]
ep3 = epsilon[3]
# Setup the number of collocation points.
Nx = N[1]
Ny = N[2]
Nz = N[3]
#Setup the number of evaluation points.
Ne1 = Ne[1]
Ne2 = Ne[2]
Ne3 = Ne[3]
# Generate node points.
xk = generate_nodes_x(Nx)
yk = generate_nodes_y(Ny)
zk = generate_nodes_z(Nz)
grid_col = ndgrid(xk, yk, zk)
xx_k = grid_col[1]
yy_k = grid_col[2]
zz_k = grid_col[3]
# Generate evaluation points
xe = linspace(xmin, xmax, Ne1)
ye = linspace(ymin, ymax, Ne2)
ze = linspace(zmin, zmax, Ne3)
grid_eval = ndgrid(xe, ye, ze)
xx_e = grid_eval[1]
yy_e = grid_eval[2]
zz_e = grid_eval[3]
## ------------------- ACTUAL INTERPOLATION -----