Team:iBowu-China/Model pythoncode simple model optimize
simple_model_optimize.py
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def model(x, a_mr, r, be_mr, a_pr, be_pr, a_leak, a_mb, k_b, be_mb, a_pb, be_pb):
N = x[0]
A = x[1]
mr = a_mr * r * N / be_mr
pr = a_pr * mr / be_pr
mb = (a_leak + a_mb * r * N * (k_b * A * pr)**2 / (1 + (k_b * A * pr)**2)) / be_mb
pb = a_pb * mb / be_pb
return pb
def readData(fileName):
x = [[],[]]
y = []
with open(fileName, "r") as f_o:
lines = f_o.read().split("\n")
for i in range(1, len(lines)):
line = lines[i]
if not line:
continue
l_s = line.split(",")
a = float(l_s[0])
n = float(l_s[1])
g1 = float(l_s[3])
g2 = float(l_s[5])
g3 = float(l_s[7])
x[0].append(n)
x[1].append(a)
#y.append((g1+g2+g3) / 3)
y.append(g1)
x[0].append(n)
x[1].append(a)
y.append(g2)
x[0].append(n)
x[1].append(a)
y.append(g3)
return np.array(x), np.array(y)
if __name__ == "__main__":
x, y = readData("train_data.csv")
p = np.random.rand(11)
bound = ([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [np.inf, 1, np.inf,np.inf,np.inf,np.inf,np.inf,np.inf,np.inf,np.inf,np.inf])
popt, popv = curve_fit(model, x, y,p, bounds=bound)
x_new = np.arange(0, 35, 1)
y_new = np.arange(0, 1500, 10)
x_new, y_new = np.meshgrid(x_new, y_new)
z_new = np.zeros_like(x_new)
for i in range(x_new.shape[1]):
for j in range(y_new.shape[0]):
z_new[j, i] = model([x_new[0, i], y_new[j, 0]], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], popt[7], popt[8], popt[9], popt[10])
print(popt)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(x_new, y_new, z_new, cmap='rainbow', label = "simulated distribution")
ax.scatter(x[0], x[1], y, c = "black", depthshade=True, label = "experimental data distribution")
ax.set_xlabel("DNA($ng/uL)$")
ax.set_ylabel("AHL($nM$)")
ax.set_zlabel("GFP($ng/ul$)")
ax.grid(False)
ax.legend()
plt.show()
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def model(x, a_mr, r, be_mr, a_pr, be_pr, a_leak, a_mb, k_b, be_mb, a_pb, be_pb):
N = x[0]
A = x[1]
mr = a_mr * r * N / be_mr
pr = a_pr * mr / be_pr
mb = (a_leak + a_mb * r * N * (k_b * A * pr)**2 / (1 + (k_b * A * pr)**2)) / be_mb
pb = a_pb * mb / be_pb
return pb
def readData(fileName):
x = [[],[]]
y = []
with open(fileName, "r") as f_o:
lines = f_o.read().split("\n")
for i in range(1, len(lines)):
line = lines[i]
if not line:
continue
l_s = line.split(",")
a = float(l_s[0])
n = float(l_s[1])
g1 = float(l_s[3])
g2 = float(l_s[5])
g3 = float(l_s[7])
x[0].append(n)
x[1].append(a)
#y.append((g1+g2+g3) / 3)
y.append(g1)
x[0].append(n)
x[1].append(a)
y.append(g2)
x[0].append(n)
x[1].append(a)
y.append(g3)
return np.array(x), np.array(y)
if __name__ == "__main__":
x, y = readData("train_data.csv")
p = np.random.rand(11)
bound = ([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [np.inf, 1, np.inf,np.inf,np.inf,np.inf,np.inf,np.inf,np.inf,np.inf,np.inf])
popt, popv = curve_fit(model, x, y,p, bounds=bound)
x_new = np.arange(0, 35, 1)
y_new = np.arange(0, 1500, 10)
x_new, y_new = np.meshgrid(x_new, y_new)
z_new = np.zeros_like(x_new)
for i in range(x_new.shape[1]):
for j in range(y_new.shape[0]):
z_new[j, i] = model([x_new[0, i], y_new[j, 0]], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], popt[7], popt[8], popt[9], popt[10])
print(popt)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(x_new, y_new, z_new, cmap='rainbow', label = "simulated distribution")
ax.scatter(x[0], x[1], y, c = "black", depthshade=True, label = "experimental data distribution")
ax.set_xlabel("DNA($ng/uL)$")
ax.set_ylabel("AHL($nM$)")
ax.set_zlabel("GFP($ng/ul$)")
ax.grid(False)
ax.legend()
plt.show()