metabolism.py

import numpy as np
import scipy.linalg as linalg
import scipy.optimize as optimize

# r0: rate constant for initial partition -> active partition
# r1: rate constant for active partition -> inactive partition
# rate_constant = np.log(2)/half_life
def single_compound_matrix(r0, r1):
	return np.array([[-r0, r0, 0], [0, -r1, r1], [0, 0, 0]]).T

# evaluates active partition only
def _eval_at_time(mat, t):
	return linalg.expm(mat*t)[1, 0]

eval_at_time = np.vectorize(_eval_at_time, excluded=[0])

def find_max(mat):
	result = optimize.minimize_scalar(lambda t: -eval_at_time(mat, t))
	return result.x, -result.fun

def find_single_compound_matrix(half_life, t_max):
	r1 = np.log(2)/half_life
	def objective(r0):
		if r0 <= 0:
			return np.inf
		return (find_max(single_compound_matrix(r0, r1))[0] - t_max)**2
	result = optimize.minimize_scalar(objective)
	return single_compound_matrix(result.x, r1)