{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solver Interface\n", "\n", "Each cobrapy solver must expose the following API. The solvers all will have their own distinct LP object types, but each can be manipulated by these functions. This API can be used directly when implementing algorithms efficiently on linear programs because it has 2 primary benefits:\n", "\n", "1. Avoid the overhead of creating and destroying LP's for each operation\n", "\n", "2. Many solver objects preserve the basis between subsequent LP's, making each subsequent LP solve faster\n", "\n", "We will walk though the API with the cglpk solver, which links the cobrapy solver API with [GLPK](http://www.gnu.org/software/glpk/)'s C API." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import cobra.test\n", "\n", "model = cobra.test.create_test_model(\"textbook\")\n", "solver = cobra.solvers.cglpk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Attributes and functions\n", "\n", "Each solver has some attributes:\n", "\n", "### solver_name\n", "\n", "The name of the solver. This is the name which will be used to select the solver in cobrapy functions." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'cglpk'" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.solver_name" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.optimize(solver=\"cglpk\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### \\_SUPPORTS_MILP\n", "\n", "The presence of this attribute tells cobrapy that the solver supports mixed-integer linear programming" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver._SUPPORTS_MILP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### solve\n", "\n", "Model.optimize is a wrapper for each solver's solve function. It takes in a cobra model and returns a solution" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.solve(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### create_problem\n", "\n", "This creates the LP object for the solver." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lp = solver.create_problem(model, objective_sense=\"maximize\")\n", "lp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### solve_problem\n", "\n", "Solve the LP object and return the solution status" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'optimal'" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.solve_problem(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### format_solution\n", "\n", "Extract a cobra.Solution object from a solved LP object" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.format_solution(lp, model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### get_objective_value\n", "\n", "Extract the objective value from a solved LP object" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.8739215069684909" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.get_objective_value(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### get_status\n", "\n", "Get the solution status of a solved LP object" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'optimal'" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.get_status(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### change_variable_objective\n", "\n", "change the objective coefficient a reaction at a particular index. This does not change any of the other objectives which have already been set. This example will double and then revert the biomass coefficient." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "12" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.reactions.index(\"Biomass_Ecoli_core\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.7478430139369818" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.change_variable_objective(lp, 12, 2)\n", "solver.solve_problem(lp)\n", "solver.get_objective_value(lp)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.8739215069684909" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.change_variable_objective(lp, 12, 1)\n", "solver.solve_problem(lp)\n", "solver.get_objective_value(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### change variable_bounds\n", "\n", "change the lower and upper bounds of a reaction at a particular index. This example will set the lower bound of the biomass to an infeasible value, then revert it." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'infeasible'" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.change_variable_bounds(lp, 12, 1000, 1000)\n", "solver.solve_problem(lp)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'optimal'" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.change_variable_bounds(lp, 12, 0, 1000)\n", "solver.solve_problem(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### change_coefficient\n", "\n", "Change a coefficient in the stoichiometric matrix. In this example, we will set the entry for ADP in the ATMP reaction to in infeasible value, then reset it." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "16" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.metabolites.index(\"atp_c\")" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "10" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.reactions.index(\"ATPM\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'infeasible'" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.change_coefficient(lp, 16, 10, -10)\n", "solver.solve_problem(lp)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'optimal'" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.change_coefficient(lp, 16, 10, -1)\n", "solver.solve_problem(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### set_parameter\n", "\n", "Set a solver parameter. Each solver will have its own particular set of unique paramters. However, some have unified names. For example, all solvers should accept \"tolerance_feasibility.\"" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "solver.set_parameter(lp, \"tolerance_feasibility\", 1e-9)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.set_parameter(lp, \"objective_sense\", \"minimize\")\n", "solver.solve_problem(lp)\n", "solver.get_objective_value(lp)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.8739215069684912" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solver.set_parameter(lp, \"objective_sense\", \"maximize\")\n", "solver.solve_problem(lp)\n", "solver.get_objective_value(lp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example with FVA\n", "\n", "Consider flux variability analysis (FVA), which requires maximizing and minimizing every reaction with the original biomass value fixed at its optimal value. If we used the cobra Model API in a naive implementation, we would do the following:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 171 ms, sys: 0 ns, total: 171 ms\n", "Wall time: 171 ms\n" ] } ], "source": [ "%%time\n", "# work on a copy of the model so the original is not changed\n", "m = model.copy()\n", "\n", "# set the lower bound on the objective to be the optimal value\n", "f = m.optimize().f\n", "for objective_reaction, coefficient in m.objective.items():\n", " objective_reaction.lower_bound = coefficient * f\n", "\n", "# now maximize and minimze every reaction to find its bounds\n", "fva_result = {}\n", "for r in m.reactions:\n", " m.change_objective(r)\n", " fva_result[r.id] = {\n", " \"maximum\": m.optimize(objective_sense=\"maximize\").f,\n", " \"minimum\": m.optimize(objective_sense=\"minimize\").f\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead, we could use the solver API to do this more efficiently. This is roughly how cobrapy implementes FVA. It keeps uses the same LP object and repeatedly maximizes and minimizes it. This allows the solver to preserve the basis, and is much faster. The speed increase is even more noticeable the larger the model gets." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 8.28 ms, sys: 25 µs, total: 8.31 ms\n", "Wall time: 8.14 ms\n" ] } ], "source": [ "%%time\n", "# create the LP object\n", "lp = solver.create_problem(model)\n", "\n", "# set the lower bound on the objective to be the optimal value\n", "solver.solve_problem(lp)\n", "f = solver.get_objective_value(lp)\n", "for objective_reaction, coefficient in model.objective.items():\n", " objective_index = model.reactions.index(objective_reaction)\n", " # old objective is no longer the objective\n", " solver.change_variable_objective(lp, objective_index, 0.)\n", " solver.change_variable_bounds(\n", " lp, objective_index, f * coefficient,\n", " objective_reaction.upper_bound)\n", "\n", "# now maximize and minimze every reaction to find its bounds\n", "fva_result = {}\n", "for index, r in enumerate(model.reactions):\n", " solver.change_variable_objective(lp, index, 1.)\n", " result = {}\n", " solver.solve_problem(lp, objective_sense=\"maximize\")\n", " result[\"maximum\"] = solver.get_objective_value(lp)\n", " solver.solve_problem(lp, objective_sense=\"minimize\")\n", " result[\"minimum\"] = solver.get_objective_value(lp)\n", " solver.change_variable_objective(lp, index, 0.)\n", " fva_result[r.id] = result" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.3" } }, "nbformat": 4, "nbformat_minor": 0 }