Double Loop in RASCaltulator can be dramatically sped up

[stschiff]

One of the critical speed-bottlenecks in RASCalculator is arguably the double loop that happens for every site of the input file:

RAStools/RASCalculator.py

Lines 134 to 135 in a9c66c6

	for Lftidx, leftPop in enumerate(LeftPops):
	for Tstidx, testPop in enumerate(TestPops):

Effectively, this means that the algorithm scales quadratically with the number of populations (or individuals, for ungrouped input).

I think we can speed this up to scale effectively linear, as long as the allele frequency is low. The idea is that within the first loop, we first ask whether there is any derived allele frequency in the first population. If not, we just move on to the next outer loop item, because we know that all internal loop items would just add zero.

So I'm proposing to change this

 for Lftidx, leftPop in enumerate(LeftPops):
        for Tstidx, testPop in enumerate(TestPops):

to this:

 for Lftidx, leftPop in enumerate(LeftPops):
    if afDict[leftPop] > 0:
        for Tstidx, testPop in enumerate(TestPops):

or something to that effect.

[TCLamnidis]

Currently sites where both afDict[leftPop] and afDict[testPop] are 0 are counted in the number of observations. Additionally, the polarisation by the outgroup population means that even if both of these (leftPop and testPop) are 0, a non-zero ras will be added to the accumulator when the afDict[outgroup] is not 0.

Would this maybe be a safer change? Or do you not think this will make a difference?

 for Lftidx, leftPop in enumerate(LeftPops):
    if afDict[leftPop] == 0 and afDict[args.outgroup] == 0:
        continue
    else:
        for Tstidx, testPop in enumerate(TestPops):

with the addition of an accumulator for mj before the continue.

[stschiff]

OK, forgot about the ougroup business. Well, this is weird what you're writing about the accumulator. Do you mean for the denominator of the RAS calculation? Shouldn't it then be

for Lftidx, leftPop in enumerate(LeftPops):
    if afDict[leftPop] == 0 and afDict[args.outgroup] == 0 or afDict[leftPop] == 1 and afDict[args.outgroup] == 1:
        continue
    else:
        for Tstidx, testPop in enumerate(TestPops):

In any case, by the way, I think the outgroup thing is a real nuisance and we introduced it only to make this tool also work for ordinary F3 statistics. I would almost remove that feature, as normally with rare variants you can safely assume that the minor allele is the "derived" allele, and the major allele is the "ancestral" allele. So I would just remove the outgroup option and generally compute afDict as the minor allele frequency.

[stschiff]

I'm currently preparing a comprehensive writeup on all of this stuff (https://github.com/Schiffels-Popgen/RAS_exploration/tree/doc_writeup_stephan), and I think we don't need an outgroup. To be discussed...