Reading datasets from files to lists

Question

I'm working on a simple program and am trying to make it faster.

import time
import numpy as np
import matplotlib.pyplot as plt
import pp

first_file = open("ac_data.dat", 'r') #real dataset includes about 20000 terms
res_file = open("res.dat", 'w')
times_file = open("times.dat", 'a')

#times = [0.000000, 0.000500, 0.001000, 0.001500, 0.002000, 0.002500, 0.003000, 0.003500, 0.004000, 0.004500, 0.005000, 0.005500, 0.006000, 0.006500, 0.007000, 0.007500, 0.008000, 0.008500, 0.009000, 0.009500]
#magn = [-13.876622, -10.014824, -16.356894, -11.639914, -13.103313, -14.335239, -12.250072, -10.727098, -8.701272, -9.632907, -9.673712, -10.541722, -14.075446, -13.097790, -12.495679, -10.322924, -14.979391, -14.895666, -11.874325, -9.287736]

times = []
magn = []

for i in first_file:
  dat = [float(j) for j in i.split()]
  times.append(dat[0])
  magn.append(dat[1])

length = len(magn) #supposed to be equal 20000 for original data

#autocorrelation function
def rxx_func(amp):
    N = len(amp)
    rxx = [0]*N

    for m in xrange(N):
        for n in xrange(N-m):
            rxx[m]+=amp[n]*amp[n+m]
    return rxx

#just prove with in-built   
def autocorr(x): 
    result = np.correlate(x, x, mode='full')
    return result[result.size/2:]

#Parallel or ordinary?
answer = int(raw_input('Non-arallel = 0, parallel = 1 '))

if answer == 0:
    print 'Non-parallel calc was started'
    start = time.time()

    rxx = rxx_func(magn)

    end = time.time()
    calc_time = end - start

    time_string = 'Non-parallel: N = %i    T = %f\n'%(length,calc_time)


else:
    print 'Parallel calc was started'
    ppservers = ()
    ncp = 4
    job_server = pp.Server(ncp, ppservers=ppservers)
    print "Starting pp with", job_server.get_ncpus(), "workers"

    arg_n = tuple(magn)
    job = job_server.submit(rxx_func, (arg_n,), (), ())

    start = time.time()

    rxx = job()
    job_server.print_stats()

    end = time.time()
    calc_time = end - start

    time_string = 'Parallel with %i CPUs: N = %i    T = %f\n'%(ncp,length,calc_time)


print (" \n Task for %i terms takes %f seconds for calc" %(length, calc_time))
print (" Max value of Autocorrelation func achieves %f" %(max(rxx)))
print (" And it'll be normalized to 50 \n")

#normalization to 50
norm_const = 50/max(rxx)

proves = autocorr(magn)

for k in xrange(length):
    rxx[k] = rxx[k]*norm_const
    proves[k] = proves[k]*norm_const

#plotting
plt.plot(times, magn)      
plt.plot(times, rxx)        
plt.plot(times, proves,'*')

plt.show()


for j in xrange(length): #saving results
    st = '%f    %f    %f\n'%(times[j], rxx[j], proves[j])
    res_file.write(st)

times_file.write(time_string) #saving calc times for comparsion

first_file.close()
res_file.close()
times_file.close()

The program reads datasets from files to lists. Typical datasets look like commented #magn and #time. Real datasets will consist of 20000 lines or more.

I've trying use the "parallel Python" package, but it runs even more slower than non-parallel code.

For example, some results:

Non-parallel: N = 20000    T = 74.530000 sec
Parallel with 2 CPUs: N = 20000    T = 80.229000 sec
Parallel with 4 CPUs: N = 20000    T = 80.594000 sec

And I can't figure out why. Maybe I don't understand how it must be used.

Answer 1

I don't do much python so I might be wrong but to my understanding, this:

job = job_server.submit(rxx_func, (arg_n,), (), ())

just submits a single job but it doesn't automatically parallelize it. If you want to process the input in parallel you need to submit n jobs each working on 1-nth of the input and then combine the results. I think your parallel execution code should look something like this:

slice_size = len(magn) / ncp;

# submit a job for each chunk
jobs = [job_server.submit(rxx_func, (magn,(i-1)*slice_size, slice_size), (), ()) for i in xrange(ncp)]

# combine the results into one list, requires #import itertools
rxx = list(itertools.chain.from_iterable([job() for job in jobs]))

You will have to change your rxx_func to accept a start index and a count which defines for how many items it is responsible:

def rxx_func(amp, start_index, count):
    N = len(amp)
    rxx = [0]*count

    for m in xrange(start_index, start_index + count - 1):
        for n in xrange(N-m):
            rxx[m]+=amp[n]*amp[n+m]
    return rxx

I'm sure there is plenty which can be optimized in the above.