modified the debugging and context to make documentation more complete so it's a good outcome

Author: adolgert

docs/make.jl

@@ -111,12 +111,14 @@ makedocs(;
             "commonrandom.md",
             "importance_skills.md",
             "hamiltonianmontecarlo.md",
-            "gen/overview.md",
-            "gen/distribution.md",
-            "gen/generative_function.md",
-            "gen/observation_likelihood.md",
-            "gen/importance_mixture.md",
-            "gen/hmc_paths.md",
+            "Gen.jl Integration" => [
+                "gen/overview.md",
+                "gen/distribution.md",
+                "gen/generative_function.md",
+                "gen/observation_likelihood.md",
+                "gen/importance_mixture.md",
+                "gen/hmc_paths.md",
+            ],
             "gen/turing_dist.md",
             "gen/survival_snippet.md",
         ],

docs/src/choosing_sampler.md

@@ -2,16 +2,6 @@
 
 How to make a sampler for events in CompetingClocks.jl.
 
-## Quick Reference
-
-| Name           | Distribution support | Best use cases                 |
-|----------------|----------------------|--------------------------------|
-| First-to-fire  | All                  | Fastest for simple simulation  |
-| Next Reaction  | All                  | Best for common random numbers |
-| First Reaction | All                  | Fastest for very few events    |
-| Direct         | Exponential-only     | Quicker likelihood calculation |
-| Petri          | All                  | Debug rare errors              |
-
 ## High-level Interface
 
 Every event in a simulation has an identifying key which can be any immutable Julia
@@ -26,6 +16,7 @@ The sampler takes a `Random.AbstractRNG` as its second argument.
 Specify features for the sampler with keyword arguments.
 These features determine the type of the sampler.
 
+ * `method=nothing`---Setting this chooses a single sampling algorithm for all clocks.
  * `step_likelihood=false`---Setting this to `true` let's you calculate the
    the likelihood of the next event and time before firing it.
  * `path_likelihood=false`---Set to true in order to calculate likelihood of a
@@ -41,9 +32,32 @@ These features determine the type of the sampler.
  * `recording=false`---This will create a vector of every enable and disable
    event for test and debug.
 
+If you want to try a particular sampler instead of having the `SamplingContext`
+pick one for you, use the `method` parameter.
+
+```julia
+sampler = SamplingContext(KeyType, Float64, rng; method=DirectMethod(:keep, :tree))
+```
+
+## Quick Reference
+
+These are all of the samplers.
+
+| Name           | Distribution support | Best use cases                 |
+|----------------|----------------------|--------------------------------|
+| [`FirstToFireMethod`](@ref)  | All                  | Fastest for simple simulation  |
+| [`NextReactionMethod`](@ref)  | All                  | Best for common random numbers |
+| [`FirstReactionMethod`](@ref) | All                  | Fastest for very few events    |
+| [`DirectMethod`](@ref)         | Exponential-only     | Quicker likelihood calculation |
+| [`RejectionMethod`](@ref)         | Exponential-only     | RSSA for large systems |
+| [`PartialPropensityMethod`](@ref)         | Exponential-only     | For reaction networks. |
+| [`PetriMethod`](@ref)          | All                  | Debug rare errors              |
+
+
 ## Hierarchical Samplers
 
-For spatial simulations or for chemical simulations that aren't well-mixed,
+For spatial simulations, chemical simulations that aren't well-mixed, or
+simulations with some fast and many slow clocks,
 it can help to split a sampler into multiple buckets so that each bucket can
 update its own list of what fires next. This package does this by adding
 sampler groups.
@@ -55,7 +69,11 @@ A sampler group is
    whether a given event belongs to this sampler.
  * An optional `method` to say what kind of sampler this group should use.
 
+This example sends events with keys like `(:infect, 3)` tto the sampler called
+`:forthright` and sends keys like `(:recover, 7)` to the sampler called
+`:sparky`.
 ```julia
+const KeyType = Tuple{Symbol,Integer}
 builder = SamplerBuilder(KeyType, Float64)
 add_group!(builder, :sparky => (x,d) -> x[1] == :recover, method=NextReaction())
 add_group!(builder, :forthright=>(x,d) -> x[1] == :infect)

docs/src/contextinterface.md

@@ -41,6 +41,10 @@ CompetingClocks.clone(::SamplingContext{K,T,Sampler,RNG,Like,CRN,Dbg}, rng::RNG)
 CompetingClocks.freeze_crn!(::SamplingContext)
 CompetingClocks.reset_crn!(::SamplingContext)
 CompetingClocks.enabled(::SamplingContext)
+CompetingClocks.enabled_history(::SamplingContext)
+CompetingClocks.disabled_history(::SamplingContext)
+CompetingClocks.EnablingEntry
+CompetingClocks.DisablingEntry
 Base.length(::SamplingContext)
 CompetingClocks.isenabled(::SamplingContext{K}, ::K) where {K}
 Base.keytype(::SamplingContext{K}) where {K}

docs/src/debugging.md

@@ -27,3 +27,14 @@ The `DebugWatcher` records every time an event is enabled or disabled.
 You wouldn't normally want to use all of this time and memory (and memory churn),
 but it can be helpful to see the trace of all events enabled and disabled,
 in addition to those that fired.
+
+```julia
+for enabling_event in enabled_history(sampler)
+    println("$(enabling_event.clock), $(enabling_event.when)")
+end
+for disabling_event in disabled_history(sampler)
+    println("$(disabling_event.clock), $(disabling_event.when)")
+end
+```
+
+Look at [`CompetingClocks.enabled_history`](@ref) and [`CompetingClocks.disabled_history`](@ref).

docs/src/index.md

@@ -37,23 +37,21 @@ or to calculate the likelihood of a sample path for statistical estimation.
 
 ## Usage
 
-The library provides you with samplers. Each sampler has the same interface. Here, a distribution is a [Distributions.ContinuousUnivariateDistribution](https://juliastats.org/Distributions.jl/stable/univariate/#Continuous-Distributions), `RNG` is a [random number generator](https://docs.julialang.org/en/v1/stdlib/Random/#Generators-(creation-and-seeding)), the `key` is some identifier (maybe an integer) for the event, and an enabling time is a zero-time for the given distribution.
+The library provides you with samplers. Each sampler has the same interface. Here, a distribution is a [Distributions.ContinuousUnivariateDistribution](https://juliastats.org/Distributions.jl/stable/univariate/#Continuous-Distributions), `RNG` is a [random number generator](https://docs.julialang.org/en/v1/stdlib/Random/#Generators-(creation-and-seeding)), the `key` is some identifier (maybe an integer, tuple, or pair) for the event, and an enabling time is a zero-time for the given distribution.
 
  * `enable!(sampler, key, distribution, enabling_time)` - to start the clock on when an event will fire next.
 
  * `disable!(sampler, key)` - to turn off an event so it can't fire.
 
- * `next(sampler)` - to ask this library who could fire next.
+ * `next(sampler)` - to ask this library who could fire next, and it return a (time, key).
 
  * `fire!(sampler, key, time)` - choose which event happens at what time.
 
-Different samplers are specialized for sampling more quickly and accurately for different applications. For instance, some applications have very few events enabled at once, while some have many. Some applications use only exponentially-distributed events, while some have a mix of distribution types. Because continuous-time discrete event systems can fire many events, the literature has focused on reducing the number of CPU instructions required to sample each event, and this library reflects that focus.
-
 ## When NOT to use Competing Clocks
 
- * **Pure exponential distributions?** JumpProcesses.jl is a complete framework for this.
- * **Need ODE coupling?** Again, it's easier to stay within SciML and JumpProcesses.jl.
- * **Want high-level frameworks?** Try Agents.jl.
+ * **Pure exponential distributions and chemical systems?** JumpProcesses.jl is a complete framework for this.
+ * **Need ODE coupling?** Again, it's easier to stay within SciML and JumpProcesses.jl. CompetingClocks.jl supports ODEs as Dirac distributions.
+ * **Want high-level frameworks?** Try [Agents.jl](https://juliadynamics.github.io/Agents.jl/stable/) or [ConcurrentSim.jl](https://juliadynamics.github.io/ConcurrentSim.jl/stable/).
 
 CompetingClocks.jl is for:
 

docs/src/integration-guide.md

@@ -49,7 +49,8 @@ performant for sampling.
 This package doesn't use Tasks or thread guards because it is so common to
 run multiple copies of a simulation, each on its own thread. If you run
 multiple simulations, it can be more correct to use random number generators
-that work in parallel, such as those in the `Random123` package.
+that work in parallel, such as those in the `Random123` package or the `jump!()`
+method for the `Xoshiro` sampler in the `Random` package.
 
 ### State and Dependencies
 
@@ -58,14 +59,14 @@ events. It can be very helpful for performance if you use information about
 relationships among events to limit updates to the samplers. There are a few
 kinds of dependency graphs that may help.
 
- * When one event fires, it changes state. Which other events could be, or are,
-   enabled because of the changed state?
- * When one event fires, it changes state. Which other events have transition
-   rates that depend on that state? The event may remain enabled, but its
+ * Dependency Graph---When one event fires, it changes state. Which other events could be
+   enabled or disabled because of the changed state?
+ * Rate Dependency Graph---When one event fires, it changes state. Which other events have transition
+   *rates* that depend on that state? The event may remain enabled, but its
    distribution of firing times may have a different parameter value.
- * When one event fires, you might know directly which other events have
+ * Event-to-event Graph---When one event fires, you might know directly which other events have
    enabling rules or transition rates that depend on that first event, without
-   needing to ask what state changed.
+   needing to ask what state changed and what events depend on the state.
 
 How you represent a dependency graph depends on the problem.
 
@@ -94,16 +95,16 @@ Split the sampler into groups so that each region is faster.
 ```julia
 builder = SamplerBuilder(EventKey, Float64)
 for region in spatial_grid
-    add_group!(builder, regin.name => (key, dist) -> key.location == region.id)
+    add_group!(builder, region.name => (key, dist) -> key.location == region.id)
 end
 ```
 
 ### Chemical Reaction Networks
 Map reactions to events.
 ```julia
 function enable_reaction!(sampler, reaction, state, params)
-    rate = calculate_propensity(reaction, state, params)
-    enable!(sampler, reaction.id, Exponential(1/rate))
+    propensity = calculate_propensity(reaction, state, params)
+    enable!(sampler, reaction.id, Exponential(1/propensity))
 end
 ```
 
@@ -123,17 +124,22 @@ sampler = SamplingContext(builder, rng)
 # After simulation, the log_prob is a Float64.
 log_prob = pathloglikelihood(sampler, end_time)
 ```
+The resulting value is the likelihood of a single trajectory of the system.
+The `pathloglikelihood` is a regular function and can be differentiated.
 
 ## Variance reduction
 Construct the sampler with the `common_random` option.
 ```julia
 builder = SamplerBuilder(KeyType, Float64; common_random=true)
 sampler = SamplingContext(builder, rng)
 ```
-Run it a bunch of times. Each run will pin more random draws. Then use
+Run the simulation a bunch of times, calling `reset!()` between runs to
+clear the sampler's memory of the previous run but save the common random
+numbers. Each run will pin more random draws. Then use
 `freeze_crn!` to stop collecting random draws.
 ```julia
 for i in 1:warm_up
+    reset!(sampler)
     results1 = run_simulation(model, sampler)
 end
 freeze_crn!(sampler)  # Lock random draws
@@ -173,7 +179,7 @@ of each copy by one over the number of copies.
 
 CompetingClocks offers a `split!` function that copies the state to a vector
 of clocks and keeps a `split_weight` as a record of how many times splitting
-was done. Take a look at the implementation of `split!()`.
+was done. Take a look at the implementation of [`CompetingClocks.split!()`](@ref).
 
 Your code might keep the sampler in the main simulation framework, in which case
 you would include the sampler in the cloning and copying process.
@@ -187,22 +193,21 @@ sampler = SamplingContext(builder, rng)
 ```
 Then run with logging on.
 ```julia
-with_logger(ConsoleLogger(stdout, Logging.Info)) do
+with_logger(ConsoleLogger(stdout, Logging.Debug)) do
     observed, importance = run_epochs(100)
 end
 ```
 
-There is also a chaos-monkey option called a `Petri` sampler.
+There is also a [chaos-monkey](https://en.wikipedia.org/wiki/Chaos_engineering) option called a `Petri` sampler.
 This sampler always advances time by the same time step `dt=1.0` by default,
 and it ignores the clock distributions. It chooses any enabled sampler at
 each step, with the goal of finding unusual sampling paths. This is surprisingly
 effective at finding logic errors in code to change state.
 
 ## Performance Considerations
 
-**Type Stability**---The ClockKey is a place where it can be very helpful to
-use different types to describe different kinds of events, but it's also a place
-where type stability makes the calculations faster.
+**Type Stability**---Use a concrete type for the
+ClockKey for more speed.
 
 **Sampler Choice**---If there are a lot of clocks enabled at the same time,
 the the sampler choice can be important. Experts in traffic simulation use
@@ -212,5 +217,3 @@ samplers about whether clock keys are reused.
 
 **Memory Allocation**---Instead of making a new sampler for each run, use
 `reset!(sampler)` between runs to clear out the clocks.
-
-In short, make a sampler choice and benchmark.

docs/src/low_level_interface.md

@@ -2,7 +2,10 @@
 # Low-level Sampler Interface
 
 The main interface to the package is the [`SamplingContext`](@ref).
-Within that class are the samplers themselves. This describes their interface.
+The `SamplingContext` contains sampling algorithms and other helper classes. This describes the interface to underlying sampling algorithms, which differ slightly:
+
+ * They don't hold a random number generator internally so it is passed as an argument.
+ * All input times are absolute, so the enabling time of a distribution is given in absolute time.
 
 ## Common Interface to Low-level Samplers
 
@@ -45,8 +48,8 @@ CompetingClocks.DirectCall{K,T,P}
 CompetingClocks.FirstReaction{K,T}
 CompetingClocks.FirstToFire{K,T}
 CompetingClocks.PSSACR{K,T}
-CompetingClocks.Petri{K,T}
 CompetingClocks.RSSA{K,T}
+CompetingClocks.Petri{K,T}
 ```
 
 ## Hierarchical Sampling

docs/src/quickstart.md

@@ -1,4 +1,4 @@
-# CompetingClocks.jl Quickstart
+# Quickstart
 
 Learn CompetingClocks.jl basics by building a minimal simulation.
 
@@ -20,8 +20,7 @@ population = Set(1:person_cnt)
 # Events will be identified by a tuple of type and individual.
 const ClockKey = Tuple{Symbol,Int}
 # Create a sampler.
-builder = SamplerBuilder(ClockKey, Float64)
-sampler = SamplingContext(builder, rng)
+sampler = SamplingContext(ClockKey, Float64, rng)
 enable!(sampler, (:birth, 0), Exponential(inv(length(population) * 2.0)))
 for mort in population
     enable!(sampler, (:death, mort), Gamma(inv(1.5)))
@@ -51,7 +50,7 @@ end
  3. Create a sampler.
  4. Enable initial events.
  5. Ask what happens next and
-    - `fire!` that event.
+    - `fire!` that event, which also disables the event.
     - Handle changes to state and enabled events.
 
 You tell it what *could* happen with `enable!()`, ask what happens `next()`,
@@ -62,20 +61,23 @@ So that you can stop a simulation at a fixed time rather than after an event.
 
 Learn More Techniques
 
- - State management: your-first-sim.md - Better state handling.
- - Choosing samplers: choosign-samplers.md - When to use which algorithm.
+ - State management: [Sample Main Loop](@ref "Sample Main Loop") - Better state handling.
+ - Choosing samplers: [Choosing a Sampler](@ref "Choosing a Sampler") - When to use which algorithm.
  - Real examples:
-   - SIR
-   - reliability
+   - [SIR Model](sir.md)
+   - [Reliability](reliability.md)
+   - [Birth-Death Process](constant_birth.md)
+   - [Gene Expression](gene_expression.md)
 
 Advanced Features
 
- - Get likelihoods: likelihood-tutorial
- - Variance reduction: variance-reduction.md
- - Build a framework: integration-guide.md
+ - [Get likelihoods](@ref "Importance Sampling for Simulation")
+ - [Variance reduction](@ref "Common Random Numbers")
+ - [Build a framework](@ref "Integration Guide")
 
 API Reference
 
- - api.md#enable - Full parameter documentation
- - api.md#next - Return value details
- - api.md#samplerbuilder - Configuration options
+ - [Build a Context](@ref "Building a Context")
+ - [`CompetingClocks.enable!`](@ref) - Parameter documenation.
+ - [`CompetingClocks.next`](@ref) - Return value details.
+ - [`CompetingClocks.SamplerBuilder`](@ref) - Configuration options.

src/context.jl

@@ -1,5 +1,5 @@
 export SamplingContext, enable!, fire!, isenabled, freeze_crn!
-export sample_from_distribution!
+export sample_from_distribution!, enabled_history, disabled_history
 
 
 mutable struct SamplingContext{K,T,Sampler<:SSA{K,T},RNG,Like,CRN,Dbg}
@@ -78,7 +78,8 @@ We keep the internal logic simple. If something is present, call it.
  * `step_likelihood` - whether you will call `steploglikelihood` before each `fire!`
  * `path_likelihood` - whether you will call `pathloglikelihood`
     at the end of a simulation run.
- * `debug` - Print log messages at the debug level.
+ * `debug` - Print log messages at the debug level. Enabling this
+   stores every enabling and disabling event so don't leave it on.
  * `recording` - Store every enable and disable for later examination.
  * `common_random` - Use common random numbers during sampling.
  * `method` - If you want a single, particular sampler, put its `SamplerSpec` here.
@@ -412,6 +413,36 @@ function enabled(ctx::SamplingContext)
 end
 
 
+"""
+    enabled_history(ctx::SamplingContext)
+
+Returns a `Vector{EnablingEntry{K,T}}` that has every time a clock was enabled.
+
+See [`CompetingClocks.EnablingEntry`](@ref).
+"""
+function enabled_history(ctx::SamplingContext)
+    if !isnothing(ctx.debug)
+        return enabled_history(ctx.debug)
+    else
+        error("In order to get history of enabling create context with `recording=true`")
+    end
+end
+
+"""
+    disabled_history(ctx::SamplingContext)
+
+Returns a `Vector{DisablingEntry{K,T}}` that has every time a clock was enabled.
+
+See [`CompetingClocks.DisablingEntry`](@ref).
+"""
+function disabled_history(ctx::SamplingContext)
+    if !isnothing(ctx.debug)
+        return disabled_history(ctx.debug)
+    else
+        error("In order to get history of disabling create context with `recording=true`")
+    end
+end
+
 """
     length(sampling)