I am Ludovico. I have been accepted to the Julia season of contribution, and during the community bonding time I started to tackle the first part of my proposal. I am building a new package called “Surrogates.jl”. My goal is to implement methods to approximate solution of Differential equations.
I started developing two methods: “Kriging” and “Response surfaces”. They both deal with radial basis functions and polynomials. As of right now, I just have the code for the one dimensional case, but the general case is not much more difficult. I am pretty sure that the code I have written is not “Julia-like” at all, so I would love to hear all the suggestions you have for it, if you do not mind.
I guess I should just post the first function, because I think that all the suggestions would work for the second one as well.
Below you can take a look at the code:
function Radial_1D(x,y,a,b,kind::String,lambda = 0)
#=
(x,y) set of nodes
(a,b) interval
kind is type of radial basis function
lambda is optional parameter with kind == multiquadric
=#
#1D Chebyshev is suggested
Chebyshev(x,k) = cos(k*acos(-1 + 2/(b-a)*(x-a)))
#Type of Radial basis function
#q is the number of polynomials in the basis
#The numbers are suggested by papers
if lambda === nothing
if kind == "linear"
q = 1
phi = z -> abs(z)
elseif kind == "cubic"
q = 2
phi = z -> abs(z)^3
elseif kind == "thinplate"
q = 2
phi = z -> z^2*log(abs(z))
else
error("Wrong type")
end
else
if kind == "multiquadric"
q = 1
phi = z -> sqrt(abs(z)^2 + lambda^2)
else
error("Wrong type")
end
end
if length(x) != length(y)
error("Data length does not match")
end
n = length(x)
#Find coefficients for both radial basis functions and polynomial terms
size = n+q
D = zeros(Float32, size, size)
d = zeros(Float32,size)
#In this array I have in the first n entries the coefficient of the radial
#basis function, in the last q term the coefficients for polynomial terms
coeff = zeros(Float32,size)
#Matrix made by 4 blocks:
#=
A | B
B^t | 0
A nxn, B nxq A,B symmetric so the matrix D is symmetric as well.
=#
for i = 1:n
d[i] = y[i]
for j = 1:n
D[i,j] = phi(x[i] - x[j])
end
for k = n+1:size
#Symmetric
D[i,k] = Chebyshev(x[i],k)
D[k,i] = D[i,k]
end
end
#Vector of size n + q containing in the first n terms the coefficients of
# the radial basis function and in the last q term the coefficient for the
# polynomials
return D\d
end
I am not sure if you are aware of the packages KrigingEstimators.jl and ScatteredInterpolation.jl which deal with the methods you are working with. You can probably take inspiration from those; alternatively, I am sure they’d be glad to accept a pull request adding any extra functionality you may need.
For any and all Differential Equations + Julia questions… a good go-to person would be @ChrisRackauckas (lead developer of DifferentialEquations.jl ecosystem). As your project develops if possible to get feedback from Chris, I am sure it would be of great benefit. Alternatively, a compilation of good people to check with are found here.
This kind of code makes life pretty hard for the compiler because it doesn’t know the value of kind until runtime, which means it can’t do any function inlining etcetera.
Instead, I would tend to suggest replacing the kind and lambda arguments with a basisfunc::AbstractBasisFunction argument, where you define functor types for the different basis functions. (Use functors because you want to both call the basis function and you want to extract other information like the degree.)
Don’t hard-code a type like Float32. Instead, e.g. you could use eltype(x) or float(eltype(x)) to decide on the type; that is, use the type of the input to decide on the type of the computation. That way the same code will work for different numeric types.
Thank you @MrUrq for sharing KrigingEstimators.jl.
@ludoro I would be happy to understand your plans for Surrogates.jl and help you get started. If all you need is a working Kriging implementation to approximate surfaces, then you are all set. The code has been tested in many cases (through GeoStats.jl) and the performance is great after many iterations to make it type stable.
As has been said, replacing kind::String and lambda with a type that can be dispatched on to avoid all the if elseif ... else statements would be nice as it will make your code more readable and easier to extend with new functionality by other packages that have their own kinds.
In places where you do find you want to specify function arguments as a string, it’s customary to instead use symbols instead ie :cubic instead of "cubic".
Finally, you’ll notice that the sort of comments you provide in your code are actually not very common in julia. Instead, such comments usually will end up in the documentation string so it can be looked up easily.
it looks like you are guaranteed that your for loops never go out of bounds by construction, ie. I think the user can’t provide an input that causes this function to go out of bounds. In this case, I think it’s advisable to put @inbounds in front of the outer for loop. This will turn off all your array bounds checking and possibly give a performance boost.
Instead of manually symmetrizing D in your for loop, you should maybe instead make D a Symmetric type. This could give a performance boost in D \ d
Welcome to Julia
You might find Documentation · The Julia Language valuable. It describes how to document functions.
I think most of the performance improvement ideas are already mentioned.
Just a small thing that isn’t Julia specific. I find it most readable and it improves the performance very very very slightly if you handle exceptions at the beginning of the function so I would put this at the very beginning of the function.
if length(x) != length(y)
error("Data length does not match")
end
if you handle exceptions at the beginning of the function so I would put this at the very beginning of the function.