format SciML Style

.JuliaFormatter.toml

@@ -0,0 +1 @@
+style = "sciml"

.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'

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)

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",
-]
+]

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
 
 """

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,21 +30,21 @@ 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
 
     # 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