added 3d plot option

Author: Philipp Schneider

conda-recipe/meta.yaml

@@ -1,4 +1,4 @@
-{% set version = '0.38' %}
+{% set version = '0.33' %}
 
 package:
   name: straindesign

doc/plot.png

@@ -8,13 +8,15 @@
 from straindesign.names import *
 from typing import Dict, Tuple
 from pandas import DataFrame
-from numpy import floor, sign, mod, nan, isnan, unique, inf, isinf, zeros, full, linspace, prod
+from numpy import floor, sign, mod, nan, isnan, unique, inf, isinf, full, linspace, \
+                  prod, array, mean, flip, ceil, floor
+from numpy.linalg import matrix_rank
 from os import cpu_count
 from contextlib import redirect_stdout, redirect_stderr
 from io import StringIO
 import matplotlib.pyplot as plt
-from scipy.spatial import ConvexHull   
-import mpl_toolkits.mplot3d as a3
+from scipy.spatial import ConvexHull, Delaunay
+from matplotlib.cm import ScalarMappable, get_cmap
 
 from straindesign.parse_constr import linexpr2mat, linexprdict2str
 
@@ -496,13 +498,13 @@ def yopt(model,**kwargs):
     else:
         status = INFEASIBLE
 
-def phase_plane(model, axes, **kwargs):
+def plot_flux_space(model, axes, **kwargs):
     reaction_ids = model.reactions.list_attr("id")
 
     if CONSTRAINTS in kwargs: 
         kwargs[CONSTRAINTS] = parse_constraints(kwargs[CONSTRAINTS],reaction_ids)
     else:
-        kwargs[CONSTRAINTS] = None
+        kwargs[CONSTRAINTS] = []
 
     if SOLVER not in kwargs:
         kwargs[SOLVER] = None
@@ -513,7 +515,7 @@ def phase_plane(model, axes, **kwargs):
     else:
         show = True    
 
-    axes = list(axes) # cast to list of lists
+    axes = [list(ax) if not isinstance(ax,str) else [ax] for ax in axes] # cast to list of lists
     num_axes = len(axes)
     if num_axes not in [2,3]:
         raise Exception('Please define 2 or 3 axes as a list of tuples [ax1, ax2, (optional) ax3] with ax1 = (den,num).\n'+\
@@ -529,28 +531,54 @@ def phase_plane(model, axes, **kwargs):
 
     ax_name = ["" for _ in range(num_axes)]
     ax_limits = [(nan,nan) for _ in range(num_axes)]
+    ax_type = ["" for _ in range(num_axes)]
     for i,ax in enumerate(axes):
-        ax = parse_linexpr(ax,reaction_ids)[0]
-        ax_name[i] = linexprdict2str(ax)
-        sol_min = fba(model,obj=ax,constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='minimize')
-        sol_max = fba(model,obj=ax,constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='maximize')
-        # abort if any of the fluxes are unbounded or undefined
-        unbnd = [i+1 for i,v in enumerate([sol_min,sol_max]) if v.status == UNBOUNDED]
-        if any(unbnd):
-            raise Exception('One of the specified reactions is unbounded or undefined. Phase plane cannot be generated.')
-        ax_limits[i] = [min((0,sol_min.objective_value)),max((0,sol_max.objective_value))]
-        axes[i] = ax
+        if len(ax) == 1:
+            ax_type[i] = 'rate'
+        elif len(ax) == 2:
+            ax_type[i] = 'yield'
+        if ax_type[i] == 'rate':
+            ax[0] = parse_linexpr(ax[0],reaction_ids)[0]
+            ax_name[i] = linexprdict2str(ax[0])
+            sol_min = fba(model,obj=ax[0],constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='minimize')
+            sol_max = fba(model,obj=ax[0],constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='maximize')
+            # abort if any of the fluxes are unbounded or undefined
+            inval = [i+1 for i,v in enumerate([sol_min,sol_max]) if v.status == UNBOUNDED or v.status == INFEASIBLE]
+            if any(inval):
+                raise Exception('One of the specified reactions is unbounded or problem is infeasible. Plot cannot be generated.')
+            ax_limits[i] = [min((0,sol_min.objective_value)),max((0,sol_max.objective_value))]
+        elif ax_type[i] == 'yield':
+            ax[0] = parse_linexpr(ax[0],reaction_ids)[0]
+            ax[1] = parse_linexpr(ax[1],reaction_ids)[0]
+            ax_name[i] = '('+linexprdict2str(ax[0])+') / ('+linexprdict2str(ax[1])+')'
+            sol_min = yopt(model,obj_num=ax[0],obj_den=ax[1],constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='minimize')
+            sol_max = yopt(model,obj_num=ax[0],obj_den=ax[1],constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='maximize')
+            # abort if any of the yields are unbounded or undefined
+            inval = [i+1 for i,v in enumerate([sol_min,sol_max]) if v.status == UNBOUNDED or v.status == INFEASIBLE]
+            if any(inval):
+                raise Exception('One of the specified yields is unbounded or undefined or problem is infeasible. Plot cannot be generated.')
+            ax_limits[i] = [min((0,ceil_dec(sol_min.objective_value,9))),max((0,floor_dec(sol_max.objective_value,9)))]
 
     # compute points
     x_space = linspace(ax_limits[0][0], ax_limits[0][1], num=points)
     lb = full(points, nan)
     ub = full(points, nan)
     for i,x in enumerate(x_space):
-        constr = [axes[0],'=',x]
-        sol_vmin = fba(model,constraints=constr,obj=axes[1],obj_sense='minimize')
-        lb[i] = sol_vmin.objective_value
-        sol_vmax = fba(model,constraints=constr,obj=axes[1],obj_sense='maximize')
-        ub[i] = sol_vmax.objective_value
+        constr = kwargs[CONSTRAINTS].copy()
+        if ax_type[0] == 'rate':
+            constr += [[axes[0][0],'=',x]]
+        elif ax_type[0] == 'yield':
+            constr += [[{**axes[0][0], **{k:-v*x for k,v in axes[0][1].items()}},'=',0]]
+        if ax_type[1] == 'rate':
+            sol_vmin = fba(model,constraints=constr,obj=axes[1][0],obj_sense='minimize')
+            lb[i] = ceil_dec(sol_vmin.objective_value,9)
+            sol_vmax = fba(model,constraints=constr,obj=axes[1][0],obj_sense='maximize')
+            ub[i] = floor_dec(sol_vmax.objective_value,9)
+        elif ax_type[1] == 'yield':
+            sol_vmin = yopt(model,constraints=constr,obj_num=axes[1][0],obj_den=axes[1][1],obj_sense='minimize')
+            lb[i] = ceil_dec(sol_vmin.objective_value,9)
+            sol_vmax = yopt(model,constraints=constr,obj_num=axes[1][0],obj_den=axes[1][1],obj_sense='maximize')
+            ub[i] = floor_dec(sol_vmax.objective_value,9)
 
     if num_axes == 2:
         x = [v for v in x_space] + [v for v in reversed(x_space)]
@@ -559,6 +587,7 @@ def phase_plane(model, axes, **kwargs):
             x.extend([x_space[0], x_space[0]])
             y.extend([lb[0],      ub[0]])
         plot1 = plt.fill(x, y)
+        # plot1 = plt.plot(x, y)
         plot1[0].axes.set_xlabel(ax_name[0])
         plot1[0].axes.set_ylabel(ax_name[1])
         plot1[0].axes.set_xlim(ax_limits[0][0]*1.05,ax_limits[0][1]*1.05)
@@ -570,22 +599,87 @@ def phase_plane(model, axes, **kwargs):
     elif num_axes == 3:
         max_diff_y = max([abs(l-u) for l,u in zip(lb,ub)])
         datapoints = []
-        for x,l,u in zip(x_space,lb,ub):
+        datapoints_top = []
+        datapoints_bottom = []
+        for i,(x,l,u) in enumerate(zip(x_space,lb,ub)):
             if l-u != 0:
                 y_space = linspace(l,u,int(-(-points // (max_diff_y/abs(l-u)))))
             else:
                 y_space = [0.0]
-            for y in y_space:
-                constr = [[axes[0],'=',x], [axes[1],'=',y]]
-                sol_vmin = fba(model,constraints=constr,obj=axes[2],obj_sense='minimize')
-                datapoints += [[x,y,sol_vmin.objective_value]]
-                sol_vmax = fba(model,constraints=constr,obj=axes[2],obj_sense='maximize')
-                datapoints += [[x,y,sol_vmax.objective_value]]
-        hull = ConvexHull(datapoints)
+            datapoints_top += [[]]
+            datapoints_bottom += [[]]
+            for j,y in enumerate(y_space):
+                constr = kwargs[CONSTRAINTS].copy()
+                if ax_type[0] == 'rate':
+                    constr += [[axes[0][0],'=',x]]
+                elif ax_type[0] == 'yield':
+                    constr += [[{**axes[0][0], **{k:-v*x for k,v in axes[0][1].items()}},'=',0]]
+                if ax_type[1] == 'rate':
+                    constr += [[axes[1][0],'=',y]]
+                elif ax_type[1] == 'yield':
+                    constr += [[{**axes[1][0], **{k:-v*y for k,v in axes[1][1].items()}},'=',0]]
+                if ax_type[2] == 'rate':
+                    sol_vmin = fba(model,constraints=constr,obj=axes[2][0],obj_sense='minimize')
+                    sol_vmax = fba(model,constraints=constr,obj=axes[2][0],obj_sense='maximize')
+                elif ax_type[2] == 'yield':
+                    sol_vmin = yopt(model,constraints=constr,obj_num=axes[2][0],obj_den=axes[2][1],obj_sense='minimize')
+                    sol_vmax = yopt(model,constraints=constr,obj_num=axes[2][0],obj_den=axes[2][1],obj_sense='maximize')
+                datapoints_top[i] += [len(datapoints)]
+                datapoints += [[x,y,floor_dec(sol_vmax.objective_value,9)]]
+                datapoints_bottom[i] += [len(datapoints)]
+                datapoints += [[x,y,ceil_dec(sol_vmin.objective_value,9)]]
+        
+        # Construct Denaunay triangles for plotting from all 6 perspectives
+        triang = []
+        # triangles top
+        for i in range(len(datapoints_top)-1): 
+            temp_points = datapoints_top[i]+datapoints_top[i+1]
+            pts = [[datapoints[idx_p][0],datapoints[idx_p][1]] for idx_p in temp_points]
+            if matrix_rank(pts) > 1:
+                triang_temp = Delaunay(pts).simplices
+                triang += [[temp_points[idx] for idx in p] for p in triang_temp]
+        # triangles bottom
+        for i in range(len(datapoints_bottom)-1): 
+            temp_points = datapoints_bottom[i]+datapoints_bottom[i+1]
+            pts = [[datapoints[idx_p][0],datapoints[idx_p][1]] for idx_p in temp_points]
+            if matrix_rank(pts) > 1:
+                triang_temp = Delaunay(pts).simplices
+                triang += [[temp_points[idx] for idx in flip(p)] for p in triang_temp]
+        # triangles front
+        for i in range(len(x_space)-1): 
+            temp_points = [datapoints_top[i][0],datapoints_top[i+1][0],datapoints_bottom[i][0],datapoints_bottom[i+1][0]]
+            pts = [[datapoints[idx_p][0],datapoints[idx_p][2]] for idx_p in temp_points]
+            if matrix_rank(pts) > 1:
+                triang_temp = Delaunay(pts).simplices
+                triang += [[temp_points[idx] for idx in flip(p)] for p in triang_temp]
+        # triangles back
+        for i in range(len(x_space)-1): 
+            temp_points = [datapoints_top[i][-1],datapoints_top[i+1][-1],datapoints_bottom[i][-1],datapoints_bottom[i+1][-1]]
+            pts = [[datapoints[idx_p][0],datapoints[idx_p][2]] for idx_p in temp_points]
+            if matrix_rank(pts) > 1:
+                triang_temp = Delaunay(pts).simplices
+                triang += [[temp_points[idx] for idx in flip(p)] for p in triang_temp]
+        # triangles left
+        for i in range(len(datapoints_top[0])-1): 
+            temp_points = [datapoints_top[0][i],datapoints_top[0][i+1],datapoints_bottom[0][i],datapoints_bottom[0][i+1]]
+            pts = [[datapoints[idx_p][1],datapoints[idx_p][2]] for idx_p in temp_points]
+            if matrix_rank(pts) > 1:
+                triang_temp = Delaunay(pts).simplices
+                triang += [[temp_points[idx] for idx in p] for p in triang_temp]
+        # triangles right
+        for i in range(len(datapoints_top[-1])-1): 
+            temp_points = [datapoints_top[-1][i],datapoints_top[-1][i+1],datapoints_bottom[-1][i],datapoints_bottom[-1][i+1]]
+            pts = [[datapoints[idx_p][1],datapoints[idx_p][2]] for idx_p in temp_points]
+            if matrix_rank(pts) > 1:
+                triang_temp = Delaunay(pts).simplices
+                triang += [[temp_points[idx] for idx in p] for p in triang_temp]
+
+        # hull = ConvexHull(datapoints)
         x = [d[0] for d in datapoints]
         y = [d[1] for d in datapoints]
         z = [d[2] for d in datapoints]
-        ax = a3.Axes3D(plt.figure())
+        # ax = a3.Axes3D(plt.figure())
+        ax = plt.figure().add_subplot(projection='3d')
         ax.dist=10
         ax.azim=30
         ax.elev=10
@@ -595,112 +689,27 @@ def phase_plane(model, axes, **kwargs):
         ax.set_xlabel(ax_name[0])
         ax.set_ylabel(ax_name[1])
         ax.set_zlabel(ax_name[2])
-        plot1 = ax.plot_trisurf(x,y,z,triangles=hull.simplices,linewidth=0.2, antialiased=True, alpha=0.67)
+        # set higher color value for 'larger' values in all dimensions
+        # use middle value of triangles instead of means and then norm 
+        colors = [nan for _ in triang]
+        for i,t in enumerate(triang):
+            p = [datapoints[i] for i in t]
+            interior_point = [0.0, 0.0, 0.0]
+            for d in range(3):
+                v = [q[d] for q in p]
+                interior_point[d] = mean([max(v),min(v)])/max(abs(array(ax_limits[d])))
+            colors[i] = sum(interior_point) # norm(interior_point)
+        lw = min([1,6.0/len(triang)])
+        plot1 = ax.plot_trisurf(x,y,z,triangles=triang,linewidth=lw,edgecolors='black', antialiased=True, alpha=0.90) #  array=colors, cmap=plt.cm.winter
+        colors = colors/max(colors)
+        colors = get_cmap("Spectral")(colors)
+        plot1.set_fc(colors)
         if show:
             plt.show()
         return plot1
-    
-def yield_space(model, axes, **kwargs):
-    reaction_ids = model.reactions.list_attr("id")
-    
-    if CONSTRAINTS in kwargs: 
-        kwargs[CONSTRAINTS] = parse_constraints(kwargs[CONSTRAINTS],reaction_ids)
-    else:
-        kwargs[CONSTRAINTS] = None
-
-    if SOLVER not in kwargs:
-        kwargs[SOLVER] = None
-    solver = select_solver(kwargs[SOLVER],model)
-        
-    if 'show' in kwargs:
-            show = kwargs['show']
-    else:
-        show = True    
-    
-    axes = [list(ax) for ax in axes] # cast to list of lists
-    num_axes = len(axes)
-    if num_axes not in [2,3]:
-        raise Exception('Please define 2 or 3 axes as a list of tuples [ax1, ax2, (optional) ax3] with ax1 = (den,num).\n'+\
-                        '"den" and "num" being linear expressions.')
-            
-    if 'points' in kwargs:
-        points = kwargs['points']
-    else:
-        if num_axes == 2:
-            points = 40
-        else:
-            points = 25
-    
-    ax_name = ["" for _ in range(num_axes)]
-    ax_limits = [(nan,nan) for _ in range(num_axes)]
-    for i,ax in enumerate(axes):
-        ax[0] = parse_linexpr(ax[0],reaction_ids)[0]
-        ax[1] = parse_linexpr(ax[1],reaction_ids)[0]
-        ax_name[i] = '('+linexprdict2str(ax[0])+') / ('+linexprdict2str(ax[1])+')'
-        sol_min = yopt(model,obj_num=ax[0],obj_den=ax[1],constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='minimize')
-        sol_max = yopt(model,obj_num=ax[0],obj_den=ax[1],constraints=kwargs[CONSTRAINTS],solver=solver,obj_sense='maximize')
-        # abort if any of the yields are unbounded or undefined
-        unbnd = [i+1 for i,v in enumerate([sol_min,sol_max]) if v.status == UNBOUNDED]
-        if any(unbnd):
-            raise Exception('One of the specified yields is unbounded or undefined. Yield space cannot be generated.')
-        ax_limits[i] = [min((0,sol_min.objective_value)),max((0,sol_max.objective_value))]
-
-    # compute points
-    x_space = linspace(ax_limits[0][0], ax_limits[0][1], num=points)
-    lb = full(points, nan)
-    ub = full(points, nan)
-    for i,x in enumerate(x_space):
-        constr = [{**axes[0][0], **{k:-v*x for k,v in axes[0][1].items()}},'=',0]
-        sol_vmin = yopt(model,constraints=constr,obj_num=axes[1][0],obj_den=axes[1][1],obj_sense='minimize')
-        lb[i] = sol_vmin.objective_value
-        sol_vmax = yopt(model,constraints=constr,obj_num=axes[1][0],obj_den=axes[1][1],obj_sense='maximize')
-        ub[i] = sol_vmax.objective_value
 
-    if num_axes == 2:
-        x = [v for v in x_space] + [v for v in reversed(x_space)]
-        y = [v for v in lb] + [v for v in reversed(ub)]
-        if lb[0] != ub[0]:
-            x.extend([x_space[0], x_space[0]])
-            y.extend([lb[0],      ub[0]])
-        plot1 = plt.fill(x, y)
-        plot1[0].axes.set_xlabel(ax_name[0])
-        plot1[0].axes.set_ylabel(ax_name[1])
-        plot1[0].axes.set_xlim(ax_limits[0][0]*1.05,ax_limits[0][1]*1.05)
-        plot1[0].axes.set_ylim(ax_limits[1][0]*1.05,ax_limits[1][1]*1.05)
-        if show:
-            plt.show()
-        return plot1
+def ceil_dec(v,n):
+    return ceil(v*(10**n))/(10**n)
 
-    elif num_axes == 3:
-        max_diff_y = max([abs(l-u) for l,u in zip(lb,ub)])
-        datapoints = []
-        for x,l,u in zip(x_space,lb,ub):
-            if l-u != 0:
-                y_space = linspace(l,u,int(-(-points // (max_diff_y/abs(l-u)))))
-            else:
-                y_space = [0.0]
-            for y in y_space:
-                constr = [[{**axes[0][0], **{k:-v*x for k,v in axes[0][1].items()}},'=',0],\
-                          [{**axes[1][0], **{k:-v*y for k,v in axes[1][1].items()}},'=',0]]
-                sol_vmin = yopt(model,constraints=constr,obj_num=axes[2][0],obj_den=axes[2][1],obj_sense='minimize')
-                datapoints += [[x,y,sol_vmin.objective_value]]
-                sol_vmax = yopt(model,constraints=constr,obj_num=axes[2][0],obj_den=axes[2][1],obj_sense='maximize')
-                datapoints += [[x,y,sol_vmax.objective_value]]
-        hull = ConvexHull(datapoints)
-        x = [d[0] for d in datapoints]
-        y = [d[1] for d in datapoints]
-        z = [d[2] for d in datapoints]
-        ax = a3.Axes3D(plt.figure())
-        ax.dist=10
-        ax.azim=30
-        ax.elev=10
-        ax.set_xlim(ax_limits[0])
-        ax.set_ylim(ax_limits[1])
-        ax.set_zlim(ax_limits[2])
-        ax.set_xlabel(ax_name[0])
-        ax.set_ylabel(ax_name[1])
-        ax.set_zlabel(ax_name[2])
-        plot1 = ax.plot_trisurf(x,y,z,triangles=hull.simplices,linewidth=0.2, antialiased=True, alpha=0.67)
-        if show:
-            plt.show()
-        return plot1
+def floor_dec(v,n):
+    return floor(v*(10**n))/(10**n)