# Code review for first serious piece of Julia code
**URL:** https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577
**Category:** New to Julia
**Tags:** review
**Created:** [May 24, 2019, 3:57pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577 "2019-05-24T15:57:11Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@ludoro](https://discourse.julialang.org/u/ludoro)
#### Post date: [May 24, 2019, 3:57pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/1 "2019-05-24T15:57:11Z")
Hello everyone!
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:
```julia
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
```
---
### Author:  [@MrUrq](https://discourse.julialang.org/u/MrUrq)
#### Post date: [May 24, 2019, 4:22pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/2 "2019-05-24T16:22:53Z")
Hi Ludoro,
I am not sure if you are aware of the packages [KrigingEstimators.jl](https://github.com/juliohm/KrigingEstimators.jl) and [ScatteredInterpolation.jl](https://github.com/eljungsk/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.
---
### Author:  [@Nosferican](https://discourse.julialang.org/u/Nosferican)
#### Post date: [May 24, 2019, 6:29pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/3 "2019-05-24T18:29:06Z")
> [@ludoro](#):
>
> Differential equations
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](http://docs.juliadiffeq.org/latest/#Acknowledgements-1).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 24, 2019, 11:51pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/4 "2019-05-24T23:51:48Z")
> [@ludoro](#):
>
> `if kind == "linear"; q = 1; phi = z -> abs(z); elseif kind == "cubic"l q = 2; phi = z -> abs(z)^3; `
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.)
> [@ludoro](#):
>
> `D = zeros(Float32, size, size)`
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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 25, 2019, 12:34am UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/5 "2019-05-25T00:34:15Z")
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](https://github.com/juliohm/GeoStats.jl)) and the performance is great after many iterations to make it type stable.
---
### Author:  [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [May 25, 2019, 2:14am UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/6 "2019-05-25T02:14:20Z")
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 `kind`s.
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.
---
### Author:  [@ludoro](https://discourse.julialang.org/u/ludoro)
#### Post date: [May 25, 2019, 12:55pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/7 "2019-05-25T12:55:38Z")
For sure I will take a look at that package and get back to you! Thanks a lot!
---
### Author:  [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [May 25, 2019, 3:19pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/8 "2019-05-25T15:19:23Z")
A couple of other random things:
- 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`
---
### Author:  [@Wikunia](https://discourse.julialang.org/u/Wikunia)
#### Post date: [May 25, 2019, 7:21pm UTC](https://discourse.julialang.org/t/code-review-for-first-serious-piece-of-julia-code/24577/9 "2019-05-25T19:21:18Z")
Welcome to Julia 😉
You might find [Documentation · The Julia Language](https://docs.julialang.org/en/v1/manual/documentation/index.html) 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.
```julia
if length(x) != length(y)
error("Data length does not match")
end
```