fixed numerical problem in BFGS optimization

stamatak · Aug 5, 2015 · 9a02452 · 9a02452
1 parent 4319e00
commit 9a02452
Show file tree

Hide file tree

Showing 7 changed files with 81 additions and 73 deletions.
diff --git a/Makefile.SSE3.PTHREADS.gcc b/Makefile.SSE3.PTHREADS.gcc
@@ -3,7 +3,7 @@
 
 CC = gcc 
 
-CFLAGS = -D_USE_PTHREADS -D__SIM_SSE3 -O2 -D_GNU_SOURCE -msse3 -fomit-frame-pointer -funroll-loops  #-Wall -pedantic -Wunused-parameter -Wredundant-decls  -Wreturn-type  -Wswitch-default -Wunused-value -Wimplicit  -Wimplicit-function-declaration  -Wimplicit-int -Wimport  -Wunused  -Wunused-function  -Wunused-label -Wno-int-to-pointer-cast -Wbad-function-cast  -Wmissing-declarations -Wmissing-prototypes  -Wnested-externs  -Wold-style-definition -Wstrict-prototypes  -Wdeclaration-after-statement -Wpointer-sign -Wextra -Wredundant-decls -Wunused -Wunused-function -Wunused-parameter -Wunused-value  -Wunused-variable -Wformat  -Wformat-nonliteral -Wparentheses -Wsequence-point -Wuninitialized -Wundef -Wbad-function-cast
+CFLAGS = -D_USE_PTHREADS -D__SIM_SSE3 -D_GNU_SOURCE -msse3 #-O2 -fomit-frame-pointer -funroll-loops  #-Wall -pedantic -Wunused-parameter -Wredundant-decls  -Wreturn-type  -Wswitch-default -Wunused-value -Wimplicit  -Wimplicit-function-declaration  -Wimplicit-int -Wimport  -Wunused  -Wunused-function  -Wunused-label -Wno-int-to-pointer-cast -Wbad-function-cast  -Wmissing-declarations -Wmissing-prototypes  -Wnested-externs  -Wold-style-definition -Wstrict-prototypes  -Wdeclaration-after-statement -Wpointer-sign -Wextra -Wredundant-decls -Wunused -Wunused-function -Wunused-parameter -Wunused-value  -Wunused-variable -Wformat  -Wformat-nonliteral -Wparentheses -Wsequence-point -Wuninitialized -Wundef -Wbad-function-cast
 
 
 LIBRARIES = -lm -pthread 

diff --git a/Makefile.SSE3.gcc b/Makefile.SSE3.gcc
@@ -3,7 +3,7 @@
 
 CC = gcc 
 
-CFLAGS = -D__SIM_SSE3  -msse3 -D_GNU_SOURCE -O2 -fomit-frame-pointer -funroll-loops  #-pedantic -Wall  -Wunused-parameter -Wredundant-decls  -Wreturn-type  -Wswitch-default -Wunused-value -Wimplicit  -Wimplicit-function-declaration  -Wimplicit-int -Wimport  -Wunused  -Wunused-function  -Wunused-label -Wno-int-to-pointer-cast -Wbad-function-cast  -Wmissing-declarations -Wmissing-prototypes  -Wnested-externs  -Wold-style-definition -Wstrict-prototypes -Wpointer-sign -Wextra -Wredundant-decls -Wunused -Wunused-function -Wunused-parameter -Wunused-value  -Wunused-variable -Wformat  -Wformat-nonliteral -Wparentheses -Wsequence-point -Wuninitialized -Wundef -Wbad-function-cast
+CFLAGS = -D__SIM_SSE3   -msse3 -D_GNU_SOURCE #-O2 -fomit-frame-pointer -funroll-loops  #-pedantic -Wall  -Wunused-parameter -Wredundant-decls  -Wreturn-type  -Wswitch-default -Wunused-value -Wimplicit  -Wimplicit-function-declaration  -Wimplicit-int -Wimport  -Wunused  -Wunused-function  -Wunused-label -Wno-int-to-pointer-cast -Wbad-function-cast  -Wmissing-declarations -Wmissing-prototypes  -Wnested-externs  -Wold-style-definition -Wstrict-prototypes -Wpointer-sign -Wextra -Wredundant-decls -Wunused -Wunused-function -Wunused-parameter -Wunused-value  -Wunused-variable -Wformat  -Wformat-nonliteral -Wparentheses -Wsequence-point -Wuninitialized -Wundef -Wbad-function-cast
 
 LIBRARIES = -lm
 

diff --git a/axml.c b/axml.c
@@ -7991,6 +7991,16 @@ void printModelParams(tree *tr, analdef *adef)
 		int 
 		  k;
 
+		printBothOpen("\nLG4X rates: ");
+		for(k = 0; k < 4; k++)
+		  printBothOpen("%f ", tr->partitionData[model].gammaRates[k]);		 
+
+		printBothOpen("\n\nLG4X weights: ");
+		for(k = 0; k < 4; k++)
+		  printBothOpen("%f ", tr->partitionData[model].weights[k]);	
+
+		printBothOpen("\n\n");
+
 		for(k = 0; k < 4; k++)
 		  {
 		    printBothOpen("LGM %d\n", k);
@@ -8764,7 +8774,7 @@ static void copyLG4(tree *localTree, tree *tr, int model, const partitionLengths
     {
       int 
 	k;
-
+          
       for(k = 0; k < 4; k++)
 	{
 	   memcpy(localTree->partitionData[model].EIGN_LG4[k],        tr->partitionData[model].EIGN_LG4[k],        pl->eignLength * sizeof(double));
@@ -8881,8 +8891,11 @@ static void execFunction(tree *tr, tree *localTree, int tid, int n)
 	  dlnLdlz[NUM_BRANCHES],
 	  d2lnLdlz2[NUM_BRANCHES];
 
-	memcpy(localTree->coreLZ,   tr->coreLZ,   sizeof(double) *  localTree->numBranches);
-	memcpy(localTree->executeModel, tr->executeModel, sizeof(boolean) * localTree->NumberOfModels);
+	if(tid > 0)
+	  {
+	    memcpy(localTree->coreLZ,   tr->coreLZ,   sizeof(double) *  localTree->numBranches);
+	    memcpy(localTree->executeModel, tr->executeModel, sizeof(boolean) * localTree->NumberOfModels);
+	  }
 
 	execCore(localTree, dlnLdlz, d2lnLdlz2);
 

diff --git a/axml.h b/axml.h
@@ -1210,6 +1210,9 @@ typedef struct
 
 /****************************** FUNCTIONS ****************************************************/
 
+extern void optimizeRatesBFGS(tree *tr);
+extern void setRateModel(tree *tr, int model, double rate, int position);
+
 extern void ascertainmentBiasSequence(unsigned char tip[32], int numStates, int dataType, int nodeNumber, int *ascMissingVector);
 
 extern void computePlacementBias(tree *tr, analdef *adef);

diff --git a/eigen.c b/eigen.c
@@ -31,7 +31,9 @@
 
 #include <math.h>
 #include <assert.h>
-
+#include "axml.h"
+#include <string.h>
+#include <stdlib.h>
 
 static void mytred2(double **a, const int n, double *d, double *e)
 {
@@ -198,3 +200,5 @@ void makeEigen(double **_a, const int n, double *d, double *e)
   mytred2(_a, n, d, e);
   mytqli(d, e, n, _a);
 }
+
+
diff --git a/newviewGenericSpecial.c b/newviewGenericSpecial.c
@@ -3414,7 +3414,7 @@ static void newviewGTRGAMMA(int tipCase,
 		 if(useFastScaling)
 		   addScale += wgt[i];
 		 else
-		   ex3[i]  += 1;
+		   ex3[i]  += 1;	     
 	       }
 	     else
 	       {
@@ -3594,7 +3594,7 @@ static void newviewGTRGAMMA(int tipCase,
 	     if(useFastScaling)
 	       addScale += wgt[i];
 	     else
-	       ex3[i]  += 1;
+	       ex3[i]  += 1;	     
 	   }
 	 else
 	   {