Stackoverflow Error in plot of 3D data point

Movicmatics · May 19, 2023, 3:59am

Hello everyone, I am trying to run my code but I am encountering an error, I would be glad if help can be render.
The Error mesage is:
Screenshot 2023-05-19 125724

ERROR: StackOverflowError:
Stacktrace:
  [1] RecipesPipeline.Volume(v::Array{Float64, 3}, x_extents::Array{Float64, 3}, y_extents::Array{Float64, 3}, z_extents::Tuple{Float64, Float64}) (repeats 65054 times)
    @ RecipesPipeline C:\Users\morak\.julia\packages\RecipesPipeline\BGM3l\src\utils.jl:109
  [2] RecipesPipeline.Volume(v::Array{Float64, 3}, x_extents::Array{Float64, 3}, y_extents::Array{Float64, 3})
    @ RecipesPipeline C:\Users\morak\.julia\packages\RecipesPipeline\BGM3l\src\utils.jl:109
  [3] macro expansion
    @ C:\Users\morak\.julia\packages\RecipesPipeline\BGM3l\src\user_recipe.jl:210 [inlined]
  [4] apply_recipe(::AbstractDict{Symbol, Any}, ::AbstractArray{<:Union{Missing, Number}, 3}, ::Any, ::Any)
    @ RecipesPipeline C:\Users\morak\.julia\packages\RecipesBase\BRe07\src\RecipesBase.jl:300
  [5] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline C:\Users\morak\.julia\packages\RecipesPipeline\BGM3l\src\user_recipe.jl:38
  [6] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline C:\Users\morak\.julia\packages\RecipesPipeline\BGM3l\src\RecipesPipeline.jl:72
  [7] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots C:\Users\morak\.julia\packages\Plots\esM5q\src\plot.jl:223
  [8] plot(::Any, ::Vararg{Any}; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
    @ Plots C:\Users\morak\.julia\packages\Plots\esM5q\src\plot.jl:102
  [9] scatter(::Any, ::Vararg{Any}; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
    @ Plots C:\Users\morak\.julia\packages\RecipesBase\BRe07\src\RecipesBase.jl:427
 [10] eval
    @ .\boot.jl:370 [inlined]
 [11] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1864
 [12] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base .\essentials.jl:816
 [13] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base .\essentials.jl:813
 [14] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\eval.jl:233
 [15] (::VSCodeServer.var"#66#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\eval.jl:157
 [16] withpath(f::VSCodeServer.var"#66#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\repl.jl:249
 [17] (::VSCodeServer.var"#65#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\eval.jl:155
 [18] hideprompt(f::VSCodeServer.var"#65#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})  
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [19] (::VSCodeServer.var"#64#68"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\eval.jl:126
 [20] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:514
 [21] with_logger
    @ .\logging.jl:626 [inlined]
 [22] (::VSCodeServer.var"#63#67"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\eval.jl:225
 [23] #invokelatest#2
    @ .\essentials.jl:816 [inlined]
 [24] invokelatest(::Any)
    @ Base .\essentials.jl:813
 [25] macro expansion
    @ c:\Users\morak\.vscode\extensions\julialang.language-julia-1.47.2\scripts\packages\VSCodeServer\src\eval.jl:34 [inlined]
 [26] (::VSCodeServer.var"#61#62")()
    @ VSCodeServer .\task.jl:514`




This is the code that i am trying to run



'''using LinearAlgebra
using SpecialFunctions
using Plots
pythonplot() 
using QuadGK
using FastGaussQuadrature


# STEP 1 parameters, coordinate of X(n_1, n_2)
N_1 = 50
N_2 = 50
N_3 = 50
R_1 = 10
R_2 = 10
R_3 = 10
P = 10
N_phi = 90
N_r = 30
N_theta = 60
a_1 = 1
a_2 = 1.3
a_3 = 1.5

# Meshgrid
x_1 = range(-0.5 * R_1, stop = 0.5 * R_1, length = N_1)
x_2 = range(-0.5 * R_2, stop = 0.5 * R_2, length = N_2)
x_3 = range(-0.5 * R_3, stop = 0.5 * R_3, length = N_3)

X1 = similar(x_1, N_1, N_2, N_3)
X2 = similar(x_2, N_1, N_2, N_3)
X3 = similar(x_3, N_1, N_2, N_3)

for i in 1:N_1
    for j in 1:N_2
        for k in 1:N_3
            X1[i, j, k] = x_1[i]
            X2[i, j, k] = x_2[j]
            X3[i, j, k] = x_3[k]
        end
    end
end
#X1, X2, X3 = [collect(it) for it in Iterators.product(x_1, x_2, x_3)]
println("step 1 is done")

# STEP 2: Evaluate Rho at X(n_1, n_2, n_3), Beta, g(X)
B = 1 / a_1 + 1 / a_2 + 1 / a_3
g = zeros(N_1, N_2, N_3)
for i in 1:N_1
    for j in 1:N_2
        for k in 1:N_3
            g[i, j, k] = (x_1[i]^2 / a_1) + (x_2[j]^2 / a_2) + (x_3[k]^2 / a_3)
        end
    end
end


# Computation of Rho
Rho = (-2 * B .+ 4 * ((x_1 .^ 2) ./ a_1 .+ (x_2 .^ 2) ./ a_2 .+ (x_3 .^ 2) ./ a_3)) .* exp.(-g)
println("Step 2 is done")

# Transformation, Spherical grid plot
X_h, W_h = gauss(N_r)        # Nodes

r = 0.5 * P * (X_h .+ 1)   # transformation for node
w = 0.5 * P * W_h         # transformation for weight

K1 = zeros(N_phi, N_r)
K2 = zeros(N_phi, N_r)
K3 = zeros(N_phi, N_r)

h_phi = 2 * π / N_phi       # step size in phi-direction
h_theta = π / N_theta       # step size in theta-direction

phi = range(0, stop=2 * π, length=N_phi)
theta = range(0, stop=π, length=N_theta)

cs = cos.(phi)
sn = sin.(phi)

for q in 1:N_r
    for t in 1:N_theta
        K1[:, q] .= r[q] * cs
        K2[:, q] .= r[q] * sn
        K3[:, q] .= r[q] * sin(theta[t])
    end
end

display(plot(K1, K2, K3, seriestype = :scatter, aspect_ratio = :equal, title = "A polar grid in Fourier domain"))
savefig("polargrid.png")
println("step 3 is done")

# Computation of P_d(K) when d=3
p_3_k = exp.(-K1.^2 / a_1 .- K2.^2 / a_2 .- K3.^2 / a_3)

println("p_3_k done")

# Step 3: Apply oversampling factor
Oversampling_factor = 2
N_1_os = Oversampling_factor * N_1
N_2_os = Oversampling_factor * N_2
N_3_os = Oversampling_factor * N_3

R_1_os = R_1
R_2_os = R_2
R_3_os = R_3


x_1_os = range(-0.5 * R_1_os, stop = 0.5 * R_1_os, length = N_1_os)
x_2_os = range(-0.5 * R_2_os, stop = 0.5 * R_2_os, length = N_2_os)
x_3_os = range(-0.5 * R_3_os, stop = 0.5 * R_3_os, length = N_3_os)

X1_os = similar(x_1_os, N_1_os, N_2_os, N_3_os)
X2_os = similar(x_2_os, N_1_os, N_2_os, N_3_os)
X3_os = similar(x_3_os, N_1_os, N_2_os, N_3_os)

for i in 1:N_1
    for j in 1:N_2
        for k in 1:N_3
            X1_os[i, j, k] = x_1_os[i]
            X2_os[i, j, k] = x_2_os[j]
            X3_os[i, j, k] = x_3_os[k]
        end
    end
end

println("oversampling factor done")

Rho_T = zeros(ComplexF64, N_phi, N_r, N_theta)

for n1 in 1:N_1
    for n2 in 1:N_2
        for n3 in 1:N_3
            l = K1 .* x_1[n1] .+ K2 .* x_2[n2] .+ K3 .* x_3[n3]
            Rho_T .+= reshape(exp.(-1im .* l), (N_phi, N_r, 1)) .* Rho[n1, n2, n3]
        end
    end
end

Rho_T = Rho_T .* (R_1 / N_1) .* (R_2 / N_2) .* (R_3 / N_3)
println("Rho_T is done")


# Computation of K_T
K_T = 1 ./ (K1 .^ 2 .+ K2 .^ 2 .+ K3 .^ 2)
println("K_T done")

# Computation of I_1
I_1 = zeros(Complex{Float64}, N_1_os, N_2_os, N_3_os)
for i in 1:N_phi
    for q in 1:N_r
        for t in 1:N_theta
            thet = K1[i, q] .* X1_os .+ K2[i, q] .* X2_os .+ K3[i, q] .* X3_os
            I_1 .+= w[q] .* exp.(1im .* thet) .* K_T[i, q] .* Rho_T[i, q, t] .* (1 .- p_3_k[i, q]) .* h_phi
        end
    end
end
I_1 ./= (2 * π) ^ 3
println("I_1 is done")

# Computation of I_2
I_2 = zeros(Complex{Float64}, N_1, N_2, N_3)
for i in 1:N_phi
    for q in 1:N_r
        for t in 1:N_theta
            thet = K1[i, q] .* X1 .+ K2[i, q] .* X2 .+ K3[i, q] .* X3
            I_2 .+= w[q] .* exp.(1im .* thet) .* K_T[i, q] .* Rho_T[i, q, t] .* (1 .- p_3_k[i, q]) .* h_phi
        end
    end
end
I_2 ./= (2 * π) ^ 3
println("I_1 is done")

C = zeros(ComplexF64, N_1, N_2, N_3)
for i in 1:N_phi
    for q in 1:N_r
        for t in 1:N_theta
            if i <= size(C, 1) && q <= size(C, 2) && t <= size(C, 3)
                C[i, q, t] = I_1[i, q, t] + I_2[i, q, t]
            end
        end
    end
end
println("C is done")

#Computations of analytic solution C_a
C_a = - exp.(-g)
println("C_a is done")


Plots.scatter(
    X1,
    X2,
    X3,
    color = C_a,
    marker = (:auto, 0.8, :viridis),
    xlabel = "X1",
    ylabel = "X2",
    zlabel = "X3",
    title = "Graph of C_a"
)'''

ufechner7 · May 19, 2023, 5:08am

Can you tell us the values of X1, X2 and X3?

This works fine for me:

using Plots
X1=rand(10)
X2=rand(10)
X3=rand(10)
Plots.scatter(X1,X2,X3)

Movicmatics · May 19, 2023, 5:38am

# Meshgrid
x_1 = range(-0.5 * R_1, stop = 0.5 * R_1, length = N_1)
x_2 = range(-0.5 * R_2, stop = 0.5 * R_2, length = N_2)
x_3 = range(-0.5 * R_3, stop = 0.5 * R_3, length = N_3)

X1 = similar(x_1, N_1, N_2, N_3)
X2 = similar(x_2, N_1, N_2, N_3)
X3 = similar(x_3, N_1, N_2, N_3)[quote="ufechner7, post:2, topic:99093, full:true"]
Can you tell us the values of X1, X2 and X3?

This works fine for me:

using Plots
X1=rand(10)
X2=rand(10)
X3=rand(10)
Plots.scatter(X1,X2,X3)

[/quote]

i created a meshgrid with the values of X1,X2,X3

ufechner7 · May 19, 2023, 7:23am

Your example is still not self-contained. Please read: make it easier to help you , in particular point four.

Movicmatics · May 19, 2023, 8:23am

i am trying to make a plot with multiple axis

using LinearAlgebra
using SpecialFunctions
using Plots
pythonplot() 
using QuadGK
using FastGaussQuadrature


# STEP 1 parameters, coordinate of X(n_1, n_2)
N_1 = 10
N_2 = 10
N_3 = 10
R_1 = 5
R_2 = 5
R_3 = 5
P = 10
N_phi = 90
N_r = 30
N_theta = 60
a_1 = 1
a_2 = 1.3
a_3 = 1.5

# Meshgrid
x_1 = range(-0.5 * R_1, stop = 0.5 * R_1, length = N_1)
x_2 = range(-0.5 * R_2, stop = 0.5 * R_2, length = N_2)
x_3 = range(-0.5 * R_3, stop = 0.5 * R_3, length = N_3)

X1 = similar(x_1, N_1, N_2, N_3)
X2 = similar(x_2, N_1, N_2, N_3)
X3 = similar(x_3, N_1, N_2, N_3)

for i in 1:N_1
    for j in 1:N_2
        for k in 1:N_3
            X1[i, j, k] = x_1[i]
            X2[i, j, k] = x_2[j]
            X3[i, j, k] = x_3[k]
        end
    end
end
#X1, X2, X3 = [collect(it) for it in Iterators.product(x_1, x_2, x_3)]
println("step 1 is done")

# STEP 2: Evaluate Rho at X(n_1, n_2, n_3), Beta, g(X)
B = 1 / a_1 + 1 / a_2 + 1 / a_3
g = zeros(N_1, N_2, N_3)
for i in 1:N_1
    for j in 1:N_2
        for k in 1:N_3
            g[i, j, k] = (x_1[i]^2 / a_1) + (x_2[j]^2 / a_2) + (x_3[k]^2 / a_3)
        end
    end
end


# Computation of Rho
Rho = (-2 * B .+ 4 * ((x_1 .^ 2) ./ a_1 .+ (x_2 .^ 2) ./ a_2 .+ (x_3 .^ 2) ./ a_3)) .* exp.(-g)
println("Step 2 is done")

# Transformation, Spherical grid plot
X_h, W_h = gauss(N_r)        # Nodes

r = 0.5 * P * (X_h .+ 1)   # transformation for node
w = 0.5 * P * W_h         # transformation for weight

K1 = zeros(N_phi, N_r, N_theta)
K2 = zeros(N_phi, N_r, N_theta)
K3 = zeros(N_phi, N_r, N_theta)

h_phi = 2 * π / N_phi       # step size in phi-direction
h_theta = π / N_theta       # step size in theta-direction

phi = range(0, stop=2 * π, length=N_phi)
theta = range(0, stop=π, length=N_theta)

cs = cos.(phi)
sn = sin.(phi)

for q in 1:N_r
    for t in 1:N_theta
        for p in 1:N_phi
            K1[p, q, t] = r[q] * cs[p] * sin(theta[t])
            K2[p, q, t] = r[q] * sn[p] * sin(theta[t])
            K3[p, q, t] = r[q] * cos(theta[t])
        end
    end
end

display(plot(K1, K2, K3, seriestype = :scatter, aspect_ratio = :equal, title = "A polar grid in Fourier domain"))

cormullion · May 19, 2023, 12:16pm

might be a good idea to report this on the plots.jl repo on github. I don’t think it should do a StackOverflow.

ufechner7 · May 19, 2023, 12:51pm

I can reproduce your problem.

K1, K2 and K3 have the type:

Array{Float64, 3}

and a size of

julia> size(K1)
(90, 30, 60)

Does this make sense?

As far as I understand the plot function expects three one dimensional vectors for creating a 3D scatter plot. Why do you pass three three dimensional arrays?

Movicmatics · May 22, 2023, 12:20am

Okay
Thanks

Movicmatics · May 22, 2023, 12:25am

yes you are right but in this case i want to plot 3 vector which have 3 dimensional array each