I am not sure if this issue is already known and is being fixed with the upcoming refactoring. After getting around the known problem with the ArrayIndexOutOfBoundsException when using the -Strace option, I realized that my output trees looked very strange. The branch lengths in one of my two populations were always zero. My command line is:
I started wondering if the problem could be in the way the external trace data is read and processed by msms. I only got around the ArrayIndexOutOfBoundsException when I defined the trace up to time 0, i.e. the last entry in the trace was always 0 0 1 0 1. It seemed that msms takes the frequency(A) of 0 instead of 1 at time 0 for the second population, which causes all lineages from that population to coalesce immediately. I then started looking through the code and found the following two issues:
at.mabs.simulators.ReaderSelection: line 79-89
double[] fs = lastPastward.frequencys;
FrequencyState state = new FrequencyState(fs.length, fs.length);
boolean fixation=true;
for (int i = 0; i < fs.length; i++) {
state.setFrequency(i, 1, fs[i]);
state.setFrequency(i, 0, 1 - fs[i]);
if(fs[i]>0 && fs[i]<1){
fixation=false;
}
}
But fs is in the format freq(a)_pop1, freq(A)_pop1, freq(a)_pop2, freq(A)_pop2. If I understand this correctly, then it should rather be something like this:
double[] fs = lastPastward.frequencys;
FrequencyState state = new FrequencyState(fs.length/2, fs.length/2);
boolean fixation=true;
for (int i = 0; i < fs.length/2; i++) {
state.setFrequency(i, 1, fs[2_i+1]);
state.setFrequency(i, 0, 1 - fs[2_i+1]);
if(fs[2_i+1]>0 && fs[2_i+1]<1){
fixation=false;
}
}
Similar problem here:
at.mabs.simulators.ReaderSelection: line 162-175
int index = 1;
double[] entryData = e.frequencys;
double[] nextEntryData = next.frequencys;
for (int d = 0; d < demeCount; d++) {
double total = 0;
double f = entryData[index] * (1 - t) + nextEntryData[index] * t;
index++;
total += f;
array[d] = f;
}
Line 169 should read index+=2 instead.
I hope that I am not totally wrong with this, but my results make a lot more sense since I did these corrections.
I am not sure if this issue is already known and is being fixed with the upcoming refactoring. After getting around the known problem with the ArrayIndexOutOfBoundsException when using the -Strace option, I realized that my output trees looked very strange. The branch lengths in one of my two populations were always zero. My command line is:
java -Xmx256M -jar ./msms.jar 20 1 -N 10000 -t 5 -I 2 10 10 -ej 10 2 1 -SI 0.25 2 0 0 -Strace trace.txt -T
I started wondering if the problem could be in the way the external trace data is read and processed by msms. I only got around the ArrayIndexOutOfBoundsException when I defined the trace up to time 0, i.e. the last entry in the trace was always 0 0 1 0 1. It seemed that msms takes the frequency(A) of 0 instead of 1 at time 0 for the second population, which causes all lineages from that population to coalesce immediately. I then started looking through the code and found the following two issues:
at.mabs.simulators.ReaderSelection: line 79-89
double[] fs = lastPastward.frequencys; FrequencyState state = new FrequencyState(fs.length, fs.length); boolean fixation=true; for (int i = 0; i < fs.length; i++) { state.setFrequency(i, 1, fs[i]); state.setFrequency(i, 0, 1 - fs[i]); if(fs[i]>0 && fs[i]<1){ fixation=false; } }
But fs is in the format freq(a)_pop1, freq(A)_pop1, freq(a)_pop2, freq(A)_pop2. If I understand this correctly, then it should rather be something like this:
double[] fs = lastPastward.frequencys; FrequencyState state = new FrequencyState(fs.length/2, fs.length/2); boolean fixation=true; for (int i = 0; i < fs.length/2; i++) { state.setFrequency(i, 1, fs[2_i+1]); state.setFrequency(i, 0, 1 - fs[2_i+1]); if(fs[2_i+1]>0 && fs[2_i+1]<1){ fixation=false; } }
Similar problem here:
at.mabs.simulators.ReaderSelection: line 162-175
int index = 1; double[] entryData = e.frequencys; double[] nextEntryData = next.frequencys; for (int d = 0; d < demeCount; d++) { double total = 0; double f = entryData[index] * (1 - t) + nextEntryData[index] * t; index++; total += f; array[d] = f; }
Line 169 should read index+=2 instead.
I hope that I am not totally wrong with this, but my results make a lot more sense since I did these corrections.
Thanks a lot for your efforts!