Mit der Julia-Sprache habe ich eine Funktion definiert, mit der Punkte innerhalb der Kugel mit dem Radius 3.14 gleichmäßig abgetastet werden können. Dabei wird die Ablehnungsabtastung wie folgt verwendet:

function spherical_sample(N::Int64)

    # generate N points uniformly distributed inside sphere 
    # using rejection sampling:

    points = pi*(2*rand(5*N,3).-1.0)

    ind = sum(points.^2,dims=2) .<= pi^2

    ## ideally I wouldn't have to do this: 
    ind_ = dropdims(ind,dims=2)

    return points[ind_,:][1:N,:]

end

Ich habe einen Hack für die Untergruppe von Arrays gefunden:

ind = sum(points.^2,dims=2) .<= pi^2

## ideally I wouldn't have to do this: 
ind_ = dropdims(ind,dims=2)

Grundsätzlich sollte die Array-Indizierung jedoch einzeilig sein. Wie könnte ich das in Julia besser machen?

5
Aidan Rocke 7 Feb. 2020 im 19:16

4 Antworten

Beste Antwort

Das Problem ist, dass Sie einen zweidimensionalen Indexvektor erstellen. Sie können dies vermeiden, indem Sie eachrow verwenden:

ind = sum.(eachrow(points.^2)) .<= pi^2

Damit wäre Ihre vollständige Antwort:

function spherical_sample(N::Int64)
    points = pi*(2*rand(5*N,3).-1.0)
    ind = sum.(eachrow(points.^2)) .<= pi^2
    return points[ind,:][1:N,:]
end
4
David Varela 7 Feb. 2020 im 16:35

Dies beantwortet Ihre Frage nicht direkt (da Sie bereits zwei Vorschläge erhalten haben), sondern ich wollte eher darauf hinweisen, wie Sie das gesamte Verfahren anders implementieren können, wenn Sie möchten, dass es effizient ist.

Der erste Punkt besteht darin, das Generieren von 5*N Datenzeilen zu vermeiden. Das Problem besteht darin, dass es sehr wahrscheinlich nicht ausreicht, N gültige Stichproben zu generieren. Der Punkt ist, dass die Wahrscheinlichkeit einer gültigen Stichprobe in Ihrem Modell ~ 50% beträgt, sodass möglicherweise nicht genügend Punkte zur Auswahl stehen und die Auswahl von [1:N, :] einen Fehler auslöst.

Unten ist der Code, den ich verwenden würde, um dieses Problem zu vermeiden:

function spherical_sample(N::Integer) # no need to require Int64 only here
    points = 2 .* pi .* rand(N, 3) .- 1.0 # note that all operations are vectorized to avoid excessive allocations
    while N > 0 # we will run the code until we have N valid rows
        v = @view points[N, :] # use view to avoid allocating
        if sum(x -> x^2, v) <= pi^2 # sum accepts a transformation function as a first argument
            N -= 1 # row is valid - move to the previous one
        else
            rand!(v) # row is invalid - resample it in place
            @. v = 2 * pi * v - 1.0 # again - do the computation in place via broadcasting
        end
    end
    return points
end
4
Bogumił Kamiński 7 Feb. 2020 im 17:37

Dieser ist ziemlich schnell und verwendet StaticArrays. Sie können wahrscheinlich auch etwas Ähnliches mit gewöhnlichen Tupeln implementieren:

using StaticArrays

function sphsample(N)
    T = SVector{3, Float64}
    v = Vector{T}(undef, N)
    n = 1
    while n <= N
        p = rand(T) .- 0.5
        @inbounds v[n] = p .* 2π
        n += (sum(abs2, p) <= 0.25)
    end
    return v
end

Auf meinem Laptop ist es ~ 9x schneller als die Lösung mit view s.

3
DNF 7 Feb. 2020 im 19:31

Hier ist ein Einzeiler:

points[(sum(points.^2,dims=2) .<= pi^2)[:],:][1:N, :]

Beachten Sie, dass [:] eine Dimension löscht, damit BitArray für die Indizierung verwendet werden kann.

4
Przemyslaw Szufel 7 Feb. 2020 im 16:38