Mercurial > repos > bimib > cobraxy

comparison COBRAxy/utils/CBS_backend.py @ 4:41f35c2f0c7b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author luca_milaz
date Wed, 18 Sep 2024 10:59:10 +0000
parents
children
comparison
equal deleted inserted replaced
3:1f3ac6fd9867 4:41f35c2f0c7b
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
182 val_max=np.max([df_fva.loc[reaction,"minimum"],df_fva.loc[reaction,"maximum"]])
183
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
200 return coefficients_df