Take numbers form two intervals in concentric spheres in Julia

Question

I am trying to take numbers from two intervals in Julia. The problem is the following, I am trying to create concentric spheres and I need to generate vectors of dimension equal to 15 filled with numbers taken from each circle. The code is:

rmax = 5
ra = fill(0.0,1,rmax)

for i=1:rmax-1
    ra[:,i].=rad/i
    ra[:,rmax].= 0
end
for i=1:3
    ptset = Any[]
    for j=1:200
        yt= 0
        yt= rand(Truncated(Normal(0, 1), -ra[i], ra[i] ))
        if -ra[(i+1)] < yt <= -ra[i] || ra[(i+1)] <= yt < ra[i]
            push!(ptset,yt)
            if length(ptset) == 15
                break
            end
        end
        
    end
end

Here, I am trying to generate spheres with uniform random numbers inside of each one; In this case, yt is only part of the construction of the numbers inside the sphere.

I would like to generate points in a sphere with radius r0 (ra[:,4] for this case), then points distributed from the edge of the first sphere to the second one wit radius r1 (here ra[:,3]) and so on.

In order to do that, I try to take elements that fulfill one of the two conditions -ra[(i+1)] < yt <= -ra[i] or ra[(i+1)] <= yt < ra[i], i.e. I would like to generate a vector with positive and negative numbers. I used the operator || but it seems to take only the positive part. I am new in Julia and I am not sure how to take the elements from both parts of the interval. Does anyone has a hit on how to do it?. Thanks in advance

You have a couple of weird constructions there. Can you try to explain what distribution you want to achieve in plain words, without Julia concepts mixed in? — phipsgabler, Jun 25 '20 at 15:52
@phipsgabler Hi, this is just some part of the code, actually, I am trying to generate spheres with uniform random numbers inside, here `yt` is only part of the construction, but first I was trying to implement the vector's construction to see if this part is working. So, I want to generate points in a spere with radious r0, then pint distributed form the egde of the firs sphere to the second one wit radious r1, and so on. — mors, Jun 25 '20 at 16:00
So, to sample points uniformly distributed in successive [shells](https://en.wikipedia.org/wiki/Spherical_shell), each described by radii `r[i]` and `r[i+1]`? — phipsgabler, Jun 25 '20 at 16:11
Yes, I would like to have points uniformly distributes between the radii, I would like to have vectors of dimension 15 with numbes from each shell but with both negative and positive part, like `-r0<=yt<=r0` for the first, -`r0 — mors, Jun 25 '20 at 16:25
OK, I'll try, but you'll have to wait a bit until [this question](https://stackoverflow.com/q/62580541/1346276) is resolved :) — phipsgabler, Jun 25 '20 at 17:01
Thank you, I really appreciate it. I will aldotry to clarity the problem. — mors, Jun 25 '20 at 17:30
Yeah, please add the relevent information to the question body. — phipsgabler, Jun 25 '20 at 19:13

score 1 · Accepted Answer · answered Jun 25 '20 at 19:34

I hope I understood you correctly. First, we need to be able to sample uniformly from an N-dimensional shell with radii r0 and r1:

using Random
using LinearAlgebra: normalize

struct Shell{N}
    r0::Float64
    r1::Float64
end

Base.eltype(::Type{<:Shell}) = Vector{Float64}

function Random.rand(rng::Random.AbstractRNG, d::Random.SamplerTrivial{Shell{N}}) where {N}
    shell = d[]
    Δ = shell.r1 - shell.r0
    θ = normalize(randn(N))   # uniformly distributed N-dimensional direction of length 1
    r = shell.r0 .* θ  # scale to a point on the interior of the shell
    return r .+ Δ .* θ .* .√rand(N)  # add a uniformly random segment between r0 and r1
end

(See here for more info about hooking into Random. You could equally implement a new Distribution, but that's not really necessary.)

Most importantly, a truncated normal will not result in a uniform distribution, but neither will adding a uniform scaling into the right direction: see here for why the square root is necessary (and I hope I got it right; you should check the math once more).

Then we can just create a sequence of shell samples with nested radii:

rmax = 5
rad = 10.0
ra = range(0, rad, length=rmax)

ptset = [rand(Shell{2}(ra[i], ra[i+1]), 15) for i = 1:(rmax - 1)]

(This part I wasn't really sure about, but the point should be clear.)

Thanks! it is exactly what I needed. Now I just need to find how to change the output array form `ptset` from `15-element Array{Array{Float64,1},1}:` to `15-element Array{Float64,1}:` because I need to send the output to a function. Also, I need to learn how to construct properly a function as you did. Thank you again. — mors, Jun 26 '20 at 19:36

Take numbers form two intervals in concentric spheres in Julia

1 Answers1