Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

stochastic degradation toy model #884

Merged
merged 56 commits into from
Aug 31, 2019
Merged

stochastic degradation toy model #884

merged 56 commits into from
Aug 31, 2019

Conversation

danielfridman98
Copy link
Contributor

I've coded up the stochastic degradation toy model (Radek Erban et al., 2007) which will be a useful toy model for ABC. I've uploaded the Jupyter notebook file into examples under the name toy-model-stochastic-degradation (Erban et al., 2007)

@codecov
Copy link

codecov bot commented Jul 30, 2019

Codecov Report

Merging #884 into master will not change coverage.
The diff coverage is 100%.

Impacted file tree graph

@@          Coverage Diff           @@
##           master   #884    +/-   ##
======================================
  Coverage     100%   100%            
======================================
  Files          53     56     +3     
  Lines        5157   5749   +592     
======================================
+ Hits         5157   5749   +592
Impacted Files Coverage Δ
pints/toy/_stochastic_degradation_model.py 100% <100%> (ø)
pints/toy/__init__.py 100% <100%> (ø) ⬆️
pints/_mcmc/__init__.py 100% <0%> (ø) ⬆️
pints/toy/_fitzhugh_nagumo_model.py 100% <0%> (ø) ⬆️
pints/toy/_annulus.py 100% <0%> (ø) ⬆️
pints/_log_likelihoods.py 100% <0%> (ø) ⬆️
pints/_mcmc/_slice_doubling.py 100% <0%> (ø)
pints/_mcmc/_slice_stepout.py 100% <0%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f0fa8be...4787385. Read the comment docs.

@danielfridman98
Copy link
Contributor Author

danielfridman98 commented Jul 31, 2019

Just updated my stochastic degradation toy model, including all of Chon's comments as well as additional changes. I've written up an example notebook as well for the stochastic degradation model.

Next step is to write the unit tests

Copy link
Collaborator

@ben18785 ben18785 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice work Daniel -- great start! Don't be discouraged by the number of comments by Chon and me. Creating good software takes a lot of time for everyone. I, in particular, am crap at the software-writing side and people always have tons of suggestions as to how to improve things!

A few things:

  • I think you'll want to remove the _stochastic_decay script.
  • See my note about adding mean and variance functions for this toy model.
  • Chat with Chon and perhaps Michael about how to include a likelihood for this model -- we may want to create a separate StochasticDegradationLogLikelihood class, I suspect.
  • In terms of the example notebook, I think you should add/change the following:
    • In the description, you need to make a distinction between the model (A -> phi at rate k) and how we simulate from the model (using the Gillespie algorithm - which is an exact stochastic simulation algorithm).
    • For the model, I would cite Erban; for the algorithm, cite the original Gillespie paper.
    • No need for caps on the Stochast Degradation Model in the title or descriptive text.

Finally, thinking a bit ahead here. There are ABC methods (for example, lazy ABC) that work with approximate simulation algorithms. The idea behind these is that you first simulate from a computationally cheap algorithm then check for 'not rejecting'. If you do not reject, you then simulate from the computationally more intense non-approximate algorithm.

Gillespie (which is what is coded currently) is an exact SSA (stochastic simulation algorithm). Tau leaping is an approximate SSA, which is much faster than Gillespie to simulate from. As such, I think we probably want an approx_simulate method (or a number of such methods) for this and other stochastic reaction equations. To cope with this, I suspect we want to create an abstract parent class for all stochastic reaction methods that has this extra functionality. Happy to chat this over Skype if easier to convey!

pints/toy/_stochastic_decay_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_decay_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_decay_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_decay_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_decay_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_degradation_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_degradation_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_degradation_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_degradation_model.py Outdated Show resolved Hide resolved
pints/toy/_stochastic_degradation_model.py Outdated Show resolved Hide resolved
Copy link
Member

@MichaelClerx MichaelClerx left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @danielfridman98 and @ben18785 !
This is looking really good, but there are still two issues to deal with:

  1. Please check all your maths in the docstrings. The syntax is explained here https://github.com/pints-team/pints/blob/master/CONTRIBUTING.md#using-maths-in-documentation and just below that it shows you how to test if you've done it right. I've had a look myself so think it's fine now, but please double-check I didn't change the mathematics

  2. Please update the unit tests

I've also changed "concentration" to "molecule count" everywhere, as that's what's being simulated

docs/source/toy/index.rst Outdated Show resolved Hide resolved
examples/toy-model-stochastic-degradation.ipynb Outdated Show resolved Hide resolved
times = np.linspace(0, 100, 101)
values = model.simulate(parameters, times)

# Test output of Gillespie algorithm
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These tests need to go!
A unit test should test that a class does what its supposed to do, but not how it actually does it!
So you never check internals / private properties of a class, you only check that the methods return what they ought to return. (It gets a bit trickier if the methods update some internal state of the object, but that's not the case here!)
(And don't solve this by adding a method to return these properties either ;-) )

So all your tests should be tests of the values returned

Copy link
Collaborator

@ben18785 ben18785 Aug 31, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've kept these in because I asked for them -- it's very easy to make a mistake with the interpolation function for the stochastic models, so would like to keep them please.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What would happen to values if it went wrong? Can we check for that?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not without knowing the internals, no.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually I've figured a work-round. Give me a few mins

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, no I haven't -- can we keep this as is please?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why don't you seed the random generator, print the output, check it manually, and then hardcode those values in the test?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've created separate functions that do the interpolation and return the raw values (and added tests for these). Let us know what you think @MichaelClerx

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still don't quite understand why that's needed. Will users benefit from having these methods as part of the API? If so, then great, if not then this is not the way to go.
But if we can't tell whether its working from the output of simulate(), then why do we care?

(I found this article https://enterprisecraftsmanship.com/2017/10/23/unit-testing-private-methods/ that says the stuff I'm saying a bit better, even if it doesn't give that many arguments)

pints/tests/test_toy_stochastic_degradation_model.py Outdated Show resolved Hide resolved
Copy link
Member

@MichaelClerx MichaelClerx left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's leave this discussion for now, I'll trust you that the newly exposed methods are useful for the users :-)

@MichaelClerx MichaelClerx dismissed stale reviews from chonlei and ben18785 August 31, 2019 18:26

outdated

@MichaelClerx MichaelClerx merged commit cd46a42 into master Aug 31, 2019
@MichaelClerx MichaelClerx deleted the toy_models-daniel branch August 31, 2019 18:27
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants