diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..453925c --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1 @@ +style = "sciml" \ No newline at end of file diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 0000000..2a3517a --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,42 @@ +name: format-check + +on: + push: + branches: + - 'master' + - 'release-' + tags: '*' + pull_request: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v1 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/docs/make.jl b/docs/make.jl index 7f42f6f..60897d0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,18 +2,14 @@ using Documenter, PoissonRandom include("pages.jl") -makedocs( - sitename="PoissonRandom.jl", - authors="Chris Rackauckas", - modules=[PoissonRandom], - clean=true,doctest=false, - format = Documenter.HTML(analytics = "UA-90474609-3", - assets = ["assets/favicon.ico"], - canonical="https://poissonrandom.sciml.ai/stable/"), - pages=pages -) +makedocs(sitename = "PoissonRandom.jl", + authors = "Chris Rackauckas", + modules = [PoissonRandom], + clean = true, doctest = false, + format = Documenter.HTML(analytics = "UA-90474609-3", + assets = ["assets/favicon.ico"], + canonical = "https://poissonrandom.sciml.ai/stable/"), + pages = pages) -deploydocs( - repo = "github.com/SciML/PoissonRandom.jl.git"; - push_preview = true -) +deploydocs(repo = "github.com/SciML/PoissonRandom.jl.git"; + push_preview = true) diff --git a/docs/pages.jl b/docs/pages.jl index e1e30f6..9a91b48 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,6 +1,6 @@ # Put in a separate page so it can be used by SciMLDocs.jl -pages=[ +pages = [ "Home" => "index.md", "pois_rand.md", -] \ No newline at end of file +] diff --git a/src/PoissonRandom.jl b/src/PoissonRandom.jl index 61bccc5..7f840b0 100644 --- a/src/PoissonRandom.jl +++ b/src/PoissonRandom.jl @@ -25,70 +25,70 @@ end # ad_rand(λ) = ad_rand(Random.GLOBAL_RNG, λ) function ad_rand(rng::AbstractRNG, λ) - s = sqrt(λ) - d = 6.0*λ^2 - L = floor(Int,λ-1.1484) - # Step N - G = λ + s*randn(rng) - - if G >= 0.0 - K = floor(Int,G) - # Step I - if K >= L - return K - end - - # Step S - U = rand(rng) - if d*U >= (λ-K)^3 - return K - end - - # Step P - px,py,fx,fy = procf(λ,K,s) - - # Step Q - if fy*(1-U) <= py*exp(px-fx) - return K - end - end - - while true - # Step E - E = randexp(rng) - U = 2.0*rand(rng)-1.0 - T = 1.8+copysign(E,U) - if T <= -0.6744 - continue - end - - K = floor(Int,λ + s*T) - px,py,fx,fy = procf(λ,K,s) - c = 0.1069/λ - - # Step H - @fastmath if c*abs(U) <= py*exp(px+E)-fy*exp(fx+E) - return K - end - end + s = sqrt(λ) + d = 6.0 * λ^2 + L = floor(Int, λ - 1.1484) + # Step N + G = λ + s * randn(rng) + + if G >= 0.0 + K = floor(Int, G) + # Step I + if K >= L + return K + end + + # Step S + U = rand(rng) + if d * U >= (λ - K)^3 + return K + end + + # Step P + px, py, fx, fy = procf(λ, K, s) + + # Step Q + if fy * (1 - U) <= py * exp(px - fx) + return K + end + end + + while true + # Step E + E = randexp(rng) + U = 2.0 * rand(rng) - 1.0 + T = 1.8 + copysign(E, U) + if T <= -0.6744 + continue + end + + K = floor(Int, λ + s * T) + px, py, fx, fy = procf(λ, K, s) + c = 0.1069 / λ + + # Step H + @fastmath if c * abs(U) <= py * exp(px + E) - fy * exp(fx + E) + return K + end + end end # log(1+x)-x # accurate ~2ulps for -0.227 < x < 0.315 function log1pmx_kernel(x::Float64) - r = x/(x+2.0) - t = r*r + r = x / (x + 2.0) + t = r * r w = @evalpoly(t, - 6.66666666666666667e-1, # 2/3 - 4.00000000000000000e-1, # 2/5 - 2.85714285714285714e-1, # 2/7 - 2.22222222222222222e-1, # 2/9 - 1.81818181818181818e-1, # 2/11 - 1.53846153846153846e-1, # 2/13 - 1.33333333333333333e-1, # 2/15 - 1.17647058823529412e-1) # 2/17 - hxsq = 0.5*x*x - r*(hxsq+w*t)-hxsq + 6.66666666666666667e-1, # 2/3 + 4.00000000000000000e-1, # 2/5 + 2.85714285714285714e-1, # 2/7 + 2.22222222222222222e-1, # 2/9 + 1.81818181818181818e-1, # 2/11 + 1.53846153846153846e-1, # 2/13 + 1.33333333333333333e-1, # 2/15 + 1.17647058823529412e-1) # 2/17 + hxsq = 0.5 * x * x + r * (hxsq + w * t) - hxsq end # use naive calculation or range reduction outside kernel range. @@ -97,48 +97,48 @@ function log1pmx(x::Float64) if !(-0.7 < x < 0.9) return log1p(x) - x elseif x > 0.315 - u = (x-0.5)/1.5 - return log1pmx_kernel(u) - 9.45348918918356180e-2 - 0.5*u + u = (x - 0.5) / 1.5 + return log1pmx_kernel(u) - 9.45348918918356180e-2 - 0.5 * u elseif x > -0.227 return log1pmx_kernel(x) elseif x > -0.4 - u = (x+0.25)/0.75 - return log1pmx_kernel(u) - 3.76820724517809274e-2 + 0.25*u + u = (x + 0.25) / 0.75 + return log1pmx_kernel(u) - 3.76820724517809274e-2 + 0.25 * u elseif x > -0.6 - u = (x+0.5)*2.0 - return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5*u + u = (x + 0.5) * 2.0 + return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5 * u else - u = (x+0.625)/0.375 - return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625*u + u = (x + 0.625) / 0.375 + return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625 * u end end # Procedure F function procf(λ, K::Int, s::Float64) # can be pre-computed, but does not seem to affect performance - ω = 0.3989422804014327/s - b1 = 0.041666666666666664/λ - b2 = 0.3*b1*b1 - c3 = 0.14285714285714285*b1*b2 - c2 = b2 - 15.0*c3 - c1 = b1 - 6.0*b2 + 45.0*c3 - c0 = 1.0 - b1 + 3.0*b2 - 15.0*c3 + ω = 0.3989422804014327 / s + b1 = 0.041666666666666664 / λ + b2 = 0.3 * b1 * b1 + c3 = 0.14285714285714285 * b1 * b2 + c2 = b2 - 15.0 * c3 + c1 = b1 - 6.0 * b2 + 45.0 * c3 + c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3 if K < 10 px = -float(λ) - py = λ^K/factorial(K) + py = λ^K / factorial(K) else - δ = 0.08333333333333333/K - δ -= 4.8*δ^3 - V = (λ-K)/K - px = K*log1pmx(V) - δ # avoids need for table - py = 0.3989422804014327/sqrt(K) + δ = 0.08333333333333333 / K + δ -= 4.8 * δ^3 + V = (λ - K) / K + px = K * log1pmx(V) - δ # avoids need for table + py = 0.3989422804014327 / sqrt(K) end - X = (K-λ+0.5)/s + X = (K - λ + 0.5) / s X2 = X^2 - fx = -0.5*X2 # missing negation in pseudo-algorithm, but appears in fortran code. - fy = ω*(((c3*X2+c2)*X2+c1)*X2+c0) - return px,py,fx,fy + fx = -0.5 * X2 # missing negation in pseudo-algorithm, but appears in fortran code. + fy = ω * (((c3 * X2 + c2) * X2 + c1) * X2 + c0) + return px, py, fx, fy end """ diff --git a/test/runtests.jl b/test/runtests.jl index 0a652d5..6d71025 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,8 +7,8 @@ n_tsamples = 10^5 function test_samples(rand_func, distr::Distributions.DiscreteUnivariateDistribution, n::Int; # number of samples to generate - q::Float64=1.0e-8, # confidence interval, 1 - q as confidence - verbose::Bool=false) # show intermediate info (for debugging) + q::Float64 = 1.0e-8, # confidence interval, 1 - q as confidence + verbose::Bool = false) # show intermediate info (for debugging) # The basic idea # ------------------ @@ -30,8 +30,8 @@ function test_samples(rand_func, vmin = minimum(distr) vmax = maximum(distr) - rmin = floor(Int,quantile(distr, 0.00001))::Int - rmax = floor(Int,quantile(distr, 0.99999))::Int + rmin = floor(Int, quantile(distr, 0.00001))::Int + rmax = floor(Int, quantile(distr, 0.99999))::Int m = rmax - rmin + 1 # length of the range p0 = Distributions.pdf.((distr,), rmin:rmax) # reference probability masses @assert length(p0) == m @@ -39,12 +39,12 @@ function test_samples(rand_func, # determine confidence intervals for counts: # with probability q, the count will be out of this interval. # - clb = Vector{Int}(undef,m) - cub = Vector{Int}(undef,m) - for i = 1:m + clb = Vector{Int}(undef, m) + cub = Vector{Int}(undef, m) + for i in 1:m bp = Distributions.Binomial(n, p0[i]) - clb[i] = floor(Int,quantile(bp, q/2)) - cub[i] = ceil(Int,Distributions.cquantile(bp, q/2)) + clb[i] = floor(Int, quantile(bp, q / 2)) + cub[i] = ceil(Int, Distributions.cquantile(bp, q / 2)) @assert cub[i] >= clb[i] end @@ -54,7 +54,7 @@ function test_samples(rand_func, # scan samples and get counts cnts = zeros(Int, m) - for i = 1:n + for i in 1:n @inbounds si = samples[i] if rmin <= si <= rmax cnts[si - rmin + 1] += 1 @@ -65,7 +65,7 @@ function test_samples(rand_func, end # check the counts - for i = 1:m + for i in 1:m verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts[i])") clb[i] <= cnts[i] <= cub[i] || error("The counts are out of the confidence interval.") @@ -75,15 +75,15 @@ end println("testing count random sampler") for λ in [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0] - test_samples(PoissonRandom.count_rand, Distributions.Poisson(λ), n_tsamples) + test_samples(PoissonRandom.count_rand, Distributions.Poisson(λ), n_tsamples) end println("testing ad random sampler") for λ in [5.0, 10.0, 15.0, 20.0, 30.0] - test_samples(PoissonRandom.ad_rand, Distributions.Poisson(λ), n_tsamples) + test_samples(PoissonRandom.ad_rand, Distributions.Poisson(λ), n_tsamples) end println("testing mixed random sampler") for λ in [5.0, 10.0, 15.0, 20.0, 30.0] - test_samples(pois_rand, Distributions.Poisson(λ), n_tsamples) + test_samples(pois_rand, Distributions.Poisson(λ), n_tsamples) end