Global pass on docstrings

Author: guilgautier

Project.toml

@@ -4,15 +4,10 @@ authors = ["Guillaume Gautier <[email protected]> and contributors"]
 version = "0.1.0"
 
 [deps]
-Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
-LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

src/PartialRejectionSampling.jl

@@ -1,5 +1,7 @@
 module PartialRejectionSampling
 
+# Imports
+
 using LinearAlgebra
 const LA = LinearAlgebra
 
@@ -11,16 +13,44 @@ const LG = LightGraphs
 using SimpleWeightedGraphs
 const SWG = SimpleWeightedGraphs
 
-using Plots, GraphPlot, Colors
-
 using Distances
 
-using LazySets
-const LS = LazySets
-
 using SpecialFunctions
 const SF = SpecialFunctions
 
+# Exports
+
+const PRS = PartialRejectionSampling
+export PRS
+
+export AbstractWindow
+export AbstractSpatialWindow,
+    RectangleWindow,
+    SquareWindow,
+    BallWindow
+export AbstractDiscreteWindow,
+    GraphNode
+
+export AbstractPointProcess
+export AbstractSpatialPointProcess,
+    HomogeneousPoissonPointProcess,
+    HardCorePointProcess,
+    StraussPointProcess
+export AbstractGraphPointProcess,
+    Ising,
+    HardCoreGraph,
+    RootedSpanningForest,
+    SinkFreeGraph
+export PatternFreeString
+
+export generate_sample,
+    generate_sample_prs,
+    generate_sample_grid_prs,
+    generate_sample_dcftp,         # For Spatial point processes
+    generate_sample_gibbs_perfect  # For the Ising model
+
+# Code inclusions
+
 include("common.jl")
 
 include("utils.jl")
@@ -35,17 +65,13 @@ include("poisson.jl")
 include("strauss.jl")
 include("hard_core_spatial.jl")
 
-include("hard_core_graph.jl")
-include("ising.jl")
-
 # Graph point processes
+include("ising.jl")
+include("hard_core_graph.jl")
 include("rooted_spanning_forest.jl")
 include("sink_free_graph.jl")
 
 # Misc
 include("pattern_free_string.jl")
 
-# Display
-include("display.jl")
-
 end

src/common.jl

@@ -9,3 +9,4 @@ dimension(pp::AbstractSpatialPointProcess) = dimension(window(pp))
 abstract type AbstractGraphPointProcess{T} <: AbstractPointProcess{T} end
 
 function generate_sample end
+function generate_sample_prs end

src/display.jl

@@ -1,18 +1,60 @@
-## Dominated Coupling From The Past (dCFTP)
-#  - [Kendall \& Moller's orginal formulation of dCFTP](https://www.researchgate.net/publication/  2821877_Perfect_Metropolis-Hastings_simulation_of_locally_stable_point_processes)
-#  - Huber, Perfect Simulation
-#  - [Kendall's notes on perfect simulation](https://warwick.ac.uk/fac/sci/statistics/staff/  academic-research/kendall/personal/ppt/428.pdf)
-#
-# Requirement: the target spatial point process must have the following methods
-#  - function papangelou_conditional_intensity end
-#  - function upper_bound_papangelou_conditional_intensity end
-#  - function window end
+"""
+Implementation of dominated Coupling From The Past (dCFTP)
+- [KeMo99](@cite) and [KeMo00](@cite) orginal formulation of dCFTP
+- [Hub16](@cite)
+- [Kendall's notes on perfect simulation](https://warwick.ac.uk/fac/sci/statistics/staff/  academic-research/kendall/personal/ppt/428.pdf)
+"""
 
-function generate_sample_dcftp(
+@doc raw"""
+    papangelou_conditional_intensity(pp, x, X)
+
+Compute the [Papangelou conditional intensity](https://en.wikipedia.org/wiki/Point_process#Papangelou_intensity_function) of the point process `pp`, as the ratio of densities ``\frac{f(X \cup x)}{f(X)}``.
+
+**See also**
+
+- Section 6.1.1 [MoWa04](@cite)
+- [wiki](https://en.wikipedia.org/wiki/Point_process#Papangelou_intensity_function)
+"""
+function papangelou_conditional_intensity end
+
+"""
+    upper_bound_papangelou_conditional_intensity
+
+Compute an upper bound of [`papangelou_conditional_intensity`](@ref)
+
+**See also**
+
+- Equation 2.1 [KeMo99](@cite)
+"""
+function upper_bound_papangelou_conditional_intensity end
+
+function isrepulsive end
+function isattractive end
+
+"""
+    generate_sample_dcftp(
         pp::AbstractSpatialPointProcess{T};
-        n₀::Int=1,
         win::Union{Nothing,AbstractWindow}=nothing,
+        n₀::Int=1,
         rng=-1
+    )::Vector{T} where {T}
+
+Generate an exact sample from a spatial point process `pp` on window optional window `win` using dominated coupling from the past.
+
+- Default window is `window(pp)=pp.window`
+- Initial coalescence check performed after `n₀` steps.
+- Seed or random number generator is passed via `rng`.
+
+**See also**
+
+    - [KeMo99](@cite), [KeMo00](@cite)
+    - Section 11.2.6 [MoWa04](@cite)
+"""
+function generate_sample_dcftp(
+    pp::AbstractSpatialPointProcess{T};
+    win::Union{Nothing,AbstractWindow}=nothing,
+    n₀::Int=1,
+    rng=-1
 )::Vector{T} where {T}
 
     @assert n₀ >= 1
@@ -23,29 +65,48 @@ function generate_sample_dcftp(
     birth_rate = β * volume(window_)
 
     # Dominating process
-    hp = HomogeneousPoissonPointProcess(β, window_)
-    D = Set{T}(eachcol(generate_sample(hp; rng=rng)))
+    hppp = HomogeneousPoissonPointProcess(β, window_)
+    D = Set{T}(eachcol(generate_sample(hppp; rng=rng)))
 
     M = Float64[]  # Marking process
     R = T[]        # Recording process
 
     steps = -1:-1:-n₀
     while true
-        backward_update!(D, M, R, steps, birth_rate, window_; rng=rng)
+        backward_extend!(D, M, R, steps, birth_rate, window_; rng=rng)
         coupling, L = forward_coupling(D, M, R, pp, β)
         coupling && return collect(L)
         steps = (steps.stop-1):-1:(2*steps.stop)
     end
 end
 
-function backward_update!(
+"""
+    backward_extend!(
         D::Set{T},          # Dominating process
         M::Vector{Float64}, # Marking process
         R::Vector{T},       # Recording process
         steps::StepRange,   # Number of backward steps
         birth_rate::Real,
         win::AbstractWindow;
         rng=-1
+    ) where {T}
+
+Sample from the dominating birth-death process backwards in time according to `steps`.
+The marks `M` and the points `R` which were added (uniform mark) / deleted (mark=0) along the run are recorded (`pushfirst!`).
+The final state of the dominating process is updated to `D`.
+
+**See also**
+
+- SBDevolve() [KeMo99](@cite)
+"""
+function backward_extend!(
+    D::Set{T},          # Dominating process
+    M::Vector{Float64}, # Marking process
+    R::Vector{T},       # Recording process
+    steps::StepRange,   # Number of backward steps
+    birth_rate::Real,
+    win::AbstractWindow;
+    rng=-1
 ) where {T}
     rng = getRNG(rng)
     for _ in steps
@@ -66,6 +127,21 @@ function backward_update!(
     end
 end
 
+"""
+    forward_coupling(
+        D::Set{T},           # Dominating process
+        M::Vector{Float64},  # Marking process
+        R::Vector{T},        # Recording process
+        pp::AbstractSpatialPointProcess{T},
+        β::Real              # Upper bound on papangelou conditional intensity
+    ) where {T}
+
+Check if coalescence occured between the lower and upper bounding processes and return the state of the lower bounding process at time 0.
+
+**See also**
+
+- SBDadd() [KeMo99](@cite)
+"""
 function forward_coupling(
         D::Set{T},           # Dominating process
         M::Vector{Float64},  # Marking process