Module: YPetri::Core::Timed::Gillespie

Defined in:
lib/y_petri/core/timed/gillespie.rb

Overview

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.

Instance Method Summary collapse

Instance Method Details

#choose_from_discrete_distribution(distribution) ⇒ Object

Given a discrete probability distributions, this function makes a random choice of a category.



61
62
63
64
65
66
67
 
# 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

#choose_TS_transition(propensities) ⇒ Object

Chooses the transition to fire.



71
72
73
 
# File 'lib/y_petri/core/timed/gillespie.rb', line 71

def choose_TS_transition( propensities )
  transitions.fetch choose_from_discrete_distribution( propensities )
end

#fire!(transition) ⇒ Object

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.



79
80
81
82
83
84
85
86
 
# 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) ⇒ Object

Computes Δ for the period of Δt.



51
52
53
54
55
56
 
# 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) ⇒ Object

Step.



44
45
46
47
 
# File 'lib/y_petri/core/timed/gillespie.rb', line 44

def gillespie_step! propensities
  t = choose_TS_transition( propensities )
  fire! t
end

#rng ⇒ Object

Returns a random number generator, only created once.



14
15
16
 
# File 'lib/y_petri/core/timed/gillespie.rb', line 14

def rng
  @rng ||= ::Random
end

#step!(Δt = simulation.step) ⇒ Object

Makes a stochastic number of Gillespie steps necessary to span the period Δt.



20
21
22
23
24
25
26
27
28
29
30
31
32
33
 
# 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) ⇒ Object

This method updates next firing time given propensities.



37
38
39
40
 
# 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