module YPetri::Core::Timed::Gillespie
Plain Gillespie
algorithm.
The distinguishing quality of Gillespie
method is, that it does not work from from Δt towards Δstate. Instead, it makes a random choice of the transition to fire (weighted by the transition propensities) and a random choice of Δt. Both next transition to fire, and the size of Δt to slice off from the time axis are thus stochastically determined.
Public Instance Methods
Chooses the transition to fire.
# File lib/y_petri/core/timed/gillespie.rb, line 71 def choose_TS_transition( propensities ) transitions.fetch choose_from_discrete_distribution( propensities ) end
Given a discrete probability distributions, this function makes a random choice of a category.
# File lib/y_petri/core/timed/gillespie.rb, line 61 def choose_from_discrete_distribution( distribution ) sum = rng.rand * distribution.reduce( :+ ) distribution.index do |p| sum -= p sum <= 0 end end
Fires a transitions. More precisely, performs a single transition event with certain stoichiometry, adding / subtracting the number of quanta to / from the codomain places as indicated by the stoichiometry.
# File lib/y_petri/core/timed/gillespie.rb, line 79 def fire!( transition ) cd, sto = transition.codomain, transition.stoichiometry mv = simulation.m_vector cd.zip( sto ).each { |pl, coeff| mv.set( pl, mv.fetch( pl ) + pl.quantum * coeff ) } @gillespie_time = @next_gillespie_time end
Computes Δ for the period of Δt.
# File lib/y_petri/core/timed/gillespie.rb, line 51 def gillespie_delta_time( propensities ) sum = Σ propensities # mean_period = 1 / sum # TODO: This line seem to be useless # Exponential distribution Distribution::Exponential.p_value( rng.rand, sum ) end
Step.
# File lib/y_petri/core/timed/gillespie.rb, line 44 def gillespie_step! propensities t = choose_TS_transition( propensities ) fire! t end
Returns a random number generator, only created once.
# File lib/y_petri/core/timed/gillespie.rb, line 14 def rng @rng ||= ::Random end
Makes a stochastic number of Gillespie
steps necessary to span the period Δt.
# File lib/y_petri/core/timed/gillespie.rb, line 20 def step! Δt=simulation.step @gillespie_time = curr_time = simulation.time target_time = curr_time + Δt propensities = propensity_vector_TS.column_to_a update_next_gillespie_time( propensities ) until ( @next_gillespie_time > target_time ) gillespie_step! propensities simulation.recorder.alert! propensities = propensity_vector_TS.column_to_a update_next_gillespie_time( propensities ) end simulation.increment_time! Δt print '.' end
This method updates next firing time given propensities.
# File lib/y_petri/core/timed/gillespie.rb, line 37 def update_next_gillespie_time( propensities ) @next_gillespie_time = @gillespie_time + gillespie_delta_time( propensities ) end