Gillespie started major rebase

[zierenberg]

A lot of things have changed, including some API.

It all started with implementing a proper Gillespie wrapper around KineticMonteCarlo, which included developing some nice event handler. After talking with George, I remove all modular substructure and cleaned up API.

@fmikulasch please check if the current API works for you and can be integrated into your workflow.

[fmikulasch]

The last commit seems to lack some of the routines. Can you fix this? Then I will go over it.

[zierenberg]

Which routines are you missing. I rather cleaned the package and changed the API, so routines may have a new name ...

For the inhomogeneous pp for example you should call next(Inhomogeneous Poisson(), rng, rates,...)

On Tue, Mar 3, 2020, 13:25 fmikulasch notifications@github.com wrote:

The last commit seems to lack some of the routines. Can you fix this? Then I will go over it.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/zierenberg/MonteCarloX.jl/pull/33?email_source=notifications&email_token=AGBMN5ZNHQL6KTALX3Q62JTRFTZK3A5CNFSM4K766MJ2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOENTJEWA#issuecomment-593924696, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGBMN52ADKHHQ3MXKZC4AR3RFTZK3ANCNFSM4K766MJQ .

[fmikulasch]

Sorry I wasn't clear. It seems the last commit removed almost all files in src/ or am I seeing this wrongly?

[zierenberg]

Sorry, I renamed files to match with typical package convention in julia and ... forgot to add the new files to git. Hence local tests passed and I did not notice. Now, you should be able to see them.

[codecov-io]

Codecov Report

Merging #33 into master will decrease coverage by 15.74%. The diff coverage is 64.33%.

@@             Coverage Diff             @@
##           master      #33       +/-   ##
===========================================
- Coverage   80.21%   64.46%   -15.75%     
===========================================
  Files           9        9               
  Lines         187      273       +86     
===========================================
+ Hits          150      176       +26     
- Misses         37       97       +60     

Impacted Files	Coverage Δ
src/MonteCarloX.jl	100.00% <ø> (+50.00%)	:arrow_up:
src/gillespie.jl	0.00% <0.00%> (ø)
src/kinetic_monte_carlo.jl	30.30% <30.30%> (ø)
src/poisson_process.jl	42.42% <42.42%> (ø)
src/importance_sampling.jl	72.00% <72.00%> (ø)
src/event_handler.jl	81.66% <81.66%> (ø)
src/utils.jl	93.75% <93.75%> (ø)
src/cluster_wolff.jl	100.00% <100.00%> (ø)
src/reweighting.jl	100.00% <100.00%> (ø)
... and 7 more

---

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 1300e47...e93a54e. Read the comment docs.

[fmikulasch]

Hi, sorry took me a while to get to it. I tried to get the Gillespie algorithm running but I am unsure how it is intended to be used. Is the evolution to be stored through the event handler? If so I find it a bit inconvenient to have to define my own event handler before I can use it. I think the user should have more of a "plug and play" experience. Could you give me an example of how you think it should be used?

[zierenberg]

Hi, there is a set of event handlers implemented. I have a running program to show how these can be used. This will be added to examples at some point. If you do not want to use this feature, you can always go back to kineticMC (which is the proper algorithm where you can pass a list of weights). The Gillespie wrapper was (naming convention wise) an artifact from me not knowing which me to give the child. This we can still debate.

However, right now the question was if all the changes I did are still okay with YOUR usage of the framework ;). Namely the Inhomogeneous poisson process.

Cheers,

On Mon, Mar 9, 2020, 22:43 fmikulasch notifications@github.com wrote:

Hi, sorry took me a while to get to it. I tried to get the Gillespie algorithm running but I am unsure how it is intended to be used. Is the evolution to be stored through the event handler? If so I find it a bit inconvenient to have to define my own event handler before I can use it. I think the user should have more of a "plug and play" experience. Could you give me an example of how you think it should be used?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/zierenberg/MonteCarloX.jl/pull/33?email_source=notifications&email_token=AGBMN56DJ5LFA2R6N3D6MHTRGVPG5A5CNFSM4K766MJ2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEOJFV7A#issuecomment-596794108, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGBMN55ASYL5W63KDUINEGDRGVPG5ANCNFSM4K766MJQ .

[zierenberg]

Suggestion:

I move the current Gillespie advance!() function to kinetic_monte_carlo.jl (because it belongs there) and we write an additional wrapper for Gillespie, which essentially implements the original context of Gillespie (namely for chemical reactions) to be consistent with the wikipedia entry stating:

Gillespie developed two different, but equivalent formulations; the direct method and the first reaction method. Below is a summary of the steps to run the algorithm (mathematics omitted): 1) Initialization: Initialize the number of molecules in the system, reaction constants, and random number generators. 2) Monte Carlo step: Generate random numbers to determine the next reaction to occur as well as the time interval. The probability of a given reaction to be chosen is proportional to the number of substrate molecules, the time interval is exponentially distributed with mean 1 / R T O T {\displaystyle 1/R{\mathrm {TOT} }} {\displaystyle 1/R{\mathrm {TOT} }}.

3. Update: Increase the time by the randomly generated time in Step 2. Update the molecule count based on the reaction that occurred.
4. Iterate: Go back to Step 2 unless the number of reactants is zero or the simulation time has been exceeded. see more here

We could then start with

• 1) initialize (which returns an event handler)
• 2) advance (which essentially does 2-4)
• 3) step (which is 2)

[fmikulasch]

I found that it is difficult to debug the kineticMC when you mistakenly hand it negative or zero rates - it doesn't produce an error and the next sampled timestep is infinty or negative. We could think about sacrificing performance and warn the user? Otherwise I approve :) Edit: I think we could also merge some of the search functions used to select events...

[zierenberg]

I found that it is difficult to debug the kineticMC when you mistakenly hand it negative or zero rates - it doesn't produce an error and the next sampled timestep is infinty or negative. We could think about sacrificing performance and warn the user? Otherwise I approve :) Edit: I think we could also merge some of the search functions used to select events...

Good point. I would argue that negative timestep would be warning enough :). Maybe this can be checked without much effort via

if not(sum(rates)>0) 
  ...
end

The problem is rather: "where would we set the bound?". This could already fail for sum(rates) about 1e-10 if all rates are zero but the sum has cumulated floating point errors (actually thinking of which ... is this a thing we should look out for? In terms that every so and so often the rates.sum should be recomputed....)

[zierenberg]

@fmikulasch : I think for now, we could merge this into master to keep on working on details in different (smaller) development steps. What do you think? I would be curious to see if you can merge the pull request.