Commit 3508cb21 authored by anna's avatar anna
Browse files

HermiteGF interpolation v0.1

parents
#=The BSD license:
——————————————
Copyright (c) 2017 Anna Yurova, Max-Planck-Institut für Plasmaphysik.
All rights reserved.
Redistribution and use in source and binary forms, with or without modifi-
cation, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
* Neither the name of the Max-Planck-Institut für Plasmaphysik nor the
names of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ”AS
IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABIL-
ITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
IN NO EVENT SHALL MAX-PLANCK-INSTITUT FÜR PLASMAPHYSIK
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EX-
EMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLI-
GENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
——————————————
In addition, we ask you to cite the following reference in scientific publica-
tions which contain results obtained with this software and developments:
A. Yurova, K. Kormann
“Stable evaluation of Gaussian radial basis functions using Hermite polyno-
mials”
=#
function evaluate_hermite(xk, N, epsilon, gamma)
# 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.
# This function should be used both for computation of the
# collocation and evaluation matrices.
# Initialize the result matrix
result = zeros(Int(length(xk)), N);
# Write the values of the first two Hermite functions to initialize the
# three-term recurrence
result[:, 1] = exp(-0.5.*((gamma*xk).^2));
result[:, 2] = gamma*xk.*sqrt(2).*exp(-0.5.*((gamma*xk).^2));
# Three term recurrence for the Hermite functions with argument gamma*x
for i = 3:N
result[:, i] = sqrt(2.0/(i-1)).*(gamma*xk).*result[:, i-1] - sqrt((i-2.0)/(i-1)).*result[:, i-2];
end
# Scaling the Hermite functions with the exponential factor
# exp(-epsilon^2 x^2 + (gamma x)^2/2)
for i=1:N
result[:, i] = (pi^(1/4))*result[:, i].*exp((xk.^2).*(gamma*gamma*0.5 - epsilon^2));
end
return result;
end
#=The BSD license:
——————————————
Copyright (c) 2017 Anna Yurova, Max-Planck-Institut für Plasmaphysik.
All rights reserved.
Redistribution and use in source and binary forms, with or without modifi-
cation, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
* Neither the name of the Max-Planck-Institut für Plasmaphysik nor the
names of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ”AS
IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABIL-
ITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
IN NO EVENT SHALL MAX-PLANCK-INSTITUT FÜR PLASMAPHYSIK
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EX-
EMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLI-
GENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
——————————————
In addition, we ask you to cite the following reference in scientific publica-
tions which contain results obtained with this software and developments:
A. Yurova, K. Kormann
“Stable evaluation of Gaussian radial basis functions using Hermite polyno-
mials”
=#
###
# In this file the functions for the serial computing of the interpolated values 's'
# in 1-4D are implemented. Here the tensor representation of 's' is used (see Sec. 4.1 of the paper).
###
function evaluate_s_4D(X_all, Y_all, Z_all, W_all, F)
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
Nw = size(W_all, 2);
Ne4 = size(W_all, 1);
# Initialize the result tensor
s = zeros(Ne1, Ne2, Ne3, Ne4);
# Computing s = (W_all \otimes Z_all \otimes Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
@inbounds for col_dim_4 = 1:Nw
@inbounds for eval_dim_4 = 1:Ne4
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3, eval_dim_4] = s[eval_dim_1,eval_dim_2, eval_dim_3, eval_dim_4] + X_all[eval_dim_1, col_dim_1]*Y_all[eval_dim_2, col_dim_2]*Z_all[eval_dim_3, col_dim_3]*W_all[eval_dim_4, col_dim_4]*F[col_dim_1, col_dim_2, col_dim_3, col_dim_4];
end
end
end
end
end
end
end
end
return s;
end
function evaluate_s_3D(X_all, Y_all, Z_all, F)
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
# Initialize the result tensor
s = zeros(Ne1, Ne2, Ne3);
# Computing s = (Z_all \otimes Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3] = s[eval_dim_1,eval_dim_2, eval_dim_3] + X_all[eval_dim_1, col_dim_1]*Y_all[eval_dim_2, col_dim_2]*Z_all[eval_dim_3 , col_dim_3]*F[col_dim_1, col_dim_2, col_dim_3];
end
end
end
end
end
end
return s;
end
function evaluate_s_2D(X_all, Y_all, F)
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
# Initialize the result tensor
s = zeros(Ne1, Ne2);
# Computing s = (Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2] = s[eval_dim_1,eval_dim_2] + X_all[eval_dim_1, col_dim_1]*Y_all[eval_dim_2, col_dim_2]*F[col_dim_1, col_dim_2];
end
end
end
end
return s;
end
function evaluate_s_1D(X_all, F)
# Extract the number of collocation and evaluation points
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
# Initialize the result tensor
s = zeros(Ne1);
# Computing s = X_all * vec(F) (see Sec. 4.1 of the paper)
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1] = s[eval_dim_1] + X_all[eval_dim_1, col_dim_1]*F[col_dim_1];
end
end
return s;
end
#=The BSD license:
——————————————
Copyright (c) 2017 Anna Yurova, Max-Planck-Institut für Plasmaphysik.
All rights reserved.
Redistribution and use in source and binary forms, with or without modifi-
cation, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
* Neither the name of the Max-Planck-Institut für Plasmaphysik nor the
names of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ”AS
IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABIL-
ITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
IN NO EVENT SHALL MAX-PLANCK-INSTITUT FÜR PLASMAPHYSIK
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EX-
EMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLI-
GENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
——————————————
In addition, we ask you to cite the following reference in scientific publica-
tions which contain results obtained with this software and developments:
A. Yurova, K. Kormann
“Stable evaluation of Gaussian radial basis functions using Hermite polyno-
mials”
=#
# The structure of the parallelization is taken from https://docs.julialang.org/en/latest/manual/parallel-computing
###
# In this file the functions for the parallel computing of the interpolated values 's'
# in 4 and 5D are implemented. Here the tensor representation of 's' is used (see Sec. 4.1 of the paper)
###
function evaluate_s_4D_exec!(X_all, Y_all, Z_all, W_all, F, s)
# This function is what we call from the rest of the code.
# It distributes the work among processes and collects the results.
@sync begin
for p in procs(s)
@show p
@async remotecall_wait(evaluate_s_4D_parallel_chunk!, p, X_all, Y_all, Z_all, W_all, F, s)
end
end
return s;
end
function evaluate_s_5D_exec!(X_all, Y_all, Z_all, W_all, V_all, F, s)
# This function is what we call from the rest of the code.
# It distributes the work among processes and collects the results.
@sync begin
for p in procs(s)
@show p
@async remotecall_wait(evaluate_s_5D_parallel_chunk!, p, X_all, Y_all, Z_all, W_all, V_all, F, s)
end
end
return s;
end
@everywhere function myrange_4D(q::SharedArray)
# This function retuns the (irange,jrange) indexes assigned to a worker.
# In the case of the HermiteGF code, q corresponds to the tensor 's' containing
# the values of the interpolant at the evaluation points
idx = indexpids(q)
if idx == 0
# This worker is not assigned a piece
return 1:0
end
# Split the number of points in 4-th dimension into <number of threads>
# equal chuncks
nchunks = length(procs(q))
splits = [round(Int, s) for s in linspace(0,size(q,4),nchunks+1)]
return splits[idx]+1:splits[idx+1]
end
@everywhere function myrange_5D(q::SharedArray)
# This function retuns the (irange,jrange) indexes assigned to a worker
# In the case of the HermiteGF code, q corresponds to the tensor 's' containing
# the values of the interpolant at the evaluation points
idx = indexpids(q)
if idx == 0
# This worker is not assigned a piece
return 1:0
end
# Split the number of points in 5-th dimension into <number of threads>
# equal chuncks
nchunks = length(procs(q))
splits = [round(Int, s) for s in linspace(0,size(q,5),nchunks+1)]
return splits[idx]+1:splits[idx+1]
end
@everywhere function evaluate_s_4D_parallel!(X_all, Y_all, Z_all, W_all, F, s, irange)
# This is the kernel where the actual multiplication happens.
# Every thread takes care of the slices within 'irange' of the evaluation points in 4-th dimension
#@show (irange) # display what irange each worker is taking care of
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
Nw = size(W_all, 2);
Ne4 = size(W_all, 1);
# Computing s(:, :, :, irange) = (W_all(irange, :) \otimes Z_all \otimes Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
for eval_dim_4 in irange
@inbounds for col_dim_4 = 1:Nw
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3, eval_dim_4] = s[eval_dim_1,eval_dim_2, eval_dim_3, eval_dim_4] + X_all[eval_dim_1, col_dim_1]*Y_all[eval_dim_2, col_dim_2]*Z_all[eval_dim_3, col_dim_3]*W_all[eval_dim_4, col_dim_4]*F[col_dim_1, col_dim_2, col_dim_3, col_dim_4];
end
end
end
end
end
end
end
end
return s;
end
@everywhere function evaluate_s_5D_parallel!(X_all, Y_all, Z_all, W_all, V_all, F, s, irange)
# This is the kernel where the actual multiplication happens.
# Every thread takes care of the slices within the 'irange' of the evaluation points in 5-th dimension
#@show (irange) # display what irange each worker is taking care of
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
Nw = size(W_all, 2);
Ne4 = size(W_all, 1);
Nv = size(V_all, 2);
Ne5 = size(V_all, 1);
# Computing s(:, :, :, :, irange) = (V_all(irange, :) \otimes W_all \otimes Z_all \otimes Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
for eval_dim_5 in irange
@inbounds for col_dim_5 = 1:Nv
@inbounds for col_dim_4 = 1:Nw
@inbounds for eval_dim_4 = 1:Ne4
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3, eval_dim_4, eval_dim_5] = s[eval_dim_1,eval_dim_2, eval_dim_3, eval_dim_4, eval_dim_5] + X_all[eval_dim_1, col_dim_1]*Y_all[eval_dim_2, col_dim_2]*Z_all[eval_dim_3, col_dim_3]*W_all[eval_dim_4, col_dim_4]*V_all[eval_dim_5, col_dim_5]*F[col_dim_1, col_dim_2, col_dim_3, col_dim_4, col_dim_5];
end
end
end
end
end
end
end
end
end
end
return s;
end
# This is an interface simplifying the call of the function from the rest of the code
@everywhere evaluate_s_4D_parallel_chunk!(X_all, Y_all, Z_all, W_all, F, s) = evaluate_s_4D_parallel!(X_all, Y_all, Z_all, W_all, F, s, myrange_4D(s))
@everywhere evaluate_s_5D_parallel_chunk!(X_all, Y_all, Z_all, W_all, V_all, F, s) = evaluate_s_5D_parallel!(X_all, Y_all, Z_all, W_all, V_all, F, s, myrange_5D(s))
#=The BSD license:
——————————————
Copyright (c) 2017 Anna Yurova, Max-Planck-Institut für Plasmaphysik.
All rights reserved.
Redistribution and use in source and binary forms, with or without modifi-
cation, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
* Neither the name of the Max-Planck-Institut für Plasmaphysik nor the
names of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ”AS
IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABIL-
ITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
IN NO EVENT SHALL MAX-PLANCK-INSTITUT FÜR PLASMAPHYSIK
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EX-
EMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLI-
GENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
——————————————
In addition, we ask you to cite the following reference in scientific publica-
tions which contain results obtained with this software and developments:
A. Yurova, K. Kormann
“Stable evaluation of Gaussian radial basis functions using Hermite polyno-
mials”
=#
###
# In this file the some helper functions are implemented. In particular,
# - evaluation of the radial basis
# - generation of chebyshev and uniform nodes (1D)
# - trapezoidal rule 1-5D for the UNIFORM grid
###
function evaluate_radial(vc, ve, N, epsilon)
# Evaluate radial basis
result = zeros(Int(length(ve)), Int(length(vc)));
for i = 1:Int(length(ve)) # loop over evaluation points
for i_c = 1:Int(length(vc)) # loop over center points
result[i, i_c] = exp((-epsilon^2).*(ve[i] - vc[i_c]).^2);
end
end
return result;
end
function generate_chebyshev_nodes(N, xmin, xmax)
# Generate chebyshev nodes on the interval [xmin xmax]
theta1 = (pi:-pi/(N-1):-1e-14);
xk = ((cos(theta1)+1)*(xmax-xmin))/2 + xmin;
return xk
end
function generate_uniform_nodes(N, xmin, xmax)
# Generate uniform nodes on the interval [xmin xmax]
nodes = linspace(-1, 1, N);
xk = ((nodes + 1) * (xmax - xmin)) /2 + xmin;
return xk
end
# Trapezoidal rule for the uniform grid for 1-5D
function trapz_1D(values, dx)
return (sum(values) - values[1]/2 - values[end]/2)*dx;
end
function trapz_2D(values, dx, dy)
temp = zeros(size(values, 1));
for i = 1:size(values, 1)
temp[i] = trapz_1D(values[i, :], dy);
end
return trapz_1D(temp, dx);
end
function trapz_3D(values, dx, dy, dz)
temp = zeros(size(values, 1));
for i = 1:size(values, 1)
temp[i] = trapz_2D(values[i, :, :], dy, dz);
end
return trapz_1D(temp, dx);
end
function trapz_4D(values, dx, dy, dz, dw)
temp = zeros(size(values, 1));
for i = 1:size(values, 1)
temp[i] = trapz_3D(values[i, :, :, :], dy, dz, dw);
end
return trapz_1D(temp, dx);
end
function trapz_5D(values, dx, dy, dz, dw, dv)
temp = zeros(size(values, 1));
for i = 1:size(values, 1)
temp[i] = trapz_4D(values[i, :, :, :, :], dy, dz, dw, dv);
end
return trapz_1D(temp, dx);
end
#=The BSD license:
——————————————
Copyright (c) 2017 Anna Yurova, Max-Planck-Institut für Plasmaphysik.
All rights reserved.
Redistribution and use in source and binary forms, with or without modifi-
cation, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
* Neither the name of the Max-Planck-Institut für Plasmaphysik nor the
names of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ”AS
IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABIL-
ITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
IN NO EVENT SHALL MAX-PLANCK-INSTITUT FÜR PLASMAPHYSIK
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EX-
EMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLI-
GENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
——————————————
In addition, we ask you to cite the following reference in scientific publica-
tions which contain results obtained with this software and developments:
A. Yurova, K. Kormann
“Stable evaluation of Gaussian radial basis functions using Hermite polyno-
mials”
=#
###
# Include all necessary files for the interpolation. This has to be ran every time
# when Julia is started.
###
if !isdefined(:interpolate_1D)
include("helper_functions.jl")
include("ndgrid.jl")
include("evaluate_tensor_product.jl")
include("evaluate_tensor_product_parallel.jl")
include("evaluate_hermite.jl")
include("interpolation.jl")
else
print("You already included all necessary functions! If you changed the code of some function, recompile this function separately!")
end
This diff is collapsed.
Instructions for Linux.
In order to install Julia on your computer, perform the following steps:
- Download julia from https://julialang.org/downloads/ and unzip it.
- Note that you need to make sure curl and cmake are installed. On Ubuntu:
sudo apt-get install curl
sudo apt-get install cmake
- Download atom from https://atom.io/ and install it (On Ubuntu the package manager can be used).
- Install uber-juno through installation manager in atom.
- Set the Julia path to the Julia binary that was installed in the first step (Use settings -> packages -> Julia -> setting).
- Start Julia.
- If you are installing PyCall and PyPlot for the first time, just do ENV["PYTHON"]="" followed by Pkg.build("PyCall") before running Pkg.add("PyPlot").
Before running simulations, always run the file "init.jl". It is necessary to do it every time you start Julia, but not for every simulation.
In order to run parallel simulations, it is necessary to start Julia with appropriate amount of processes. Command line example:
julia -p 16 -L $HOME/hermiteGF/Julia/init.jl $HOME/hermiteGF/Julia/test_dependence_on_N.jl
# This file is a part of Julia. License is MIT: https://julialang.org/license
ndgrid(v::AbstractVector) = copy(v)
function ndgrid{T}(v1::AbstractVector{T}, v2::AbstractVector{T})
m, n = length(v1), length(v2)
v1 = reshape(v1, m, 1)
v2 = reshape(v2, 1, n)
(repmat(v1, 1, n), repmat(v2, m, 1))
end
function ndgrid_fill(a, v, s, snext)
for j = 1:length(a)
a[j] = v[div(rem(j-1, snext), s)+1]
end
end
function ndgrid{T}(vs::AbstractVector{T}...)
n = length(vs)
sz = map(length, vs)
out = ntuple(i->Array{T}(sz), n)
s = 1
for i=1:n
a = out[i]::Array
v = vs[i]
snext = s*size(a,i)