## 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.040.09 0.81 0.09 0.04 0.04 0.05 0.16 0.22 0.16 0.050.04 0.09 0.04 0.04 0.04 0.16 0.41 0.56 0.41 0.160.04 0.04 0.04 0.04 0.04 0.22 0.56 1.0 0.56 0.220.04 0.04 0.04 0.04 0.04 0.16 0.41 0.56 0.41 0.160.04 0.04 0.04 0.04 0.04 0.05 0.16 0.22 0.16 0.050.04 0.07 0.09 0.07 0.04 0.04 0.04 0.04 0.04 0.040.07 0.17 0.28 0.17 0.07 0.04 0.04 0.04 0.04 0.040.09 0.28 0.81 0.28 0.09 0.04 0.04 0.04 0.04 0.040.07 0.17 0.28 0.17 0.07 0.04 0.04 0.04 0.04 0.04`

 `import matplotlib.pyplot as pltimport numpy as npimport HM_Utils_fracdef init_data(N=30): A = np.zeros(N**2) A.shape = (N,N) return Adef 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.vectorizedef 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**4N = 40ax,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.04A = scale(A)if N <= 10: HM_Utils_frac.show_array(A) HM_Utils_frac.make_heatmap(A,ax)plt.savefig('example.png')`