|
293
|
1 from swiglpk import *
|
|
|
2 import random
|
|
|
3 import pandas as pd
|
|
|
4 import numpy as np
|
|
|
5 import cobra as cb
|
|
|
6
|
|
|
7 # Initialize LP problem
|
|
|
8 def initialize_lp_problem(S):
|
|
|
9
|
|
|
10 len_vector=len(S.keys())
|
|
|
11 values=list(S.values())
|
|
|
12 indexes=list(S.keys())
|
|
|
13 ia = intArray(len_vector+1);
|
|
|
14 ja = intArray(len_vector+1);
|
|
|
15 ar = doubleArray(len_vector+1);
|
|
|
16
|
|
|
17 i=0
|
|
|
18 ind_row=[indexes[i][0]+1 for i in range(0, len(values) )]
|
|
|
19 ind_col=[indexes[i][1]+1 for i in range(0, len(values) )]
|
|
|
20 for i in range(1, len(values) + 1):
|
|
|
21 ia[i]=ind_row[i-1]
|
|
|
22 ja[i]=ind_col[i-1]
|
|
|
23 ar[i] = values[i-1]
|
|
|
24
|
|
|
25 nrows=S.shape[0]
|
|
|
26 ncol=S.shape[1]
|
|
|
27
|
|
|
28 return len_vector, values, indexes, ia, ja, ar, nrows, ncol
|
|
|
29
|
|
|
30
|
|
|
31
|
|
|
32 # Solve LP problem from the structure of the metabolic model
|
|
|
33 def create_and_solve_lp_problem(lb,ub,nrows, ncol, len_vector, ia, ja, ar,
|
|
|
34 obj_coefs,reactions,return_lp=False):
|
|
|
35
|
|
|
36
|
|
|
37 lp = glp_create_prob();
|
|
|
38 glp_set_prob_name(lp, "sample");
|
|
|
39 glp_set_obj_dir(lp, GLP_MAX);
|
|
|
40 glp_add_rows(lp, nrows);
|
|
|
41 eps = 1e-16
|
|
|
42 for i in range(nrows):
|
|
|
43 glp_set_row_name(lp, i+1, "constrain_"+str(i+1));
|
|
|
44 glp_set_row_bnds(lp, i+1, GLP_FX, 0.0, 0.0);
|
|
|
45 glp_add_cols(lp, ncol);
|
|
|
46 for i in range(ncol):
|
|
|
47 glp_set_col_name(lp, i+1, "flux_"+str(i+1));
|
|
|
48 glp_set_col_bnds(lp, i+1, GLP_DB,lb[i]-eps,ub[i]+eps);
|
|
|
49 glp_load_matrix(lp, len_vector, ia, ja, ar);
|
|
|
50
|
|
|
51 try:
|
|
|
52 fluxes,Z=solve_lp_problem(lp,obj_coefs,reactions)
|
|
|
53 if return_lp:
|
|
|
54 return fluxes,Z,lp
|
|
|
55 else:
|
|
|
56 glp_delete_prob(lp);
|
|
|
57 return fluxes,Z
|
|
|
58 except Exception as e:
|
|
|
59 glp_delete_prob(lp)
|
|
|
60 raise Exception(e)
|
|
|
61
|
|
|
62
|
|
|
63 # Solve LP problem from the structure of the metabolic model
|
|
|
64 def solve_lp_problem(lp,obj_coefs,reactions):
|
|
|
65
|
|
|
66 # Set the coefficients of the objective function
|
|
|
67 i=1
|
|
|
68 for ind_coef in obj_coefs:
|
|
|
69 glp_set_obj_coef(lp, i, ind_coef);
|
|
|
70 i+=1
|
|
|
71
|
|
|
72 # Initialize the parameters
|
|
|
73 params=glp_smcp()
|
|
|
74 params.presolve=GLP_ON
|
|
|
75 params.msg_lev = GLP_MSG_ALL
|
|
|
76 params.tm_lim=4000
|
|
|
77 glp_init_smcp(params)
|
|
|
78
|
|
|
79 # Solve the problem
|
|
|
80 glp_scale_prob(lp,GLP_SF_AUTO)
|
|
|
81
|
|
|
82 value=glp_simplex(lp, params)
|
|
|
83
|
|
|
84 Z = glp_get_obj_val(lp);
|
|
|
85
|
|
|
86 if value == 0:
|
|
|
87 fluxes = []
|
|
|
88 for i in range(len(reactions)): fluxes.append(glp_get_col_prim(lp, i+1))
|
|
|
89 return fluxes,Z
|
|
|
90 else:
|
|
|
91 raise Exception("error in LP problem. Problem:",str(value))
|
|
|
92
|
|
|
93
|
|
|
94 # Create LP structure
|
|
|
95 def create_lp_structure(model):
|
|
|
96
|
|
|
97 reactions=[el.id for el in model.reactions]
|
|
|
98 coefs_obj=[reaction.objective_coefficient for reaction in model.reactions]
|
|
|
99
|
|
|
100 # Lower and upper bounds
|
|
|
101 lb=[reaction.lower_bound for reaction in model.reactions]
|
|
|
102 ub=[reaction.upper_bound for reaction in model.reactions]
|
|
|
103
|
|
|
104 # Create S matrix
|
|
|
105 S=cb.util.create_stoichiometric_matrix(model,array_type="dok")
|
|
|
106
|
|
|
107 return S,lb,ub,coefs_obj,reactions
|
|
|
108
|
|
|
109 # CBS sampling interface
|
|
|
110 def randomObjectiveFunctionSampling(model, nsample, coefficients_df, df_sample):
|
|
|
111
|
|
|
112 S,lb,ub,coefs_obj,reactions = create_lp_structure(model)
|
|
|
113 len_vector, values, indexes, ia, ja, ar, nrow, ncol = initialize_lp_problem(S)
|
|
|
114
|
|
|
115 for i in range(nsample):
|
|
|
116
|
|
|
117 coefs_obj=coefficients_df.iloc[:,i].values
|
|
|
118
|
|
|
119 if coefs_obj[-1]==1: #minimize
|
|
|
120 coefs_obj= coefs_obj[0:-1] * -1
|
|
|
121 else:
|
|
|
122 coefs_obj=coefs_obj[0:-1]
|
|
|
123
|
|
|
124 fluxes,Z = create_and_solve_lp_problem(lb,ub, nrow, ncol, len_vector,
|
|
|
125 ia, ja, ar, coefs_obj,reactions,return_lp=False)
|
|
|
126 df_sample.loc[i] = fluxes
|
|
|
127 pass
|
|
|
128
|
|
|
129 def randomObjectiveFunctionSampling_cobrapy(model, nsample, coefficients_df, df_sample):
|
|
|
130
|
|
|
131 for i in range(nsample):
|
|
|
132
|
|
|
133 dict_coeff={}
|
|
|
134 if(coefficients_df.iloc[-1][i]==1):
|
|
|
135 type_problem = -1 #minimize
|
|
|
136 else:
|
|
|
137 type_problem = 1
|
|
|
138
|
|
|
139 for rxn in [reaction.id for reaction in model.reactions]:
|
|
|
140 dict_coeff[model.reactions.get_by_id(rxn)] = coefficients_df.loc[rxn][i] * type_problem
|
|
|
141
|
|
|
142 model.objective = dict_coeff
|
|
|
143 solution = model.optimize().fluxes
|
|
|
144 for rxn, flux in solution.items():
|
|
|
145 df_sample.loc[i][rxn] = flux
|
|
|
146
|
|
|
147 pass
|
|
|
148
|
|
|
149 # Create random coefficients for CBS
|
|
|
150 def randomObjectiveFunction(model, n_samples, df_fva, seed=0):
|
|
|
151
|
|
|
152
|
|
|
153 #reactions = model.reactions
|
|
|
154 reactions = [reaction.id for reaction in model.reactions]
|
|
|
155 cont=seed
|
|
|
156 list_ex=reactions.copy()
|
|
|
157 list_ex.append("type_of_problem")
|
|
|
158 coefficients_df = pd.DataFrame(index=list_ex,columns=[str(i) for i in range(n_samples)])
|
|
|
159
|
|
|
160 for i in range(0, n_samples):
|
|
|
161
|
|
|
162 cont=cont+1
|
|
|
163 random.seed(cont)
|
|
|
164
|
|
|
165 # Genera un numero casuale tra 0 e 1
|
|
|
166 threshold = random.random() #coefficiente tra 0 e 1
|
|
|
167
|
|
|
168 for reaction in reactions:
|
|
|
169
|
|
|
170 cont=cont+1
|
|
|
171 random.seed(cont)
|
|
|
172
|
|
|
173 val=random.random()
|
|
|
174
|
|
|
175 if val>threshold:
|
|
|
176
|
|
|
177 cont=cont+1
|
|
|
178 random.seed(cont)
|
|
|
179
|
|
|
180 c=2*random.random()-1 #coefficiente tra -1 e 1
|
|
|
181
|
|
325
|
182 val_max = np.max([abs(df_fva.loc[reaction,"minimum"]), abs(df_fva.loc[reaction,"maximum"])])
|
|
|
183
|
|
293
|
184 if val_max!=0: #solo se la fva รจ diversa da zero
|
|
|
185 coefficients_df.loc[reaction,str(i)] = c/val_max #divido per la fva
|
|
|
186 else:
|
|
|
187 coefficients_df.loc[reaction,str(i)] = 0
|
|
|
188
|
|
|
189 else:
|
|
|
190 coefficients_df.loc[reaction,str(i)] = 0
|
|
|
191
|
|
|
192 cont=cont+1
|
|
|
193 random.seed(cont)
|
|
|
194
|
|
|
195 if random.random()<0.5:
|
|
|
196 coefficients_df.loc["type_of_problem",str(i)] = 0 #maximize
|
|
|
197 else:
|
|
|
198 coefficients_df.loc["type_of_problem",str(i)] = 1 #minimize
|
|
|
199
|
|
325
|
200 return coefficients_df
|
