cMCS_enumerator_test.py

# exec(open('optlang_mcs_test5.py').read())
#%%
import cobra
# import cobra.util.array
from optlang_enumerator.cMCS_enumerator import *
from optlang_enumerator.mcs_computation import *
import time
import numpy
import os
import sys
import efmtool_link.efmtool_intern as efmtool_intern
import efmtool_link.efmtool4cobra as efmtool4cobra
import pickle
import optlang_enumerator
import optlang_enumerator.cobra_cnapy
from importlib import reload
#%%
reload(optlang_enumerator.cMCS_enumerator)
from optlang_enumerator.cMCS_enumerator import *
reload(optlang_enumerator.mcs_computation)
from optlang_enumerator.mcs_computation import *

#%%
ex = optlang_enumerator.cobra_cnapy.CNApyModel.read_sbml_model(r"tests/metatool_example_no_ext.xml")
# ex.solver = 'cplex'
stdf = cobra.util.array.create_stoichiometric_matrix(ex, array_type='DataFrame')
rev = [r.reversibility for r in ex.reactions]
reac_id = stdf.columns.tolist()
reac_id_symbols = get_reac_id_symbols(reac_id)
# target = [(equations_to_matrix(ex, ["-1 Pyk", "-1 Pck"]), [-1, -1])] # -Pyk does not work
target = [[("Pyk", ">=", 1), ("Pck", ">=", 1)]]
target = [relations2leq_matrix(parse_relations(t, reac_id_symbols=reac_id_symbols), reac_id) for t in target]
#target.append(target[0]) # duplicate target
flux_expr= [r.flux_expression for r in ex.reactions] # !! lose validity when the solver is changed !!
target_constraints= get_leq_constraints(ex, target)
kn = efmtool_intern.null_rat_efmtool(stdf.values)
#res = cobra.flux_analysis.single_reaction_deletion(ex, processes=1) # no interactive multiprocessing on Windows
    
# %%
info = dict()
e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, stdf.values, rev, target, 
                                        bigM= 100, threshold=0.1, split_reversible_v=True, irrev_geq=True)
# e.model.objective = e.minimize_sum_over_z
# e.write_lp_file('testM')
mcs,_ = e.enumerate_mcs(max_mcs_size=5, info=info)
print(info)

e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, stdf.values, rev, target,
                                        threshold=0.1, split_reversible_v=False, irrev_geq=False, kn=kn) #, ref_set=mcs)
# e.model.objective = e.minimize_sum_over_z
# e.write_lp_file('testI')
#e.model.configuration.verbosity = 3
#e.model.configuration.presolve = 'auto' # presolve remains off on the CPLEX side
#e.model.configuration.presolve = True # presolve on, this is CPLEX default
#e.model.configuration.lp_method = 'auto' # sets lpmethod back to automatic on the CPLEX side
#e.model.problem.parameters.reset() # works (CPLEX specific)
# without reset() optlang switches presolve off and fixes lpmethod
mcs2,_ = e.enumerate_mcs(max_mcs_size=5, info=info)
print(info)

print(len(set(mcs).intersection(set(mcs2))))
print(set(mcs) == set(mcs2))
all(check_mcs(ex, target[0], mcs, optlang.interface.INFEASIBLE))

# %%
# e.model.problem.solution.MIP.get_best_objective()
# e.model.problem.solution.MIP.get_mip_relative_gap()
e = ConstrainedMinimalCutSetsEnumerator(optlang.gurobi_interface, stdf.values, rev, target, 
                                        bigM= 100, threshold=0.1, split_reversible_v=True, irrev_geq=True)
# e.model.problem.parameters.mip.tolerances.mipgap.set(0.99)
e.model.problem.Params.MipGap = 0.99
info = dict()
mcs3,_ = e.enumerate_mcs(max_mcs_size=5, enum_method=3, model=ex, targets=target, info=info)
print(info)
print(e.evs_sz_lb)
print(set(mcs) == set(mcs3))

#%%
e = ConstrainedMinimalCutSetsEnumerator(optlang.gurobi_interface, stdf.values, rev, target, 
                                        bigM= 100, threshold=0.1, split_reversible_v=True, irrev_geq=True)
mcs2,_ = e.enumerate_mcs(max_mcs_size=5, enum_method=2, max_mcs_num=20, info=info)
print(info)
print(set(mcs) == set(mcs2))
# %%
e.model.configuration.verbosity = 3
e.model.problem.Params.PoolSearchMode = 2
e.model.problem.Params.PoolSolutions = 100
e.evs_sz.ub = 1
e.evs_sz_lb = 1
e.model.problem.optimize()
z_vars = [e.model.problem.getVarByName(z.name) for z in e.z_vars]
print(e.model.problem.getAttr("Xn", z_vars))
e.model.problem.Params.SolutionNumber = 1
e.model.problem.getAttr("Xn", z_vars)
print(e.model.problem.status)

# %%
e = ConstrainedMinimalCutSetsEnumerator(optlang.gurobi_interface, stdf.values, rev, target, 
                                        bigM= 100, threshold=0.1, split_reversible_v=True, irrev_geq=True)
e.model.objective = e.minimize_sum_over_z
e.model.configuration.verbosity = 3
z_vars = [e.model.problem.getVarByName(z.name) for z in e.z_vars]
cut_set_cb = GUROBImakeMCSCallback(z_vars, ex, target) #, redundant_constraints=False)
e.model.problem.Params.PoolSearchMode = 2
e.model.problem.Params.LazyConstraints = 1 # must be activated explicitly
def call_back_func(model, where): # encapsulating with functools.partial not accepted by Gurobi
    cut_set_cb.invoke(model, where) # passing this directly not accepted by Gurobi
e.model.problem.optimize(call_back_func)
len(cut_set_cb.minimal_cut_sets)

#%%
e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, stdf.values, rev, target, 
                                        bigM= 100, threshold=0.1, split_reversible_v=True, irrev_geq=True)
info = dict()
e.model.configuration._iocp.mip_gap = 0.99
mcs3,_ = e.enumerate_mcs(max_mcs_size=5, enum_method=3, model=ex, targets=target, info=info)
print(info)
print(e.evs_sz_lb)
print(set(mcs) == set(mcs3))

#%%
# flux_lb= -1000*numpy.ones(stdf.values.shape[1])
# flux_lb[numpy.logical_not(rev)] = 0
# flux_ub= 1000*numpy.ones(stdf.values.shape[1])
flux_lb = numpy.array([ex.reactions[i].lower_bound for i in range(len(ex.reactions))])
flux_ub = numpy.array([ex.reactions[i].upper_bound for i in range(len(ex.reactions))])
cobra_config = cobra.Configuration()
flux_lb[numpy.isinf(flux_lb)] = cobra_config.lower_bound
flux_ub[numpy.isinf(flux_ub)] = cobra_config.upper_bound
desired = [(equations_to_matrix(ex, ["-1 AspCon"]), [-1], flux_lb, flux_ub)]
# desired[0] = list(desired[0])
# desired[0][2] = numpy.array([   0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0.
# ,0.,0.,0.,1.,0., -999., -999., -999., -999.,1.,1.,0.])
# desired[0][3] = numpy.array([1000.,1000., 999., 999., 999., 500., 500.,1000., 999.,1000.,
#  1000.,1000., 499.5,999., 999.,1000., 999., 999., 999.,1000.,
#  1000.,1000.,1000., 999. ])
# cuts = numpy.array([ True,True,True,True,True,True,True, False,True,True,True,True
# ,True,True,True, False,True,True,True,True,True, False, False,True])
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, stdf.values, rev, target,
                                        threshold=0.1, split_reversible_v=False, irrev_geq=False,
                                        desired=desired)
# e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, stdf.values, rev, target,
#                                         threshold=0.1, split_reversible_v=True, irrev_geq=True,
#                                         desired=desired, SOS_constraints=True)
# e.model.objective = e.minimize_sum_over_z
# e.write_lp_file('testA')
mcs3,_ = e.enumerate_mcs() #max_mcs_size=5)
print(len(mcs3))
print(all(check_mcs(ex, target[0], mcs3, optlang.interface.INFEASIBLE)))
print(all(check_mcs(ex, desired[0], mcs3, optlang.interface.OPTIMAL)))
#%%
with ex as exb: # cobrapy FVA raises error with unbounded reactions
    exb.solver = 'cplex'
    for r in exb.reactions:
        if r.lower_bound == -numpy.inf:
            r.lower_bound = -1000
        if r.upper_bound == numpy.inf:
            r.upper_bound = 1000
    # mcs4,_ = compute_mcs(exb, target, desired=desired, enum_method=1, max_mcs_size=20)
    mcs4,_ = compute_mcs(exb, target, enum_method=2, max_mcs_size=20, desired=desired, network_compression=True)
print(len(mcs4))
set(mcs3) == set(mcs4)
#%%
# subset_compression = efmtool_intern.CompressionMethod[:]([efmtool_intern.CompressionMethod.CoupledZero, efmtool_intern.CompressionMethod.CoupledCombine, efmtool_intern.CompressionMethod.CoupledContradicting])
rd, subT, comprec = efmtool_intern.compress_rat_efmtool(stdf.values, rev, remove_cr=True,
            compression_method= efmtool_intern.subset_compression) #[0:2]
# rd = stdf.values
# subT = numpy.eye(rd.shape[1])
rev_rd = numpy.logical_not(numpy.any(subT[numpy.logical_not(rev), :], axis=0))
#%%
# exc, subT = efmtool4cobra.compress_model(ex)
exc = ex.copy()
subT = efmtool4cobra.compress_model_sympy(exc)
rd = cobra.util.array.create_stoichiometric_matrix(exc, array_type='dok') #, dtype=numpy.object) # does actually work with rationals
rev_rd = [r.reversibility for r in exc.reactions]

#%%
target_rd = [(T@subT, t) for T, t in target]
# e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, rd, rev_rd, target_rd, 
#                                         bigM= 100, threshold=0.1, split_reversible_v=True, irrev_geq=True)
# e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, rd, rev_rd, target_rd, 
#                                         bigM= 100, threshold=0.1, split_reversible_v=True, kn=efmtool_intern.null_rat_efmtool(rd))
e = ConstrainedMinimalCutSetsEnumerator(optlang.coinor_cbc_interface, rd, rev_rd, target_rd, 
                                        bigM= 100, threshold=0.1, split_reversible_v=True, kn=efmtool_intern.null_rat_efmtool(rd))
# e.model.objective = e.minimize_sum_over_z
# e.model.configuration._iocp.mip_gap = 0.99
rd_mcs,_ = e.enumerate_mcs(max_mcs_size=5, enum_method=1, model=exc, targets=target_rd)
# rd_mcs,_ = e.enumerate_mcs(max_mcs_size=5)
print(len(rd_mcs))
#set(expand_mcs(rd_mcs, subT)) == set(map(lambda x: tuple(numpy.where(x)[0]), mcs2))
set(expand_mcs(rd_mcs, subT)) == set(mcs2)

#%%
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, stdf.values, rev, target,
                                        threshold=0.1, split_reversible_v=False, irrev_geq=False) #, ref_set=mcs)
# e.model.configuration.verbosity = 3
# e.model.objective = e.minimize_sum_over_z
mcs2p,_ = e.enumerate_mcs(max_mcs_size=5, enum_method=2) #, max_mcs_num=1)
print(set(mcs) == set(mcs2p))
# e.model.problem.solution.get_objective_value() == 0

#%%
cuts = numpy.full(24, True, dtype=bool)
cuts[0] = False
cuts[23] = False
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, stdf.values, rev, target, 
                                        cuts=cuts, threshold=0.1, split_reversible_v=True, irrev_geq=True)
e.model.objective = e.minimize_sum_over_z
# e.write_lp_file('testI2')
t = time.time()
mcs3,_ = e.enumerate_mcs(max_mcs_size=5)
print(time.time() - t)

e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, stdf.values, rev, target, 
                                        cuts=cuts, bigM=100, threshold=0.1, split_reversible_v=True, irrev_geq=False)
e.model.objective = e.minimize_sum_over_z
# e.write_lp_file('testI2')
t = time.time()
mcs4,_ = e.enumerate_mcs(max_mcs_size=5)
print(time.time() - t)

print(set(mcs3) == set([m for m in mcs if 0 not in m and 23 not in m]))
print(set(mcs4) == set(mcs3))

#%%
ecc2 = cobra.io.read_sbml_model(r"../cnapy-projects/ECC2comp/model.sbml")
# ecc2.solver = 'glpk'
ecc2_stdf = cobra.util.array.create_stoichiometric_matrix(ecc2, array_type='DataFrame')
cuts= numpy.full(ecc2_stdf.shape[1], True, dtype=bool) # results do not agree when exchange reactions can be cut, problem with tiny fluxes (and M too small)
# for r in ecc2.boundary:
#     cuts[ecc2_stdf.columns.get_loc(r.id)] = False
# for r in range(len(ecc2.reactions)):
#     if ecc2.reactions[r].boundary:
#         cuts[r] = False
cuts = numpy.array([not r.boundary for r in ecc2.reactions])
reac_id = ecc2_stdf.columns.tolist()
reac_id_symbols = get_reac_id_symbols(reac_id)
# ecc2_mue_target = [[("Growth", ">=", 0.01), ("GlcUp", "<=", 10), ("AcUp", "<=", 10), ("GlycUp", "<=", 10), ("SuccUp", "<=", 10)]]
ecc2_mue_target = [[("Growth", ">=", 0.01)]]
ecc2_mue_target = [relations2leq_matrix(parse_relations(t, reac_id_symbols=reac_id_symbols), reac_id) for t in ecc2_mue_target]
bounds_mat, bounds_rhs = reaction_bounds_to_leq_matrix(ecc2)
ecc2_mue_target = [(scipy.sparse.vstack((t[0], bounds_mat), format='csr'), numpy.hstack((t[1], bounds_rhs))) for t in ecc2_mue_target]
ecc2_mue_target_constraints= get_leq_constraints(ecc2, ecc2_mue_target)
# %%
e = ConstrainedMinimalCutSetsEnumerator(optlang.gurobi_interface, ecc2_stdf.values, [r.reversibility for r in ecc2.reactions], ecc2_mue_target,
                                        cuts=cuts, threshold=0.1, split_reversible_v=True, irrev_geq=True, bigM=1000)#, SOS_constraints=True)
# %%
# e.model.objective = e.minimize_sum_over_z
# e.model.configuration.tolerances.feasibility = 1e-9
# e.model.configuration.tolerances.optimality = 1e-9
# e.model.configuration.tolerances.integrality = 1e-10
#e.write_lp_file('testI')
#e.model.configuration.verbosity = 3
e.evs_sz_lb = 1 
info = dict()
ecc2_mcs,_ = e.enumerate_mcs(max_mcs_size=3, enum_method=2, model=ecc2, targets=ecc2_mue_target, info=info)
# ecc2_mcs,_ = e.enumerate_mcs(max_mcs_size=3, enum_method=1, info=info)
print(len(ecc2_mcs), info) # in this example SOS are somewhat slower than indicators
# %%
e.evs_sz.lb = 1
e.evs_sz.ub = 3
e.model.objective = e.minimize_sum_over_z
e.model.configuration.verbosity = 3
z_vars = [e.model.problem.getVarByName(z.name) for z in e.z_vars]
cut_set_cb = GUROBImakeMCSCallback(z_vars, ecc2, ecc2_mue_target) #, redundant_constraints=False)
e.model.problem.Params.PoolSearchMode = 2
e.model.problem.Params.LazyConstraints = 1 # must be activated explicitly
def call_back_func(model, where): # encapsulating with functools.partial not accepted by Gurobi
    cut_set_cb.invoke(model, where) # passing this directly not accepted by Gurobi
e.model.problem.optimize(call_back_func)

# %%
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, ecc2_stdf.values, [r.reversibility for r in ecc2.reactions], ecc2_mue_target,
                                        cuts=cuts, threshold=0.1, split_reversible_v=True, irrev_geq=True)#, SOS_constraints=True)
e.model.problem.parameters.mip.pool.intensity.set(4)
e.model.problem.parameters.mip.strategy.search.set(1) # traditional branch-and-cut search
z_idx = e.model.problem.variables.get_indices([z.name for z in e.z_vars])

cut_set_cb = MakeMinimalCutSetCallback(z_idx, ecc2, ecc2_mue_target) #, redundant_constraints=False)
e.model.problem.set_callback(cut_set_cb, cplex.callbacks.Context.id.candidate)
e.model.problem.parameters.mip.limits.populate.set(e.model.problem.parameters.mip.pool.capacity.get())
e.model.configuration.verbosity = 3
e.model.problem.parameters.emphasis.mip.set(1) # integer feasibility
e.model.objective = e.minimize_sum_over_z
e.evs_sz.lb = 1
e.evs_sz.ub = 3
e.model.problem.populate_solution_pool()
print(e.model.problem.solution.get_status_string(), cut_set_cb.candidate_count,
      len(set(cut_set_cb.minimal_cut_sets)))

# %% MCS guessing heuristics simulation
correct = []
known = []
unclear = []
not_cs = []
for i in range(1, len(cut_set_cb.minimal_cut_sets)-1): # assume the i-th MCS has been found
    for j in range(0, i): # try heuristic with all previous MCS
        common = cut_set_cb.minimal_cut_sets[i].intersection(cut_set_cb.minimal_cut_sets[j])
        if len(common) > 0:
            guess = (cut_set_cb.minimal_cut_sets[i].union(cut_set_cb.minimal_cut_sets[j]) - common) # or iteratively take out only 1 of the common?
            # check that guess is not a superset of any other known MCS
            if check_mcs(ecc2, ecc2_mue_target[0], [guess], optlang.interface.INFEASIBLE)[0]:
                guess = frozenset(make_minimal_cut_set(ecc2, list(guess), ecc2_mue_target_constraints))
                if guess in cut_set_cb.minimal_cut_sets[i+1:]: # check if the guess is contained in the future MCS
                    correct.append(guess)
                elif guess in cut_set_cb.minimal_cut_sets[:i]:
                    known.append(guess)
                else:
                    unclear.append(guess) # are really larger MCS
            else:
                not_cs.append(guess)
print(len(correct), len(set(correct)), len(unclear), len(set(unclear)), len(known))
# unclear = [sorted(tuple(m)) for m in set(unclear)]
# is_mcs = check_mcs(ecc2, ecc2_mue_target[0], unclear, optlang.interface.INFEASIBLE)

# %%
# ecc2_mcs,_ = compute_mcs(ecc2, ecc2_mue_target, [], cuts, 2, 3, 1000, 100)
ecc2_mcs,_ = compute_mcs(ecc2, ecc2_mue_target, cuts=cuts, enum_method=3, max_mcs_size=3, network_compression=True)
print(len(ecc2_mcs))
ecc2_mcs_rxns= [tuple(ecc2_stdf.columns[r] for r in mcs) for mcs in ecc2_mcs]
print(ecc2_mcs_rxns)
ecc2_mcsF,_ = compute_mcs(ecc2, ecc2_mue_target, cuts=cuts, enum_method=2, max_mcs_size=3, network_compression=False)
print(all(check_mcs(ecc2, ecc2_mue_target[0], ecc2_mcs, optlang.interface.INFEASIBLE)))
print(set(ecc2_mcs) == set(ecc2_mcsF))

# %%
with open("ecc2_mcs.pkl","wb") as f:
    pickle.dump(set(ecc2_mcs), f)

# %%
# ecc2.solver = 'glpk_exact'
# ecc2.solver = 'glpk'
# ecc2.solver = 'cplex'
with ecc2 as model:
    model.objective = model.problem.Objective(0)
    model.solver = 'glpk_exact'
    model_target_constraints= get_leq_constraints(model, ecc2_mue_target)
    model.add_cons_vars(model_target_constraints[0])
    res = cobra.flux_analysis.single_reaction_deletion(model, reaction_list=ecc2_stdf.columns[cuts], processes=1) # no interactive multiprocessing on Windows
del model
single_cuts= list(map(tuple, res.index[res.values[:, 1] == optlang.interface.INFEASIBLE]))
print(set(single_cuts) == set(mcs for mcs in ecc2_mcs_rxns if len(mcs) == 1))
print(set(single_cuts) - set(mcs for mcs in ecc2_mcs_rxns if len(mcs) == 1))
check_mcs(ecc2, ecc2_mue_target[0], single_cuts, optlang.interface.INFEASIBLE)

# %% FVA with copied model
model = ecc2.copy() 
model.objective = model.problem.Objective(0)
fva_tol = 1e-9 # with CPLEX 1e-8 leads to removal of EX_adp_c, 1e-9 keeps EX_adp_c
model.tolerance = fva_tol # prevent essential EX_meoh_ex from being blocked, sets solver feasibility/optimality tolerances
# model.solver.configuration.tolerances.feasibility = 1e-9
# model.solver = 'glpk' #'glpk_exact' # appears to make problems for context management
# model_mue_target_constraints= get_leq_constraints(model, ecc2_mue_target)
# model.add_cons_vars(model_mue_target_constraints[0])
fva_res = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0, processes=1) # no interactive multiprocessing on Windows
print(fva_res.loc['EX_adp_c',:])
print(fva_res.loc['EX_AMOB_ex'])
blocked = []
blocked_rxns = []
for i in range(fva_res.values.shape[0]):
    # if fva_res.values[i, 0] == 0 and fva_res.values[i, 1] == 0:
    if fva_res.values[i, 0] >= -fva_tol and fva_res.values[i, 1] <= fva_tol:
        blocked.append(i)
        blocked_rxns.append(fva_res.index[i])
print(blocked_rxns)

# %% FVA in context
#model = ecc2.copy() # copy model because switching solver in context sometimes gives an error (?!?)
ecc2.solver = 'coinor_cbc'
ecc2.solver.problem.opt_tol = 1e-9
with ecc2 as model:
    model.objective = model.problem.Objective(0)
    fva_tol = 1e-9 # with CPLEX 1e-8 leads to removal of EX_adp_c, 1e-9 keeps EX_adp_c
    model.tolerance = fva_tol # prevent essential EX_meoh_ex from being blocked, sets solver feasibility/optimality tolerances
    # model.solver.configuration.tolerances.feasibility = 1e-9
    # model.solver = 'glpk' #'glpk_exact' # appears to make problems for context management
    # model_mue_target_constraints= get_leq_constraints(model, ecc2_mue_target)
    # model.add_cons_vars(model_mue_target_constraints[0])
    fva_res = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0, processes=1) # no interactive multiprocessing on Windows
print(fva_res.loc['EX_adp_c',:])
print(fva_res.loc['EX_AMOB_ex']) # !!! with GLPK significantly different compared to the copied model, OK with CPLEX
print(fva_res.loc['EX_meoh_ex',:])
blocked = []
blocked_rxns = []
for i in range(fva_res.values.shape[0]):
    # if fva_res.values[i, 0] == 0 and fva_res.values[i, 1] == 0:
    if fva_res.values[i, 0] >= -fva_tol and fva_res.values[i, 1] <= fva_tol:
        blocked.append(i)
        blocked_rxns.append(fva_res.index[i])
print(blocked_rxns)

# %% FVA with try/finally model restoration
model = ecc2
previous_tolerance = model.tolerance
previous_objective = model.objective
try:
    model.objective = model.problem.Objective(0)
    fva_tol = 1e-9 # with CPLEX 1e-8 leads to removal of EX_adp_c, 1e-9 keeps EX_adp_c
    model.tolerance = fva_tol # prevent essential EX_meoh_ex from being blocked, sets solver feasibility/optimality tolerances
    fva_res = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0, processes=1) # no interactive multiprocessing on Windows
finally:
    model.tolerance = previous_tolerance
    model.objective = previous_objective
print(fva_res.loc['EX_adp_c',:])
print(fva_res.loc['EX_AMOB_ex']) # !!! also incorrect ?!?
blocked = []
blocked_rxns = []
for i in range(fva_res.values.shape[0]):
    # if fva_res.values[i, 0] == 0 and fva_res.values[i, 1] == 0:
    if fva_res.values[i, 0] >= -fva_tol and fva_res.values[i, 1] <= fva_tol:
        blocked.append(i)
        blocked_rxns.append(fva_res.index[i])
print(blocked_rxns)

#%% only glpk_exact recognizes EX_adp_c as essential, but knocking it out is
# not a problem with CPLEX as its flux is below the minimal tolerance
# however compression with rationals will not work properly when EX_adp_c is removed
# because then growth will appear to be blocked
print(ecc2.slim_optimize())
for i in blocked:
    with ecc2 as model:
        model_mue_target_constraints= get_leq_constraints(model, ecc2_mue_target)
        model.add_cons_vars(model_mue_target_constraints[0])
        model.reactions[i].knock_out()
        print(model.reactions[i].id, model.slim_optimize())

# %% try FASTCC
model = ecc2.copy() # copy model because switching solver in context gives an error (?!?)
model.tolerance= 1e-9
model.objective = model.problem.Objective(0)
model.solver = 'glpk_exact' # appears to use CPLEX anyway?
# modelcc = cobra.flux_analysis.fastcc(ecc2, zero_cutoff=1e-9) # gives some falsely blocked reactions
modelcc = cobra.flux_analysis.fastcc(model, flux_threshold=10, zero_cutoff=1e-9) # still identifies EX_adp_c as blocked
del model
set(ecc2.reactions.list_attr("id")) - set(modelcc.reactions.list_attr("id"))

#%% FVA mit verschiedenen solvern machen um zu sehen ob {'EX_adp_c'} ein cut set ist (boundary müssen cuts sein)
ecc2.solver = 'glpk_exact'
with ecc2 as model:
    model.objective = model.problem.Objective(0)
    #print(model.objective.expression)
    flux_expr= [r.flux_expression for r in model.reactions] # !! lose validity when the solver is changed !!
    rexpr = matrix_row_expressions(ecc2_mue_target[0][0], flux_expr)
    model.add_cons_vars(leq_constraints(model.problem.Constraint, rexpr, ecc2_mue_target[0][1]))
    res = cobra.flux_analysis.flux_variability_analysis(model, [ecc2.reactions.get_by_id(r[0]) for r in single_cuts], fraction_of_optimum=0, processes=1) # no interactive multiprocessing on Windows
print(res.loc['EX_adp_c',:])
# CPLEX: 
# minimum    2.042810e-13
# maximum    9.059420e-14
# Name: EX_adp_c, dtype: float64
# glpk_exact
# minimum   -1.117471e-09
# maximum   -4.705152e-12
# Name: EX_adp_c, dtype: float64
numpy.linalg.matrix_rank(ecc2_stdf.values) # no conservation relations

#%%
e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, ecc2_stdf.values, [r.reversibility for r in ecc2.reactions], ecc2_mue_target,
                                        cuts=cuts, bigM=1000, threshold=0.1, split_reversible_v=True, irrev_geq=True)
# e.model.objective = e.minimize_sum_over_z
#e.write_lp_file('testM')
#e.model.configuration.verbosity = 3
e.model.configuration._iocp.mip_gap = 0.99
ecc2_mcsB,_ = e.enumerate_mcs(max_mcs_size=3, enum_method=3, model=ecc2, targets=ecc2_mue_target)
print(set(ecc2_mcs) == set(ecc2_mcsB), len(ecc2_mcsB))

# %% sympy and efmtool conversions to rational are not necessarily the same, nsimplify appears to match
# v = 1/3
# v = 0.1
v = 0.2
r1 = sympy.Rational(v)
r2 = sympy.nsimplify(v, rational=True) # same as sympy.Rational with rational_conversion='exact'
print(r1, r2, abs(float(r1-r2)))
efmtool4cobra.jRatMat2sympyRatMat(efmtool_intern.numpy_mat2jBigIntegerRationalMatrix(numpy.array([[v]]))).values()

#%%
# not for enum_method 3 because this requires as compressed model for the minimality checks
rev = [r.reversibility for r in ecc2.reactions]
# rd, subT = efmtool_intern.compress_rat_efmtool(ecc2_stdf.values, rev, remove_cr=True, # if CR are not removed problem with enumeration
#             compression_method=efmtool_intern.subset_compression)[0:2]
rd, subT = efmtool_intern.compress_rat_efmtool(ecc2_stdf.values, rev, remove_cr=True, # if CR are not removed problem with enumeration
            compression_method=efmtool_intern.subset_compression, remove_rxns=blocked)[0:2] # OK when EX_adp_c is removed
rev_rd = numpy.logical_not(numpy.any(subT[numpy.logical_not(rev), :], axis=0))
kn = efmtool_intern.null_rat_efmtool(rd)
print(rd.shape)
print(numpy.linalg.matrix_rank(rd))
print(kn.shape)
print(numpy.max(abs(rd@kn)))

#%%
# compression_tolerance = numpy.finfo(float).eps # does not work
compression_tolerance = 1e-10 # OK
# compression_tolerance = 1e-9 # misses one CR
# ecc2c, subT = efmtool4cobra.compress_model(ecc2)
ecc2c, subT = efmtool4cobra.compress_model(ecc2, remove_rxns=blocked_rxns, tolerance=compression_tolerance) # tolerance=0 here also OK
rev_rd = [r.reversibility for r in ecc2c.reactions]
# rd = cobra.util.array.create_stoichiometric_matrix(ecc2c, array_type='dok')
# bc = efmtool_intern.basic_columns_rat(rd.transpose().toarray(), tolerance=1e-12) # needs non-zero tolerance or...
# rd = rd[numpy.sort(bc), :] # ...it misses one CR
efmtool4cobra.remove_conservation_relations(ecc2c, tolerance=compression_tolerance)
rd = cobra.util.array.create_stoichiometric_matrix(ecc2c, array_type='dok')
print(numpy.sort(list(map(abs, rd.values())))[[0, -1]])
kn = efmtool_intern.null_rat_efmtool(rd)
print(rd.shape)
print(numpy.linalg.matrix_rank(rd.toarray()))
print(kn.shape)
print(numpy.max(abs(rd@kn)))

#%%
ecc2c = ecc2.copy()
# subT = efmtool4cobra.compress_model_sympy(ecc2)
subT = efmtool4cobra.compress_model_sympy(ecc2c, remove_rxns=blocked_rxns)
print(len(ecc2c.metabolites), len(ecc2c.reactions))
rev_rd = [r.reversibility for r in ecc2c.reactions]
efmtool4cobra.remove_conservation_relations_sympy(ecc2c)
reduced = cobra.util.array.create_stoichiometric_matrix(ecc2c, array_type='dok', dtype=numpy.object)
# reduced, bc = efmtool4cobra.remove_conservation_relations_sympy(ecc2c, return_reduced_only=True)[0:2]
# rdb = scipy.sparse.dok_matrix((reduced.shape[0], reduced.shape[1]))
# for (r, c), v in reduced.items():
#     rdb[r, c] = float(v)
# rd = rdb
rd = efmtool4cobra.dokRatMat2lilFloatMat(reduced)
# for m in [ecc2c.metabolites[i].id for i in set(range(len(ecc2c.metabolites))) - set(bc)]:
#     ecc2c.metabolites.get_by_id(m).remove_from_model()
#%%
# kn = efmtool_intern.null_rat_efmtool(rd)
import ch.javasoft.smx.ops.Gauss as Gauss
jkn = Gauss.getRationalInstance().nullspace(efmtool4cobra.sympyRatMat2jRatMat(reduced))
kn = efmtool_intern.jpypeArrayOfArrays2numpy_mat(jkn.getDoubleRows())
print(rd.shape)
print(numpy.linalg.matrix_rank(rd.toarray())) # toarray() makes it a full matrix, otherwise weird result with scipy sparse
print(kn.shape) # !! is incomplete when tolerance=0
print(numpy.max(abs(rd@kn)))
#%% cannot get scipy sparse matrix multiplication to work, probably because the matrices have object dtype
import sympy.matrices.sparse
srd = sympy.matrices.sparse.MutableSparseMatrix(reduced.shape[0], reduced.shape[1], reduced.items())
# srd = sympy.matrices.SparseMatrix(reduced.shape[0], reduced.shape[1], dict(reduced.items()))
# skn = srd.nullspace() # extremely slow
skn = efmtool4cobra.jRatMat2sympyRatMat(jkn)
sympy.matrices.sparse.MutableSparseMatrix(skn.shape[0], skn.shape[1], skn.items())
max(abs(srd@skn))
#%%
#ecc2.solver = 'glpk_exact'
#ecc2c.solver = 'glpk_exact'
sol = ecc2.optimize()
growth_idx = ecc2.reactions.index('Growth')
growth_idx_c = numpy.where(subT[growth_idx, :])[0][0]
ecc2c.objective = ecc2c.reactions[growth_idx_c]
sol_c = ecc2c.optimize()
print(sol, sol_c, abs(sol.objective_value - sol_c.objective_value))
# # %%
# for carbon_uptake in ['AcUp', 'GlycUp', 'SuccUp', 'GlcUp']:
#     idx = ecc2.reactions.index(carbon_uptake)
#     idx_c = numpy.where(subT[idx, :])[0][0]
#     print(sol.fluxes[carbon_uptake], sol_c.fluxes[ecc2c.reactions[idx_c].id], subT[idx, idx_c])
# # at least the bound on GlcUp is lost

# #%%
# print(ecc2_stdf.values.shape) # has 40 extra exchange reactions for the 40 external metabolites in ECC2comp
# print(ecc2_stdf.columns[0])
# ecc2_stdf.columns[numpy.where(subT[:, 0])]

#%%
target_rd = [(T@subT, t) for T, t in ecc2_mue_target]

# e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, rd, rev_rd, target_rd, 
#                                         threshold=0.1, split_reversible_v=True, irrev_geq=True,
#                                         cuts=numpy.any(subT[cuts, :], axis=0))
# e.model.problem.parameters.mip.tolerances.mipgap.set(0.99)

e = ConstrainedMinimalCutSetsEnumerator(optlang.coinor_cbc_interface, rd, rev_rd, target_rd, 
                                        threshold=0.1, bigM=1000, split_reversible_v=True, irrev_geq=True,
                                        cuts=numpy.any(subT[cuts, :], axis=0))#, kn=kn) #efmtool_intern.null_rat_efmtool(rd))
#e.model.problem.max_mip_gap = 0.99 # does not appear to make the solver stop early
e.model.problem.max_solutions = 1 # stops with feasible solutions
e.model.problem.threads = -1 # default does not appear to use multi-threading
#e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, rd, rev_rd, target_rd, 
#                                        threshold=0.1, bigM=1000, split_reversible_v=True, #irrev_geq=True,
#                                        cuts=numpy.any(subT[cuts, :], axis=0), kn=kn) #efmtool_intern.null_rat_efmtool(rd))
#e.model.configuration._iocp.mip_gap = 0.99

# here a subset is repressible when one ot its reactions is repressible
# e.model.objective = e.minimize_sum_over_z
info = dict()
# rd_mcs,_ = e.enumerate_mcs(max_mcs_size=3, info=info)
# rd_mcs,_ = e.enumerate_mcs(max_mcs_size=3, enum_method=3, model=ecc2c, targets=target_rd, info=info)
# with ecc2c as tm:
    # tm.solver = 'glpk' #'glpk_exact' # actually optlang runs regular GLPK first and only if the results is optimal glpk_exact
    # tm.solver.configuration.verbosity = 3
rd_mcs,_ = e.enumerate_mcs(max_mcs_size=3, enum_method=3, model=ecc2c, targets=target_rd, info=info)
print(info)
print(len(rd_mcs))
xsubT= subT.copy()
xsubT[numpy.logical_not(cuts), :] = 0 # only expand to reactions that are repressible 
xmcs = expand_mcs(rd_mcs, xsubT)
set(xmcs) == set(ecc2_mcs)

#%%
from swiglpk import *
test_model = ecc2c.copy()
test_model.solver = 'glpk_exact'
target_constraints= get_leq_constraints(test_model, target_rd)
print(make_minimal_cut_set(test_model, (9,38,64), target_constraints)) # should give (38, 64), looks like GLPK has warm start issues
test_model.add_cons_vars(target_constraints[0])
# test_model.reactions[9].knock_out()
test_model.reactions[38].knock_out()
test_model.reactions[64].knock_out()
return_value = glp_exact(test_model.solver.problem, test_model.solver.configuration._smcp)
glpk_status = glp_get_status(test_model.solver.problem)
print(return_value, glpk_status)
print(test_model.solver._run_glp_exact())
test_model2 = test_model.copy() # with copy
return_value = glp_exact(test_model2.solver.problem, test_model2.solver.configuration._smcp)
glpk_status = glp_get_status(test_model2.solver.problem)
print(return_value, glpk_status)
print(test_model2.solver._run_glp_exact())
test_model3 = glp_create_prob()
glp_copy_prob(test_model3, test_model.solver.problem, GLP_OFF)
# glp_std_basis(test_model3) # without this the old basis is still used
glp_adv_basis(test_model3, 0) # without this the old basis is still used
return_value = glp_exact(test_model3, test_model.solver.configuration._smcp)
glpk_status = glp_get_status(test_model3)
print(return_value, glpk_status)

# %%
from swiglpk import *
test_model = ecc2c.copy()
# test_model = ecc2.copy()
test_model.solver = 'glpk_exact' # actually optlang runs regular GLPK first and only if the results is optimal glpk_exact
test_model.objective = test_model.reactions.EX_Biomass
print(test_model.objective)
print(test_model.solver._run_glp_exact())
ov = test_model.optimize()
print(ov)
# print(test_model.solver._run_glp_exact())
return_value = glp_exact(test_model.solver.problem, test_model.solver.configuration._smcp)
glpk_status = glp_get_status(test_model.solver.problem)
print(return_value, glpk_status)
target_constraints= get_leq_constraints(test_model, target_rd)
print(make_minimal_cut_set(test_model, (9,14,49), target_constraints)) # looks like GLPK has warm start issues
test_model.add_cons_vars(target_constraints[0])
# print(test_model.solver._run_glp_exact()) # -> undefined
return_value = glp_exact(test_model.solver.problem, test_model.solver.configuration._smcp)
glpk_status = glp_get_status(test_model.solver.problem)
print(return_value, glpk_status)
print(test_model.solver._run_glp_exact())
test_model.reactions[14].knock_out()
test_model.reactions[49].knock_out()
print(test_model.slim_optimize(), test_model.solver.status)

#%% FVA in reduced system

# %%
num_reac = len(ecc2c.reactions)
j_max = numpy.full((num_reac, num_reac), numpy.nan)
j_min = numpy.full((num_reac, num_reac), numpy.nan)
j_min_max = numpy.full((num_reac, num_reac), numpy.nan)
for i in range(num_reac):
    with ecc2c as mKOi:
        mKOi.reactions[i].knock_out()
        for j in range(i+1, num_reac): # hier doch über alle Reaktionen (ausser i)
            mKOi.objective = mKOi.reactions[j]
            mKOi.objective.direction = "max"
            j_max[i, j] = mKOi.slim_optimize()
            j_min_max[i, j] = mKOi.slim_optimize()
            mKOi.objective.direction = "min"
            j_min[i, j] = mKOi.slim_optimize()
            j_min_max[j, i] = mKOi.slim_optimize()
            if j_min_max[i, j] == j_min_max[j, i] and j_min_max[j, i] == 0:
                print(i, j)

#%%
iJO1366 = cobra.io.read_sbml_model(r"../cnapy-projects/iJO1366/model.sbml")
# iJO1366.solver = 'glpk_exact'
# %% full FVA
fva_tol = 1e-9
model = iJO1366.copy()
# model.solver = 'cplex' # 14.1 s, 878 blocked
# model.solver = 'coinor_cbc' # 128.3 s, 880 blocked -> probably does not warm start
model.solver = 'glpk' # 26.6 seconds, 879 blocked,  would be extremely slow with glpk_exact 
# model.objective = model.problem.Objective(0)
model.tolerance = fva_tol
start_time = time.monotonic()
fva_res = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0, processes=1) # no interactive multiprocessing on Windows
# fva_res = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0.99, processes=1) # fraction 1.0 may lead to too tight optimum constraint
print(time.monotonic() - start_time) # 20 seconds in iJO1366 without remove_rxns
del model
blocked_rxns = []
for i in range(fva_res.values.shape[0]):
    # if fva_res.values[i, 0] == 0 and fva_res.values[i, 1] == 0:
    if fva_res.values[i, 0] >= -fva_tol and fva_res.values[i, 1] <= fva_tol:
        blocked_rxns.append(fva_res.index[i])
print(len(blocked_rxns))
# %% problematic result close to zero near the solver tolerance
zero_threshold = fva_tol/100
nonzero_threshold = fva_tol*100
abs_fva_res = numpy.abs(fva_res.values)
idx = numpy.where(numpy.logical_or(numpy.logical_and(abs_fva_res[:, 0] < nonzero_threshold, abs_fva_res[:, 0] >= zero_threshold),
        numpy.logical_and(abs_fva_res[:, 1] < nonzero_threshold, abs_fva_res[:, 1] >= zero_threshold)))[0]
print(len(idx))
# %%
model = iJO1366.copy()
model.solver = 'glpk_exact'
# model.objective = model.problem.Objective(0)
model.tolerance = fva_tol
fva_res_gx = cobra.flux_analysis.flux_variability_analysis(model, reaction_list=[model.reactions[i] for i in idx], fraction_of_optimum=0.99, processes=1)
numpy.max(numpy.abs(fva_res_gx.values - fva_res.values[idx, :])) # no meaningful accuracy gain here

#%%
iJO1366c = iJO1366.copy()
# subT = efmtool4cobra.compress_model_sympy(iJO1366c)
subT = efmtool4cobra.compress_model_sympy(iJO1366c, blocked_rxns)
# %% save with rational stoichiometric coefficients
f = open("iJO1366c_subT.pkl","wb")
pickle.dump((cobra.io.model_to_dict(iJO1366c), subT),f)
f.close()
#%% loading back appears to work correctly
with open("iJO1366c_subT.pkl", "rb") as f:
    m, subT = pickle.load(f)
    iJO1366c = cobra.io.model_from_dict(m)
iJO1366c.solver = "glpk_exact" # solver type appears not to be restored

#%%
sol = iJO1366.optimize()
sol_c = iJO1366c.optimize()
print(sol, sol_c, abs(sol.objective_value - sol_c.objective_value))

#%%
rev_rd = [r.reversibility for r in iJO1366c.reactions]
efmtool4cobra.remove_conservation_relations_sympy(iJO1366c)
reduced = cobra.util.array.create_stoichiometric_matrix(iJO1366c, array_type='dok', dtype=numpy.object)
rd = efmtool4cobra.dokRatMat2lilFloatMat(reduced)
# reduced, bc = efmtool4cobra.remove_conservation_relations_sympy(iJO1366c)[0:2]
# rdb = scipy.sparse.dok_matrix((reduced.shape[0], reduced.shape[1]))
# for (r, c), v in reduced.items():
#     rdb[r, c] = float(v)
# rd = rdb
# for m in [iJO1366c.metabolites[i].id for i in set(range(len(iJO1366c.metabolites))) - set(bc)]:
#     iJO1366c.metabolites.get_by_id(m).remove_from_model()
sol_c = iJO1366c.optimize()
print(sol, sol_c, abs(sol.objective_value - sol_c.objective_value))

#%%
# kn = efmtool_intern.null_rat_efmtool(rd)
import ch.javasoft.smx.ops.Gauss as Gauss
jkn = Gauss.getRationalInstance().nullspace(efmtool4cobra.sympyRatMat2jRatMat(reduced))
kn = efmtool_intern.jpypeArrayOfArrays2numpy_mat(jkn.getDoubleRows())
print(rd.shape)
print(numpy.linalg.matrix_rank(rd.toarray())) # toarray() makes it a full matrix, otherwise weird result with scipy sparse
print(kn.shape) # !! is incomplete when tolerance=0
print(numpy.max(abs(rd@kn)))

#%%
iJO1366_stdf = cobra.util.array.create_stoichiometric_matrix(iJO1366, array_type='DataFrame')
cuts= numpy.full(iJO1366_stdf.shape[1], True, dtype=bool)
for r in iJO1366.boundary:
    cuts[iJO1366_stdf.columns.get_loc(r.id)] = False
#cuts = None
# iJO1366_mue_target = [(equations_to_matrix(iJO1366, 
#                     ["-1 BIOMASS_Ec_iJO1366_core_53p95M", "-1 EX_glc__D_e"]), [-0.01, 10])]
# iJO1366_mue_target_constraints= get_leq_constraints(iJO1366, iJO1366_mue_target)
reac_id = iJO1366.reactions.list_attr("id")
reac_id_symbols = get_reac_id_symbols(reac_id)
iJO1366_mue_target = [[("BIOMASS_Ec_iJO1366_core_53p95M", ">=", 0.01)]]
iJO1366_mue_target = [relations2leq_matrix(parse_relations(t, reac_id_symbols=reac_id_symbols), reac_id) for t in iJO1366_mue_target]
bounds_mat, bounds_rhs = reaction_bounds_to_leq_matrix(iJO1366)
iJO1366_mue_targetB = [(scipy.sparse.vstack((t[0], bounds_mat), format='csr'), numpy.hstack((t[1], bounds_rhs))) for t in iJO1366_mue_target]
iJO1366_mue_target_constraintsB= get_leq_constraints(iJO1366, iJO1366_mue_targetB)
# %% 
iJO1366_mcs,_ = compute_mcs(iJO1366, iJO1366_mue_target, cuts=cuts, enum_method=2, max_mcs_size=1, network_compression=True)
print(len(iJO1366_mcs))
# %% 
iJO1366_mcs,_ = compute_mcs(iJO1366, iJO1366_mue_target, cuts=cuts, enum_method=3, max_mcs_size=1, network_compression=True)
print(len(iJO1366_mcs))
# %%
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, iJO1366_stdf.values, [r.reversibility for r in iJO1366.reactions], iJO1366_mue_targetB,
                                        cuts=cuts, threshold=0.1, split_reversible_v=False, irrev_geq=False)
# e.model.objective = e.minimize_sum_over_z
#e.write_lp_file('testI')
# e.model.configuration.verbosity = 3
iJO1366_mcs,_ = e.enumerate_mcs(max_mcs_size=1, enum_method=2)
iJO1366_mcs_rxns = [tuple(reac_id[r] for r in mcs) for mcs in iJO1366_mcs]
print(len(iJO1366_mcs)) # 271
print(all(check_mcs(iJO1366, iJO1366_mue_targetB[0], iJO1366_mcs, optlang.interface.INFEASIBLE)))
# !!!!!!                                      ^^^ really need the implicit bounds in the target
# otherwise ATPM as KO would be missed in this example because its >= 3.15 requirement would vanisch through the KO in the check
#%% misses some cuts with CPLEX, OK with glpk_exact
candidates = [r for r, c in zip(iJO1366.reactions.list_attr("id"), cuts) if c]
with iJO1366 as model:
    model.objective = model.problem.Objective(0)
    # model.solver = 'glpk_exact' # again problems with context
    model_target_constraints= get_leq_constraints(model, iJO1366_mue_targetB)
    model.add_cons_vars(model_target_constraints[0])
    resD = cobra.flux_analysis.single_reaction_deletion(model, reaction_list=candidates,
                                                         processes=1) # no interactive multiprocessing on Windows
    single_cuts= list(map(tuple, resD.index[resD.values[:, 1] == optlang.interface.INFEASIBLE]))
    print(len(single_cuts)) # 255 ??
    print(set(iJO1366_mcs_rxns) - set(single_cuts))
# print(set(single_cuts) == set(mcs for mcs in iJO1366_mcs_rxns if len(mcs) == 1))
#%% wird auch infeasible wenn die nächste cell vorher ausgeführt wird?!?!??!
# muss was mit dem internen solver status zu tun haben da das Modell sich nicht ändert (str(cobra.io.model_to_dict(iJO1366) bleibt identisch)
# auch infeasible wenn der solver nach dem Laden auf 'glpk_exact' umgestellt wird
with iJO1366 as model:
    model.objective = model.problem.Objective(0)
    model_target_constraints= get_leq_constraints(model, iJO1366_mue_targetB)
    model.add_cons_vars(model_target_constraints[0])
    model.reactions.DBTS.knock_out()
    res = model.optimize()
    print(res, res.status)
#%% 
check_mcs(iJO1366, iJO1366_mue_targetB[0], [(iJO1366.reactions.index('DBTS'),)], optlang.interface.INFEASIBLE)

#%% 
model = iJO1366
model.objective = model.problem.Objective(0)
model_target_constraints= get_leq_constraints(model, iJO1366_mue_targetB)
model.add_cons_vars(model_target_constraints[0])
model.reactions.DBTS.knock_out()
res = model.optimize()
print(res, res.status)
#%% misses some cuts
candidates = [r for r, c in zip(iJO1366.reactions.list_attr("id"), cuts) if c]
model = iJO1366
model.objective = model.problem.Objective(0)
# model.solver = 'glpk_exact' # again problems with context
model_target_constraints= get_leq_constraints(model, iJO1366_mue_targetB)
model.add_cons_vars(model_target_constraints[0])
resD = cobra.flux_analysis.single_reaction_deletion(model, reaction_list=candidates,
                                                        processes=1) # no interactive multiprocessing on Windows
single_cuts= list(map(tuple, resD.index[resD.values[:, 1] == optlang.interface.INFEASIBLE]))
print(len(single_cuts)) # 255 ??
print(set(iJO1366_mcs_rxns) - set(single_cuts))
# print(set(single_cuts) == set(mcs for mcs in iJO1366_mcs_rxns if len(mcs) == 1))

#%% OK
candidates = [(c,) for c in numpy.where(cuts)[0]]
res = check_mcs(iJO1366, iJO1366_mue_targetB[0], candidates, optlang.interface.INFEASIBLE)
print(sum(res == True))
m = set(iJO1366_mcs) - set([c for c, b in zip(candidates, res) if b])
print(len(m))

# %%
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, iJO1366_stdf.values, [r.reversibility for r in iJO1366.reactions], iJO1366_mue_target,
                                        cuts=cuts, bigM=10000, threshold=0.0001, split_reversible_v=False, irrev_geq=False)
# e.model.objective = e.minimize_sum_over_z
#e.write_lp_file('testM')
# e.model.configuration.verbosity = 3
e.model.configuration.tolerances.feasibility = 1e-9
e.model.configuration.tolerances.optimality = 1e-9
e.model.configuration.tolerances.integrality = 1e-10
iJO1366_mcsB,_ = e.enumerate_mcs(max_mcs_size=1)
print(set(iJO1366_mcs) == set(iJO1366_mcsB), len(iJO1366_mcsB))

# %%
target_rd = [(T@subT, t) for T, t in iJO1366_mue_targetB]
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, rd, rev_rd, target_rd, 
                                        threshold=0.1, split_reversible_v=True, irrev_geq=True,
                                        cuts=numpy.any(subT[cuts, :], axis=0))
info = dict()
rd_mcs,_ = e.enumerate_mcs(max_mcs_size=1, enum_method=2, info=info)
print(info)
print(len(rd_mcs))
xsubT= subT.copy()
xsubT[numpy.logical_not(cuts), :] = 0 # only expand to reactions that are repressible 
xmcs = expand_mcs(rd_mcs, xsubT)
set(xmcs) == set(iJO1366_mcs)

#%% OK, a bit slower than CPLEX populate
start_time = time.monotonic()
candidates = [(c,) for c in numpy.where(numpy.any(subT[cuts, :], axis=0))[0]]
res = check_mcs(iJO1366c, target_rd[0], candidates, optlang.interface.INFEASIBLE)
print(time.monotonic() - start_time)
print(sum(res == True))
m = set(rd_mcs) - set([c for c, b in zip(candidates, res) if b])
print(len(m))

#
# %% does not work with bigM out of the box
# e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, rd, rev_rd, target_rd, 
e = ConstrainedMinimalCutSetsEnumerator(optlang.cplex_interface, rd, rev_rd, target_rd, 
                                        threshold=0.1, bigM=10000, split_reversible_v=True, irrev_geq=True,
                                        cuts=numpy.any(subT[cuts, :], axis=0))
rd_mcs,_ = e.enumerate_mcs(max_mcs_size=1, enum_method=1)

#%% GLPK struggles with numerical instability 
e = ConstrainedMinimalCutSetsEnumerator(optlang.glpk_interface, iJO1366_stdf.values, [r.reversibility for r in iJO1366.reactions], iJO1366_mue_target,
                                        cuts=cuts, bigM=1000, threshold=0.001, split_reversible_v=False, irrev_geq=False)
e.model.objective = e.minimize_sum_over_z
#e.write_lp_file('testM')
e.model.configuration.verbosity = 3
e.model.configuration.tolerances.feasibility = 1e-8
e.model.configuration._iocp.tol_obj= 1e-8
e.model.configuration._iocp.tol_int= 1e-10
iJO1366_mcsB,_ = e.enumerate_mcs(max_mcs_size=1)
print(set(iJO1366_mcs) == set(iJO1366_mcsB), len(iJO1366_mcsB))


# %%
print(len(iJO1366.genes), sum(iJO1366.genes.list_attr('functional')))
# %%
ast = cobra.core.gene.parse_gpr(iJO1366.reactions[436].gene_reaction_rule)
cobra.core.gene.eval_gpr(ast[0], {'b0351'})
cobra.core.gene.eval_gpr(ast[0], ast[1])
#sympy.sympify('b0351 | b1241')
#sympy.logic.boolalg.to_dnf(sympy.sympify('(b4213) & (b0351 | b1241)'))
# %%

# %%
iJO1366 = cobra.io.read_sbml_model(r"..\cnapy-projects\iJO1366\model.sbml")

# %% see how much memory resources GLPK needs...
from swiglpk import *
test_model = iJO1366.copy()
# test_model = ecc2.copy()
test_model.solver = 'glpk_exact'
test_model.objective = test_model.problem.Objective(0.0)
print(test_model.objective)
start_time = time.monotonic()
#print(test_model.solver._run_glp_exact()) # to make the basis
print(test_model.slim_optimize()) # to make the basis
print("Initial simplex", time.monotonic() - start_time) # iJO1366 glp_exact: 56 s, slim_optimize with glp_simplex: 0.1 s
num_models = 1000 # about a 1.1 GB python process (10000 for ECC2); 2.3 GB for 1000 iJO1366
model_collection = [None] * num_models
return_value = [None] * num_models
glpk_status = [None] * num_models
start_time = time.monotonic()
for i in range(num_models):
    model_collection[i] = glp_create_prob()
    glp_copy_prob(model_collection[i], test_model.solver.problem, GLP_OFF)
print("Copied models in", time.monotonic() - start_time)
start_time = time.monotonic()
for i in range(num_models):
    return_value[i] = glp_simplex(model_collection[i], test_model.solver.configuration._smcp)
    return_value[i] = glp_exact(model_collection[i], test_model.solver.configuration._smcp)
    glpk_status[i] = glp_get_status(model_collection[i])
print(time.monotonic() - start_time) # 9 seconds with ECC2, 20 with iJO1366


# %%