Merge pull request #62 from adolgert/docs/read-through-text

Docs/read through text

docs/make.jl

@@ -54,7 +54,7 @@ for (source, target) in adliterate
             chunk_source = joinpath(example_base, chunk)
             chunk_target = joinpath(example_base, "literate", chunk)
             @info "Moving $chunk_source to $chunk_target"
-            mv(chunk_source, chunk_target)
+            mv(chunk_source, chunk_target; force=true)
         end
     end
 end
@@ -91,7 +91,11 @@ makedocs(;
             "SIR Model" => "sir.md",
             "Reliability" => "reliability.md"
         ],
-        "Reference" => "reference.md"
+        "Reference" => [
+            "interface.md",
+            "reference.md",
+            "algorithms.md"
+        ]
     ],
 )
 

docs/src/algorithms.md

@@ -0,0 +1,18 @@
+```@meta
+CurrentModule = Fleck
+```
+
+# Algorithms
+
+Many samplers depend on data structures to allow efficient querying of clocks ordered with respect to some value, usually the firing time. These types and methods implement them for Fleck.
+
+```@docs
+CumSumPrefixSearch
+BinaryTreePrefixSearch
+KeyedKeepPrefixSearch
+KeyedRemovalPrefixSearch
+choose
+setindex!
+rand
+set_multiple!
+```

docs/src/distributions.jl

@@ -7,13 +7,13 @@ using LaTeXStrings #hide
 using QuadGK #hide
 using ForwardDiff #hide
 
-# Fleck simulates generalized semi-Markov processes (GSMP). Every event in a generalized semi-Markov process is chosen as the result of a competion among clocks to see which fires next.
+# Fleck is a sampler for generalized semi-Markov processes (GSMP). Every event in a generalized semi-Markov process is chosen as the result of a competion among clocks to see which fires next.
 # 
 # In a *process-oriented* simulation (like [SimJulia](https://simjuliajl.readthedocs.io/en/stable/welcome.html)), control flow is based on tasks. Each task performs some action on the state, rolls the dice, and sets a wake-up time. It might wake up as expected and possible execute code, or it might be interrupted by another task's actions. In contrast, an *event-oriented* simulation using Fleck will create a set of possible next events, assign a probability distribution for *when* each can happen, and the timing of which happens first determines *which* next event happens. Let's look at how a probability distribution describes the time for an event to happen and then how they compete in Fleck.
 # 
 # ## Distributions in Time
 # 
-# Let's say you have a cold. You know you aren't going to recover immediately, but, as days go by, you're more and more sure you'll recover soon. You can read this graph below as the probability to recover given that you have not yet recovered.
+# Let's say you have a cold. You know you aren't going to recover immediately, but, as days go by, you're more and more sure you'll recover soon. This graph below shows recovery as a *hazard rate, which is the probability, per unit time, given that the event has not yet happened.*
 
 const Gamma51 = Gamma(5.0, 1.0) #hide
 survival = ccdf #hide
@@ -29,10 +29,10 @@ title!(p, "Hazard of a Gamma Distribution") #hide
 DisplayAs.PNG(p) #hide
 
 # 
-# It starts at zero, meaning there's no way you'll recover as soon as you're sick. It gets more likely over time that you're at the tail end of being sick. This is called a *hazard rate,* and the hazard rate shown is that of a Gamma distribution, commonly used to describe the rate of recovery for a population of individuals who are sick.
+# This hazard rate starts at zero, meaning there's no way you'll recover when you're first sick. It gets more likely over time that you're at the tail end of being sick. The hazard rate shown is that of a Gamma distribution, commonly used to describe the rate of recovery for a population of individuals who are sick.
 # 
-# If, instead, you want to see the number of people who recover on any given day, that is called a probability distribution function.
-
+# If, instead, you want to see the number of people who recover on any given day, that is called a probability distribution function (pdf), which is a much more common way to display a distribution in time.
+#
 x = 0:0.01:15 #hide
 y = pdf.(Gamma51, x) #hide
 p = plot(x, y, label="Gamma(5,1)") #hide
@@ -42,7 +42,7 @@ title!(p, "PDF of a Gamma Distribution") #hide
 DisplayAs.PNG(p) #hide
 
 # 
-# Where the hazard rate is an instantaneous quantity at a point in time, the probability distribution function (pdf) integrates over all possible future times. If we call the hazard rate $\lambda(t)$ and call the pdf $f$, we get this relationship.
+# Where the hazard rate is an instantaneous quantity at a point in time, the probability distribution function (pdf) integrates over all possible future times. If we call the hazard rate $\lambda(t)$ and call the pdf $f(t)$, we get this relationship.
 # 
 # ```math
 # f(t) = \lambda(t) e^{-\int_0^t \lambda(s)ds}
@@ -55,10 +55,10 @@ DisplayAs.PNG(p) #hide
 t=5.0 #hide
 x1 = 5:0.01:15 #hide
 y1 = shiftpdf.(Gamma51, t, x1) #hide
-sgp = plot(x1, y1, label="Shifted Gamma(5,1)") #hide
+sgp = plot(x1, y1, label="New PDF for Restart") #hide
 x = 0:0.01:15 #hide
 y = pdf.(Gamma51, x) #hide
-plot!(sgp, x, y, label="Gamma(5,1)") #hide
+plot!(sgp, x, y, label="Original PDF") #hide
 vline!(sgp, [t], label="New start time") #hide
 xlabel!(sgp, "Time") #hide
 ylabel!(sgp, "Probability Density Function") #hide
@@ -72,7 +72,7 @@ DisplayAs.PNG(sgp) #hide
 # f(t;t>t_0) = \lambda(t) e^{-\int_{t_0}^t \lambda(s)ds}
 # ```
 # 
-# In some sense, the underlying nature of when a particular event will happen is described more by the hazard rate, whereas the distribution function tells us about ensembles of events.
+# The hazard rate describes a flow of probability, whereas the distribution function tells us about ensembles of events.
 # 
 # The hazard rate is related to the well known cumulative distribution function (CDF) by an integral. The CDF tells us what is the overall probability the event occured some time in the interval ``[t_0,t_1)``.
 # 
@@ -83,7 +83,7 @@ DisplayAs.PNG(sgp) #hide
 # Equally important for simulation is the survival function (sometimes called the complementary cumulative distribution function), which is the probability the event will not occur until after ``t_1``.
 # 
 # ```math
-# S(t_1) = 1 - F(t_1;t_1) = e^{-\int_{t_0}^{t_1} \lambda(s) ds}
+# S(t_1) = 1 - F(t_1;0) = e^{-\int_{0}^{t_1} \lambda(s) ds}
 # ```
 
 y1 = cdf.(Gamma51, x) #hide
@@ -95,6 +95,8 @@ ylabel!(p, "Probability") #hide
 title!(p, "Cumulative distribution and Survival functions \nof a Gamma Distribution") #hide
 DisplayAs.PNG(p) #hide
 
+# For our example, survival is the chance the cold lasts longer than the given time.
+
 # ## Competition
 # 
 # ### Individual Distributions
@@ -128,15 +130,15 @@ DisplayAs.PNG(indp) #hide
 # 
 # The separate hazard rates are what we put into the simulation. Given their competition, the hazard rates will remain unchanged, but the pdfs will change.
 # 
-# ### Total Probability
+# ### Marginal Probability
 # 
 # Each of the three clock distributions above corresponds to a unique event ``E_i``, which has a probability that it will be the first to fire. We calculate this probability by marginalizing over the other events, which ends up being an integral over the distribution, multiplied by the survivals of the other events.
 # 
 # ```math
 # P[E_i] = \int_0^\infty f_i(t) \prod_{j\ne i} S_j(t) dt
 # ```
 # 
-# That gives the graph on the left.
+# That gives the chart on the left, where the sum of all ``P[E_i]`` is one.
 # 
 function marginal_probability_event(dists, which, t) #hide
     core_matrix = pdf(dists[which], t) #hide
@@ -147,7 +149,7 @@ function marginal_probability_event(dists, which, t) #hide
     end #hide
     core_matrix #hide
 end #hide
-
+#hide
 total_probability(which, lim) = quadgk(t -> marginal_probability_event(com_dists, which, t), 0, lim, rtol=1e-3)[1] #hide
 tp = total_probability.([1, 2, 3], Inf) #hide
 @show tp #hide
@@ -162,19 +164,27 @@ bp = plot(com_x, com_pdfs, labels=com_labels, title="PDFs of Competing Clocks")
 condp = plot(bb, bp, size=(1000, 400)) #hide
 DisplayAs.PNG(condp) #hide
 # 
-# The graph on the right shows the conditional distribution in time for each event, given that it was the one that fired, so it is $P[t_i | E_i]$. You can see that these distributions don't match the distributions for the individual events. They are modified by competition.
+# The graph on the right shows the conditional distribution in time for each event, given that it was the one that fired, so it is $P[t_i | E_i]$. In the language of semi-Markov processes, these distributions are called *holding times.* You can see that these distributions don't match the distributions for the individual events. They are modified by competition.
 # 
 # 
-# ### Total Time
+# ### Marginal Time
 # 
+# What if we split the marginal and conditional the other way? Instead of marginalizing the probability of which event fires, start with a marginal of the probability for when any event fires. One way to calculate this is to say that the hazard rate for any event to fire is the sum of the hazard rates.
+# ```math
+# \lambda_m(t) = \sum_i\lambda_i(t)
+# ```
+# From here, we know the pdf for the first firing time.
+#
+# ```math
+# f(t) = \lambda_m(t)\exp\left(-\int_0^t\lambda_m(s)ds\right)
+# ```
+#
 # From the graph above, if we pick a time, $t_1=10$, we can read from the graph three hazard rates, $(\lambda_1(t_1),\lambda_2(t_1), \lambda_3(t_1))$. Each hazard rate is the rate, per unit time, of that event. We know that, if the simulation makes it to $t=10$ without any event happening, the conditional probability for any one of those events is the ratio of hazard rates.
 # 
 # ```math
 # P[E_1|t=t_1] = \frac{\lambda_1(t_1)}{\lambda_1(t_1)+\lambda_2(t_1)+\lambda_3(t_1)}
 # ```
-# 
-# If all distributions are Exponential, then they all have a constant hazard rate, and that conditional probability stays the same over time. If any distributions are non-Exponential, then that conditional probability changes.
-# 
+# Now we can plot, on the left, the pdf for who fires first and, on the right, the probability of which event fires, given the firing time.
 function total_pdf(t) #hide
     factor1 = zero(Float64) #hide
     for dist in com_dists #hide
@@ -187,22 +197,22 @@ function total_pdf(t) #hide
     return factor1 * factor2 #hide
 end #hide
 a = plot(com_x, total_pdf.(com_x), legend=:none, title="PDF for First to Fire") #hide
-
+#hide
 function com_cond_prob(i, t) #hide
     hazard(com_dists[i], t) / sum(hazard.(com_dists, t)) #hide
 end #hide
 com_rel_prob = stack([com_cond_prob.(i, com_x) for i in eachindex(com_dists)]) #hide
 b = plot(com_x, com_rel_prob, labels=com_labels, ylabel="Conditional Probability", #hide
     title="Which Event Fires Given Firing Time") #hide
 DisplayAs.PNG(plot(a, b, size=(1000, 400))) #hide
-
 # 
 # On the left of this graph is the pdf for the first event of the three to fire. We can see this as a marginal $P[t]$ and then the right-hand graph as the conditional $P[E_i|t]$.
 # 
 # ## Specification of a Simulation
-# 
-# We simulated using competing clocks when we measure using survival analysis. One is the inverse of the other.
-# 
-# A classic example of survival analysis is a hospital study where participants enter a trial for a drug, and doctors observe the date at which each patient recovers. How would you draw a probability distribution function for such a trial? You would draw a histogram where each week shows a count of how many recovered. What about the hazard rate for such a trial? Here, you would look only at the number of participants remaining in the trial and ask what fraction recover each week. That's a rate of recovery given that participants have not yet recovered.
-# 
-# There are many disciplines that use survival analysis to measure rates of processes. Epidemiology measures rates of contagion, death, and recovery. Demography measures time to have children or move to another country. Reliability studies the time for a part to break and be fixed. Drug discovery measures the time to approve a chemical for the next trial stage. For traffic flow, it's the time to the next intersection. For package delivery, it's the time to the next package scan. All of these do use, or can use, survival analysis to observe a system to simulate.
+#
+# If we imagine a drug trial, where patients can recover, die, or exit the trial for some other reason, there are three mutually-exclusive events, like the example above. If we pick the recovery event and plot its distribution in time, which of the above plots will we see? This will be a holding time. It won't be the pdf that represents the rate of recovery in the absence of competing events. However, given observations of competing events, it is possible to calculate back to the original hazard rates using survival analysis.
+#
+# [Survival analysis](https://en.wikipedia.org/wiki/Survival_analysis) uses observations of event times and event cancellations to estimate hazard rates for each event. It helps you tease apart the effects of competition to see the underlying probability per unit time that any event would fire, given that it has not yet fired.
+#
+# Simulation is the opposite of survival analysis. It allows you to take rules for how any event would fire, in the absence of competition, and to place it in a more complicated environment where competition happens. When you specify a continuous-time simulation, it isn't specified with the pdfs of holding times but with the pdfs of rates derived from survival analysis.
+

docs/src/guide.md

@@ -4,10 +4,6 @@ A discrete event simulation tracks the effects of individual events in time. We
 
 We write discrete event simulations by defining the state of a system and events that can act on that state. The state of a system can be the count of each chemical species for a chemical reaction system. It can be different numbers of chickens at different weights in different chicken coops for a farming simulation. It can be the age of each widget in an engine for a reliability simulation. These states are called the *physical state* of the system and vary widely in what they represent.
 
-Even though nearly every simulation represents a unique system, it is often possible to find enough in common among systems that we can make libraries to help create, simulate, and observe changes of state. For instance, a large class of physical states are counts of things, whether they are molecules, vehicles, or chickens. These states can be represented in code as vectors of integers, and there is a long history of models represented by vectors of integers. Most of this work focuses on using chemical simulation codes to represent not only chemicals but also other systems like vehicles. An older, more general, form of this simulation is the [vector addition system](https://en.wikipedia.org/wiki/Vector_addition_system). However, states don't need to be restricted to vectors of integers.
-
-We usually want to simulate systems that have real-valued observations such as ages, temperatures, and distances. These continuous parts of the state of the system mean that the system rarely returns to exactly the same state. There are now an infinite number of possible states of the system, which changes how we would think of the statistics of the system, but it doesn't much change how we handle simulation, because simulation is still driven by events.
-
 An *event* is an observable change of state at a time. It is a function where the argument is the current state of the system and the return value is the new state of the system. We say an event occurs, happens, or fires when we apply that function to the state. How that state changes is up to the model and how it defines the state, but what about *when* that state changes?
 
 Given the current state of the system, there is a set of possible next events which could happen. We call the possible next events *enabled.* Think of the moment in time just after the simulation has begun or an event has just happened. At this moment, there may be multiple enabled events. Which one is next depends on how often that event happens. If our simulation includes both the radioactive decay of an unstable element and decay of iron, the unstable element will usually decay first. We describe the rates of these events using probability distributions in time. Each event has an associated continuous univariate distribution where the variable is time.
@@ -26,6 +22,8 @@ All these sets are visualized in the figure below, reproduced from "Stochastic p
 
 ![](assets/ClockUpdate.png)
 
+Fleck is responsible for knowing ``E^{*}``, but it is the main event loop of the simulation using Fleck that decides which clocks to disable and enable.
+
 ## References
 
 [Haas:2006]	PJ. Haas, “[Stochastic petri nets: Modelling, stability, simulation](https://link.springer.com/book/10.1007/b97265),” Springer Science & Business Media, 2006.

docs/src/interface.md

@@ -0,0 +1,25 @@
+```@meta
+CurrentModule = Fleck
+```
+
+# Interface
+
+These are methods which are defined for any samplers subtyping `<:SSA`, the abstract sampler type.
+
+## Use a Sampler
+
+```@docs
+enable!
+disable!
+next
+reset!
+```
+
+## Query a Sampler
+
+```@docs
+getindex
+keys
+keytype
+length
+```

docs/src/reference.md

@@ -2,57 +2,26 @@
 CurrentModule = Fleck
 ```
 
-# Reference
-
- * [Interface](@ref)
- * [Samplers](@ref)
- * [Algorithms](@ref)
-
-
-## Interface
-
-These are methods which are defined for any samplers subtyping `<:SSA`, the abstract sampler type.
-
-### Use a Sampler
-
-```@docs
-enable!
-disable!
-next
-reset!
-```
-
-### Query a Sampler
-
-```@docs
-getindex
-keys
-keytype
-length
-```
-
-## Samplers
+# Samplers
 
 The choice of sampler determines specific algorithms that are used to sample, update, and disable clocks. Helpers also exist that are useful for logging, utilizing common random numbers, and hierarchical sampling.
 
-### Sampler Supertype
+## Sampler Supertype
 
 ```@docs
 SSA
 ```
 
-### Sampler Types
+## Sampler Types
 
 ```@docs
 FirstReaction
 FirstToFire
 DirectCall
 CombinedNextReaction
-consume_survival
-sampling_space
 ```
 
-### Sampling Helpers
+## Sampling Helpers
 
 ```@docs
 CommonRandomRecorder
@@ -64,19 +33,6 @@ SingleSampler
 ChatReaction
 DebugWatcher
 TrackWatcher
-```
-
-## Algorithms
-
-Many samplers depend on data structures to allow efficient querying of clocks ordered with respect to some value, usually the firing time. These types and methods implement them for Fleck.
-
-```@docs
-CumSumPrefixSearch
-BinaryTreePrefixSearch
-KeyedKeepPrefixSearch
-KeyedRemovalPrefixSearch
-choose
-setindex!
-rand
-set_multiple!
+consume_survival
+sampling_space
 ```