1

So i was trying this and I find it really unfeasible. I am not that aware about smart ways to do the following. Can somebody help ? Also inputs of lists are quite big.

This task was to create an image from the values generated by me. center_star contains list of [x,y] pairs which are centers of various point like objects.

1800 value represents that image to be generated is of 1800x1800 pixel.

Sigma variable has value 2 by default.

final=[[0]*1800]*1800
for i in range(len(center_stars)):
    xi=center_stars[i][0]
    yi=center_stars[i][1]
    print(i)
    for j in range(1800):
        for k in range(1800):
            final[j][k]+=gauss_amplitude[i]*(math.e**((-1*((xi-j)**2+(yi-k)**2))/2*sigma*sigma))

Is there a smarter way to save time using some of the numpy operation and execute this piece of code in less time?

Murtuza Vadharia
  • 126
  • 10
  • As a minor fix, save expressions that are repeatedly reused. `2*sigma*sigma` while small is recalculated every time through the loop. You should save this expression outside the loop over i. Similarly, in the loop over k, you recompute `(xi-j)**2` every iteration even though j isn't changing. – Tyberius Apr 29 '20 at 18:52
  • 1
    This looks like something to vectorize with `numpy`. – Prune Apr 29 '20 at 18:53
  • Some things you can check/read about: https://stackoverflow.com/questions/11442191/parallelizing-a-numpy-vector-operation, https://stackoverflow.com/questions/13068760/parallelise-python-loop-with-numpy-arrays-and-shared-memory, https://stackoverflow.com/questions/51350640/parallel-for-loop-over-numpy-matrix – tevemadar Apr 29 '20 at 18:53
  • All the 1800 elements of list `final` are the same lists!, changing `final[j][k]` in loop will change **ALL** `final[:][k]`s. Probably not what you intend. Definitely use numpy. – Ehsan Apr 29 '20 at 19:35
  • Btw, are you aware of that `a=[[0]*2]*2;a[0][0]=1;a[1][0]` displays `1`? The `*` creates multiple references to the same *object*, which is not a problem for an immutable number in 1D, but it is for a list for the 2nd dimension. Use list comprehensions instead, like `b=[[0 for i in range(2)] for j in range(2)];b[0][0]=1;b[1][0]`, which displays `0`. – tevemadar May 07 '20 at 13:47

3 Answers3

0

Something like that:

import math
import numpy as np
N = 1800
final = np.empty((N, N))
final1 = np.empty((N, N))
j = np.arange(N)
k = np.arange(N)
jj, kk = np.meshgrid(j, k)
sigma = 2.
s = 0.5 / (sigma * sigma)
for i in range(len(center_stars)):
    xi = center_stars[i][0]
    yi = center_stars[i][1]
    final += gauss_amplitude[i] * np.exp(- ((xi - jj.T)**2 + (yi - kk.T)**2) * s)
# Code below is for comparison:
    for j in range(N):
        for k in range(N):
            final1[j][k]+=gauss_amplitude[i] * (math.e** (-((xi-j)**2+(yi-k)**2)/(2*sigma*sigma)))

Besides, I assume you missed brackets around 2*sigma*sigma

Alexey Mints
  • 408
  • 3
  • 12
  • I tried this, and tried testing on center_star list with around 1400 elements it took around 250 seconds while for 2300 elements it took around 410 seconds. Is it possible to still optmize it ? My data-list is quite large. – Murtuza Vadharia Apr 30 '20 at 01:01
  • You should obviously remove the last 4 lines from my example - they are just for comparison. – Alexey Mints Apr 30 '20 at 06:15
  • Also, if you need yet more speed, take a look at scipy.ndimage.convolve – Alexey Mints Apr 30 '20 at 06:19
  • ```final += gauss_amplitude[i] * np.exp(- ((xi - jj.T)**2 + (yi - kk.T)**2) * s)``` I changed this into ```final[x_min:x_max,y_min:y_max] += gauss_amplitude[i] * np.exp(- ((xs - jj.T[x_min:x_max,y_min:y_max]) ** 2 + (ys - kk.T[x_min:x_max,y_min:y_max]) ** 2) * s)```. This way I dont't have to do operation on such big meshgrid. I am calculating the x_max, x_min, y_max and y_min using width given by the user. – Murtuza Vadharia Apr 30 '20 at 12:24
0

you can try to compact your code like this:

Gauss=lambda i,j,k,xi,yi:gauss_amplitude[i]*(math.e**((-((xi-j)**2+(yi-k)**2))/(2*sigma*sigma)))
final=[[Gauss(i,j,k,x[0],x[1]) for j in range(1800) for k in range(1800)] for i,x in enumerate(center_starts)]
Sho_Lee
  • 149
  • 4
  • This only decreases the line of code. It still loops 1800*1800 times even in this list comprehension form. Its way slower than numpy. – Murtuza Vadharia Apr 30 '20 at 00:59
0

If your sigmas are all the same, you can achieve this with out any loop by using scipy.signal.convolve2d.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve2d
from scipy.stats import multivariate_normal

sigma = 3
width = 200   # smaller than yours so you can see the single pixels
n_stars = 50

# draw some random stars
star_row = np.random.randint(0, width, n_stars)
star_col = np.random.randint(0, width, n_stars)
star_amplitude = np.random.normal(50, 10, n_stars)

# assign amplitudes to center pixel of stars
amplitudes = np.zeros((width, width))
amplitudes[star_row, star_col] = star_amplitude


# create 2d gaussian kernel
row = col = np.arange(-4 * sigma, 4 * sigma + 1)
grid = np.stack(np.meshgrid(row, col)).T
kernel = multivariate_normal(
    [0, 0],
    [[sigma**2, 0], [0, sigma**2]]
).pdf(grid)
kernel /= kernel.sum()


# convolve with 2d gaussian
final = convolve2d(amplitudes, kernel, mode='same')


fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(9, 3))

img = ax1.imshow(amplitudes)
fig.colorbar(img, ax=ax1)
ax1.set_label('Before Convolution')

img = ax2.imshow(kernel)
fig.colorbar(img, ax=ax2)
ax2.set_label('Convolution Kernel')

img = ax3.imshow(final)
fig.colorbar(img, ax=ax3)
ax3.set_label('After Convolution')

fig.tight_layout()
fig.savefig('conv2d.png', dpi=300)

Result:

result

If the sigmas differ, you get a way with a single loop over the possible sigmas.

MaxNoe
  • 14,470
  • 3
  • 41
  • 46