title: Speeding up analysis on arrays in numpy
body: <p>I have a python code, which imports 4 column txt file with numbers
  first three columns are x,y,z coordinate and fourth column is a density at that coordinate.</p>
  
  <p>below is the code that reads, converts to ndarray, Fourier transform that field, calculate the distance from origin (k=(0,0,0)) and a transformed coordinate, and taking average, and plot them.
  Thanks to pandas (python library for data analysis) and python FFT, loading 256^3 rows and Fourier transform them are very fast and done in few seconds.</p>
  
  <p>However, transforming the loaded txt to numpy ndarray, calculating the average density (average values of each coordinate), and calculating distance from the origin (k=(0,0,0)) take very long time.</p>
  
  <p>I think the problem is np.around part at the end but I can't figure out the way to optimize it.</p>
  
  <p>I have a resource of 32 core machine.</p>
  
  <p>Could someone teach me how to speed up, make it a multiprocess code, or something like that so that those can be done very quickly? Thanks.</p>
  
  <p>(If you are cosmologist and ever need this code, you may use it but please contact me, if you can. Thanks)</p>
  
  <pre><code>from __future__ import division
  import numpy as np
  
  ngridx = 128
  ngridy = 128    
  ngridz = 128
  
  maxK = max(ngridx,ngridy,ngridz)
  
  #making input file
  f = np.zeros((ngridx*ngridy*ngridz,4))
  
  i = 0
  for i in np.arange(len(f)):
      f[i][0] = int(i/(ngridy*ngridz))
      f[i][1] = int((i/ngridz))%ngridy
      f[i][2] = int(i%ngridz)
      f[i][3] = np.random.rand(1)
      if i%1000000 ==0:
          print i
  #This takes forever
  #end making input file
  
  #Thanks to Mike,
  a = f[:,3].reshape(ngridx,ngridy,ngridz)
  
  avg =np.sum(f[:,3])/len(f)
  a /= avg
  p = np.fft.fftn(a)
  #This part is much much faster than before (Original Post).
  
  #Keeping track of corresponding wavenumbers (k_x, k_y,k_z) for each element in p
  #This is just a convension on fourier transformation so you can ignore this part
  kValx = np.fft.fftfreq( ngridx , (1.0 / ngridx ) )
  kValy = np.fft.fftfreq( ngridy , (1.0 / ngridy ) )
  kValz = np.fft.fftfreq( ngridz , (1.0 / ngridz ) )
  kx = np.zeros((ngridx,ngridy,ngridz))
  ky = np.zeros((ngridx,ngridy,ngridz))
  kz = np.zeros((ngridx,ngridy,ngridz))
  rangecolx = np.arange(ngridx)
  rangecoly = np.arange(ngridy)
  rangecolz = np.arange(ngridz)
  for row in np.arange(ngridx):
      for column in np.arange(ngridy):
          for height in np.arange(ngridz):
              kx[row][column][height] = (kValx[row])
              ky[row][column][height] = (kValy[column])
              kz[row][column][height] = (kValz[height])
      if row%10 == 0:
          print row
  print 'wavenumber generate complete!'
  
  #Calculating the average powerspectrum in terms of fixed K (Distance from origin to a point in fourier space)
  #by taking the spherical shell of thickness 1 and averaging out the values inside it.
  #I am sure that this process can be optimised somehow, but I gave up.
  
  qlen = maxK/2 #Nyquist frequency
  q = np.zeros(((qlen),4),dtype=complex)
  #q is a four column array with length maxK/2.
  #q[:,0] is integer wavenumber (K, which is the distance from the origin = sqrt(kx^2+ky^2+kz^2))
  #q[:,1] is the sum of square of the fourier transformed value 
  #q[:,2] is the sum of the fourier transformed value, 
  #and q[:,3] is the total number of samples with K=q[:,0]
  
  for i in  np.arange(len(q)):
      q[i][0] = i
  i = 0
  for i in np.arange(len(p)):
      for r in np.arange(len(p[0])):
          for s in np.arange(len(p[0,0])):
              K = np.around(np.sqrt(kx[i,r,s]**2+ky[i,r,s]**2+kz[i,r,s]**2))
              if K &lt; qlen:
                  q[K][1]=q[K][1]+np.abs(p[i,r,s])**2
                  q[K][2]=q[K][2]+p[i,r,s]
                  q[K][3]=q[K][3]+1   
      if i%10 ==0:
          print 'i = ',i,' !'
  print q
  </code></pre>