The best way of writing loops in numpy is not writing loops and instead using vectorized operations. For example:
c = 0
for i in range(len(a)):
c += a[i] + b[i]
becomes
c = np.sum(a + b, axis=0)
For a and b with a shape of (100000, 100) this takes 0.344 seconds in the first variant, and 0.062 seconds in the second.
In the case presented in your question the following does what you want:
sol['ems'][naxis:] = numpy.ravel(
numpy.repeat(
data['ems'][naxis:,ispec,numpy.newaxis] * ems_mod,
nplanes,
axis=1
),
order='F'
)
This could be further optimized with some tricks, but that would reduce clarity and is probably premature optimization because:
plain python: 0.064 seconds
numpy: 0.002 seconds
The solution works as follows:
Your original version contains jp = naxis + ip which merely skips the first naxis elements [naxis:] selects all but the first naxis elements. Your inner loop repeats the value of data[jp,ispec] for nplanes times and writes it to multiple locations ip3d = jp + npoints_per_plane * ipl which is equivalent to a flattened 2D array offset by naxis. Therefore a second dimension is added via numpy.newaxis to the (previously 1D) data['ems'][naxis:, ispec], the values are repeated nplanes times along this new dimension via numpy.repeat. The resulting 2D array is then flattened again via numpy.ravel (in Fortran order, i.e., with the lowest axis having the smallest stride) and written to the appropriate subarray of sol['ems']. If the target array was actually 2D, the repeat could be skipped by using automatic array broadcasting.
If you run into a situation where you cannot avoid using loops, you could use Cython (which supports efficient buffer views on numpy arrays).
