-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcriteria_helpers.py
192 lines (151 loc) · 6.7 KB
/
criteria_helpers.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
from models import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def run_counter(alpha1 = 34.73, alpha2 = 49.36 , alpha3 = 32.73, alpha4 = 49.54, delta1 = 1.93, delta2 = 0.69, Kd = 10.44, n = 4.35):
"""
TESTING
"""
clock_enable = 0 # clock enable - 1 is on, 0 is off
set_number = 2
# simulation parameters
t_end = 200
N = 1000
# # model parameters
# alpha1 = 34.73 # protein_production 0.1 to 50
# alpha2 = 49.36 # protein_production 0.1 to 50
# alpha3 = 32.73 # protein_production 0.1 to 50
# alpha4 = 49.54 # protein_production 0.1 to 50
# delta1 = 1.93 # protein_degradation 0.001 to 50
# delta2 = 0.69 # protein_degradation 0.001 to 50
# Kd = 10.44 # Kd 0.01 to 250
# n = 4.35 # hill 1 to 5
inhibition_rate = 200
params_ff = (alpha1, alpha2, alpha3, alpha4, delta1, delta2, Kd, n)
# three-bit counter with external clock
# a1, not_a1, q1, not_q1, a2, not_a2, q2, not_q2, a3, not_a3, q3, not_q3
# Y0 = np.array([0]*12) # initial state
Y0 = np.array([0] * 12) # initial state
T = np.linspace(0, t_end, N) # vector of timesteps
# numerical interation
Y = odeint(counter_model_3, Y0, T, args=(params_ff, clock_enable * inhibition_rate, set_number))
Y_reshaped = np.split(Y, Y.shape[1], 1)
Q1 = Y_reshaped[2]
not_Q1 = Y_reshaped[3]
Q2 = Y_reshaped[6]
not_Q2 = Y_reshaped[7]
Q3 = Y_reshaped[10]
not_Q3 = Y_reshaped[11]
clk = modulated_clock(T, clock_enable * inhibition_rate, Kd, n)
# plotting the results for current params
# plt.subplot(1, 2, 1)
# plt.plot(T, Q1, label='q1')
# plt.plot(T, Q2, label='q2')
# plt.plot(T, Q3, label='q3')
# plt.plot(T, clk, '--', linewidth=2, label="CLK", color='black', alpha=0.25)
# plt.legend()
# MY CODE:
def normalize(np_array):
"""
Normalizes a numpy array in range 0 to 1.
:param np_array: np_array to normalize
:return: normalized np_array
"""
return (np_array - np_array.min()) / (np_array.max() - np_array.min())
normalized_clk = normalize(clk)
normalized_T = normalize(T)
normalized_Q1 = normalize(Q1)
normalized_Q2 = normalize(Q2)
normalized_Q3 = normalize(Q3)
# plotting normalized graph
# plt.subplot(1, 2, 2)
# plt.plot(normalized_T, normalized_Q1, label='q1')
# plt.plot(normalized_T, normalized_Q2, label='q2')
# plt.plot(normalized_T, normalized_Q3, label='q3')
# plt.plot(normalized_T, normalized_clk, '--', linewidth=2, label="CLK", color='black', alpha=0.25)
# plt.tight_layout()
# plt.show()
return normalized_Q1, normalized_Q2, normalized_Q3, normalized_clk, normalized_T
# print(normalized_T)
def get_ideal_response(period, T, clk, mode, threshold):
"""
Function returns values for an ideal electrical response.
:param period: user defined values to repeat ie. [0,1] or [0,0,1,1]
:param T: timestep (x value)
:param clk: clock (sin) value (y value)
:param mode: trigger mode for clock 1 - rising edge, 2 - lowering edge, 3 - positive state, 4 - negative state
:param threshold: the value at which we evaluate the output
:return: array of ideal values of size T
"""
# At threshold value of the clock the output should change immediately
current_value = 0
ideal_values = []
prev_clk_value = clk[0]
if mode == 1:
# Positive edge
for timestep in range(len(T)):
if clk[timestep] > threshold >= prev_clk_value:
current_value += 1
if current_value >= len(period):
current_value = 0
ideal_values.append(period[current_value])
prev_clk_value = clk[timestep]
elif mode == 2:
# Negative edge
for timestep in range(len(T)):
if clk[timestep] <= threshold < prev_clk_value:
current_value += 1
if current_value >= len(period):
current_value = 0
ideal_values.append(current_value)
prev_clk_value = clk[timestep]
# These modes can be made with rising and lowering edge
elif mode == 3:
for timestep in range(len(T)):
if clk[timestep] < 0.005:
if changed == 1:
current_value = not current_value
changed = 0
else:
changed = 1
ideal_values.append(current_value)
elif mode == 4:
for timestep in range(len(T)):
if clk[timestep] > 0.995:
if changed == 1:
current_value = not current_value
changed = 0
else:
changed = 1
ideal_values.append(current_value)
else:
print("Error mode is incorrect")
return
return np.array(ideal_values)
# # plotting normalized graph
# plt.plot(normalized_T, normalized_Q1, label='q1')
# plt.plot(normalized_T, ideal_response_q1, label='q1 ideal')
# # plt.plot(normalized_T, normalized_Q2, label='q2')
# # plt.plot(normalized_T, ideal_response_q2, label='q2 ideal')
# # plt.plot(normalized_T, normalized_Q3, label='q3')
# # plt.plot(normalized_T, ideal_response_q3, label='q3 ideal')
# plt.plot(normalized_T, normalized_clk, '--', linewidth=2, label="CLK", color='black', alpha=0.25)
# plt.show()
from sklearn.metrics import mean_squared_error
def fitness_function_calc(normalized_Q1, ideal_response_q1, normalized_Q2, ideal_response_q2, normalized_Q3, ideal_response_q3):
# Calculate Mean Squared Error for all three outputs, if one signal has significantly higher error, it will be penalized
mse_q1 = mean_squared_error(normalized_Q1, ideal_response_q1)
mse_q2 = mean_squared_error(normalized_Q2, ideal_response_q2)
mse_q3 = mean_squared_error(normalized_Q3, ideal_response_q3)
# Calculate the fitness
fitness = mse_q1**2 + mse_q2**2 + mse_q3**2
return -fitness
# Define your fitness function
def fitness_function(params):
alpha1, alpha2, alpha3, alpha4, delta1, delta2, Kd, n = params
normalized_Q1, normalized_Q2, normalized_Q3, normalized_clk, normalized_T = run_counter(alpha1, alpha2, alpha3, alpha4, delta1, delta2, Kd, n)
# Your fitness function here, e.g., the sum of Q3 values
ideal_response_q1 = get_ideal_response([0, 1], normalized_T, normalized_clk, 1, 0.01)
ideal_response_q2 = get_ideal_response([1, 1, 0, 0], normalized_T, normalized_clk, 1, 0.01)
ideal_response_q3 = get_ideal_response([0, 0, 1, 1, 1, 1, 0, 0], normalized_T, normalized_clk, 1, 0.01)
fitness = fitness_function_calc(normalized_Q1, ideal_response_q1, normalized_Q2, ideal_response_q2, normalized_Q3, ideal_response_q3)
return -fitness # Minimize, so we negate the fitness