|
381
|
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_ERR
|
|
|
76 params.tm_lim=4000
|
|
|
77 glp_init_smcp(params)
|
|
|
78
|
|
|
79 glp_term_out(GLP_OFF)
|
|
|
80
|
|
|
81 try:
|
|
|
82
|
|
|
83 # Solve the problem
|
|
|
84 glp_scale_prob(lp,GLP_SF_AUTO)
|
|
|
85
|
|
|
86 value=glp_simplex(lp, params)
|
|
|
87
|
|
|
88 Z = glp_get_obj_val(lp);
|
|
|
89
|
|
|
90 if value == 0:
|
|
|
91 fluxes = []
|
|
|
92 for i in range(len(reactions)): fluxes.append(glp_get_col_prim(lp, i+1))
|
|
|
93 return fluxes,Z
|
|
|
94 else:
|
|
|
95 raise Exception("error in LP problem. Problem:",str(value))
|
|
|
96 except Exception as e:
|
|
|
97 # Re-enable terminal output for error reporting
|
|
|
98 glp_term_out(GLP_ON)
|
|
|
99 raise Exception(e)
|
|
|
100 finally:
|
|
|
101 # Re-enable terminal output after solving
|
|
|
102 glp_term_out(GLP_ON)
|
|
|
103
|
|
|
104 # Create LP structure
|
|
|
105 def create_lp_structure(model):
|
|
|
106
|
|
|
107 reactions=[el.id for el in model.reactions]
|
|
|
108 coefs_obj=[reaction.objective_coefficient for reaction in model.reactions]
|
|
|
109
|
|
|
110 # Lower and upper bounds
|
|
|
111 lb=[reaction.lower_bound for reaction in model.reactions]
|
|
|
112 ub=[reaction.upper_bound for reaction in model.reactions]
|
|
|
113
|
|
|
114 # Create S matrix
|
|
|
115 S=cb.util.create_stoichiometric_matrix(model,array_type="dok")
|
|
|
116
|
|
|
117 return S,lb,ub,coefs_obj,reactions
|
|
|
118
|
|
|
119 # CBS sampling interface
|
|
|
120 def randomObjectiveFunctionSampling(model, nsample, coefficients_df, df_sample):
|
|
|
121
|
|
|
122 S,lb,ub,coefs_obj,reactions = create_lp_structure(model)
|
|
|
123 len_vector, values, indexes, ia, ja, ar, nrow, ncol = initialize_lp_problem(S)
|
|
|
124
|
|
|
125 for i in range(nsample):
|
|
|
126
|
|
|
127 coefs_obj=coefficients_df.iloc[:,i].values
|
|
|
128
|
|
|
129 if coefs_obj[-1]==1: #minimize
|
|
|
130 coefs_obj= coefs_obj[0:-1] * -1
|
|
|
131 else:
|
|
|
132 coefs_obj=coefs_obj[0:-1]
|
|
|
133
|
|
|
134 fluxes,Z = create_and_solve_lp_problem(lb,ub, nrow, ncol, len_vector,
|
|
|
135 ia, ja, ar, coefs_obj,reactions,return_lp=False)
|
|
|
136 df_sample.loc[i] = fluxes
|
|
|
137 pass
|
|
|
138
|
|
|
139 def randomObjectiveFunctionSampling_cobrapy(model, nsample, coefficients_df, df_sample):
|
|
|
140
|
|
|
141 for i in range(nsample):
|
|
|
142
|
|
|
143 dict_coeff={}
|
|
|
144 if(coefficients_df.iloc[-1][i]==1):
|
|
|
145 type_problem = -1 #minimize
|
|
|
146 else:
|
|
|
147 type_problem = 1
|
|
|
148
|
|
|
149 for rxn in [reaction.id for reaction in model.reactions]:
|
|
|
150 dict_coeff[model.reactions.get_by_id(rxn)] = coefficients_df.loc[rxn][i] * type_problem
|
|
|
151
|
|
|
152 model.objective = dict_coeff
|
|
|
153 solution = model.optimize().fluxes
|
|
|
154 for rxn, flux in solution.items():
|
|
|
155 df_sample.loc[i][rxn] = flux
|
|
|
156
|
|
|
157 pass
|
|
|
158
|
|
|
159 # Create random coefficients for CBS
|
|
|
160 def randomObjectiveFunction(model, n_samples, df_fva, seed=0):
|
|
|
161
|
|
|
162
|
|
|
163 #reactions = model.reactions
|
|
|
164 reactions = [reaction.id for reaction in model.reactions]
|
|
|
165 cont=seed
|
|
|
166 list_ex=reactions.copy()
|
|
|
167 list_ex.append("type_of_problem")
|
|
|
168 coefficients_df = pd.DataFrame(index=list_ex,columns=[str(i) for i in range(n_samples)])
|
|
|
169
|
|
|
170 for i in range(0, n_samples):
|
|
|
171
|
|
|
172 cont=cont+1
|
|
|
173 random.seed(cont)
|
|
|
174
|
|
|
175 # Genera un numero casuale tra 0 e 1
|
|
|
176 threshold = random.random() #coefficiente tra 0 e 1
|
|
|
177
|
|
|
178 for reaction in reactions:
|
|
|
179
|
|
|
180 cont=cont+1
|
|
|
181 random.seed(cont)
|
|
|
182
|
|
|
183 val=random.random()
|
|
|
184
|
|
|
185 if val>threshold:
|
|
|
186
|
|
|
187 cont=cont+1
|
|
|
188 random.seed(cont)
|
|
|
189
|
|
|
190 c=2*random.random()-1 #coefficiente tra -1 e 1
|
|
|
191
|
|
|
192 val_max = np.max([abs(df_fva.loc[reaction,"minimum"]), abs(df_fva.loc[reaction,"maximum"])])
|
|
|
193
|
|
|
194 if val_max!=0: #solo se la fva รจ diversa da zero
|
|
|
195 coefficients_df.loc[reaction,str(i)] = c/val_max #divido per la fva
|
|
|
196 else:
|
|
|
197 coefficients_df.loc[reaction,str(i)] = 0
|
|
|
198
|
|
|
199 else:
|
|
|
200 coefficients_df.loc[reaction,str(i)] = 0
|
|
|
201
|
|
|
202 cont=cont+1
|
|
|
203 random.seed(cont)
|
|
|
204
|
|
|
205 if random.random()<0.5:
|
|
|
206 coefficients_df.loc["type_of_problem",str(i)] = 0 #maximize
|
|
|
207 else:
|
|
|
208 coefficients_df.loc["type_of_problem",str(i)] = 1 #minimize
|
|
|
209
|
|
|
210 return coefficients_df
|
