Piano variations on Dark Eyes (with sheet music)

A few weeks ago I made a video of swing variations on the famous russian/gipsy tune Dark Eyes.

As some people have been asking for the sheet music, here it is !
Les yeux noirs / Dark Eyes

Like my previous one it is under a Creative Common licence, which means you can do pretty much everything of it but you must cite A. de Lamarmotte as the original author. Don’t hesistate to tell me if you see any obvious mistakes !

And here is the source-code of the file. Enjoy !

% OPEN WITH LILYPOND

\version "2.12.3"

%{
*
*  L E S   Y E U X  N O I R S
*
* Version 0 - 31/03/2012
*
*
* The full text of the licence can be found here:
*
* In a nutshell, YOU ARE FREE to share (copy, distribute, transmit),
* to remix or adapt the work, and to make commercial use of the work.
*
* BUT I ask you to cite "A. de la Marmotte" as the original author
* if you do any of the things above, just so that I can see what this
* little piece is becoming.
*
* Ideas for improvements are most welcome !
*
%}

title = " Les yeux noirs - Dark eyes - Очи чёрные"
subtitle = "Variations by A. de Lamarmotte"
subsubtitle = "On the famous Russian Gipsy song"
arranger="Grenoble, France, August 2012"
composer = "This is a free work ( licence CC BY 3.0)"
}

% In what follows, one line represents one bar.

global = {
\key d \minor
\time 4/4
}

right = \new Voice \with {
}{
\global
\partial 4. {gis'8 a'4}
<g' bes'>4. <g' a'>4. r4
r2. \acciaccatura{gis'16}<e' a'>4
<f' bes'>4. <d' a'>4. r4
r4 r8 d'8 f' a' d''4
d''2 \grace{e'16[ f' fis']} g'8 d''4 cis''8~cis''4.
gis'8 a'8 f''8 cis''4
e''8-^ f''16 e'' d''8 a' gis' g' f' d'
f' g' gis' a' f''4 g''4
f''\prall d''16 f'' g'' a'' bes''4 d''8 d'''4
\ottava #1
<f''' bes'''>4.-^ <g''' c''''>4-^ <f''' bes'''>
<e''' gis'''>8 <f''' a'''>8 a'' <cis''' e'''> <d''' f'''> \ottava #0  f''

<gis'' cis'''> <a'' d'''> d'' <e'' gis''> <f'' a''> a'

<cis'' e''> <d'' f''> f' gis'
\times 2/3 {a' bis' a'} g' f' e' d' cis' bes
\change Staff = "left"
\stemUp
a  g f e d cis bes, a,
<<
{ r4 <f a b>4. <f a b>4 <f a b>8-^~<f a b>4} \\
{d1 r4}
>>
\change Staff = "right"
\stemNeutral
a'8 gis' a' c'' b'4
bes'4 a'8 gis' a' c'' b' bes'4
a'8 gis' a' f'' e'' cis'' a'
\times 2/3 {g' a' g'} e' g'  \times 2/3 {f' g' f'} d' f'
\times 2/3 {e' f' e'} cis' e'  d' a c' b
bes d' f' a' d' f' fis' g'
bes' d'' e'' g' bes' d'' e'' g''
<<
{ a''1 r2} \\
{r4 c''8 f'' \times 2/3 {e'' f'' e''} a' e''
\times 2/3 {d'' e'' d''} a' }
>>
f''4-^ g''4-^
\times 2/3 { f''8-^ g'' f''} bes' bes'' \times 2/3 { e''-^ f'' e''} bes' bes''
\times 2/3 { f''-^ g'' f''} bes' <e'' bes''>4 <e'' bes'' e'''>4.
\times 2/3 { e''8-^ f'' e''} a' a'' \times 2/3 { d''-^ e'' d''} a' a''
\times 2/3 { e''-^ f'' e''} a' <d'' a''>4 <d'' a''>4 a'8
\times 2/3 { d'' e'' d''} cis'' c'' b' bes' a' g'
e' g' bes' cis'' e'' g'' bes'' cis'''
d'''4 <d'' f'' a'' b''> r4
\clef bass
<d f a b >
\clef treble
r4 <d' d''>4 <e' e''> <f' f''>
<g'' fis'>8 g' bes' d'' e'' d'' bes' g'
\times 2/3 { f' g' f'} e' dis' e' g' bes' gis'
a' cis'' f'' cis'' e'' d'' a' f'
\times 2/3 { e' f' e'} \times 2/3 { d' f' d'}
\times 2/3 { cis' f' cis'} \times 2/3 { b f' b}
bes d' f' bes' d'' bes' g' fis'
f' e' <g' c''>2 <g' bes'>4
<f' a'>4-.  <e'' gis''>8(<f'' a''>4) a'4.
\grace{bes'16[ b'16]} c''4 <gis'' b''>8(<c''' a''>4) a'4.
bes'8 d'' f'' fis'' g'' f'' d'' bes'
\times 2/3 { a' bes' a'} gis' a' cis'' e'' g'' bes''
a''4 <d'' f'' d'''>8 <d'' f'' d'''>4 <d'' f'' cis'''>8 <d'' f'' cis'''>4
<d'' f'' c'''>8 <d'' f'' c'''>4 <d'' f'' b''>4 <d'' f'' b''>8 <d'' f'' bes''>4
\times 2/3 { a''8 bes'' a''} g'' f'' e'' d'' cis'' bes'
a' g' f' e' d' cis' bes a
\time 3/4
\acciaccatura{<d' f a>8}<d' f a>4-^
\acciaccatura{<d' e gis>8}<d' e gis>4-^
\acciaccatura{<d' e g>8}<d' e g>4-^
\time 4/4
<<
{\arpeggioArrowDown <d' e f a>4\arpeggio r2.  } \\
{r4 \times 12/20 {e16[ f a ais b e' f' a' ais' b' e'' f'' a'' ais'' b''
\ottava #1
e''' f''' a''' ais''' b'''] }}>>
<d''' f''' a''' d''''>1

\bar "|."

}

left = \new Voice \with {
}{
\global
\partial 4. {r4.}
<e g bes d'>4 <e g bes d'>4 <e g bes d'>4 <e g bes d'>4
<e g a cis'>4 <e g a cis'>4 <e g a cis'>4 <e g a cis'>4
< d f a b>4 < d f a b>4 < d f a b>4 < d f a b>4
< d f a b>4 < d f a b>4 < d f a b>4 < d f a b>4
<e g bes d'>4 <e g bes d'>4 <e g bes d'>4 <e g bes d'>4
<e g a cis'>4 <e g a cis'>4 <e g a cis'>4 <e g a cis'>4
< d f a b>4 r2.
r1
<g bes d' f'>2 r2
r2
\clef treble
<gis'' cis'''>2
<a'' d'''>4 r2.
r1
<a cis' e'>4 r2.
\clef bass
\stemDown
<a, e,>4 r2.
d,2 a,,4. d,4.
\stemNeutral
r4 f2
e4 <g bes d'> d <g bes d'>
cis <e g a> a, <e g a>
d <f a b> a, <f a b>
d <f a b> a, <f a b>
c <g bes c'> d <g bes c'>
e <g bes c'> c <g bes c'>
f <a c' d'> e <a c' d'>
d <a c' d'> c <a c' d'>
<g, g> <bes d' f'> <bes, bes> <bes d' e'>
<c c'> <bes d' f'> <cis cis'> <cis' e' g'>
<d d'> <d' f' a'> <cis cis'> <d' f'>
<c c'> <d' f'> <b, b> <d' f'>
<bes, bes> <e g a> a, <e g a>
b, <e g a> cis <e g a>
d4.\sustainOn
\clef treble
<d' f' a' b'>4 <d' f' a' b'>4.
\clef bass
\ottava #-1
<d,, d,>2  \sustainOff
\ottava #0
r2
<g, g>4 <bes d' e'> d, <bes d' e'>
a, <g bes d'> cis <e g a>
d <f a b> e <f a b>
f <d a b> a, <f a b>
g, <f bes d'> d <f bes d'>
c <g bes c'> e <g bes c'>
<f a c' d'>4 r2 <f a c' d'>4
<fis a c' d'>4 r2 <d fis a c'>4
g4 <bes d' e'> d4 <g bes d'>
cis <e  a g> a, <e  a g>
<d, d>-^ <f a b> <d, d>-^ <dis, dis>-^
<e, e>-^ <f, f>-^ <fis, fis>-^ <d, d>-^
<g g,>4 <bes d' e'> d,4 <g bes d>
a, <e  a g> cis <e  a g>

\time 3/4
\acciaccatura{<d f>8}<d f>4-^
\acciaccatura{<cis e>8}<cis e>4-^
\acciaccatura{<c ees>8}<c ees>4-^
\time 4/4
<d, a, d>1\arpeggio
\clef treble
<d'' a''>2
\clef bass
<d d,>2
\bar "|."

}
#(set-global-staff-size 18)

\score {
\new PianoStaff \with {
instrumentName = "Piano"
} <<
\new Staff = "right" \with {
midiInstrument = "acoustic grand"
} \right
\new Staff = "left" \with {
midiInstrument = "acoustic grand"
} { \clef bass \left }
>>
\layout { }
\midi {
\context {
\Score
tempoWholesPerMinute = #(ly:make-moment 100 4)
}
}
}



Some like it hot

Here is an example of what I believe to be good communications skills in data representation. The Scientific Visualization Studio of the NASA Goddard Space Flight Center made this video showing 130 years of global warming:

As many of its viewers, this sure frightened me, so I went to look for the source. The map shows differences in temperature with the 1950’s. According to the NASA, the world average temperature has increased of 0.51 celcius degrees between 1951 and 2011. This surprised me a little, I would have thought it was more after watching the video ! Then it struck me that the projection of the world map (known as Mercator) gives an exagerated importance to the poles, which partly explains the amount of red in the picture. You can see how grotesquely big the antartic appears. For the artic (which isn’t reported on the map as it is not a land), it is more complicated to see, but think that Greenland is actually four times smaller than Brasil. If you have a look at the world maps you can find on the net you will see that some other solutions exist, and that this projection was a deliberate choice of the vizualization studio.

Now, how big is the distorsion ? I wrote a small script to analyse the 2011 map and found out that the average temperature increase as it appears on the picture is of +0.81°C, which is about 60% higher than the 0.51°C announced on the NASA’s website !

My methods

In case you are interested, I did the analysis with Python. I first took a screenshot of the 2011 frame in high definition, then I extracted a one-pixel high rectangle of the color bar. Finally I painted in black the zones I didn’t want to be considered (country borders, year, colorbar):

Then I used the following Python script (which requires Scipy and Pylab):

from scipy import ndimage
from scipy.interpolate import NearestNDInterpolator
from pylab import *

# load the color bar and the map

# create a function that will attribute a value to a color by looking
# at the nearest color in the colorbar

scaleValues = linspace(-2,2,len(colorScale)) # scale of the colorbar
colorInterpolator = NearestNDInterpolator(colorScale, scaleValues)
color2value = lambda color :  colorInterpolator([color])[0]

# compute the corresponding value and the norm of the colors of all
# points of the map

values = apply_along_axis(color2value, 2,map2011)
norms = apply_along_axis(norm, 2,map2011)

# compute the mean of the values on the non-black pixels

print values[norms != 0].mean()


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()



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[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


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

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[g1][g2]['sign'] = sign

import pickle

# Save the graph as EcoliNetwork.gn

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

#### 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()



Interception of a linear trajectory with constant speed

Today I have been asked about a very common interception problem, which was originally about pirates, but let us wrap it into a different story :

Alice just spotted a white rabbit urging to its rabbit hole ! Given the coordinates of the positions $A$ , $B$, $H$ of Alice, the bunny and the hole, as well as the respective speeds $S_a$ and $S_b$ of Alice and the rabbit, say whether Alice can catch the Rabbit before it disappears, and give the time and place of the fastest possible interception !

I guess that I am not the first one to solve this problem. However I couldn’t find any simple solution on the internet. So here is one, proving how helpful trigonometry can be when it comes to catching rabbits. The nice thing about it is that, while it relies on trigonometry, it doesn’t actually require to compute any trigonometrical function !

First it is obvious that, since we are looking for the fastest catch, Alice’s trajectory is a straight line. Let us call $C$ and $t_C$ the location and the time of the catch. We have then $BC = S_b t_C$ and $AC = S_a t_C$. In our problem, finding the length $BC$ would tell us whether Alice can catch the rabbit before is reaches the rabbit hole (case $BC < BH$), and would immediately lead to both the location and time of the interception :

$C = B + \dfrac{BC}{BH}\overrightarrow{BH}$

$t_C = BC/S_b$

To express BC using the coordinates of the points, let us apply the famous Law of Sines to the triangle $ABC$ :

$\dfrac{\sin \alpha}{BC} = \dfrac{\sin \beta}{AC} = \dfrac{\sin \gamma}{AB}$

$BC = \dfrac {\sin \alpha}{\sin \gamma} AB = \dfrac {\sin \alpha}{\sin \gamma} \sqrt{(x_B-x_A)^2+(y_B-y_A)^2}$

Now all we have to do is to express $\sin \alpha$ and $\sin \gamma$ in function of the given data. To do so we will first compute $\sin \beta$, then we will express $\sin \alpha$ with $\sin \beta$, and finally $\sin \gamma$ using the two first sines. The value of $\sin \beta$ can be computed from the points coordinates as follows:

$\sin \beta = \dfrac{det(\overrightarrow{BA},\overrightarrow{BH})}{ BA * BH } = \dfrac{(x_A - x_B)(y_H-y_B) - (y_A - y_B)(x_H-x_B)}{\sqrt{(x_B-x_A)^2+(y_B-y_A)^2} \sqrt{(x_B-x_H)^2+(y_B-y_H)^2}}$

Then we use the Law of Sines again, to compute $\sin \alpha$:

$\sin \alpha = \dfrac{BC}{AC} \sin \beta = \dfrac{S_b t_C}{S_a t_C} \sin \beta = \dfrac{S_b}{S_a} \sin \beta$

This only makes sense, of course, if $\frac{S_a}{S_b} | \sin \beta | \leq 1$. If this is not the case we can deduce that Alice will never catch the rabbit, which solves the problem.
Finally we use the fact that the angles of a triangle sum to $\pi$ to compute $\sin \gamma$:

$\sin \gamma = \sin (\pi - \alpha - \beta) = \sin (\alpha + \beta) = \sin \alpha \cos \beta + \cos \alpha \sin \beta$

$\sin \gamma = \sin \alpha \sqrt{1 - \sin^2 \beta} + \sin \beta \sqrt{1 - \sin^2 \alpha}$

Here we are ! Below is a script implementing this technique using Python's pylab module (but any other langage would do !).


from pylab import *

def interception(A, B, H, Sa, Sb):
''' Returns None if A cannot catch B before B reaches H.
Otherwise it returns the time and position of the interception.
See attached documentation for justifications. '''

sin_beta = det(array((A-B,H-B))) / ( norm(A-B) * norm(H-B) )

sin_alpha = (Sb / Sa) * sin_beta

if abs(sin_alpha) > 1 :

print "B moves too fast to be ever caught !"
return None

else:

sin_gamma = ( sin_alpha * sqrt(1 - sin_beta**2)
+ sin_beta * sqrt(1 - sin_alpha**2) )

BC = norm(B-A) * (sin_alpha / sin_gamma)

if BC > norm(H-A):

print "B reaches H before being caught by A !"
return None

else:

t = BC / Sb
C = B + BC * (H-B)/norm(H-B)

print "A intercepted B !"
return t,C

##### EXAMPLE

# Define the constants

A = array(( 1.0 , 5.0 ))
B = array(( 4.0 , 1.0 ))
H =  array(( 6.0 , 7.0 ))
Sa = 1.1
Sb = 1.0

# Find the intersection

t,C = interception(A, B, H, Sa, Sb)

# Plot the results

scatter(*zip(A,B,H,C), s=100, color='r')

for label, point in zip(['A','B','H','C'], [A,B,H,C]):
annotate( label, xy = point, xytext = (-10, 10),
textcoords = 'offset points', fontsize = 24)

annotate("", xy=H, xytext=B, xycoords='data', textcoords='data',size=20,
arrowprops=dict(arrowstyle="simple",connectionstyle="arc3"))

annotate("", xy=C, xytext=A, xycoords='data', textcoords='data',size=20,
arrowprops=dict(arrowstyle="simple",connectionstyle="arc3"))

title("A intercepts B in C", fontsize = 24)

show()


And here is the result :

Symbolic matrix differentiation with Sympy !

After a few days spent computing strange Jacobians for my Ph.D. thesis, I figured out that my computer could actually do most of the computations for me : all I needed was an automatic matrix differentiator, i.e. an algorithm that would tell me how a function or matrices $F(X,Y,...)$ varies when the matrices $(X,Y...)$ are changed into $(X + \partial X, Y + \partial Y,...)$ where the $(\partial X,\partial Y,...)$ are small. The very useful Matrix Cookbook gives for instance

$\partial (X^T) = (\partial X)^T$

$\partial (XY) = (\partial X)Y + X (\partial Y)$

$\partial(X^{-1}) = - X^{-1}(\partial X) X^{-1}$

… And so on. We can see that it is not completely unlike classic differentiation, but the non-commutativity slightly complicates things, and the computed formulae can be monstruous : imagine differentiating by hand the hat matrix $H = X \left(X^\top \Sigma^{-1} X\right)^{-1} X^\top \Sigma^{-1}$, where $\Sigma$ is a constant 😯

Maybe I didn’t look properly, but I didn’t find any program to automatize the rules described above, so I went on building something around Sympy, a young yet sophisticated symbolic mathematics module for the Python programmation language, whose authors claim that it is easily extendable. I put that claim to the test and indeed I found it pretty easy to work my way through.

In what follows I explain what I did, in tutorial-like fashion, with the hope that this can help understanding how one can tweak Sympy to meet one’s own needs. But I am no Sympy expert, and I take any tips and comments in the comments section below ! If you are in a hurry, you’ll find the complete source code of the solution at the bottom of this blog.

First, a matrix is a symbol, which is lucky because sympy already has a Symbol class. You only have to specify that these symbols are not commutative:


def matrices(names):
''' Call with  A,B,C = matrix('A B C') '''
return symbols(names,commutative=False)



Matrix-specific transformations

All the regular operations like addition, substraction, multiplication… are already implemented in the Symbol class, so we can focus on matrix-specific transformations, like the operator $\partial$ and the inverse function $X \mapsto X^{-1}$ :

d = Function("d",commutative=False)
inv = Function("inv",commutative=False)


That was tough ! Now for the transpose function $X \mapsto X^T$ we want the general rules

$(X+Y...)^T = X^T + Y^T ...$
$(XY...)^T = ...Y^TX^T$

to automatically apply  when transposing an expression. Here is the code I wrote for that. To understand it you need  to know that when I call some Sympy Function F over a bunch of arguments X,Y,Z, the result is an object of class “F” with a list of arguments (.args)  equal to [X,Y,Z]. Moreover successive additions or multiplications are flattened into one single object: for instance
$X^2 + 2*X*Y + 3*Z + 5$
will be an object of class “Add” with arguments
$[X^2,2*X*Y,3*Z,5]$,
while
$X^2*Y*X*Z$
will be of class “Mul” with arguments
$[X^2,Y,X,Z]$.

class t(Function):
''' The transposition, with special rules
t(A+B) = t(A) + t(B) and t(AB) = t(B)t(A) '''
is_commutative = False
def __new__(cls,arg):
return Add(*[t(A) for A in arg.args])
elif arg.is_Mul:
L = len(arg.args)
return Mul(*[t(arg.args[L-i-1]) for i in range(L)])
else:
return Function.__new__(cls,arg)


That’s enough transformations for now, let us explain to sympy how to differentiate a matrix expression !

Matrix differentiation

We write all the just enough differentiation rules of the Matrix Cookbook into one dictionnary:

MATRIX_DIFF_RULES = {
# e =expression, s = a list of symbols respsect to which
# we want to differentiate

Symbol : lambda e,s : d(e) if (e in s) else 0,
Mul :  lambda e,s : Mul(matDiff(e.args[0],s),Mul(*e.args[1:]))
+  Mul(e.args[0],matDiff(Mul(*e.args[1:]),s)) ,
t :   lambda e,s : t( matDiff(e.args[0],s) ),
inv : lambda e,s : - e * matDiff(e.args[0],s) * e
}


and apply them recursively to the expression we want to treat :

def matDiff(expr,symbols):
if expr.__class__ in MATRIX_DIFF_RULES.keys():
return  MATRIX_DIFF_RULES[expr.__class__](expr,symbols)
else:
return 0


Simple as pie, and we are done ! Let’s play around with our new toy:

X,S = matrices("X S")
H= X*inv(t(X)*inv(S)*X)*t(X)*inv(S)
print  mdiff(H,X)

>>> X*(inv(t(X)*inv(S)*X)*t(d(X))*inv(S) - inv(t(X)*inv(S)*X)*(t(X)*inv(S)*d(X) + t(d(X))*inv(S)*X)*inv(t(X)*inv(S)*X)*t(X)*inv(S)) + d(X)*inv(t(X)*inv(S)*X)*t(X)*inv(S)


It works ! But we must concede that it is barely readable… so let us put some style in these expressions !

Cosmetics

You don’t need latex to have nice-looking formulae. Here is what we want to say to sympy:

• Don’t write inv(X), write X¯¹
• Don’t write t(X), write X’
• Don’t write d(X), write ∂X
• In a general way, don’t put parenthesis if the transformation applies to one symbol only.

The programmers of Sympy have made it all easy to customize the default printing method. You just write a printing method with complements for the functions of our own. First we need to include this line at the top of our source file to explain that we are going to use strange characters:

# -*- coding: utf-8 -*-


Then we write our own printer class that herits from the default one but has a few more methods :

class matStrPrinter(StrPrinter):
''' Nice printing for console mode : X¯¹, X', ∂X '''

def _print_inv(self, expr):
if expr.args[0].is_Symbol:
return  self._print(expr.args[0]) +'¯¹'
else:
return '(' +  self._print(expr.args[0]) + ')¯¹'

def _print_t(self, expr):
return  self._print(expr.args[0]) +"'"

def _print_d(self, expr):
if expr.args[0].is_Symbol:
return '∂'+  self._print(expr.args[0])
else:
return '∂('+  self._print(expr.args[0]) +')'


A this point the people who wrote the Sympy documentation want us to directly replace the default printer with this command:

Basic.__str__ = lambda self: matStrPrinter().doprint(self)


It indeed seems to be the only way to ensure that all the expressions in the recursion will print well with the standard “print” method, but it is a little to definitive for me (what if I have several customized printers to use in the same script ?). Moreover I don’t like the ‘*’ between the matrices and I didn’t find any smart way to put them away. For all these reasons I wrote my own matrix printing command:

def matPrint(m):
mem = Basic.__str__
Basic.__str__ = lambda self: matStrPrinter().doprint(self)
print str(m).replace('*','')
Basic.__str__ = mem


Now let’s try again:

matPrint(  matDiff(H,X) )

>>> X((X'S¯¹X)¯¹∂X'S¯¹ - (X'S¯¹X)¯¹(X'S¯¹∂X + ∂X'S¯¹X)(X'S¯¹X)¯¹X'S¯¹) + ∂X(X'S¯¹X)¯¹X'S¯¹


Way better ! Notice, however that the first X is a factor or two terms of a sum, one of these being in turn a factor involving a sum. We can break these imbricated expressions into a sum of products (which is simpler to analyze) by calling twice Sympy’s ‘expand’ function :

matPrint(  expand ( expand ( matDiff(H,X) ) ) )

>>> X(X'S¯¹X)¯¹∂X'S¯¹ + ∂X(X'S¯¹X)¯¹X'S¯¹ - X(X'S¯¹X)¯¹X'S¯¹∂X(X'S¯¹X)¯¹X'S¯¹ - X(X'S¯¹X)¯¹∂X'S¯¹X(X'S¯¹X)¯¹X'S¯¹


That’s nice enough for the console ! Now what if I want to report this fundamental result to the scientific community ? Let’s see how to generate the LaTeX code to embed the formula in a document.

Automatic LaTeX code generation

One of the nice features of Sympy is its ability to generate LaTeX code (and even display compiled formulae). For matrix computations we would like Sympy to follow these rules:

• $X^{-1}$ writes X^{-1}
• $X^{T}$ writes X^{T}
• $\partial X$ writes \partial X
• In a general way, don’t put parenthesis if the transformation applies to one symbol only.

Like in the previous section we create a new ‘printer’ class with these features, and then write a function inpired by Sympy’s latex method (which is simply called latex() ):

class matLatPrinter(LatexPrinter):
''' Printing instructions for latex : X^{-1},  X^T, \partial X '''

def _print_inv(self, expr):
if expr.args[0].is_Symbol:
return self._print(expr.args[0]) +'^{-1}'
else:
return '(' + self._print(expr.args[0]) + ')^{-1}'
def _print_t(self, expr):
return  self._print(expr.args[0]) +'^T'

def _print_d(self, expr):
if expr.args[0].is_Symbol:
return '\partial '+ self._print(expr.args[0])
else:
return '\partial ('+ self._print(expr.args[0]) +')'

def matLatex(expr, profile=None, **kargs):
if profile is not None:
profile.update(kargs)
else:
profile = kargs
return matLatPrinter(profile).doprint(expr)



We now try it on the derivative of H:

print matLatex( matDiff(H,X) )

>>> $X \left((X^T S^{-1} X)^{-1} \partial X^T S^{-1} - (X^T S^{-1} X)^{-1} \left(X^T S^{-1} \partial X + \partial X^T S^{-1} X\right) (X^T S^{-1} X)^{-1} X^T S^{-1}\right) + \partial X (X^T S^{-1} X)^{-1} X^T S^{-1}$


which once compiled yields
$X \left((X^T S^{-1} X)^{-1} \partial X^T S^{-1} - (X^T S^{-1} X)^{-1} \left(X^T S^{-1} \partial X + \partial X^T S^{-1} X\right) (X^T S^{-1} X)^{-1} X^T S^{-1}\right) + \partial X (X^T S^{-1} X)^{-1} X^T S^{-1}$

Yeah !

Source Code

I hope this helped. If you just came here to get a matrix differentiator, here is the code with the example. Have fun 🙂


# Declaration

# -*- coding: utf-8 -*-

#----------------------------------------------------------------------
#
# FUNCTIONS FOR THE AUTOMATIC DIFFERENTIATION  OF MATRICES WITH SYMPY
#
#----------------------------------------------------------------------

from sympy import *
from sympy.printing.str import StrPrinter
from sympy.printing.latex import LatexPrinter

#####  M  E  T  H  O  D  S

def matrices(names):
''' Call with  A,B,C = matrix('A B C') '''
return symbols(names,commutative=False)

# Transformations

d = Function("d",commutative=False)
inv = Function("inv",commutative=False)

class t(Function):
''' The transposition, with special rules
t(A+B) = t(A) + t(B) and t(AB) = t(B)t(A) '''
is_commutative = False
def __new__(cls,arg):
return Add(*[t(A) for A in arg.args])
elif arg.is_Mul:
L = len(arg.args)
return Mul(*[t(arg.args[L-i-1]) for i in range(L)])
else:
return Function.__new__(cls,arg)

# Differentiation

MATRIX_DIFF_RULES = {

# e =expression, s = a list of symbols respsect to which
# we want to differentiate

Symbol : lambda e,s : d(e) if (e in s) else 0,
Mul :  lambda e,s : Mul(matDiff(e.args[0],s),Mul(*e.args[1:]))
+  Mul(e.args[0],matDiff(Mul(*e.args[1:]),s)) ,
t :   lambda e,s : t( matDiff(e.args[0],s) ),
inv : lambda e,s : - e * matDiff(e.args[0],s) * e
}

def matDiff(expr,symbols):
if expr.__class__ in MATRIX_DIFF_RULES.keys():
return  MATRIX_DIFF_RULES[expr.__class__](expr,symbols)
else:
return 0

#####  C  O  S  M  E  T  I  C  S

# Console mode

class matStrPrinter(StrPrinter):
''' Nice printing for console mode : X¯¹, X', ∂X '''

def _print_inv(self, expr):
if expr.args[0].is_Symbol:
return  self._print(expr.args[0]) +'¯¹'
else:
return '(' +  self._print(expr.args[0]) + ')¯¹'

def _print_t(self, expr):
return  self._print(expr.args[0]) +"'"

def _print_d(self, expr):
if expr.args[0].is_Symbol:
return '∂'+  self._print(expr.args[0])
else:
return '∂('+  self._print(expr.args[0]) +')'

def matPrint(m):
mem = Basic.__str__
Basic.__str__ = lambda self: matStrPrinter().doprint(self)
print str(m).replace('*','')
Basic.__str__ = mem

# Latex mode

class matLatPrinter(LatexPrinter):
''' Printing instructions for latex : X^{-1},  X^T, \partial X '''

def _print_inv(self, expr):
if expr.args[0].is_Symbol:
return self._print(expr.args[0]) +'^{-1}'
else:
return '(' + self._print(expr.args[0]) + ')^{-1}'
def _print_t(self, expr):
return  self._print(expr.args[0]) +'^T'

def _print_d(self, expr):
if expr.args[0].is_Symbol:
return '\partial '+ self._print(expr.args[0])
else:
return '\partial ('+ self._print(expr.args[0]) +')'

def matLatex(expr, profile=None, **kargs):
if profile is not None:
profile.update(kargs)
else:
profile = kargs
return matLatPrinter(profile).doprint(expr)

#####    T  E  S  T  S

X,S = matrices("X S")
H= X*inv(t(X)*inv(S)*X)*t(X)*inv(S)

matPrint(  expand( expand( matDiff(H,X) ) ) )

print matLatex( matDiff(H,X) )