Initial commit

.gitignore

+
+docs/build

.travis.yml

+before_install:
+  - sudo apt-get install -y curl cmake
+
+addons:
+  apt:
+    update: true
+
+language: julia
+
+os: linux
+julia: 1.0
+
+before_install:
+
+notifications:
+  email: false
+
+after_success:
+  - julia -e 'using Pkg; cd(Pkg.dir("HermiteGF")); Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'
+  - julia -e 'using Pkg; cd(Pkg.dir("HermiteGF")); Pkg.add("Documenter"); include(joinpath("docs", "make.jl"))';

LICENSE

+BSD 2-Clause License
+
+Copyright (c) 2018, Julia codes for Vlasov equation
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, 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.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 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 NEGLIGENCE 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”

Project.toml

+authors = ["Anna Yurova <anna.yurova@ipp.mpg.de>"]
+name = "HermiteGF"
+uuid = "6d312152-d853-11e8-033d-978d5c49d094"
+version = "0.1.0"
+
+[deps]

README.md

+# HermiteGF.jl
+
+Stable Gaussian radial basis function interpolation based on HermiteGF expansion
+
+HermiteGF stabilization code for the RBF interpolation. 
+
+Author: Anna Yurova
+
+This is an implementation of the method described in the paper
+
+"STABLE EVALUATION OF GAUSSIAN RADIAL BASIS FUNCTIONS USING HERMITE POLYNOMIALS"
+
+by Anna Yurova and Katharina Kormann.
+
+https://arxiv.org/abs/1709.02164
+
+- HermiteGF-tensor in 1-5D. The implementation of 4-5D cases is parallel. 5D tests have to be run on a cluster to be finished in a reasonable time.
+
+In order to install Julia on your computer, perform the following steps:
+
+- Download julia from https://julialang.org/downloads/ and unzip it.  Note: The plotting currently only works with Julia v0.5.
+- 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
+
+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”*

docs/make.jl

+push!(LOAD_PATH,"../src/")
+
+using Documenter, Example
+
+makedocs(sitename="HermiteGF")

src/HermiteGF.jl

+module HermiteGF
+
+  export 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")
+
+end # module

src/evaluate_hermite.jl

+#=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

src/evaluate_tensor_product.jl

+#=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

src/evaluate_tensor_product_parallel.jl

+#=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))

src/helper_functions.jl

+#=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);