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

choose_TS_transition( propensities ) click to toggle source

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
choose_from_discrete_distribution( distribution ) click to toggle source

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
fire!( transition ) click to toggle source

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
gillespie_delta_time( propensities ) click to toggle source

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
gillespie_step!(propensities) click to toggle source

Step.

# File lib/y_petri/core/timed/gillespie.rb, line 44
def gillespie_step! propensities
  t = choose_TS_transition( propensities )
  fire! t
end
rng() click to toggle source

Returns a random number generator, only created once.

# File lib/y_petri/core/timed/gillespie.rb, line 14
def rng
  @rng ||= ::Random
end
step!(Δt=simulation.step) click to toggle source

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
update_next_gillespie_time( propensities ) click to toggle source

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