"""Implements optimization and model problems."""frommicom.dualityimportfast_dualfrommicom.utilimport(_format_min_growth,_apply_min_growth,check_modification,)frommicom.loggerimportloggerfrommicom.solutionimportsolvefromoptlang.symbolicsimportZerofromfunctoolsimportpartial
[docs]defadd_dualized_optcom(community,min_growth):"""Add dual Optcom variables and constraints to a community. Uses the original formulation of OptCom and solves the following multi-objective problem:: maximize community_growth s.t. maximize growth_rate_i for all i s.t. Sv_i = 0 lb_i >= v_i >= ub_i Notes ----- This method will only find one arbitrary solution from the Pareto front. There may exist several other optimal solutions. Arguments --------- community : micom.Community The community to modify. min_growth : positive float or array-like object. The minimum growth rate for each individual in the community. Either a single value applied to all individuals or one value for each. """logger.info("adding dual optcom to %s"%community.id)check_modification(community)min_growth=_format_min_growth(min_growth,community.taxa)prob=community.solver.interface# Temporarily subtitute objective with sum of individual objectives# for correct dual variablesold_obj=community.objectivecommunity.objective=Zeroforspincommunity.taxa:taxa_obj=community.constraints["objective_"+sp]community.objective+=taxa_obj.expression_apply_min_growth(community,min_growth)dual_coefs=fast_dual(community)logger.info("adding expressions for %d taxa"%len(community.taxa))forspincommunity.taxa:primal_const=community.constraints["objective_"+sp]coefs=primal_const.get_linear_coefficients(primal_const.variables)coefs.update({dual_var:-coeffordual_var,coefindual_coefs.items()ifspindual_var.name})obj_constraint=prob.Constraint(Zero,lb=0,ub=0,name="optcom_suboptimality_"+sp)community.add_cons_vars([obj_constraint])community.solver.update()obj_constraint.set_linear_coefficients(coefs)community.objective=old_objcommunity.modification="dual optcom"logger.info("finished adding dual optcom to %s"%community.id)
[docs]defadd_moma_optcom(community,min_growth,linear=False):"""Add a dualized MOMA version of OptCom. Solves a MOMA (minimization of metabolic adjustment) formulation of OptCom given by:: minimize cooperativity_cost s.t. maximize community_objective s.t. Sv = 0 lb >= v >= ub where community_cost = sum (growth_rate - max_growth)**2 if linear=False or community_cost = sum |growth_rate - max_growth| if linear=True Arguments --------- community : micom.Community The community to modify. min_growth : positive float or array-like object. The minimum growth rate for each individual in the community. Either a single value applied to all individuals or one value for each. linear : boolean Whether to use a non-linear (sum of squares) or linear version of the cooperativity cost. If set to False requires a QP-capable solver. """logger.info("adding dual %s moma to %s"%("linear"iflinearelse"quadratic",community.id))check_modification(community)min_growth=_format_min_growth(min_growth,community.taxa)prob=community.solver.interfaceold_obj=community.objectivecoefs=old_obj.get_linear_coefficients(old_obj.variables)# Get maximum individual growth ratesmax_gcs=community.optimize_all(progress=False)_apply_min_growth(community,min_growth)dual_coefs=fast_dual(community)coefs.update({v:-coefforv,coefindual_coefs.items()})obj_constraint=prob.Constraint(Zero,lb=0,ub=0,name="optcom_suboptimality")community.add_cons_vars([obj_constraint])community.solver.update()obj_constraint.set_linear_coefficients(coefs)obj_expr=Zerologger.info("adding expressions for %d taxa"%len(community.taxa))forspincommunity.taxa:v=prob.Variable("gc_constant_"+sp,lb=max_gcs[sp],ub=max_gcs[sp])community.add_cons_vars([v])taxa_obj=community.constraints["objective_"+sp]ex=v-taxa_obj.expressionifnotlinear:ex=ex**2obj_expr+=ex.expand()community.objective=prob.Objective(obj_expr,direction="min")community.modification="moma optcom"logger.info("finished dual moma to %s"%community.id)
[docs]defoptcom(community,strategy,min_growth,fluxes,pfba):"""Run OptCom for the community. OptCom methods are a group of optimization procedures to find community solutions that provide a tradeoff between the cooperative community growth and the egoistic growth of each individual [#p1]_. `micom` provides several strategies that can be used to find optimal solutions: - "moma": Minimization of metabolic adjustment. Simultaneously optimizes the community objective (maximize) and the cooperativity cost (minimize). This method finds an exact maximum but doubles the number of required variables, thus being slow. - "lmoma": The same as "moma" only with a linear representation of the cooperativity cost (absolute value). - "original": Solves the multi-objective problem described in [#p1]_. Here, the community growth rate is maximized simultanously with all individual growth rates. Note that there are usually many Pareto-optimal solutions to this problem and the method will only give one solution. This is also the slowest method. Parameters ---------- community : micom.Community The community to optimize. strategy : str The strategy used to solve the OptCom formulation. Defaults to "lagrangian" which gives a decent tradeoff between speed and correctness. min_growth : float or array-like The minimal growth rate required for each individual. May be a single value or an array-like object with the same length as there are individuals. fluxes : boolean Whether to return the fluxes as well. pfba : boolean Whether to obtain fluxes by parsimonious FBA rather than "classical" FBA. Returns ------- micom.CommunitySolution The solution of the optimization. If fluxes==False will only contain the objective value, community growth rate and individual growth rates. References ---------- .. [#p1] OptCom: a multi-level optimization framework for the metabolic modeling and analysis of microbial communities. Zomorrodi AR, Maranas CD. PLoS Comput Biol. 2012 Feb;8(2):e1002363. doi: 10.1371/journal.pcbi.1002363, PMID: 22319433 """ifstrategynotin_methods:raiseValueError("strategy must be one of {}!".format(",".join(_methods)))funcs=_methods[strategy]withcommunityascom:# Add needed variables etc.funcs[0](com,min_growth)returnsolve(community,fluxes=fluxes,pfba=pfba)