Merge branch 'devel' into FVA_infeasibility

Conflicts:
	cameo/solver_based_model.py

Vagrantfile

@@ -0,0 +1,65 @@
+# -*- mode: ruby -*-
+# vi: set ft=ruby :
+
+$script = <<SCRIPT
+
+sudo apt-get update
+sudo apt-get install -y make git libglpk-dev swig glpk-utils swig python-dev python-numpy python-scipy python-matplotlib
+sudo apt-get install -y python-pip
+sudo pip install python-libsbml-experimental cython -f https://dl.dropboxusercontent.com/u/22461024/pypi.drosophi.la/index.html --no-index
+sudo pip install ipython[notebook]
+sudo pip install nose
+cd /vagrant
+sudo pip install -r requirements.txt
+sudo pip install escher --pre
+sudo python setup.py develop
+# get rid of the stupid openpyxl error message
+pip uninstall openpyxl
+pip install openpyxl==1.8.6
+pip install nose rednose coverage
+sudo pip install redis
+sudo python setup.py develop
+
+SCRIPT
+
+# Vagrantfile API/syntax version. Don't touch unless you know what you're doing!
+VAGRANTFILE_API_VERSION = "2"
+
+Vagrant.configure(VAGRANTFILE_API_VERSION) do |config|
+
+  config.vm.box = "ubuntu/trusty64"
+
+  config.vm.provider "virtualbox" do |v|
+  host = RbConfig::CONFIG['host_os']
+
+  # Give VM 1/4 system memory & access to all cpu cores on the host
+  if host =~ /darwin/
+    cpus = `sysctl -n hw.ncpu`.to_i
+    # sysctl returns Bytes and we need to convert to MB
+    mem = `sysctl -n hw.memsize`.to_i / 1024 / 1024 / 4
+  elsif host =~ /linux/
+    cpus = `nproc`.to_i
+    # meminfo shows KB and we need to convert to MB
+    mem = `grep 'MemTotal' /proc/meminfo | sed -e 's/MemTotal://' -e 's/ kB//'`.to_i / 1024 / 4
+  else # sorry Windows folks, I can't help you
+    cpus = 2
+    mem = 1024
+  end
+
+  v.customize ["modifyvm", :id, "--memory", mem]
+  v.customize ["modifyvm", :id, "--cpus", cpus]
+end
+
+  # Disable automatic box update checking. If you disable this, then
+  # boxes will only be checked for updates when the user runs
+  # `vagrant box outdated`. This is not recommended.
+  # config.vm.box_check_update = false
+
+  # Create a forwarded port mapping which allows access to a specific port
+  # within the machine from a port on the host machine. In the example below,
+  # accessing "localhost:8080" will access port 80 on the guest machine.
+  config.vm.network "forwarded_port", guest: 8888, host: 8888
+
+  config.vm.provision :shell, :inline => $script
+
+end

cameo/solver_based_model.py

@@ -90,7 +90,109 @@ def to_solver_based_model(cobrapy_model, solver_interface=optlang):
     return solver_based_model
 
 
-class LazySolution(object):
+class SolutionBase(object):
+
+    def __init__(self, model, *args, **kwargs):
+        super(SolutionBase, self).__init__(*args, **kwargs)
+        self.model = model
+
+    def __str__(self):
+        """A pandas DataFrame representation of the solution.
+
+        Returns
+        -------
+        pandas.DataFrame
+        """
+        return str(self.to_frame())
+
+    def as_cobrapy_solution(self):
+        """Convert into a cobrapy Solution.
+
+        Returns
+        -------
+        cobra.core.Solution.Solution
+        """
+        return Solution(self.f, x=self.x,
+                        x_dict=self.x_dict, y=self.y, y_dict=self.y_dict,
+                        the_solver=None, the_time=0, status=self.status)
+
+    def to_frame(self):
+        """Return the solution as a pandas DataFrame.
+
+        Returns
+        -------
+        pandas.DataFrame
+        """
+        return DataFrame({'fluxes': Series(self.x_dict), 'reduced_costs': Series(self.y_dict)})
+
+
+    def get_primal_by_id(self, reaction_id):
+        """Return a flux/primal value for a reaction.
+
+        Parameters
+        ----------
+        reaction_id : str
+            A reaction ID.
+        """
+        return self.x_dict[reaction_id]
+
+
+class Solution(SolutionBase):
+    """This class mimicks the cobrapy Solution class.
+
+    Attributes
+    ----------
+    primal_dict : dict
+        A dictionary of flux values.
+    dual_dict : dict
+        A dictionary of reduced costs.
+
+    Notes
+    -----
+    See also documentation for cobra.core.Solution.Solution for an extensive list of inherited attributes.
+    """
+
+    def __init__(self, model, *args, **kwargs):
+        """
+        Parameters
+        ----------
+        model : SolverBasedModel
+        """
+        super(Solution, self).__init__(model, *args, **kwargs)
+        self.f = model.solver.objective.value
+        self.x = []
+        for reaction in model.reactions:
+            primal = reaction.variable.primal
+            if reaction.reversibility:
+                primal -= reaction.reverse_variable.primal
+            self.x.append(primal)
+        self.y = []
+        for reaction in model.reactions:
+            dual = reaction.variable.dual
+            if reaction.reversibility:
+                dual -= reaction.reverse_variable.dual
+            self.y.append(dual)
+        self.status = model.solver.status
+        self._reaction_ids = [r.id for r in self.model.reactions]
+
+    @property
+    def x_dict(self):
+        return dict(zip(self._reaction_ids, self.x))
+
+    @property
+    def y_dict(self):
+        return dict(zip(self._reaction_ids, self.y))
+
+    @property
+    def primal_dict(self):
+        return self.x_dict
+
+    @property
+    def dual_dict(self):
+        return self.y_dict
+
+
+class LazySolution(SolutionBase):
     """This class implements a lazy evaluating version of the cobrapy Solution class.
 
     Attributes
@@ -106,16 +208,36 @@ class LazySolution(object):
     See also documentation for cobra.core.Solution.Solution for an extensive list of inherited attributes.
     """
 
-    def __init__(self, model):
+    def __init__(self, model, *args, **kwargs):
         """
         Parameters
         ----------
         model : SolverBasedModel
         """
-        self.model = model
-        self._time_stamp = self.model._timestamp_last_optimization
+        super(LazySolution, self).__init__(model, *args, **kwargs)
+        if self.model._timestamp_last_optimization is not None:
+            self._time_stamp = self.model._timestamp_last_optimization
+        else:
+            self._time_stamp = time.time()
         self._f = None
 
+    def _check_freshness(self):
+        """Raises an exceptions if the solution might have become invalid due to re-optimization of the attached model.
+
+        Raises
+        ------
+        UndefinedSolution
+            If solution has become invalid.
+        """
+        if self._time_stamp != self.model._timestamp_last_optimization:
+            timestamp_formatter = lambda timestamp: datetime.datetime.fromtimestamp(timestamp).strftime(
+                "%Y-%m-%d %H:%M:%S:%f")
+            raise UndefinedSolution(
+                'The solution (captured around %s) has become invalid as the model has been re-optimized recently (%s).' % (
+                    timestamp_formatter(self._time_stamp),
+                    timestamp_formatter(self.model._timestamp_last_optimization))
+            )
+
     @property
     def f(self):
         self._check_freshness()
@@ -173,63 +295,6 @@ def primal_dict(self):
     def dual_dict(self):
         return self.y_dict
 
-    def __str__(self):
-        """A pandas DataFrame representation of the solution.
-
-        Returns
-        -------
-        pandas.DataFrame
-        """
-        return str(self.to_frame())
-
-    def _check_freshness(self):
-        """Raises an exceptions if the solution might have become invalid due to re-optimization of the attached model.
-
-        Raises
-        ------
-        UndefinedSolution
-            If solution has become invalid.
-        """
-        if self._time_stamp != self.model._timestamp_last_optimization:
-            timestamp_formatter = lambda timestamp: datetime.datetime.fromtimestamp(timestamp).strftime(
-                "%Y-%m-%d %H:%M:%S:%f")
-            raise UndefinedSolution(
-                'The solution (captured around %s) has become invalid as the model has been re-optimized recently (%s).' % (
-                    timestamp_formatter(self._time_stamp),
-                    timestamp_formatter(self.model._timestamp_last_optimization))
-            )
-
-    def as_cobrapy_solution(self):
-        """Convert into a cobrapy Solution.
-
-        Returns
-        -------
-        cobra.core.Solution.Solution
-        """
-        return Solution(self.f, x=self.x,
-                        x_dict=self.x_dict, y=self.y, y_dict=self.y_dict,
-                        the_solver=None, the_time=0, status=self.status)
-
-    def to_frame(self):
-        """Return the solution as a pandas DataFrame.
-
-        Returns
-        -------
-        pandas.DataFrame
-        """
-        return DataFrame({'fluxes': Series(self.x_dict), 'reduced_costs': Series(self.y_dict)})
-
-
-    def get_primal_by_id(self, reaction_id):
-        """Return a flux/primal value for a reaction.
-
-        Parameters
-        ----------
-        reaction_id : str
-            A reaction ID.
-        """
-        return self.x_dict[reaction_id]
-
 
 class Reaction(OriginalReaction):
     """This class extends the cobrapy Reaction class to work with SolverBasedModel.
@@ -260,6 +325,17 @@ def clone(cls, reaction, model=None):
             new_reaction._model = None
         if model is not None:
             new_reaction._model = model
+        for gene in new_reaction.genes:
+            # print gene._reaction
+            # for r in list(gene._reaction):
+            #     print r, type(r)
+            gene._reaction.remove(reaction)
+            # print gene._reaction
+            gene._reaction.add(new_reaction)
+            # print gene._reaction
+        for metabolite in new_reaction.metabolites:
+            metabolite._reaction.remove(reaction)
+            metabolite._reaction.add(new_reaction)
         return new_reaction
 
     def __init__(self, name=None):
@@ -275,7 +351,7 @@ def __init__(self, name=None):
         self._objective_coefficient = 0.
 
     def __str__(self):
-        return self.build_reaction_string()
+        return ''.join((self.id, ": ", self.build_reaction_string()))
 
     @property
     def variable(self):
@@ -652,7 +728,52 @@ def add_demand(self, metabolite, prefix="DM_"):
         self.add_reactions([demand_reaction])
         return demand_reaction
 
-    def optimize(self, new_objective=None, objective_sense='maximize', solution_type=LazySolution, **kwargs):
+    def add_ratio_constraint(self, reaction1, reaction2, ratio, prefix='ratio_constraint'):
+        """Adds a ratio constraint (reaction1/reaction2 = ratio) to the model.
+
+        Parameters
+        ----------
+        reaction1 : str or Reaction
+            A reaction or a reaction ID.
+        reaction2 : str or Reaction
+            A reaction or a reaction ID.
+        ratio : float
+            The ratio in reaction1/reaction2 = ratio
+        prefix : str
+            The prefix that will be added to the constraint ID.
+
+        Returns
+        -------
+        optlang.Constraint
+            The constraint name will be composed of `prefix`
+            and the two reaction IDs (e.g. 'ratio_constraint_reaction1_reaction2').
+
+        Examples
+        --------
+        model.add_ratio_constraint('r1', 'r2', 0.5)
+        print model.solver.constraints['ratio_constraint_r1_r2']
+        > ratio_constraint: ratio_constraint_r1_r2: 0 <= -0.5*r1 + 1.0*PGI <= 0
+        """
+        if isinstance(reaction1, types.StringType):
+            reaction1 = self.reactions.get_by_id(reaction1)
+        if isinstance(reaction2, types.StringType):
+            reaction2 = self.reactions.get_by_id(reaction2)
+
+        if reaction1.reverse_variable is not None:
+            term1 = reaction1.variable - reaction1.reverse_variable
+        else:
+            term1 = reaction1.variable
+
+        if reaction2.reverse_variable is not None:
+            term2 = reaction2.variable - reaction2.reverse_variable
+        else:
+            term2 = reaction2.variable
+
+        ratio_constraint = self.solver.interface.Constraint(term1 - ratio * term2, lb=0, ub=0, name='ratio_constraint_'+reaction1.id+'_'+reaction2.id)
+        self.solver._add_constraint(ratio_constraint, sloppy=True)
+        return ratio_constraint
+
+    def optimize(self, new_objective=None, objective_sense=None, solution_type=Solution, **kwargs):
         """OptlangBasedModel implementation of optimize. Returns lazy solution object. Exists for compatibility reasons. Uses model.solve() instead."""
         if new_objective is None or new_objective == 0:
             pass
@@ -694,18 +815,24 @@ def optimize(self, new_objective=None, objective_sense='maximize', solution_type
             "%Y-%m-%d %H:%M:%S:%f")
         self._timestamp_last_optimization = time.time()
         # logger.debug('self._timestamp_last_optimization ' + timestamp_formatter(self._timestamp_last_optimization))
-        self.solver.optimize()
+        if objective_sense is not None:
+            original_direction = self.objective.direction
+            self.objective.direction = {'minimize': 'min', 'maximize': 'max'}[objective_sense]
+            self.solver.optimize()
+            self.objective.direction = original_direction
+        else:
+            self.solver.optimize()
         solution = solution_type(self)
         # logger.debug('solution = solution_type(self) ' + timestamp_formatter(solution._time_stamp))
         self.solution = solution
         return solution
 
     def solve(self, *args, **kwargs):
         """Optimize model."""
-        solution = self.optimize(*args, **kwargs)
+        solution = self.optimize(solution_type=LazySolution, *args, **kwargs)
         if solution.status is not 'optimal':
             self.solver.configuration.presolve = True
-            solution = self.optimize(*args, **kwargs)
+            solution = self.optimize(solution_type=LazySolution, *args, **kwargs)
             self.solver.configuration.presolve = False
             if solution.status is not 'optimal':
                 status = solution.status