* отредактировал 6/17/10Python/Biomolecular Physics - попытка кодирования простого стохастического моделирования системы, демонстрирующей условное поведение!
Я пытаюсь понять, как улучшить мой код (сделайте его более питоническим). Кроме того, я заинтересован в написании более интуитивных «условностей», которые описывали бы сценарии, которые являются обычным явлением в биохимии. Условные критерии в приведенной ниже программе, которые я объяснил в ответе № 2, но меня не устраивает код - он работает отлично, но не является очевидным и нелегко реализовать для более сложных условных сценариев. Идеи приветствуются. Комментарии/критика приветствуются. Первый опыт публикации @ stackoverflow - прокомментируйте этикет, если необходимо.
код создает список значений, которые являются решением для следующего упражнения:
«В языке программирования по вашему выбору, реализация первой реакции Джиллеспи алгоритм для изучения временной ход реакции А --- > В, в котором переход от А к В может иметь место только в том случае, если присутствует другое соединение С, и где С динамически взаимно преобразуется с D, как смоделировано в Петри-сети ниже. Предположим, что существует 100 молекул A, 1 от C, и нет B или D, присутствующих в начале реакции. Установите kAB до 0,1 с-1 и оба kCD и kDC до 1,0 с-1. Имитируйте поведение системы за 100 секунд. "
def sim():
# Set the rate constants for all transitions
kAB = 0.1
kCD = 1.0
kDC = 1.0
# Set up the initial state
A = 100
B = 0
C = 1
D = 0
# Set the start and end times
t = 0.0
tEnd = 100.0
print "Time\t", "Transition\t", "A\t", "B\t", "C\t", "D"
# Compute the first interval
transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)
# Loop until the end time is exceded or no transition can fire any more
while t <= tEnd and transition >= 0:
print t, '\t', transition, '\t', A, '\t', B, '\t', C, '\t', D
t += interval
if transition == 0:
A -= 1
B += 1
if transition == 1:
C -= 1
D += 1
if transition == 2:
C += 1
D -= 1
transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)
def transitionData(A, B, C, D, kAB, kCD, kDC):
""" Returns nTransition, the number of the firing transition (0: A->B,
1: C->D, 2: D->C), and interval, the interval between the time of
the previous transition and that of the current one. """
RAB = kAB * A * C
RCD = kCD * C
RDC = kDC * D
dt = [-1.0, -1.0, -1.0]
if RAB > 0.0:
dt[0] = -math.log(1.0 - random.random())/RAB
if RCD > 0.0:
dt[1] = -math.log(1.0 - random.random())/RCD
if RDC > 0.0:
dt[2] = -math.log(1.0 - random.random())/RDC
interval = 1e36
transition = -1
for n in range(len(dt)):
if dt[n] > 0.0 and dt[n] < interval:
interval = dt[n]
transition = n
return transition, interval
if __name__ == '__main__':
sim()
Трудно разобрать ваш фактический вопрос из задыхающегося потока потока сознания. Не могли бы вы обобщить, каков ваш фактический вопрос на самом деле, чтобы мы могли на него ответить? –
Кажется, что бит монте-карло это неверно. Если 'dt [2]' меньше 1e36, тогда переход всегда будет равен 2. Я признаю, что я не совсем понимаю такое моделирование, но эти части кажутся мне подозрительными. Похоже, что 'dt' следует использовать каким-то образом в качестве весов для случайного определения того, какой процесс будет проходить. Существует случайность, но похоже, что это происходит неправильно. –
Спасибо, ребята. Мне еще предстоит посмотреть на Монте Карло, Джастин, но я прокомментирую, как только я доберусь до него. С. Лотт, ты прав, я плохо разбираюсь в этом. Я немного отредактирую свой пост. Ха-ха, я только что узнал, что означает «передача возвращаемых значений вызывающему» - это была та часть, в которой я был незнаком. Мне удалось уйти с некоторыми простыми симами, определив функции с массивами до сих пор * вздох *. Но да, суть моего поста - WTF? Я думаю, что я на правильном пути. Сегодня я расскажу об этом позже. – DocDubya