# Building a report of DNA sites and chunks

Here is the slow part of my code:

def Prepare_WCfst_file():
content_WCfst = "iteration\tlocus\tFromSite\tToSite\tWCFst\n"
nbchunks = nbsites/nbsitesPerChunk
if nbchunks%1 != 0:
print("\n\nThe number of sites per chunk for WCFst is not a divisor of the number of sites per locus.\nProcess is aborted!\n\n")
sys.exit(1)
nbchunks = int(nbchunks)
nbsitesPerChunk_minus1 = nbsitesPerChunk - 1
for iteration in xrange(nbiterations):
for locus in xrange(nbloci):
FromSite = 0
ToSite = nbsitesPerChunk_minus1
for chunk in xrange(nbchunks):
content_WCfst += str(iteration) + "\t" + str(locus) + "\t" + str(FromSite) + "\t" + str(ToSite) + "\n"
FromSite = ToSite + 1
ToSite = FromSite + nbsitesPerChunk_minus1
return content_WCfst


Typically, nbiterations and nbloci take values between 1 and 100. nbchunks however can take values of the order of 10^6 to 10^9. This function will be called about 100 to 10,000 times in total.

Suggested data to benchmark

nbloci = 10^6
nbiterations = 20
nbloci = 2
nbsitesPerChunk = 10


What the function does

Basically, what this piece of code does is to create a long string content_WCfst (which will afterward be written on a file). The file needs to contain 4 columns: iteration, locus, FromSite and ToSite. The difference between the columns FromSite and ToSite is always nbsitesPerChunk_minus1. By the end the number of lines is nbiterations * nbloci * nbchunks.

My thoughts about what to improve

I think the part that can be improved is the three loops or eventually only the most internal loop (for chunk in xrange(nbchunks):)

Note that I've tried to replace

str(iteration) + "\t" + str(locus) + "\t" + str(FromSite) + "\t" + str(ToSite) + "\n"


by

"\t".join(str(iteration), str(locus), str(FromSite), str(ToSite), "\n")


but it was even slower.

Note also that the function doesn't receive any arguments and therefore, all parameters are defined globally outside this function. I don't expect it to be an issue as they are not called over and over again.

I am hoping that a solution can exist so that the looping occurs in C, using NumPy eventually.

• Use snake_case for most names. Whilst prepare_WCfst_file is allowed, Prepare_WCfst_file is not.
• Pass varibles, don't depend on ones out of scope.
• Lines should be 79 characters max. The exception is comments at 72.

You have a very good style otherwise!

You are using archaic string formatting.

content_WCfst += str(iteration) + "\t" + str(locus) + "\t" + str(FromSite) + "\t" + str(ToSite) + "\n"


Use format instead.

content_WCfst += "{!s}\t{!s}\t{!s}\t{!s}\n".format(iteration, locus, FromSite, ToSite)


for iteration in xrange(nbiterations):
for locus in xrange(nbloci):
for chunk in xrange(nbchunks):


I think the above is the offending factor. This is as it is $O(n^3)$.

But not only that you are appending to what can be a massive list. The append, may seem small early on, but later use a lot of cpu.

To fix this you can change the function to a generator. The thing xrange is.

Pros

• Fast.
• On demand output.
• Low memory usage.
• Can be used like a list in most cases.

Cons

• If you bin the output, is worse than a list.
• Can't be duplicated.
• Can't be index or sliced. [1] or [:2] do not work.

It's simple.

content_WCfst += "{!s}\t{!s}\t{!s}\t{!s}\n".format(iteration, locus, FromSite, ToSite)
# To
yield "{!s}\t{!s}\t{!s}\t{!s}\n".format(iteration, locus, FromSite, ToSite)


If you can't use a generator, then sets are also good. However they can only store hashable data, and one of it. And don't keep order.

>>> set(['a', 'c', 'a', 'b'])
{'a', 'b', 'c'}


To implement this:

content_WCfst += "{!s}\t{!s}\t{!s}\t{!s}\n".format(iteration, locus, FromSite, ToSite)
# To


This is if content_WCfst is also a set.

I can't think of a way to solve the $O(n^3)$ part. However adding to lists and strings are very slow compared to sets and generators.

• .format is faster than string sum and a generator will avoid double memory use as he wants to save it to a fle, nice! – Caridorc Jul 9 '15 at 17:23
nbchunks = nbsites/nbsitesPerChunk
if nbchunks%1 != 0:
print("\n\nThe number of sites per chunk for WCFst is not a divisor of the number of sites per locus.\nProcess is aborted!\n\n")
sys.exit(1)
nbchunks = int(nbchunks)


This chunk of code is not good because:

nbchunks = nbsites / nbsitesPerChunk


In python-2 is already int for sure, in Python-3 you can write:

nbchunks = nbsites // nbsitesPerChunk


and it will be int for sure, so you can remove

nbchunks = int(nbchunks)


to avoid clutter.

The line:

if nbchunks%1 != 0

• (minor) should be spaced like if nbchunks % 1 != 0
• (SERIOUS) is always false as the remainder between a natural number and one is always zero!

Also the error message:

print("\n\nThe number of sites per chunk for WCFst is not a divisor of the number of sites per locus.\nProcess is aborted!\n\n")


is not connected to the if statement...

Exiting explicit like sys.exit(1) should be avoided, you should instead:

raise AdequateException(message)


Strings in Python are immutable, so string concatenation is expensive. If you need to concatenate strings, then build a list and join() the elements at the end. Alternatively, yield each line from a generator.

Better yet, use a CSV writer with Tab as the delimiter, since that is just the kind of output that you are generating.