Simulation random points on a circle

General Usage
1

Hi,

I am currently working on simulating random points on a circle based on a non-uniform distribution. Initially, I successfully generated these random points using a set of code lines. However, some of the simulated random points fell outside the circle when wrapping these lines within a function and calling it. Despite ensuring the correctness of the input variables for the function, my guess is that when the code executes within the function, the vector containing the probabilities associated with specific polar coordinates becomes disordered, leading to points being generated outside the circle.

What could be the issue? Sorry but the code is too long to be pasted here.

2

If the code is too long to share, please share an example of an input/output pair (or more than one). There is a good chance this will be enough to generate some working code, which might be even shorter.

1 Like
3

Here is a minimal code with comments

struct containing the domain
circ = Circle(0., 0., R)

#generate random number
Np = rand(Poisson(Npar))

#create a matrix containing the cartesian position of the points
X = Array{Float64}(undef, 0, Nd);

function simulate_model(particle::Particle, circ::Circle, Np::Int64, Rk::Float64, D::Float64, X::Matrix{Float64})

for i in 1:Np
    #this first block returns the correct results
    x, y = GenerateHit(circ, Rk);
    track = Track(x, y, Rk)
    #integral contains the probability of generate a point at radius Ļ, theta the porbability of generating a point at angle Īø along radius Ļ
    #Dd is a variable needed for other purposes (but it is calculated correctly)
    #radius is a vector of possible radii
    integral, theta, Dd, radius = distribute_D(particle, circ, track);

    #next line returns points outside the circle
    X_, Y_ = calculate_position(particle, circ, integral, theta, Dd, radius);

    #I check if the position is within a radius circ.r, if outside is wrong
    dist = sqrt.(X_[:,1].*X_[:,1] .+ X_[:,2].*X_[:,2])
    if size(dist[dist .> circ.r],1) != 0
        println("Error")
        return X_
    end

    #I concatenate the matrix of positions
    X = vcat(X,X_)
    D += Dd
end

return X

end

The problem is the function calculate_position that gives different answers, whether executed inside the function or as a single line of code. I noted that changing variables from global to local will give different results without converging into a stable solution.

4

I donā€™t think much can be learned from the code you posted as it depends on a lot of things not included.

Are you saying that it is the call to

calculate_position(particle, circ, integral, theta, Dd, radius)

that gives you a problem? Then I would maybe amend your if block to:

if size(dist[dist .> circ.r],1) != 0
            println("Error")
            return (; X_, particle, circ, integral, thea, Dd, radius)
        end

so you can inspect which input isnā€™t what you expected?

5

The input seems correct. I am mostly wondering how it is possible that if I run the function I get correct positions but if I run a function that calls calculate_position I get wrong positions.

6

This thing in Julia can probably be written as (not tested):

any(>(circ.r), dist)
1 Like
7

To quote code enter three backquotes (`) on a separate line, to start code block, and end with another line with three backquotes.

It would be much easier to work with the code for the problematic calculate_position function, or at the minimum the parameters causing the output to go outside circle.

8

How are you doing your conversion from polar to cartesian coordinates and how are you sampling your \theta values? My guess is that youā€™re either sampling \theta outside of the [0,2\pi] range, or thereā€™s a problem when converting to cartesian coordinates. I just played around with the example below and it seems to work fine:

using Distributions
using StatsPlots

function circle_shape(h, k, r)
	Īø = LinRange(0, 2 * Ļ€, 500)
	h .+ r * sin.(Īø), k .+ r * cos.(Īø)
end

r = 10
circ = circle_shape(0,0,r)
d = truncated(Normal(5), 0, 2Ļ€)
Īø = rand(d, 10)
x = r .* cos.(Īø)
y = r .* sin.(Īø)
plot(circ, aspect_ratio=1, legend=false)
scatter!(x,y)

image

9

It sounds like the issue youā€™re seeing is particular to the function definition that you havenā€™t shared, so Iā€™m not sure thereā€™s much we can do to help you debug that blindly.

I will add, though, that you might be interested in Meshes.jl which has a lot of built-in tooling for geometric definitions, parametric functions for generating points, and even random sampling tools. Hereā€™s a quick example of

using Distributions
using Meshes

# Create a 2D surface bounded by a unit circle
center = Meshes.Point(0, 0)
radius = 1.0
circle = Meshes.Ball(center, radius)

# Some distributions from Distributions.jl with supports in [0,1]
dist_r = Normal(0.5, 0.1)
dist_phi = Normal(0.25, 0.1)  # centered on pi/2

# Get a collection of N sample Meshes.Point's on the surface of the circle
# Using your geometric variable name as a function turns it into a 
#   parametric function, e.g. circle(r, phi) -> the Point at that location
N = 10
sample_rs = rand(dist_r, N)
sample_phis = rand(dist_phis, N)
samples = circle.(sample_rs, sample_phis)
1 Like