1

I'm new to julia so I would welcome some advice to improve the following function,

using SpecialFunctions

function rb(x, nu_max)

  bj = Array{Complex64}(length(x), nu_max)

  nu = 0.5 + (0:nu_max)
  # somehow dot broadcast isn't happy
  # bj .= [ besselj(_nu,_x)*sqrt(pi/2*_x) for _nu in nu, _x in x]

  bj   = [ besselj(_nu,_x)*sqrt(pi/2*_x) for _nu in nu, _x in x]

end

rb(1.0:0.1:2.0, 500)

basically, I'm not quite sure what's the recommended way to get a matrix over these two parameters (x and nu). The documentation doesn't offer much information, but I understand that the underlying fortran routine internally loops over nu, so I'd rather not do it again in the interest of performance.

Edit: I'm asked about the goal; it's to compute the Riccati-Bessel functions $j_1(x,\nu),h_1(x,\nu)$ for multiple values of $x$ and $\nu$.

I've stripped down stylistic questions from the original version to focus on this core issue.

learnjulia
  • 183
  • 10
  • I suggest to clarify what is your goal. Write down a simple mathematical formula, and then people may try answer it. – juliohm Feb 06 '18 at 04:59
  • At first glance, and in general, try to avoid temporary allocations by preallocating arrays and filling them in-place (e.g. using [dot broadcasting](https://julialang.org/blog/2017/01/moredots)). Also maybe use `@inbounds`. – carstenbauer Feb 06 '18 at 08:52
  • Also, which Julia version are you running? Bessel functions have been moved to `SpecialFunctions` for Julia 0.6 and later. – carstenbauer Feb 06 '18 at 08:57
  • @crstnbr thanks for your suggestions. I get an error when I pre-allocate the array ("Dimension mismatch") even though it is the right dimension. I'm guessing this has something to do with `.=` having a comprehension on the right hand side? – learnjulia Feb 06 '18 at 09:15

2 Answers2

1

This is a great example where you can take full advantage of broadcasting. It looks like you want the cartesian product between x and nu, where the rows are populated by the values of nu and the columns are x. This is exactly what broadcasting can do — you just need to reshape x such that it's a single row across many columns:

julia> using SpecialFunctions

julia> x = 1.0:0.1:2.0
1.0:0.1:2.0

julia> nu = 0.5 + (0:500)
0.5:1.0:500.5

 # this shows how broadcast works — these are the arguments and their location in the matrix
julia> tuple.(nu, reshape(x, 1, :))
501×11 Array{Tuple{Float64,Float64},2}:
 (0.5, 1.0)    (0.5, 1.1)    …  (0.5, 1.9)    (0.5, 2.0)
 (1.5, 1.0)    (1.5, 1.1)       (1.5, 1.9)    (1.5, 2.0)
 (2.5, 1.0)    (2.5, 1.1)       (2.5, 1.9)    (2.5, 2.0)
 (3.5, 1.0)    (3.5, 1.1)       (3.5, 1.9)    (3.5, 2.0)
 ⋮                           ⋱                ⋮
 (497.5, 1.0)  (497.5, 1.1)     (497.5, 1.9)  (497.5, 2.0)
 (498.5, 1.0)  (498.5, 1.1)     (498.5, 1.9)  (498.5, 2.0)
 (499.5, 1.0)  (499.5, 1.1)     (499.5, 1.9)  (499.5, 2.0)
 (500.5, 1.0)  (500.5, 1.1)  …  (500.5, 1.9)  (500.5, 2.0)

julia> bj = besselj.(nu,reshape(x, 1, :)).*sqrt.(pi/2*reshape(x, 1, :))
501×11 Array{Float64,2}:
 0.841471    0.891207     0.932039     …  0.9463       0.909297
 0.301169    0.356592     0.414341        0.821342     0.870796
 0.0620351   0.0813173    0.103815        0.350556     0.396896
 0.00900658  0.0130319    0.0182194       0.101174     0.121444
 ⋮                                     ⋱               ⋮
 0.0         0.0          0.0             0.0          0.0
 0.0         0.0          0.0             0.0          0.0
 0.0         0.0          0.0             0.0          0.0
 0.0         0.0          0.0          …  0.0          0.0
mbauman
  • 30,958
  • 4
  • 88
  • 123
  • Nice writeup. But OP is probably asking for something else (see my post update). OP should probably clarify. – carstenbauer Feb 06 '18 at 14:22
  • Regarding the matrix construction alone, I believe using a comprehension with two loops or broadcasting is similarily fast. – carstenbauer Feb 06 '18 at 14:32
  • 1
    thanks for the example, I thought I'd have to manually create all combinations of x and nu to take advantage of dot broadcasting, but this is much cleaner. Then again, if comprehension is just as fast I'm not sure which is to be preferred here, probably a matter of taste? Anyway, sorry the original question conflated various issues – my main concern is that Amos' besselj should return all values for 1:nu (i.e. a whole column here). Re-doing this one by one in julia will slow down by a factor of 500. – learnjulia Feb 06 '18 at 18:26
  • 2
    never mind, it's an [open issue](https://github.com/JuliaMath/SpecialFunctions.jl/issues/26) – learnjulia Feb 06 '18 at 18:39
0

Elaborating on my comment above. At first glance, and in general, try to avoid temporary allocations by preallocating arrays and filling them in-place (e.g. using dot broadcasting). Also maybe use @inbounds.

To give you an impression, after

using SpecialFunctions
x = 1.0
nu_max = 3
nu = 0.5 + (0:nu_max)
f(nu,x) = besselj.(nu,x).*sqrt.(pi/2*x)

compare (using BenchmarkTools) performance (and allocations) of

bj   = hcat([ besselj.(_nu,x).*sqrt.(pi/2*x) for _nu in nu]...)

and

f.(nu,x)

(Technically the output is not identical, you would have to use vcat above, but anyway)

UPDATE (after OP purified his code):

Ok, I think I (finally) see your real question now (sorry for that). What I said above was about optimizing your original code with respect to how it calls besselj and efficiently processes it's output (see @Matt B.'s post for the nice full broadcast solution here).

IIUC, you want to exploit the fact (I don't know and didn't check if this is actually true) in the calculation of besselj for given nu and x internally there is a summation over nu. In other words you want to use intermediate results of this internal summation to avoid redundant calculations.

As SpecialFunctions's besselj just seems to call the Fortran routine (probably here) I doubt that you can access any of this information. Unfortunately I can't help you along here (I'd probably look for a pure Julia implementation of besselj).

carstenbauer
  • 9,817
  • 1
  • 27
  • 40
  • can you show what you would do with x being a numeric vector? As I said I'm not happy with mixing comprehension and dot broadcasting, but it feels like nu should be treated differently than x for this particular function (where the fortran routine returns an array at once). – learnjulia Feb 06 '18 at 09:18
  • First, you could give `f` a single `nu` and a vector `x` and again use `f.(nu,x)`. – carstenbauer Feb 06 '18 at 09:25
  • Second, don't use comprehensions in your case. Why write two comprehensions with same inner for loop (i.e. `for _nu in nu` etc.)? Write the loop explicitly. – carstenbauer Feb 06 '18 at 09:26
  • I suggest you try rewrite your code with my comments in mind and update your post. Then we have a better basis for discussion. – carstenbauer Feb 06 '18 at 09:28
  • OK, I've stripped down the example to get to the core of the question. I agree that for any other function it would make sense to write a single for loop in my original question; however the key aspect here is that besselj loops over nu values in fortran. – learnjulia Feb 06 '18 at 10:05
  • thanks for your edit; indeed I don't understand why the julia wrapper seems to return only one element when the fortran routine returns a vector for all nus=1:numax. – learnjulia Feb 06 '18 at 18:28