Gillespie Algorithm( draft_2nd )

gillespie2.rb

require 'csv'

CSV_Output = false
#==================================================
#	Math Util Methods
#==================================================
def factorial(i)
	ans = 1
	1.upto(i) { |x| ans *= x }
	return ans
end

def permutation(n, k)
	ans = 1
	k.times do 
		ans *= n
		n -= 1
	end
	return ans
end

def combination(n, k)
	kk = (n - k) < k ? (n - k) : k
	return permutation(n, kk) / factorial(kk)
end

def average(list)
	temp = 0.0
	list.each {|x| temp += x.to_f() }
	return temp / list.length
end

def sum_array(array)
	s = nil
	array.each do |x|
		if s == nil
			s = x
		else
			s += x
		end
	end
	return s
end

def numerical_integrate(func, dx, start, finish)
	sum = 0.0

	x = start
	while x < finish 
		sum += func.call(x) * dx
		x += dx
	end
	
	return sum
end

#==================================================
#	Main Routine
#==================================================

def checkover0(array)
	array.each do |x|
		if x <= 0.0
			return false
		end
	end
	return true
end

class Reaction
	def initialize(reaction_left = nil, reaction_right = nil, probability = 0)
		@reaction_left = reaction_left
		@reaction_right = reaction_right
		@probability = probability
	end
	attr_accessor :reaction_right, :reaction_left, :probability
end


def gillespie(current, reactions, tlimit)
	def num_check(hash)
		hash.each do |k, v|
			if v <= 0 then return false end
		end
		return true
	end

	t = 0.0
	prev_t = 0.0
	csv = CSV.open("Gillespie.csv", "w")	if CSV_Output == true
	cnt = 0
	while t < tlimit && num_check(current) == true
		if cnt % 10 == 0
			printf("%f\n", t)
		else
			cnt = 0
		end

		if CSV_Output == true 
			if 0.001 < t - prev_t 
				temp = Array.new
				temp[0] = t
				temp << current['X'] << current['Y'] << current['Z'] << current['W']
				csv << temp
				prev_t = t
			end
		end

		a = Array.new(reactions.size)
		#calcurate c_u * h_u
		reactions.each_with_index do |react, idx|
			a[idx] = react.probability
			react.reaction_left.each do |specie, use|
				a[idx] *= combination( current[specie], use )
			end
		end
		a0 = sum_array(a)
		r1 = rand()
		r2 = rand()
		dt = (Math.log(1 / r1)) / a0.to_f

		r2 = r2 * a0 
		u   = -1;tmp = 0.0
		while tmp < r2 && u < reactions.size
			u += 1
			tmp += a[u]
		end

		# Ru occurance
		t += dt
		reactions[u].reaction_left.each {|s, n| current[s] -= n }
		reactions[u].reaction_right.each{|s, n| current[s] += n }
	end
	

	csv.close	if CSV_Output == true
end

if __FILE__ == $0
	#	X <-> Y
	#	2X <-> Z
	#	W + X <-> 2X
	current = {'X' => 10000, 'Y' => 10000, 'Z' => 10000, 'W' => 10000}
	react1 = Reaction.new({'X' => 1}, {'Y' => 1}, 1 )	# X -> Y
	react2 = Reaction.new({'Y' => 1}, {'X' => 1}, 2 )	# Y -> X
	react3 = Reaction.new({'X' => 2}, {'Z' => 1}, 3 )	# 2X -> Z
	react4 = Reaction.new({'Z' => 1}, {'X' => 2}, 4 )	# Z -> 2X
	react5 = Reaction.new({'X' => 1, 'W' => 1}, {'X' => 2}, 5)	# W + X -> 2X
	react6 = Reaction.new({'X' => 2}, {'X' => 1, 'W' => 1}, 6)	# 2X -> W + X
	gillespie(current, [react1, react2, react3, react4, react5, react6], 0.5)
end