Python for Bioinformatics: Peaks of density in numpy

Saturday, March 13, 2010

Peaks of density in numpy

Here is a Python script that generates peaks of user-defined density. It has two main parts. Given an existing matrix and a point x,y (in the interval 0-1), it constructs a matrix of the same shape containing the distance of each element from the designated point.

The second part uses the distance information to construct peaks of density to add to the matrix. The score function is vectorized. For each input element of a distance matrix, it computes a score based on (i) a threshold beyond which the score is zero, (ii) a height for the maximum score, and (iii) a function describing how the density changes with distance.

Nothing fancy, and a bit slow, but it does make it easy to play with the parameters and generate a density like that in the sample. At the end, we need to make sure each value is <= 1.0 and >= 0, because make_heatmap converts these directly to "gray" values to input to the Rectangle function.

Here is output for a 10 x 10 example:

0.04  0.09  0.04  0.04  0.04  0.04  0.04  0.04  0.04  0.04
0.09  0.81  0.09  0.04  0.04  0.05  0.16  0.22  0.16  0.05
0.04  0.09  0.04  0.04  0.04  0.16  0.41  0.56  0.41  0.16
0.04  0.04  0.04  0.04  0.04  0.22  0.56  1.0   0.56  0.22
0.04  0.04  0.04  0.04  0.04  0.16  0.41  0.56  0.41  0.16
0.04  0.04  0.04  0.04  0.04  0.05  0.16  0.22  0.16  0.05
0.04  0.07  0.09  0.07  0.04  0.04  0.04  0.04  0.04  0.04
0.07  0.17  0.28  0.17  0.07  0.04  0.04  0.04  0.04  0.04
0.09  0.28  0.81  0.28  0.09  0.04  0.04  0.04  0.04  0.04
0.07  0.17  0.28  0.17  0.07  0.04  0.04  0.04  0.04  0.04

import matplotlib.pyplot as plt
import numpy as np
import HM_Utils_frac

def init_data(N=30):
    A = np.zeros(N**2)
    A.shape = (N,N)
    return A

def eucl_distance(x0,y0,x,y):
    return np.sqrt((x-x0)**2+(y-y0)**2)

def distance_matrix(A,x,y):
    # matrix is 1.0 x 1.0 in "size"
    # given a matrix A with shape R,C
    # each row,col in M is a fractional pos
    R,C = A.shape
    L = [1.0*c/C for c in range(C)]
    M = list()
    for r in range(R):
        for c in range(C):
            M.append((L[c],L[r]))
    # now compute each dist from given point x,y
    M = [eucl_distance(x0,y0,x,y) for x0,y0 in M]
    M = np.array(M)
    M.shape = A.shape
    return M

@np.vectorize
def score(d,ht,thr,func):
    # given a distance from x,y compute score
    # first, check threshold
    if d > thr:  return 0.0
    d = thr - d
    # scale d to func(d) and height
    return ht* (func(d) / func(thr))

def make_peak(A,x,y,ht,thr,func=None):
    M = distance_matrix(A,x,y)
    if not func:  func = lambda x: x**1.5
    return score(M,ht,thr,func)
    
def scale(A):
    maxN = np.max(A)
    if maxN > 1.0:  A = A/(1.0*maxN)
    return A
    
def f(x):  return x**4

N = 40
ax,fig = HM_Utils_frac.init_fig()
A = init_data(N)
A += make_peak(A,x=0.7,y=0.3,ht=1.0,thr=0.3)
A += make_peak(A,x=0.2,y=0.8,ht=0.8,thr=0.4,func=f)
A += make_peak(A,x=0.1,y=0.1,ht=0.8,thr=0.2,func=f)
A = A + 0.04
A = scale(A)
if N <= 10:  HM_Utils_frac.show_array(A)     
HM_Utils_frac.make_heatmap(A,ax)
plt.savefig('example.png')