myprojects-hg-reborn

Fork

#! /usr/bin/env python
from scipy import *
import pylab  # used to read the .csv file

# this is a simple test code I use to understand how newaxis works in SciPy
# and to understand how Anne Archibald was able to re-express concisely the creation
# of k-mers in Smoluchowsky equation


kern=arange(16)
kern.shape=(4,4)
print 'kern is', kern

y=array([1,2,4,3]) 
#y=arange(4)

print 'y is',y

kyn=kern*y[:,newaxis]*y[newaxis,:]
# the position of the newaxis determines that one vector y is seen as y[i] and the other
# as y[j]
# See this also by calculating
#kyn=kern*y[newaxis,:] 
# and 
#kyn=kern*y[:,newaxis]

print 'kyn is', kyn

z=y[:,newaxis]+y[newaxis,:]
print 'z is',z


k=0
#This  is now the whole point: if you look at the code to solve smoluchowsky equation
# first you see that in the loop there is the condition:
# i+1+j+1=k+1 ---> j=k-i-1
#The matrix red implements exactly this: the first element is i and it varies in
# arange(k), whereas the 2nd one is j and it is equal to k-i-1.
# note that while k is fixed each time, arange(k) is an array!
# in the code by Anne, there is also a loop on k and finally red is summed to get the creator
# (apart from a trivial multiplication times a factor 0.5)
red=kyn[arange(k),k-arange(k)-1]

print 'red is',red

print 'sum(red) is',sum(red)  #this works OK, since fixing k=0 sets the creator =0 (I can
#not create the fundamental monomer

print 'So far so good'