diff -Nru raxml-8.2.3/axml.c raxml-8.2.4/axml.c --- raxml-8.2.3/axml.c 2015-08-12 07:47:19.000000000 +0000 +++ raxml-8.2.4/axml.c 2015-10-02 13:41:21.000000000 +0000 @@ -3360,8 +3360,7 @@ tr->partitionData[model].ascExpVector = (int *)rax_calloc(((size_t)tr->innerNodes) * ((size_t)tr->partitionData[model].states), sizeof(int)); - - tr->partitionData[model].ascMissingVector = (int *)rax_calloc((size_t)(tr->mxtips + 1), sizeof(int)); + tr->partitionData[model].ascSumBuffer = (double *)rax_malloc(((size_t)tr->partitionData[model].ascOffset) * sizeof(double)); @@ -4773,7 +4772,7 @@ printf(" [--flag-check][--auto-prot=ml|bic|aic|aicc]\n"); printf(" [--epa-keep-placements=number][--epa-accumulated-threshold=threshold]\n"); printf(" [--epa-prob-threshold=threshold]\n"); - printf(" [--JC69][--K80][--HKY85][--asc-miss=fraction]\n"); + printf(" [--JC69][--K80][--HKY85]\n"); #if (defined(_WAYNE_MPI) && defined(_USE_PTHREADS)) printf(" [--set-thread-affinity]\n"); #endif @@ -5212,9 +5211,7 @@ printf(" --HKY85 specify that all DNA partitions will evolve under the HKY85 model, this overrides all other model specifications for DNA partitions.\n"); printf("\n"); printf(" DEFAULT: Off\n"); - printf("\n"); - printf(" --asc-miss=fraction specify the fraction of missing data in the variable sites of the alignment you are trying to correct for \n"); - printf(" ascertainment bias. Experimental option, needs to be tested!\n"); + printf("\n"); #if (defined(_WAYNE_MPI) && defined(_USE_PTHREADS)) printf("\n"); printf(" --set-thread-affinity specify that thread-to-core affinity shall be set by RAxML for the hybrid MPI-PThreads version\n"); @@ -5350,10 +5347,6 @@ tr->useJC69 = FALSE; tr->useHKY85 = FALSE; - //asecryainment bias correction missing data - - tr->useAscMissing = FALSE; - #ifdef _BASTIEN tr->doBastienStuff = FALSE; #endif @@ -5375,7 +5368,7 @@ while(1) { static struct - option long_options[16] = + option long_options[15] = { {"mesquite", no_argument, &flag, 1}, {"silent", no_argument, &flag, 1}, @@ -5389,8 +5382,7 @@ {"epa-prob-threshold", required_argument, &flag, 1}, {"JC69", no_argument, &flag, 1}, {"K80", no_argument, &flag, 1}, - {"HKY85", no_argument, &flag, 1}, - {"asc-miss", required_argument, &flag, 1}, + {"HKY85", no_argument, &flag, 1}, {"set-thread-affinity", no_argument, &flag, 1}, {0, 0, 0, 0} }; @@ -5573,25 +5565,8 @@ case 12: tr->useHKY85 = TRUE; break; - break; + break; case 13: - if(sscanf(optarg, "%lf", &(tr->ascMissing)) != 1) - { - printf("\nError parsing ascertainment bias missing data fraction correction, RAxML expects a floating point value > 0.0 and < 1.0\n\n"); - errorExit(-1); - } - if(tr->ascMissing <= 0.0 || tr->ascMissing >= 1.0) - { - printf("\nError parsing ascertainment missing data fraction correction, RAxML expects a floating point value >= 0.0 and <= 1.0\n\n"); - errorExit(-1); - } -#ifdef _USE_PTHREADS - printf("\nError: ascertainment missing data fraction correction not implemented for PThreads version yet, exiting!\n\n"); - errorExit(-1); -#endif - tr->useAscMissing = TRUE; - break; - case 14: #if (defined(_WAYNE_MPI) && defined(_USE_PTHREADS)) adef->setThreadAffinity = TRUE; #else @@ -8693,9 +8668,7 @@ sizeof(double)); tr->partitionData[model].ascExpVector = (int *)rax_calloc(((size_t)tr->innerNodes) * ((size_t)tr->partitionData[model].states), - sizeof(int)); - - tr->partitionData[model].ascMissingVector = (int *)rax_calloc((size_t)(tr->mxtips + 1), sizeof(int)); + sizeof(int)); tr->partitionData[model].ascSumBuffer = (double *)rax_malloc(((size_t)tr->partitionData[model].ascOffset) * sizeof(double)); @@ -12902,32 +12875,7 @@ } -static void setupAscMissing(tree *tr, analdef *adef) - { - int - model, - *perm = (int*)rax_malloc(sizeof(int) * (size_t)(tr->mxtips + 1)); - - for(model = 0; model < tr->NumberOfModels; model++) - { - if(tr->partitionData[model].ascBias && tr->useAscMissing) - { - int - cutoff = (int)(tr->ascMissing * (double)(tr->mxtips) + 0.5), - i; - makePermutation(perm, 1, tr->mxtips, adef); - - for(i = 1; i <= cutoff; i++) - tr->partitionData[model].ascMissingVector[perm[i]] = 1; - } - } - - //TODO need to add barrier for copying ascMissingVector to all threads - //and we also need to copy the flag value! - - rax_free(perm); - } @@ -13080,7 +13028,7 @@ if(!adef->readTaxaOnly) { int - countNonSev = 0, + countNonSev = 0, countLG4 =0; assert(countAscBias == 0); @@ -13095,7 +13043,8 @@ for(i = 0; i < tr->NumberOfModels; i++) { if(!(tr->partitionData[i].dataType == AA_DATA || tr->partitionData[i].dataType == DNA_DATA)) - countNonSev++; + countNonSev++; + if(tr->partitionData[i].protModels == LG4 || tr->partitionData[i].protModels == LG4X) { @@ -13213,6 +13162,13 @@ setRateHetAndDataIncrement(tr, adef); + if(tr->rateHetModel == GAMMA_I && tr->saveMemory) + { + printf("\nError: Memory saving option \"-U\" not implemented for models with proportion\n"); + printf("of variable site estimates, reomve \"-U\" from your command line and re-run, exiting now.\n\n"); + errorExit(-1); + } + if(countAscBias > 0 && !adef->readTaxaOnly) { if(tr->ascertainmentCorrectionType == NOT_DEFINED) @@ -13238,8 +13194,6 @@ readAscFiles(tr); - setupAscMissing(tr, adef); - if(!adef->readTaxaOnly) { #ifdef _USE_PTHREADS diff -Nru raxml-8.2.3/axml.h raxml-8.2.4/axml.h --- raxml-8.2.3/axml.h 2015-08-12 07:47:19.000000000 +0000 +++ raxml-8.2.4/axml.h 2015-10-02 13:41:21.000000000 +0000 @@ -168,9 +168,9 @@ #define PointGamma(prob,alpha,beta) PointChi2(prob,2.0*(alpha))/(2.0*(beta)) #define programName "RAxML" -#define programVersion "8.2.3" -#define programVersionInt 8230 -#define programDate "August 12 2015" +#define programVersion "8.2.4" +#define programVersionInt 8240 +#define programDate "October 02 2015" #define TREE_EVALUATION 0 @@ -604,7 +604,6 @@ boolean ascBias; int ascOffset; int *ascExpVector; - int *ascMissingVector; double *ascSumBuffer; double *ascVector; double ascScaler[64]; @@ -1066,9 +1065,7 @@ boolean doSubtreeEPA; - double ascMissing; - boolean useAscMissing; - + } tree; @@ -1210,7 +1207,7 @@ -extern void ascertainmentBiasSequence(unsigned char tip[32], int numStates, int dataType, int nodeNumber, int *ascMissingVector); +extern void ascertainmentBiasSequence(unsigned char tip[32], int numStates, int dataType, int nodeNumber); extern void computePlacementBias(tree *tr, analdef *adef); diff -Nru raxml-8.2.3/debian/changelog raxml-8.2.4/debian/changelog --- raxml-8.2.3/debian/changelog 2015-09-22 19:29:31.000000000 +0000 +++ raxml-8.2.4/debian/changelog 2015-11-09 08:56:02.000000000 +0000 @@ -1,3 +1,9 @@ +raxml (8.2.4-1) unstable; urgency=medium + + * New upstream version + + -- Andreas Tille Mon, 09 Nov 2015 09:56:00 +0100 + raxml (8.2.3-1) unstable; urgency=medium * New upstream version diff -Nru raxml-8.2.3/evaluateGenericSpecial.c raxml-8.2.4/evaluateGenericSpecial.c --- raxml-8.2.3/evaluateGenericSpecial.c 2015-08-12 07:47:19.000000000 +0000 +++ raxml-8.2.4/evaluateGenericSpecial.c 2015-10-02 13:41:21.000000000 +0000 @@ -58,7 +58,7 @@ -void ascertainmentBiasSequence(unsigned char tip[32], int numStates, int dataType, int nodeNumber, int *ascMissingVector) +void ascertainmentBiasSequence(unsigned char tip[32], int numStates, int dataType, int nodeNumber) { assert(numStates <= 32 && numStates > 1); @@ -66,59 +66,33 @@ switch(dataType) { - case BINARY_DATA: - if(ascMissingVector[nodeNumber] == 1) - { - tip[0] = 3; - tip[1] = 3; - } - else - { - tip[0] = 1; - tip[1] = 2; - } + case BINARY_DATA: + tip[0] = 1; + tip[1] = 2; break; - case DNA_DATA: - if(ascMissingVector[nodeNumber] == 1) - { - //printf("M %d\n", nodeNumber); - tip[0] = 15; - tip[1] = 15; - tip[2] = 15; - tip[3] = 15; - } - else - { - //printf("D %d\n", nodeNumber); - tip[0] = 1; - tip[1] = 2; - tip[2] = 4; - tip[3] = 8; - } + case DNA_DATA: + tip[0] = 1; + tip[1] = 2; + tip[2] = 4; + tip[3] = 8; break; - case AA_DATA: - if(ascMissingVector[nodeNumber] == 1) - assert(0); - else - { - int - i; - - for(i = 0; i < numStates; i++) - tip[i] = i; - } + case AA_DATA: + { + int + i; + + for(i = 0; i < numStates; i++) + tip[i] = i; + } break; - case GENERIC_32: - if(ascMissingVector[nodeNumber] == 1) - assert(0); - else - { - int - i; - - for(i = 0; i < numStates; i++) - tip[i] = i; - } + case GENERIC_32: + { + int + i; + + for(i = 0; i < numStates; i++) + tip[i] = i; + } break; default: assert(0); @@ -293,7 +267,7 @@ double *x1, double *x2, double *tipVector, unsigned char *tipX1, int n, double *diagptable, const int numStates, - double *accumulator, double *weightVector, int dataType, int nodeNumber, int *ascMissingVector, + double *accumulator, double *weightVector, int dataType, int nodeNumber, double *goldmanAccumulator) { double @@ -315,7 +289,7 @@ unsigned char tip[32]; - ascertainmentBiasSequence(tip, numStates, dataType, nodeNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, nodeNumber); for (i = 0; i < n; i++) { @@ -401,7 +375,7 @@ double *x1, double *x2, double *tipVector, unsigned char *tipX1, int n, double *diagptable, const int numStates, - double *accumulator, double *weightVector, int dataType, int nodeNumber, int *ascMissingVector, + double *accumulator, double *weightVector, int dataType, int nodeNumber, double *goldmanAccumulator) { double @@ -427,7 +401,7 @@ unsigned char tip[32]; - ascertainmentBiasSequence(tip, numStates, dataType, nodeNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, nodeNumber); for (i = 0; i < n; i++) { @@ -3304,7 +3278,9 @@ #endif break; case GAMMA_I: - { + { + assert(!tr->saveMemory); + calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable); partitionLikelihood = evaluateGTRGAMMAINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant, @@ -3366,6 +3342,8 @@ break; case GAMMA_I: { + assert(!tr->saveMemory); + calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable); partitionLikelihood = evaluateGTRGAMMAPROTINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant, @@ -3572,13 +3550,13 @@ correction = evaluateCatAsc(ex1_asc, ex2_asc, x1_start_asc, x2_start_asc, tr->partitionData[model].tipVector, tip, ascWidth, diagptable, ascWidth, &accumulator, weightVector, tr->partitionData[model].dataType, - tipNodeNumber, tr->partitionData[model].ascMissingVector, &goldmanAccumulator); + tipNodeNumber, &goldmanAccumulator); } break; case GAMMA: correction = evaluateGammaAsc(ex1_asc, ex2_asc, x1_start_asc, x2_start_asc, tr->partitionData[model].tipVector, tip, ascWidth, diagptable, ascWidth, &accumulator, weightVector, tr->partitionData[model].dataType, - tipNodeNumber, tr->partitionData[model].ascMissingVector, &goldmanAccumulator); + tipNodeNumber, &goldmanAccumulator); break; default: assert(0); diff -Nru raxml-8.2.3/makenewzGenericSpecial.c raxml-8.2.4/makenewzGenericSpecial.c --- raxml-8.2.3/makenewzGenericSpecial.c 2015-08-12 07:47:19.000000000 +0000 +++ raxml-8.2.4/makenewzGenericSpecial.c 2015-10-02 13:41:21.000000000 +0000 @@ -2727,7 +2727,7 @@ #endif static void sumCatAsc(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector, - int n, const int numStates, int dataType, int qNumber, int rNumber, int *ascMissingVector, int maxtips) + int n, const int numStates, int dataType, int qNumber, int rNumber, int maxtips) { int i, k; double *left, *right, *sum; @@ -2740,8 +2740,8 @@ tip1[32], tip2[32]; - ascertainmentBiasSequence(tip1, numStates, dataType, qNumber, ascMissingVector); - ascertainmentBiasSequence(tip2, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip1, numStates, dataType, qNumber); + ascertainmentBiasSequence(tip2, numStates, dataType, rNumber); for(i = 0; i < n; i++) { @@ -2761,9 +2761,9 @@ tip[32]; if(rNumber <= maxtips) - ascertainmentBiasSequence(tip, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, rNumber); else - ascertainmentBiasSequence(tip, numStates, dataType, qNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, qNumber); for(i = 0; i < n; i++) { @@ -2795,7 +2795,7 @@ } static void sumGammaAsc(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector, - int n, const int numStates, int dataType, int qNumber, int rNumber, int *ascMissingVector, int maxtips) + int n, const int numStates, int dataType, int qNumber, int rNumber, int maxtips) { int i, l, k; double *left, *right, *sum; @@ -2811,8 +2811,8 @@ tip1[32], tip2[32]; - ascertainmentBiasSequence(tip1, numStates, dataType, qNumber, ascMissingVector); - ascertainmentBiasSequence(tip2, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip1, numStates, dataType, qNumber); + ascertainmentBiasSequence(tip2, numStates, dataType, rNumber); for(i = 0; i < n; i++) { @@ -2834,9 +2834,9 @@ tip[32]; if(rNumber < maxtips) - ascertainmentBiasSequence(tip, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, rNumber); else - ascertainmentBiasSequence(tip, numStates, dataType, qNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, qNumber); for(i = 0; i < n; i++) { @@ -4207,11 +4207,11 @@ { case CAT: sumCatAsc(tipCase, tr->partitionData[model].ascSumBuffer, x1_start_asc, x2_start_asc, tr->partitionData[model].tipVector, - states, states, tr->partitionData[model].dataType, pNumber, qNumber, tr->partitionData[model].ascMissingVector, tr->mxtips); + states, states, tr->partitionData[model].dataType, pNumber, qNumber, tr->mxtips); break; case GAMMA: sumGammaAsc(tipCase, tr->partitionData[model].ascSumBuffer, x1_start_asc, x2_start_asc, tr->partitionData[model].tipVector, - states, states, tr->partitionData[model].dataType, pNumber, qNumber, tr->partitionData[model].ascMissingVector, tr->mxtips); + states, states, tr->partitionData[model].dataType, pNumber, qNumber, tr->mxtips); break; default: assert(0); Binary files /tmp/7mF0WhCjvV/raxml-8.2.3/manual/NewManual.odt and /tmp/s7iTTDUSsL/raxml-8.2.4/manual/NewManual.odt differ Binary files /tmp/7mF0WhCjvV/raxml-8.2.3/manual/NewManual.pdf and /tmp/s7iTTDUSsL/raxml-8.2.4/manual/NewManual.pdf differ diff -Nru raxml-8.2.3/newviewGenericSpecial.c raxml-8.2.4/newviewGenericSpecial.c --- raxml-8.2.3/newviewGenericSpecial.c 2015-08-12 07:47:19.000000000 +0000 +++ raxml-8.2.4/newviewGenericSpecial.c 2015-10-02 13:41:21.000000000 +0000 @@ -528,7 +528,7 @@ double *x1, double *x2, double *x3, double *extEV, double *tipVector, int *ex3, const int n, double *left, double *right, - const int numStates, int dataType, const int qNumber, const int rNumber, int *ascMissingVector, int maxtips) + const int numStates, int dataType, const int qNumber, const int rNumber, int maxtips) { double *le, *ri, *v, *vl, *vr, @@ -545,8 +545,8 @@ tip1[32], tip2[32]; - ascertainmentBiasSequence(tip1, numStates, dataType, qNumber, ascMissingVector); - ascertainmentBiasSequence(tip2, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip1, numStates, dataType, qNumber); + ascertainmentBiasSequence(tip2, numStates, dataType, rNumber); assert(qNumber <= maxtips && rNumber <= maxtips); @@ -587,9 +587,9 @@ tip[32]; if(rNumber <= maxtips) - ascertainmentBiasSequence(tip, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, rNumber); else - ascertainmentBiasSequence(tip, numStates, dataType, qNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, qNumber); for (i = 0; i < n; i++) { @@ -690,7 +690,7 @@ double *x1, double *x2, double *x3, double *extEV, double *tipVector, int *ex3, const int n, double *left, double *right, - const int numStates, int dataType, const int qNumber, const int rNumber, int *ascMissingVector, int maxtips) + const int numStates, int dataType, const int qNumber, const int rNumber, int maxtips) { int @@ -713,8 +713,8 @@ tip1[32], tip2[32]; - ascertainmentBiasSequence(tip1, numStates, dataType, qNumber, ascMissingVector); - ascertainmentBiasSequence(tip2, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip1, numStates, dataType, qNumber); + ascertainmentBiasSequence(tip2, numStates, dataType, rNumber); assert(qNumber <= maxtips && rNumber <= maxtips); @@ -753,9 +753,9 @@ tip[32]; if(rNumber <= maxtips) - ascertainmentBiasSequence(tip, numStates, dataType, rNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, rNumber); else - ascertainmentBiasSequence(tip, numStates, dataType, qNumber, ascMissingVector); + ascertainmentBiasSequence(tip, numStates, dataType, qNumber); for (i = 0; i < n; i++) { @@ -8622,7 +8622,7 @@ tr->partitionData[model].EV, tr->partitionData[model].tipVector, ex3_asc, states, left, right, states, tr->partitionData[model].dataType, - tInfo->qNumber, tInfo->rNumber, tr->partitionData[model].ascMissingVector, tr->mxtips); + tInfo->qNumber, tInfo->rNumber, tr->mxtips); } break; case GAMMA: @@ -8632,7 +8632,7 @@ tr->partitionData[model].tipVector, ex3_asc, states, left, right, states, tr->partitionData[model].dataType, - tInfo->qNumber, tInfo->rNumber, tr->partitionData[model].ascMissingVector, tr->mxtips); + tInfo->qNumber, tInfo->rNumber, tr->mxtips); break; default: assert(0);