A GUI for the exploration of functions with Python / Matplotlib

Sliders can be a great tools for lazy engineers and modelists who want to play around with the parameters of a program without recompiling it a thousand times.
Most scientific languages ( python-pylab, matlab, scilab…) support some basic sliders methods but it can be quite long to create a GUI. In this blog I present a light, user-friendly function I implemented for my everyday work.

Example with a numerical function

def volume(x,y,z):
	""" Volume of a box with width x, heigth y, and depth z """
	return x*y*z

intervals = [ { 'label' :  'width',  'valmin': 1 , 'valmax': 5 },
              { 'label' :  'height',  'valmin': 1 , 'valmax': 5 },
              { 'label' :  'depth',  'valmin': 1 , 'valmax': 5 } ] 

inputExplorer(volume,intervals)

The new value of volume is automatically computed each time you move the sliders. It is also computed if you press Enter (which can be useful to get several values, when your function is non-deterministic). In case you are studying a function that is long to evaluate, you can decide that it will only be evaluated when Enter is pressed, by adding “wait_for_enter = True” to the arguments list.

Example with a graphical function

You can also decide to use a function that doesn’t return any value, but updates a plot. Let us illustrate this with the Lotka-Volterra model, in which wolves eat rabbits and die, leading to periodical fluctuations of both populations. From wikipedia:

\frac{dx}{dt} = x(\alpha - \beta y)
\frac{dy}{dt} = - y(\gamma - \delta  x)

This can be simulated easily in python with scipy’s function odeint :

import numpy as np
from scipy.integrate import odeint

# model

def model(state,t, a,b,c,d):
	x,y = state
	return [ x*(a-b*y) , -y*(c - d*x) ] # dx/dt, dy/dt

# vector of times

t_vec = np.linspace(0,10,500) # 500 time points between 0 and 10

# initial conditions

x0 = 0.5
y0 = 1

# parameters

a = 1
b = 2
c = 1
d = 3

# simulate and return the values at each time t in t_vec

result = odeint(model, [x0,y0], ts, args = (a,b,c,d) ) 

# plot

import matplotlib.pyplot as plt
plt.plot(ts,result)
plt.show()

However, trying by hand several values of x0, y0, a, b, c, d, can be fastidious. So let us wrap the solving and the plotting in a function and feed it to inputExplorer. Let me rewrite the whole code from scratch :

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

def model(state,t, a,b,c,d):
	x,y = state
	return [ x*(a-b*y) , -y*(c - d*x) ]

ts = np.linspace(0,10,500)

fig,ax = plt.subplots(1)

def plotDynamics(x0,y0,a,b,c,d):
	ax.clear()
	ax.plot(ts, odeint(model, [x0,y0], ts, args = (a,b,c,d)) )
	fig.canvas.draw()

sliders = [ { 'label' :  label,  'valmin': 1 , 'valmax': 5 }
		 for label in [ 'x0','y0','a','b','c','d' ] ]

inputExplorer(plotDynamics,sliders)

Code

Enjoy !

import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button

def inputExplorer(f, sliders_properties, wait_for_validation = False):
    """ A light GUI to manually explore and tune the outputs of 
        a function.
        slider_properties is a list of dicts (arguments for Slider )
        whose keys are in ( label, valmin, valmax, valinit=0.5, 
        valfmt='%1.2f', closedmin=True, closedmax=True, slidermin=None,
        slidermax=None, dragging=True)
        
        def volume(x,y,z):
            return x*y*z
    
        intervals = [ { 'label' :  'width',  'valmin': 1 , 'valmax': 5 },
                  { 'label' :  'height',  'valmin': 1 , 'valmax': 5 },
                  { 'label' :  'depth',  'valmin': 1 , 'valmax': 5 } ]
        inputExplorer(volume,intervals)
    """
        
    nVars = len(sliders_properties)
    slider_width = 1.0/nVars
    print slider_width
    
    # CREATE THE CANVAS
    
    figure,ax = plt.subplots(1)
    figure.canvas.set_window_title( "Inputs for '%s'"%(f.func_name) )
    
    # choose an appropriate height
    
    width,height = figure.get_size_inches()
    height = min(0.5*nVars,8)
    figure.set_size_inches(width,height,forward = True)
    
    
    # hide the axis
    ax.set_frame_on(False)
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    

    # CREATE THE SLIDERS
    
    sliders = []
    
    for i, properties in enumerate(sliders_properties):
        
        ax = plt.axes([0.1 , 0.95-0.9*(i+1)*slider_width,
                       0.8 , 0.8* slider_width])
        sliders.append( Slider(ax=ax, **properties) )
    
    
    # CREATE THE CALLBACK FUNCTIONS
    
    def on_changed(event) : 
        
        res = f(*(s.val for s in sliders))
        
        if res is not None:
            
            print res
    
    def on_key_press(event):
        
        if event.key is 'enter':
            
            on_changed(event)   
    
    figure.canvas.mpl_connect('key_press_event', on_key_press)
    
    # AUTOMATIC UPDATE ?
    
    if not wait_for_validation:
        
        for s in sliders :
            
            s.on_changed(on_changed)
    
    
    # DISPLAY THE SLIDERS
    
    plt.show()

Advertisements

The first game based on a gene regulation network ! (and it’s boring !)

Do you know the Wikigame ? It is an online game whose goal is to go from one article of Wikipedia to another by only clicking on the hypertext links of the articles (actually quite a lot of fun 🙂 ). In my previous post I showed how to use Python to explore E. coli’s network of gene interactions. Now, why not transpose the Wikigame to this network ? Here is how I would do it :

  1. Get the network
  2. Take its main connected component.
  3. Choose two genes of this component randomly
  4. Compute the shortest path between the two components
  5. Set the maximum number of moves allowed to two times the length of the shortest time
  6. Lead the player through the network, proposing at each nodes to go to any of its neighbors.
  7. End the game if the player reached its goal or if the moves limit has been reached

And here is the code :

import random

def iGame(geneNetwork):
	""" A game in the spirit of the wikiGame : the computer decides of
		two random genes and you must find a way from one to the other
		in a maximum of ten moves. """

	nb_allowed_moves = 10

	# SELECT THE BIGGEST CONNECTED COMPONENT OF THE NETWORK

	undirected_net = geneNetwork.to_undirected()
	component_nodes = nx.connected_components(undirected_net)[0]
	component = geneNetwork.subgraph(component_nodes)

	# SELECT TWO RANDOM GENES INSIDE THIS COMPONENT

	start,goal = random.sample(component,2)

	# COMPUTE THE SHORTEST PATH

	shortest_path = nx.shortest_path(undirected_net,start,goal)
	nb_allowed_moves = 2*len(shortest_path)

	# A RECURSIVE FUNCTION DEFINING THE RULES

	def recursive(gene, countDown):
		""" A recursive function to explore E. coli's network from
			gene to gene """

		if gene == goal:

			print "You win ! Congratulations !"

		elif countDown == 0 :

			print "So sorry, no moves left ! You lose !"

		else :

			regulators = component.predecessors(gene)
			regulated =  component.successors(gene)

			print '\n -- Goal : %s . %d moves left'%(goal,countDown)
			print 'Genes regulating %s : '%(gene), " ".join(regulators)
			print 'Genes regulated by %s : '%(gene), " ".join(regulated)

			playerChoice = ''

			while not playerChoice in ( regulators + regulated ):

				playerChoice = raw_input('Where to move next ?')

			recursive(playerChoice, countDown-1)

	# START THE GAME

	print "\n Welcome to the iGAME, You are a level 1 raider"
	print "located on the gene %s. Your quest, would you choose"%(start)
	print "to accept it, is to reach the gene %s in less"%(goal)
	print "than %d moves. Good luck !"%(nb_allowed_moves)

	recursive(start,nb_allowed_moves)

	# END THE GAME

	print "The shortest path was : "
	print " ".join(shortest_path)

This game ( that I called iGame as a reference to the iGEM competition) is certainly one of the dullest, unfunniest games I ever coded. It’s main drawback is that, without any background information, it comes to a random walk with no visibility in a maze with hundreds of bag-ends. Maybe it would be more fun if it could be played directly on the RegulonDB webpages, which give detailled descriptions on the role and location of each gene. At which case it would be one of the nerdiest games on the internet !

E. coli’s regulation network from RegulonDB to Python

RegulonDB provides a continuously updated network of interactions of the bacterium E. coli. In this post I will show how easy it is to explore this network using Python with its package NetworkX. In case you are in a hurry, I put the whole script with enough comments at the end of this post.

Getting the Graph

First we create an empty graph:

import networkx as nx
eColiNetwork = nx.DiGraph(name = 'E. Coli Network')

Then we fetch the RegulonDB file containing the informations on the network. To get an internet file in Python we use urllib2:

from urllib2 import urlopen
url = "http://regulondb.ccg.unam.mx/data/network_tf_gene.txt"
regulonFile = urlopen(url)

Now we will read this file line by line to fill our graph with nodes (genes) and edges (interactions). The file contains uninformative lines (comments begining with ‘#’ or empty lines) and lines of the form ‘gene1 gene2 sign …’ where the sign is ‘+’ or ‘-‘ according to whether the interaction is an activation or an inhibition.

for line in regulonFile :

	# ignore comments and empty lines

	if not line.startswith(('\n','\t','#')):

		# A line is of the form "gene1	gene2 sign ..."

		g1,g2, sign = line.split('\t')[:3]

		# standardize gene names ( all characters to lower case )

		g1,g2 = g1.lower(),g2.lower()

		# add the edge to the network

		eColiNetwork.add_edge(g1,g2)
		eColiNetwork[g1][g2]['sign'] = sign

Easy as pie with Python ! You can then save this graph for offline use :

nx.write_gpickle(eColiNetwork,'EcoliNetwork.gn') # To save
eColiNetwork = nx.read_gpickle('EcoliNetwork.gn') # To load

Study and Statistics

It is very easy to get informations on a NetworkX graph. Let us start with some basic statistics:

nx.info(eColiNetwork)
>>> Name: E. Coli Network
    Type: DiGraph
    Number of nodes: 1696
    Number of edges: 3809
    Average in degree:   2.2459
    Average out degree:   2.2459

Note that there are only 1700 genes out of the circa 4000 genes of E. coli. Moreover it doesn’t have so much edges, so we can expect it not to be fully connected, i.e. maybe you cannot always go from one gene to another in this graph. Let us check that out.

# undirected version of E. coli's network
undirected_net = eColiNetwork.to_undirected()
 
connected_components = nx.connected_components(undirected_net)

print 'Number of connected components : ', \
		len(connected_components )

print 'Length of the connected components : ', \
		[ len(c) for c in connected_components]

 

>>> Number of connected components :  22
>>> Length of the connected components :  [1603, 13, 7, 6, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 1]

Indeed we have 21 small satellite networks and one big central component. Let’s learn more about the latter:

main_component = nx.connected_component_subgraphs(undirected_net)[0]

print "Average shortest path in E. coli's main connected component : ", \
	   nx.average_shortest_path_length(main_component)

print "Diameter of E. coli's main connected component : ", \
	   nx.diameter(main_component)

 

>>>  Average shortest path in E. coli's main connected component :  3.61111539459
>>>  Diameter of E. coli's main connected component :  9

We can also get information on a specific gene, like its regulators (which are its predecessors in the graph, i.e. the genes pointing at it) and the genes that it regulates (its successors). Let’s have a look at acs and allr :

acs_regulators = eColiNetwork.predecessors('acs')
print "Genes regulating acs : ", " ".join( acs_regulators)

allr_targets = eColiNetwork.successors('allr')
print "Genes regulated by allr : ", " ".join(allr_targets)
>>> Genes regulating acs :  fis crp ihf
>>> Genes regulated by allr :  allb alla hyi gcl glxr alls ybby ybbw glxk

Sketching a Subgraph

Because it is so easy, let us draw some networks with NetworkX. Displaying the whole graph, or even its main component, would be ugly, as they are too big, so we choose the second connected component:

component = nx.connected_component_subgraphs(undirected_net)[1]

We use pylab for its plotting capabilities :

from pylab import *

figure()
nx.draw(component, node_size = 1000)
show()

And voilà !

E. coli subnetwork

E. coli subnetwork

Complete Code

Have fun 🙂


# (F)utilities to explore the network of Escherichia coli
# as given on RegulonDB


##### FETCH THE NETWORK ON THE INTERNET ################################


import networkx as nx
from urllib2 import urlopen

# Fetch E. coli's network on Regulon DB as a .txt file

url = "http://regulondb.ccg.unam.mx/data/network_tf_gene.txt"
regulonFile = urlopen(url)

# Create an (empty) directed graph

eColiNetwork = nx.DiGraph(name = 'E. Coli Network')

# Read the file and fill our network

for line in regulonFile :

	# ignore comments and empty lines

	if not line.startswith(('\n','\t','#')):

		# A line is of the form "gene1	gene2	+ "

		g1,g2, sign = line.split('\t')[:3]

		# standardize gene names ( all characters to lower case )

		g1,g2 = g1.lower(),g2.lower()

		# add the edge to the network

		eColiNetwork.add_edge(g1,g2)
		eColiNetwork[g1][g2]['sign'] = sign


########### SAVE AND LOAD ##############################################


import pickle

# Save the graph as EcoliNetwork.gn

with open('EcoliNetwork.gn', 'wb') as output:
	pickle.dump(eColiNetwork, output, pickle.HIGHEST_PROTOCOL)

# Load the graph

eColiNetwork = pickle.load(open('EcoliNetwork.gn', 'r'))


#### PRINT INFOS ON THE NETWORK ########################################


# undirected version of E. coli's network

undirected_net = eColiNetwork.to_undirected()

# Stats

print '-- General informations : \n', \
		nx.info(eColiNetwork)

connected_components = nx.connected_components(undirected_net)

print 'Number of connected components : ', \
		len(connected_components )

print 'Length of the connected components : ', \
		[ len(c) for c in connected_components]

main_component = nx.connected_component_subgraphs(undirected_net)[0]

print "-- Average shortest path in E. coli's main connected component : ", \
	   nx.average_shortest_path_length(main_component)

print "-- Diameter of E. coli's main connected component : ", \
	   nx.diameter(main_component)

# Infos on a specific gene :

print "-- Genes regulating acs : ", \
	   " ".join( eColiNetwork.predecessors('acs') )

print "-- Genes regulated by allr : ", \
	   " ".join( eColiNetwork.successors('allr') )


###########  D R A W  ##################################################


from pylab import *

figure()
component = nx.connected_component_subgraphs(undirected_net)[1]
nx.draw(component, node_size = 1000)
show()