Gillespie Algorithm(draft)

gillespie_draft.rb

require 'csv'
#==================================================
#	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

def simu(x, c, tlimit)
	cnt = 0
	t = 0.0
	a = Array.new
	# X: x[0]  Y: x[1]  Z: x[2]  W: x[3]
	csv = CSV.open("Output.csv", "w")
	tempresult = []
	while checkover0(x) == true && t <= tlimit

		#== Print Out
		cnt += 1
		if cnt % 5 == 0
			printf("%f, X: %d, Y: %d, Z: %d, W: %d \n", t, x[0], x[1], x[2], x[3])
			cnt = 0
		tempresult = []
		tempresult[0] = t
		tempresult += x
		csv << tempresult
		end
		#=============

		#	X <-> Y
		#	2X <-> Z
		#	W + X <-> 2X
		a[0] = c[0] * combination(x[0], 0)
		a[1] = c[1] * combination(x[1], 1)
		a[2] = c[2] * combination(x[0], 2)
		a[3] = c[3] * combination(x[2], 1)
		a[4] = c[4] * combination(x[0], 1) * combination(x[3], 1)
		a[5] = c[5] * combination(x[0], 2)
		a0 = sum_array(a)
	
		r1 = rand()
		r2 = rand()
		tau = (Math.log(1 / r1)) / a0.to_f
	
		r2 = r2 * a0 
		u   = -1
		tmp = 0
		while tmp < r2 && u < x.length
			u += 1
			tmp += a[u]
		end
	
		# Ru occurance
		# X: x[0]  Y: x[1]  Z: x[2]  W: x[3]
		t += tau
		
		#Reactions(HARD CODING)
		case u
		when 0
			x[0] = x[0] - 1
			x[1] = x[1] + 1
		when 1
			x[0] = x[0] + 1
			x[1] = x[1] - 1
		when 2
			x[0] = x[0] - 2
			x[2] = x[2] + 1
		when 3
			x[0] = x[0] + 2
			x[2] = x[2] - 1
		when 4
			x[0] = x[0] + 1
			x[3] = x[3] - 1
		when 5
			x[0] = x[0] - 1
			x[3] = x[3] + 1
		end
	end
	csv.close
end

if __FILE__ == $0
	simu([10000, 10000, 10000, 10000], [0.1, 0.2, 0.3, 0.4, 0.5, 0.6], 10.0)

end