Skip to content

Commit

Permalink
Fix gillespie (#5)
Browse files Browse the repository at this point in the history 
* update gillspie algorithm to use `sample`

* render
  • Loading branch information
@sbfnk
sbfnk authored Sep 12, 2024
1 parent 3dd22d8 commit 4940b94
Show file tree
Hide file tree
Showing 15 changed files with 90 additions and 68 deletions.
12 changes: 6 additions & 6 deletions 10_StochasticContinuous_practical.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,16 +54,16 @@ SIR_gillespie <- function(init_state, parms, tf) {
time <- time + rexp(n = 1, rate = sum(rates))
## check if next event is supposed to happen before end time
if (time <= tf) {
## generate cumulative sum of rates, to determine the type of the next
## event
cumulative_rates <- cumsum(rates)
## determine the type of the next
next_event <- sample(x = length(rates), size = 1, prob = rates)
## change to name
next_event <- names(rates)[next_event]
## determine type of next event
type <- runif(n = 1, min = 0, max = sum(rates))
if (type < cumulative_rates["infection"]) {
if (next_event == "infection") {
## infection
S <- S - 1
I <- I + 1
} else if (type < cumulative_rates["recovery"]) {
} else if (next_event == "recovery") {
## recovery
I <- I - 1
R <- R + 1
Expand Down
12 changes: 6 additions & 6 deletions 10_StochasticContinuous_solutions.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,16 +46,16 @@ SIR_gillespie <- function(init_state, parms, tf) {
time <- time + rexp(n = 1, rate = sum(rates))
## check if next event is supposed to happen before end time
if (time <= tf) {
## generate cumulative sum of rates, to determine the type of the next
## event
cumulative_rates <- cumsum(rates)
## determine the type of the next
next_event <- sample(x = length(rates), size = 1, prob = rates)
## change to name
next_event <- names(rates)[next_event]
## determine type of next event
type <- runif(n = 1, min = 0, max = sum(rates))
if (type < cumulative_rates["infection"]) {
if (next_event == "infection") {
## infection
S <- S - 1
I <- I + 1
} else if (type < cumulative_rates["recovery"]){
} else if (next_event == "recovery") {
## recovery
I <- I - 1
R <- R + 1
Expand Down
65 changes: 38 additions & 27 deletions docs/10_StochasticContinuous_practical.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>

<meta charset="utf-8">
<meta name="generator" content="quarto-1.4.555">
<meta name="generator" content="quarto-1.5.57">

<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">


<title>Modern Techniques in Modelling - 10. Stochastic continuous models</title>
<title>10. Stochastic continuous models – Modern Techniques in Modelling</title>
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
Expand All @@ -33,7 +33,7 @@
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
pre > code.sourceCode > span { display: inline-block; text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
Expand Down Expand Up @@ -101,7 +101,7 @@
}
}</script>

<script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
<script src="https://cdnjs.cloudflare.com/polyfill/v3/polyfill.min.js?features=es6"></script>
<script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml-full.js" type="text/javascript"></script>

<script type="text/javascript">
Expand Down Expand Up @@ -146,7 +146,7 @@
</a>
</div>
<div id="quarto-search" class="" title="Search"></div>
<button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
<button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" role="menu" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
<span class="navbar-toggler-icon"></span>
</button>
<div class="collapse navbar-collapse" id="navbarCollapse">
Expand All @@ -164,7 +164,7 @@
<span class="menu-text">Timetable</span></a>
</li>
<li class="nav-item dropdown ">
<a class="nav-link dropdown-toggle" href="#" id="nav-menu-practicals" role="button" data-bs-toggle="dropdown" aria-expanded="false">
<a class="nav-link dropdown-toggle" href="#" id="nav-menu-practicals" role="link" data-bs-toggle="dropdown" aria-expanded="false">
<span class="menu-text">Practicals</span>
</a>
<ul class="dropdown-menu" aria-labelledby="nav-menu-practicals">
Expand Down Expand Up @@ -204,7 +204,7 @@
</li>
</ul>
</div> <!-- /navcollapse -->
<div class="quarto-navbar-tools">
<div class="quarto-navbar-tools">
</div>
</div> <!-- /container-fluid -->
</nav>
Expand Down Expand Up @@ -299,16 +299,16 @@ <h3 class="anchored" data-anchor-id="setting-up-the-model-algorithm">Setting up
<span id="cb1-39"><a href="#cb1-39" aria-hidden="true" tabindex="-1"></a> time <span class="ot">&lt;-</span> time <span class="sc">+</span> <span class="fu">rexp</span>(<span class="at">n =</span> <span class="dv">1</span>, <span class="at">rate =</span> <span class="fu">sum</span>(rates))</span>
<span id="cb1-40"><a href="#cb1-40" aria-hidden="true" tabindex="-1"></a> <span class="do">## check if next event is supposed to happen before end time</span></span>
<span id="cb1-41"><a href="#cb1-41" aria-hidden="true" tabindex="-1"></a> <span class="cf">if</span> (time <span class="sc">&lt;=</span> tf) {</span>
<span id="cb1-42"><a href="#cb1-42" aria-hidden="true" tabindex="-1"></a> <span class="do">## generate cumulative sum of rates, to determine the type of the next</span></span>
<span id="cb1-43"><a href="#cb1-43" aria-hidden="true" tabindex="-1"></a> <span class="do">## event</span></span>
<span id="cb1-44"><a href="#cb1-44" aria-hidden="true" tabindex="-1"></a> cumulative_rates <span class="ot">&lt;-</span> <span class="fu">cumsum</span>(rates)</span>
<span id="cb1-45"><a href="#cb1-45" aria-hidden="true" tabindex="-1"></a> <span class="do">## determine type of next event</span></span>
<span id="cb1-46"><a href="#cb1-46" aria-hidden="true" tabindex="-1"></a> type <span class="ot">&lt;-</span> <span class="fu">runif</span>(<span class="at">n =</span> <span class="dv">1</span>, <span class="at">min =</span> <span class="dv">0</span>, <span class="at">max =</span> <span class="fu">sum</span>(rates))</span>
<span id="cb1-47"><a href="#cb1-47" aria-hidden="true" tabindex="-1"></a> <span class="cf">if</span> (type <span class="sc">&lt;</span> cumulative_rates[<span class="st">"infection"</span>]) {</span>
<span id="cb1-42"><a href="#cb1-42" aria-hidden="true" tabindex="-1"></a> <span class="do">## determine the type of the next</span></span>
<span id="cb1-43"><a href="#cb1-43" aria-hidden="true" tabindex="-1"></a> next_event <span class="ot">&lt;-</span> <span class="fu">sample</span>(<span class="at">x =</span> <span class="fu">length</span>(rates), <span class="at">size =</span> <span class="dv">1</span>, <span class="at">prob =</span> rates)</span>
<span id="cb1-44"><a href="#cb1-44" aria-hidden="true" tabindex="-1"></a> <span class="do">## change to name</span></span>
<span id="cb1-45"><a href="#cb1-45" aria-hidden="true" tabindex="-1"></a> next_event <span class="ot">&lt;-</span> <span class="fu">names</span>(rates)[next_event]</span>
<span id="cb1-46"><a href="#cb1-46" aria-hidden="true" tabindex="-1"></a> <span class="do">## determine type of next event</span></span>
<span id="cb1-47"><a href="#cb1-47" aria-hidden="true" tabindex="-1"></a> <span class="cf">if</span> (next_event <span class="sc">==</span> <span class="st">"infection"</span>) {</span>
<span id="cb1-48"><a href="#cb1-48" aria-hidden="true" tabindex="-1"></a> <span class="do">## infection</span></span>
<span id="cb1-49"><a href="#cb1-49" aria-hidden="true" tabindex="-1"></a> S <span class="ot">&lt;-</span> S <span class="sc">-</span> <span class="dv">1</span></span>
<span id="cb1-50"><a href="#cb1-50" aria-hidden="true" tabindex="-1"></a> I <span class="ot">&lt;-</span> I <span class="sc">+</span> <span class="dv">1</span></span>
<span id="cb1-51"><a href="#cb1-51" aria-hidden="true" tabindex="-1"></a> } <span class="cf">else</span> <span class="cf">if</span> (type <span class="sc">&lt;</span> cumulative_rates[<span class="st">"recovery"</span>]) {</span>
<span id="cb1-51"><a href="#cb1-51" aria-hidden="true" tabindex="-1"></a> } <span class="cf">else</span> <span class="cf">if</span> (next_event <span class="sc">==</span> <span class="st">"recovery"</span>) {</span>
<span id="cb1-52"><a href="#cb1-52" aria-hidden="true" tabindex="-1"></a> <span class="do">## recovery</span></span>
<span id="cb1-53"><a href="#cb1-53" aria-hidden="true" tabindex="-1"></a> I <span class="ot">&lt;-</span> I <span class="sc">-</span> <span class="dv">1</span></span>
<span id="cb1-54"><a href="#cb1-54" aria-hidden="true" tabindex="-1"></a> R <span class="ot">&lt;-</span> R <span class="sc">+</span> <span class="dv">1</span></span>
Expand Down Expand Up @@ -682,18 +682,7 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
}
return false;
}
const clipboard = new window.ClipboardJS('.code-copy-button', {
text: function(trigger) {
const codeEl = trigger.previousElementSibling.cloneNode(true);
for (const childEl of codeEl.children) {
if (isCodeAnnotation(childEl)) {
childEl.remove();
}
}
return codeEl.innerText;
}
});
clipboard.on('success', function(e) {
const onCopySuccess = function(e) {
// button target
const button = e.trigger;
// don't keep focus
Expand Down Expand Up @@ -725,15 +714,37 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
}, 1000);
// clear code selection
e.clearSelection();
}
const getTextToCopy = function(trigger) {
const codeEl = trigger.previousElementSibling.cloneNode(true);
for (const childEl of codeEl.children) {
if (isCodeAnnotation(childEl)) {
childEl.remove();
}
}
return codeEl.innerText;
}
const clipboard = new window.ClipboardJS('.code-copy-button:not([data-in-quarto-modal])', {
text: getTextToCopy
});
clipboard.on('success', onCopySuccess);
if (window.document.getElementById('quarto-embedded-source-code-modal')) {
// For code content inside modals, clipBoardJS needs to be initialized with a container option
// TODO: Check when it could be a function (https://github.com/zenorocha/clipboard.js/issues/860)
const clipboardModal = new window.ClipboardJS('.code-copy-button[data-in-quarto-modal]', {
text: getTextToCopy,
container: window.document.getElementById('quarto-embedded-source-code-modal')
});
clipboardModal.on('success', onCopySuccess);
}
var localhostRegex = new RegExp(/^(?:http|https):\/\/localhost\:?[0-9]*\//);
var mailtoRegex = new RegExp(/^mailto:/);
var filterRegex = new RegExp('/' + window.location.host + '/');
var isInternal = (href) => {
return filterRegex.test(href) || localhostRegex.test(href) || mailtoRegex.test(href);
}
// Inspect non-navigation links and adorn them if external
var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool)');
var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool):not(.about-link)');
for (var i=0; i<links.length; i++) {
const link = links[i];
if (!isInternal(link.href)) {
Expand Down
Binary file modified docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-5-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-6-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading

0 comments on commit 4940b94

Please sign in to comment.