diff -NaurdbB -xCVS -xconfigure$ -xconfig.h.in -xMakefile.in mrbayes-3.1.2.old/bayes.c mrbayes-3.1.2.new/bayes.c
--- mrbayes-3.1.2.old/bayes.c	Fri Dec 23 11:26:28 EST 2005
+++ mrbayes-3.1.2.new/bayes.c	Sat Jul 07 20:06:05 EDT 2007
@@ -692,6 +692,8 @@
 	for (i=0; i<MAX_CHAINS; i++)
 		chainParams.userTemps[i] = 1.0;              /* user-defined chain temperatures               */
 	chainParams.swapAdjacentOnly = NO;               /* swap only adjacent temperatures               */
+	chainParams.useSA = NO; /* whether to try simulated annealing - Allen */
+	chainParams.useAdapt = NO; /* whether to try adapting windows - Allen */
 	chainParams.printMax = 8;						 /* maximum number of chains to print to screen   */
 	chainParams.printAll = YES;						 /* whether to print heated chains                */
 	chainParams.runWithData = YES;				     /* whether to run with data                      */
@@ -916,6 +918,7 @@
 		modelParams[j].covswitchUni[0] = 0.0;
 		modelParams[j].covswitchUni[1] = 100.0;
 		modelParams[j].covswitchExp = 1.0;
+		modelParams[j].covswitchMaxPinvar = 1.0;
 		strcpy(modelParams[j].symPiPr, "Fixed");           /* prior for pi when unidentifiable states used */
 		modelParams[j].symBetaFix = -1.0;
 		modelParams[j].symBetaUni[0] = 0.0;
diff -NaurdbB -xCVS -xconfigure$ -xconfig.h.in -xMakefile.in mrbayes-3.1.2.old/command.c mrbayes-3.1.2.new/command.c
--- mrbayes-3.1.2.old/command.c	Sun Apr 02 10:47:51 EDT 2006
+++ mrbayes-3.1.2.new/command.c	Sat Jul 07 20:07:18 EDT 2007
@@ -169,7 +169,7 @@
 int				inMrbayesBlock;        /* flag for whether we are in a mrbayes block    */
 int				inSumtBlock;           /* are we in the sumt block                      */
 int				inValidCommand;        /* a useful flag set whenever you enter a cmd    */
-int  			isInAmbig, isInPoly;   /* flags whether we are within () or {}          */
+int  			isInPoly, isInAmbig;   /* flags whether we are within () or {}          */
 int				isTranslateDef;        /* is a translation block defined                */
 int				isUserTreeDefined;     /* flag indicating whether user tree is found    */
 char			logFileName[100];      /* name of the log file                          */
@@ -210,8 +210,8 @@
 int				longIntegerSize;       /* size of an unsigned integer                   */
 
 /* local (to this file) */
-char			*tokenP, token[CMD_STRING_LENGTH];
-CmdType			commands[] =
+static char		*tokenP, token[CMD_STRING_LENGTH];
+static CmdType		commands[] =
 				{
 			/*	Information on commands initialization:
 			 
@@ -224,9 +224,11 @@
 			 		7 = Expecting (2^TokenType) (expect) (PARAMETER = 4; SEMICOLON = 32; ALPHA = 16384; 
 			 		    ALPHA | QUESTIONMARK | DASH | NUMBER | ASTERISK | EXCLAMATIONMARK | PERCENT | WEIRD | SEMICOLON = 11715360;
 			 		    ALPHA | QUESTIONMARK | DASH | NUMBER | ASTERISK | EXCLAMATIONMARK | PERCENT | WEIRD | SEMICOLON | LEFTPAR | RIGHTPAR | LEFTCURL | RIGHTCURL = 112381728;
-			 		    PARAMETER | SEMICOLON = 36; NUMBER | ALPHA = 49152; ALPHA | SEMICOLON = 16416; EQUALSIGN = 8; NUMBER = 32768)
+			PARAMETER | SEMICOLON = 36; NUMBER | ALPHA = 49152; ALPHA | SEMICOLON = 16416; EQUALSIGN = 8;
+			NUMBER = 32768)
 			 		8 = Description of the command (cmdDescription)
-			 		9 = Where should the command be used (cmdUse) (IN_CMD = used from command line or mrbayes block; IN_FILE = used in data block or in tree block)
+			9 = Where should the command be used (cmdUse) (IN_CMD = used from command line or mrbayes block;
+			    IN_FILE = used in data block or in tree block)
 			 	   10 = Should the command be shown when "help" is typed (hiding).
 			 
 			  #1                 #2   #3                 #4  #5                                                                                                #6        #7                                                            #8       #9   #10
@@ -259,10 +261,10 @@
 			{ 25,            "Lset",  NO,            DoLset, 14,                                                     {28,29,30,31,32,33,34,40,51,52,53,91,92,131},        4,                "Sets the parameters of the likelihood model",  IN_CMD, SHOW },
 			{ 26,	       "Manual",  NO,          DoManual,  1,								    														{127},       36,				  "Prints a command reference to a text file",  IN_CMD, SHOW },
 			{ 27,          "Matrix", YES,          DoMatrix,  1,                                                                                             {11},112381728,                 "Defines matrix of characters in data block", IN_FILE, SHOW },
-			{ 28,            "Mcmc",  NO,            DoMcmc, 37,  {17,18,19,20,21,22,23,24,25,26,27,85,99,113,114,115,116,117,132,143,144,145,149,150,151,152,153,
-																														  154,155,156,157,158,159,160,161,167,170},    36,                   "Starts Markov chain Monte Carlo analysis",  IN_CMD, SHOW },
-			{ 29,           "Mcmcp",  NO,           DoMcmcp, 37,  {17,18,19,20,21,22,23,24,25,26,27,85,99,113,114,115,116,117,132,143,144,145,149,150,151,152,153,
-																														  154,155,156,157,158,159,160,161,167,170},        4, "Sets the parameters of a chain (without starting analysis)",  IN_CMD, SHOW },
+			{ 28,            "Mcmc",  NO,            DoMcmc, 39,  {17,18,19,20,21,22,23,24,25,26,27,85,99,113,114,115,116,117,132,143,144,145,149,150,151,152,153,
+																														  154,155,156,157,158,159,160,161,167,170,172,173},    36,                   "Starts Markov chain Monte Carlo analysis",  IN_CMD, SHOW },
+			{ 29,           "Mcmcp",  NO,           DoMcmcp, 39,  {17,18,19,20,21,22,23,24,25,26,27,85,99,113,114,115,116,117,132,143,144,145,149,150,151,152,153,
+																														  154,155,156,157,158,159,160,161,167,170,172,173},        4, "Sets the parameters of a chain (without starting analysis)",  IN_CMD, SHOW },
 			{ 30,        "Outgroup", YES,        DoOutgroup,  1,                                                                                             {79},    49152,                                     "Changes outgroup taxon",  IN_CMD, SHOW },
 			{ 31,           "Pairs", YES,           DoPairs,  1,                                                                                             {93},    32768,        "Defines nucleotide pairs (doublets) for stem models",  IN_CMD, SHOW },
 			{ 32,       "Partition",  NO,       DoPartition,  1,                                                                                             {16},        4,                              "Assigns a character partition",  IN_CMD, SHOW },
@@ -291,14 +293,16 @@
 		/* NOTE: If you add a command here, make certain to change NUMCOMMANDS (above, in this file) appropriately! */
 		    { 999,             NULL,  NO,              NULL,  0,                                                                                             {-1},       32,                                                           "",  IN_CMD, HIDE }  
 		};
-int					inDataBlock, inForeignBlock, isInterleaved, isFirstMatrixRead, isFirstInterleavedBlock, 
-					taxonCount, fromI, toJ, everyK, foundDash, foundSlash, foundFirst, isMixed, whichPartition,
-					isNegative, numPartitions, charOrdering, foundExp, foundColon, isFirstNode, nextAvailableNode,
-					pairId, firstPair;
-char				gapId, missingId, matchId, tempSetName[100];
-CmdType 			*commandPtr;
-ParmInfoPtr			paramPtr;
-TreeNode			*pPtr, *qPtr;
+int                             isMixed, numPartitions; /* also in bayes.c */
+static int			inDataBlock, inForeignBlock, isInterleaved, isFirstMatrixRead, isFirstInterleavedBlock;
+static int                      taxonCount, fromI, toJ, everyK, foundDash, foundSlash, foundFirst, whichPartition;
+static int                      isNegative, charOrdering, foundExp, foundColon, isFirstNode, nextAvailableNode;
+static int                      pairId, firstPair;
+char				gapId, missingId, matchId; /* also in bayes.c - initialized to 'z' */
+static char                     tempSetName[100];
+static CmdType 			*commandPtr;
+static ParmInfoPtr		paramPtr;
+static TreeNode			*pPtr, *qPtr;
 
 
 
@@ -604,16 +608,16 @@
 
 	if (memAllocs[ALLOC_CONSTRAINTNAMES] == YES)
 		goto errorExit;
-	constraintNames = (char *)SafeMalloc((size_t) (30 * 100 * sizeof(char)));
+	constraintNames = (char *)SafeMalloc((size_t) (MAX_NUM_CONSTRAINTS * 100 * sizeof(char)));
 	if (!constraintNames)
 		{
-		MrBayesPrint ("%s   Problem allocating constraintNames (%d)\n", spacer, 3000 * sizeof(char));
+		MrBayesPrint ("%s   Problem allocating constraintNames (%d)\n", spacer, MAX_NUM_CONSTRAINTS * 100 * sizeof(char));
 		goto errorExit;
 		}
-	for (i=0; i<3000; i++)
+	for (i=0; i<(MAX_NUM_CONSTRAINTS*100); i++)
 		{
 		constraintNames[i] = ' ';
-		if (i == (3000 - 1))
+		if (i == ((MAX_NUM_CONSTRAINTS*100) - 1))
 			constraintNames[i] = '\0';
 		}
 	memAllocs[ALLOC_CONSTRAINTNAMES] = YES;
@@ -2298,9 +2302,10 @@
 					return (ERROR);
 					}
 				}
-			else if (numDefinedConstraints > 30)
+			else if (numDefinedConstraints > MAX_NUM_CONSTRAINTS)
 				{
-				MrBayesPrint ("%s   You cannot define more than 30 constraints\n", spacer);
+				MrBayesPrint ("%s   You cannot define more than %d constraints\n", spacer,
+					      MAX_NUM_CONSTRAINTS);
 				return (ERROR);
 				}
 				
@@ -4898,9 +4903,8 @@
 		hasMissingAmbig = NO;
 		for (j=0; j<numTaxa; j++)
 			{
-			if (IsMissing (matrix[pos(j,i,numChar)], charInfo[i].charType) == YES)
-				hasMissingAmbig = YES;
-			if (IsAmbig (matrix[pos(j,i,numChar)], charInfo[i].charType) == YES)
+			if ((IsMissing (matrix[pos(j,i,numChar)], charInfo[i].charType) == YES) ||
+			    (IsAmbig (matrix[pos(j,i,numChar)], charInfo[i].charType) == YES))
 				hasMissingAmbig = YES;
 			}
 		if (hasMissingAmbig == YES)
@@ -5148,7 +5152,7 @@
 					}
 				else
 					{
-					if ((tkn[i] == ')' && isInAmbig == YES) || (tkn[i] == '}' && isInPoly == YES))
+					if ((tkn[i] == ')' && isInPoly == YES) || (tkn[i] == '}' && isInAmbig == YES))
 						{
 						isInAmbig = isInPoly = NO;
 						charCode = theAmbigChar;
@@ -5158,9 +5162,9 @@
 						matrix[pos(taxonCount-1,taxaInfo[taxonCount-1].charCount++,numChar)] = charCode;
 						theAmbigChar = 0;
 						}
-					else if ((tkn[i] == '(' && isInAmbig == YES) || (tkn[i] == '{' && isInPoly == YES))
+					else if ((tkn[i] == '(' && isInPoly == YES) || (tkn[i] == '{' && isInAmbig == YES))
 						{
-						if (isInAmbig == YES)
+						if (isInPoly == YES)
 							MrBayesPrint ("%s   Found an inappropriate \"(\"\n", spacer);
 						else
 							MrBayesPrint ("%s   Found an inappropriate \"{\"\n", spacer);
@@ -5182,18 +5186,13 @@
 						}
 					else if (tkn[i] == '{' && isInPoly == NO && isInAmbig == NO)
 						{
-						isInPoly = YES;
-						matrixHasPoly = YES;
-						theAmbigChar = 0;
-						}	
-					else if (tkn[i] == '(' && isInPoly == NO && isInAmbig == NO)
-						{
 						isInAmbig = YES;
 						theAmbigChar = 0;
 						}
 					else if (tkn[i] == '(' && isInPoly == NO && isInAmbig == NO)
 						{
-						isInAmbig = YES;
+						isInPoly = YES;
+						matrixHasPoly = YES;
 						theAmbigChar = 0;
 						}
 					else
@@ -9194,7 +9193,7 @@
 
 			MrBayesPrint ("   Shapepr          Uniform/Exponential/Fixed    %s", modelParams[i].shapePr);
 			if (!strcmp(modelParams[i].shapePr, "Uniform"))
-				MrBayesPrint ("(%1.1lf,%1.1lf)\n", modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
+				MrBayesPrint ("(%1.3lf,%1.1lf)\n", modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
 			else if (!strcmp(modelParams[i].shapePr, "Exponential"))
 				MrBayesPrint ("(%1.1lf)\n", modelParams[i].shapeExp);
 			else
@@ -9212,9 +9211,12 @@
 			else
 				MrBayesPrint ("(%1.1lf)\n", modelParams[i].pInvarFix);
 
-			MrBayesPrint ("   Covswitchpr      Uniform/Exponential/Fixed    %s", modelParams[i].covSwitchPr);
+			MrBayesPrint ("   Covswitchpr      Uniform/PinvarUniform/Exponential/Fixed %s", modelParams[i].covSwitchPr);
 			if (!strcmp(modelParams[i].covSwitchPr, "Uniform"))
 				MrBayesPrint ("(%1.1lf,%1.1lf)\n", modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1]);
+                        else if (!strcmp(modelParams[i].covSwitchPr, "PinvarUniform"))
+				MrBayesPrint ("(%1.1lf,%1.1lf,%1.2f)\n",
+					      modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1], modelParams[i].covswitchMaxPinvar);
 			else if (!strcmp(modelParams[i].covSwitchPr, "Exponential"))
 				MrBayesPrint ("(%1.1lf)\n", modelParams[i].covswitchExp);
 			else
@@ -9223,7 +9225,7 @@
 			MrBayesPrint ("   Symdirihyperpr   Uniform/Exponential/Fixed    %s", modelParams[i].symPiPr);
 			if (!strcmp(modelParams[i].symPiPr, "Uniform"))
 				MrBayesPrint ("(%1.1lf,%1.1lf)\n", modelParams[i].symBetaUni[0], modelParams[i].symBetaUni[1]);
-			else if (!strcmp(modelParams[i].covSwitchPr, "Exponential"))
+			else if (!strcmp(modelParams[i].symPiPr, "Exponential"))
 				MrBayesPrint ("(%1.1lf)\n", modelParams[i].symBetaExp);
 			else
 				{
@@ -11682,7 +11684,7 @@
 	PARAM   ( 41, "Shapepr",        DoPrsetParm,       "Uniform|Exponential|Fixed|\0");
 	PARAM   ( 42, "Ratecorrpr",     DoPrsetParm,       "Uniform|Fixed|\0");
 	PARAM   ( 43, "Pinvarpr",       DoPrsetParm,       "Uniform|Fixed|\0");
-	PARAM   ( 44, "Covswitchpr",    DoPrsetParm,       "Uniform|Exponential|Fixed|\0");
+	PARAM   ( 44, "Covswitchpr",    DoPrsetParm,       "PinvarUniform|Uniform|Exponential|Fixed|\0");
 	PARAM   ( 45, "Xxxxxxxxxx",     DoExcludeParm,     "\0");
 	PARAM   ( 46, "Xxxxxxxxxx",     DoIncludeParm,     "\0");
 	PARAM   ( 47, "Xxxxxxxxxx",     DoDeleteParm,      "\0");
@@ -11810,6 +11812,8 @@
 	PARAM   (169, "Aarevmatpr",     DoPrsetParm,       "Dirichlet|Fixed|\0");
 	PARAM   (170, "Nswaps",         DoMcmcParm,       "\0");
 	PARAM   (171, "Autooverwrite",  DoSetParm,         "Yes|No|\0");
+	PARAM   (172, "SimAnneal",             DoMcmcParm,       "Yes|No|\0");
+	PARAM   (173, "Adapt",             DoMcmcParm,       "Yes|No|\0");
 	/* NOTE: If a change is made to the parameter table, make certain you
 	         change the number of elements (now 200) in paramTable[] (global.h: may not be necessary 
 	         and at the top of this file). */
diff -NaurdbB -xCVS -xconfigure$ -xconfig.h.in -xMakefile.in mrbayes-3.1.2.old/Makefile mrbayes-3.1.2.new/Makefile
--- mrbayes-3.1.2.old/Makefile	Sat Mar 25 19:30:51 EST 2006
+++ mrbayes-3.1.2.new/Makefile	Thu Aug 03 20:25:54 EDT 2006
@@ -15,21 +15,28 @@
 
 DEBUG ?= no
 #OPTFLAGS = -O2 -march=pentium4 -mfpmath=sse -fomit-frame-pointers
-OPTFLAGS ?= -O3
+
+# below assumes CFLAGS already set in environment or set below - Allen
+OPTFLAGS ?= -TARG:platform=ip30:processor=r12000
+
+# resetting CFLAGS for debugging - Allen
+#CFLAGS = -g2 -show -TENV:X=1 -LANG:pch=ON:pch_dir=/tmp/allens -DEBUG:subscript_check=ON:verbose_runtime=ON
 
 # set compiler for the non-MPI version (mpicc will be used for the MPI
 # version regardless of this setting)
-CC = gcc
+CC = cc -mips4 -n32 -L/usr/share/rutgers/lib32 -L/usr/lib32/mips4 -L/usr/lib32 -L/lib32 -Wl,-rpath,/usr/share/rutgers/lib32:/usr/lib32/mips4:/usr/lib32:/lib32 -Wl,-multigot -woff 1110 -r12000
 
 # set to yes if you want to use the readline library (make sure you have it
 # installed on your machine)
 # defaults to yes for unix machines
-ifeq ($(strip $(ARCHITECTURE)), unix)
-   USEREADLINE ?= yes
-else
-   USEREADLINE ?= no
-endif
+#ifeq ($(strip $(ARCHITECTURE)), unix)
+#   USEREADLINE ?= yes
+#else
+#   USEREADLINE ?= no
+#endif
 
+USEREADLINE = no
+
 #
 # End of user configuration
 
@@ -63,7 +70,9 @@
 	CFLAGS += $(OPTFLAGS)
 endif
 
-CFLAGS   +=   -Wall 
+# GCC specific, guys! -Allen
+
+#CFLAGS   +=   -Wall 
 
 LIBS     +=   -lm
 
diff -NaurdbB -xCVS -xconfigure$ -xconfig.h.in -xMakefile.in mrbayes-3.1.2.old/mb.h mrbayes-3.1.2.new/mb.h
--- mrbayes-3.1.2.old/mb.h	Sun Mar 26 11:18:31 EST 2006
+++ mrbayes-3.1.2.new/mb.h	Thu Nov 22 12:46:54 EST 2007
@@ -10,8 +10,9 @@
    Makefile file; for compilers not using a project file or Makefile, use the #defines above to select version */
 
 typedef	double MrBFlt;		/* double used for parameter values and generally for floating point values */
-typedef float CLFlt;		/* single-precision float used for cond likes (CLFlt) to increase speed and reduce memory requirement */
+typedef double CLFlt;		/* single-precision float used for cond likes (CLFlt) to increase speed and reduce memory requirement */
 							/* set CLFlt to double if you want increased precision */
+typedef float CLFlt2; /* for things for which increased precision matters less - Allen */
 
 #undef FAST_LOG			/* define this to get a fast polynomial approximation to the log for 4by4 nucleotide data */
 						/* this results in speed gains of 10% to 20% for typical data sets */
@@ -27,13 +28,13 @@
 #undef ICC
 */
 
-#if defined GCC			/* gcc compiler */
-#define ALIGNEDMALLOC memalign
-#elif defined ICC		/* icc compiler */
+#if defined ICC			/* icc compiler */
 #define ALIGNEDMALLOC _mm_malloc
 #elif defined MS_VCPP   /* Visual .Net */
 #define ALIGNEDMALLOC _aligned_malloc
 #include <xmmintrin.h>
+#else  /* Added by Allen Smith */
+#define ALIGNEDMALLOC memalign
 #endif
 
 #ifdef PAUL /* for debugging */
@@ -305,14 +306,18 @@
 #define	P_BRSIGMA				16
 #define	P_GROWTH				17 /* NOTE: If you add another parameter, change NUM_LINKED */
 
-#define	MAX_NUM_DIVS			150
-#define	MAX_NUM_CHARSETS		150
+#define	MAX_NUM_DIVS			256 /* Allen - changed from 150 */
+#define	MAX_NUM_CHARSETS		400 /* Allen - changed from 150 */
 #define	MAX_NUM_PARTITIONS		30
 #define	MAX_NUM_TREES			200
-#define	MAX_NUM_DIV_BITS		150
+#define	MAX_NUM_DIV_BITS		256 /* Allen - changed from 150 */
 #define MAX_NUM_DIV_LONGS       1 + MAX_NUM_DIV_BITS / 32
 #define	MAX_NUM_TAXASETS		30
+#ifdef NOTOPOLOGY /* - Allen */
 #define	MAX_NUM_CONSTRAINTS		30
+#else
+#define	MAX_NUM_CONSTRAINTS		128 /* Allen - changed from 30 */
+#endif
 #define	MAX_CHAINS				256
 
 typedef void * VoidPtr;
@@ -543,6 +548,8 @@
 #define	OMEGA_10FEF						106
 #define	OMEGA_10FFB						107
 #define	OMEGA_10FFF						108
+/* Pinvar limits on Switchrates - Allen */
+#define SWITCH_PINVUNI						109
 
 /* typedef for a MoveFxn */
 typedef int (MoveFxn)(Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
@@ -562,7 +569,12 @@
 	MrBFlt		minimum[2];         /* minimum values for tuning params             */
 	MrBFlt		maximum[2];         /* maximum values for tuning params             */
 	int         parsimonyBased;     /* this move is based on parsimony (YES/NO)     */
+	int         doAdaptParam;       /* can use Corana et al sliding window adaptation for this move (YES/NO) */
+	int         whichAdaptParam;  /* which param to adapt - Allen */
+	MrBFlt          minAdaptCFactor;           /* min C-factor to use in adaptation - Allen */
+	int             adaptMultiplier;        /* whether the param to be adapted is a multiplier (YES/NO) - Allen */
 	int         level;              /* user level of this move                      */
+	int         covarion_wanted;    /* 1 if need covarion, -1 if no covarion allowed, 0 if doesn't care - Allen */
 	} MoveType;
 
 
@@ -574,8 +586,7 @@
 /* Note: This allows several different moves to affect the same parameter */
 /* It also allows the same move to affect different parameters as before */
 /* This is also a good place to keep the proposal probs */
-typedef struct
-	{
+typedef struct {
 	char		*name;				/* pointer to the name of the move type         */
 	char		*shortName;	        /* pointer to the short name of the move        */
 	MoveFxn		*moveFxn;			/* pointer to the move function					*/
@@ -584,12 +595,30 @@
 	MrBFlt		cumProposalProb;	/* the cumulative proposal probability			*/
 	int			*nAccepted;			/* number of accepted moves						*/
 	int			*nTried;			/* number of tried moves						*/
+  int             doAdaptParam;           /* can use Corana et al sliding window adaptation for this move (YES/NO) */
+  int             whichAdaptParam;        /* which param to adapt - Allen */
+  int             adaptMultiplier;        /* whether the param to be adapted is a multiplier (YES/NO) - Allen */
+  int             nTriedAfterAdapt;       /* number of times tried after last adapt, if any */
+  int             nAcceptedAfterAdapt;    /* ditto, for accepted */
+  MrBFlt          minAdapt;               /* minimum of adapt param */
+  MrBFlt          minAdaptMult;           /* minimum for multiplier of adapt param before do 2*log(x) */
+  MrBFlt          maxAdapt;               /* maximum of adapt param */
+  MrBFlt          maxAdaptMult;           /* maximum for multiplier of adapt param before do 2*log(x) */
+  MrBFlt          minAdaptCFactor;           /* min C-factor to use in adaptation - Allen */
+  MrBFlt          oldTuningParam;         /* value of adapt param before any adaptation */
+  MrBFlt          oldTuningParamMult;     /* value of adapt param before any adaptation before do 2*log(x) */
+  MrBFlt          minParamTried;           /* min value of adapt param tried */
+  MrBFlt          maxParamTried;           /* max value of adapt param tried */
+  MrBFlt          minParamRefused;        /* min value of adapt param not accepted */
+  MrBFlt          maxParamRefused;        /* max value of adapt param not accepted */
+  MrBFlt          minParamAccepted;       /* min value of adapt param accepted */
+  MrBFlt          maxParamAccepted;       /* max value of adapt param accepted */
 	MrBFlt		tuningParam[2];     /* tuning parameters for the move               */
-	} MCMCMove;
+} MCMCMove;
 
 typedef int (*LikeDownFxn)(TreeNode *, int, int);
 typedef int (*LikeRootFxn)(TreeNode *, int, int);
-typedef int (*LikeScalerFxn)(TreeNode *, int, int);
+typedef int (*LikeScalerFxn)(TreeNode *, int, int, int);
 typedef int (*LikeFxn)(TreeNode *, int, int, MrBFlt *, int);
 typedef int (*TiProbFxn)(TreeNode *, int, int);
 typedef int (*LikeUpFxn)(TreeNode *, int, int);
@@ -698,6 +727,7 @@
 	MrBFlt		covswitchFix[2];
 	MrBFlt		covswitchUni[2];
 	MrBFlt		covswitchExp;
+	MrBFlt		covswitchMaxPinvar; /* maximum ProbOff for covarion - Allen        */
 	char		symPiPr[100];      /* prior for pi when unidentifiable states used */
 	MrBFlt		symBetaFix;
 	MrBFlt		symBetaUni[2];
@@ -714,7 +744,7 @@
 	MrBFlt		brownScalesGammaMean;
 
 	char		topologyPr[100];   /* prior for tree topology                      */
-	int			activeConstraints[30]; 
+	int			activeConstraints[MAX_NUM_CONSTRAINTS]; 
 	int			numActiveConstraints;
 	int			numActiveLocks;
 	char		brlensPr[100];     /* prior on branch lengths                      */
@@ -787,6 +817,8 @@
 	int         diagnFreq;             /* mcmc diagnostics frequency                    */
 	int         relativeBurnin;        /* should a relative burnin be used ?            */
 	MrBFlt      burninFraction;        /* the sample fraction to discard as burnin      */
+	int         useSA; /* flag for whether to use simulated annealing - Allen */
+	int         useAdapt; /* flag for whether to use Corana et al-type adaptation - Allen */
 	int         allComps;              /* top conv diagnosis for all pairs?             */
 	MrBFlt      minPartFreq;           /* minimum partition frequency for conv diagn    */
 	MrBFlt      stopVal;               /* top conv diagn value to reach before stopping */
diff -NaurdbB -xCVS -xconfigure$ -xconfig.h.in -xMakefile.in mrbayes-3.1.2.old/mcmc.c mrbayes-3.1.2.new/mcmc.c
--- mrbayes-3.1.2.old/mcmc.c	Wed Apr 18 12:19:46 EDT 2007
+++ mrbayes-3.1.2.new/mcmc.c	Thu Nov 22 12:52:14 EST 2007
@@ -29,7 +29,18 @@
  * GNU General Public License for more details (www.gnu.org).
  * 
  */
+#ifdef __sgi
+# ifdef _XOPEN_SOURCE
+#  if ((_XOPEN_SOURCE+0) < 520)
+#    define _XOPEN_SOURCE 520
+#  endif
+# else
+#  define _XOPEN_SOURCE 520
+# endif
+#endif
 #include <stdio.h>
+/* Above added by Allen to try to make sure vsnprintf (for SaveSprintf) returns needed length
+   - may need to be limited in what applies to (needed for SGI IRIX at minimum) */
 #include <stdlib.h>
 #include <time.h>
 #include <math.h>
@@ -52,6 +63,16 @@
 #include <assert.h>
 #endif
 
+#if (defined(__GNUC__) || defined(_ISOC99_SOURCE) || defined(__sgi) || defined(_SVID3) || defined(_XOPEN4UX) || (defined(_XOPEN_SOURCE) && defined(_XOPEN_SOURCE_EXTENDED) && (_XOPEN_SOURCE_EXTENDED >= 1)))
+#define HAVE_CBRT
+#endif
+
+#if defined(HAVE_CBRT)
+#define CBRT(x) cbrt((double)(x))
+#else
+#define CBRT(x) pow((double)(x),(double)(1.0/3.0))
+#endif
+
 #if defined(WIN_VERSION) && !defined(__GNUC__)
 #define VISUAL
 #endif
@@ -112,6 +133,7 @@
 
 /* debugging compiler statements */
 #undef	DEBUG_COMPRESSDATA
+#define DEBUG_COMPRESSDATA2 /* - Allen */
 #undef	DEBUG_ADDDUMMYCHARS
 #undef	DEBUG_CREATEPARSMATRIX
 #undef	DEBUG_SETUPTERMSTATE
@@ -118,11 +140,13 @@
 #undef	DEBUG_INITCHAINCONDLIKES
 #undef	DEBUG_SETCHAINPARAMS
 #undef	DEBUG_BUILDSTARTTREE
+#undef  DEBUG_DEROOT_USER_TREE /* -Allen */
 #undef	DEBUG_RUNCHAIN
 #undef	DEBUG_NOSHORTCUTS
 #undef	DEBUG_NOSCALING
 #undef	DEBUG_LOCAL
 #undef	DEBUG_SPRCLOCK
+#undef  DEBUG_SCALER_NAN /* - Allen */
 #undef	DEBUG_TIPROBS_STD
 #undef	DEBUG_RUN_WITHOUT_DATA
 #undef	DEBUG_UNROOTED_SLIDER
@@ -132,6 +156,12 @@
 #undef  DEBUG_CSLIDER
 #undef  DEBUG_ExtSPRClock
 #undef  DEBUG_TREEHEIGHT
+#undef  DEBUG_MPI /* - Allen */
+
+
+#if !(defined(DEBUG_NOSCALING) || defined(NO_ALTSCALING))
+#define ALTSCALING /* do scaling only in response to LIKE_EPSILON - Allen */
+#endif
 
 /* local (to this file) data types */
 typedef struct pfnode
@@ -188,14 +218,14 @@
 #		endif
 int	    CondLikeRoot_NY98 (TreeNode *p, int division, int chain);
 int     CondLikeRoot_Std (TreeNode *p, int division, int chain);
-int	    CondLikeScaler_Gen (TreeNode *p, int division, int chain);
+int	    CondLikeScaler_Gen (TreeNode *p, int division, int chain, int do_full_scaling);
 #		if !defined (FAST_LOG)
-int	    CondLikeScaler_NUC4 (TreeNode *p, int division, int chain);
+int	    CondLikeScaler_NUC4 (TreeNode *p, int division, int chain, int do_full_scaling);
 #		else
-int	    CondLikeScaler_NUC4_fast (TreeNode *p, int division, int chain);
+int	    CondLikeScaler_NUC4_fast (TreeNode *p, int division, int chain, int do_full_scaling);
 #		endif
-int	    CondLikeScaler_NY98 (TreeNode *p, int division, int chain);
-int     CondLikeScaler_Std (TreeNode *p, int division, int chain);
+int	    CondLikeScaler_NY98 (TreeNode *p, int division, int chain, int do_full_scaling);
+int     CondLikeScaler_Std (TreeNode *p, int division, int chain, int do_full_scaling);
 int     CondLikeUp_Bin (TreeNode *p, int division, int chain);
 int     CondLikeUp_Gen (TreeNode *p, int division, int chain);
 int     CondLikeUp_NUC4 (TreeNode *p, int division, int chain);
@@ -303,6 +333,7 @@
 int     Move_NNI (Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
 int		Move_NNI_Hetero (Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
 int		Move_NodeSlider (Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
+int		Move_NodeSlider2 (Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
 int		Move_NodeSliderClock (Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
 int     Move_Omega (Param *param, int chain, long *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
 int     Move_Omega_M (Param *param, int chain, long int *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
@@ -356,7 +387,7 @@
 int     PrintTermState (void);
 void	PrintTiProbs (MrBFlt *tP, MrBFlt *bs, int nStates);
 int     PrintTopConvInfo (void);
-void    PrintToScreen (int curGen, time_t endingT, time_t startingT);
+void    PrintToScreen (int curGen, int movesAccepted, time_t endingT, time_t startingT);
 /* void    PrintToScreen (int curGen, clock_t endingT, clock_t startingT);*/
 int     PrintTree (int curGen, Tree *tree);
 int     ProcessStdChars (void);
@@ -372,7 +403,9 @@
 int     RemoveTreeSamples (int from, int to);
 int     ReopenMBPrintFiles (void);
 int     RequestAbortRun(void);
+#if !defined(ALTSCALING)
 int     ResetScalers (void);
+#endif
 int     RunChain (long int *seed);
 int     SaveSprintf(char **target, int *targetLen, char *fmt, ...);
 int     SetAARates (void);
@@ -457,7 +490,7 @@
 long int		*compMatrix;		         /* compressed character matrix					 */
 long int		*parsMatrix = NULL;		         /* parsimony (bitset) matrix for terminals  	 */
 long int		*parsSets = NULL;					 /* parsimony (bitset) matrix for int nodes 	 */
-CLFlt			*numSitesOfPat;		         /* no. sites of each pattern					 */
+CLFlt2			*numSitesOfPat;		         /* no. sites of each pattern					 */
 CLFlt			*termCondLikes = NULL;		         /* cond likes for terminals                     */
 CLFlt			*chainCondLikes;	         /* cond likes for chains						 */
 int				condLikeRowSize;	         /* row size of cond like matrices				 */
@@ -488,8 +521,8 @@
 CLFlt			**treeScaler;		         /* pointers to tree scalers for each chain		 */
 CLFlt			*nodeScalerSpace;	         /* space holding cond like node scalers         */
 CLFlt			**nodeScaler;		         /* pointers to cond like scalers for each chain */
-CLFlt			*tiProbSpace;		         /* space holding tiProbs						 */
-CLFlt			**tiProbs;			         /* pointers to tiProbs for each chain			 */
+CLFlt2			*tiProbSpace;		         /* space holding tiProbs						 */
+CLFlt2			**tiProbs;			         /* pointers to tiProbs for each chain			 */
 int				cijkRowSize;		         /* row size of cijk information                 */
 MrBFlt			*cijkSpace;			         /* space holding cijk information               */
 MrBFlt			**cijks;			         /* pointers to cijk for each chain              */
@@ -548,6 +581,45 @@
 int				highestLocalRunId;			 /* highest local run Id                         */
 #endif
 
+/* debugging scaler problems - Allen */
+
+#if defined(ALTSCALING)
+int had_like_epsilon = 0;
+int did_scaling = 0;
+#endif
+
+#define SCALER_WARN_THRESH 40.0
+#define SCALER_WARN_THRESH2 39.0
+ /* 16 =~ -log(1/2**23) */
+#define SCALER_WARN_THRESH3 16.0
+#define CL_WARN_THRESH 1.000001
+#define CL_WARN_THRESH2 2.0
+#define CL_WARN_THRESH3 20.0
+
+#define LOG_CL_WARN_THRESH 0.000001
+
+#define MAX_TIMES_CL_WARN 200
+int num_times_cl_warn = 0; /* count of num times warned */
+
+#define MAX_TIMES_SCALER_WARN 200
+int num_times_scaler_warn = 0;
+
+#define MAX_TIMES_SCALER_WARN2 200
+int num_times_scaler_warn2 = 0;
+
+#define MAX_TIMES_SCALER_WARN3 20
+int num_times_scaler_warn3 = 0; /* this is reset w/each likelihood print */
+
+#define MAX_TIMES_SCALER_WARN4 200
+int num_times_scaler_warn4 = 0;
+
+int num_times_scaler_adjusted = 0;
+/* only incremented if lnScaler would be above LIKE_EPSILON, scaler above CL_WARN_THRESH */
+int num_times_scaler_adjusted_big = 0;
+
+int num_times_scaler_adjusted2 = 0;
+/* only incremented if lnScaler would be above log(CL_WARN_THRESH) */
+int num_times_scaler_adjusted2_big = 0;
 
 /*-----------------------------------------------------------------------
 |
@@ -562,7 +634,7 @@
 				newRowSize, numDummyChars, newColumn, newChar, oldColumn, oldChar, 
 				isCompat, *tempChar, numIncompatibleChars;
 	long		*tempMatrix;
-	CLFlt		*tempSitesOfPat;
+	CLFlt2		*tempSitesOfPat;
 	ModelInfo	*m;
 	ModelParams	*mp;
 	CharInfo	cinfo;
@@ -627,7 +699,7 @@
 
 	/* allocate space for new data */
 	tempMatrix = (long *) calloc (numLocalTaxa * newRowSize, sizeof(long));
-	tempSitesOfPat = (CLFlt *) calloc (numCompressedChars, sizeof(CLFlt));
+	tempSitesOfPat = (CLFlt2 *) calloc (numCompressedChars, sizeof(CLFlt2));
 	tempChar = (int *) calloc (compMatrixRowSize, sizeof(int));
 	if (!tempMatrix || !tempSitesOfPat || !tempChar)
 		{
@@ -1138,7 +1210,10 @@
 	/* the number will be corrected further down for unbalanced scenarios */
 	numChainsForProc = (int) (chainParams.numChains * chainParams.numRuns / num_procs);
 
+#       if defined(DEBUG_MPI)
+	MrBayesPrint("%s   Attempting swap of %d and %d\n",spacer,swapA,swapB);
 #	endif
+#	endif
 
 	/* are we using character reweighting? */
 	reweightingChars = NO;
@@ -1365,6 +1440,9 @@
 					{
 					return (ERROR);
 					}
+#       if defined(DEBUG_MPI)
+				MrBayesPrint ("%s   Waiting on swap communication (A)\n",spacer);
+#       endif
 				ierror = MPI_Waitall (2, request, status);
 				if (ierror != MPI_SUCCESS)
 					{
@@ -1456,6 +1534,10 @@
 					{
 					return (ERROR);
 					}
+#       if defined(DEBUG_MPI)
+
+				MrBayesPrint ("%s   Waiting on swap communication (B)\n",spacer);
+#       endif
 				ierror = MPI_Waitall (2, request, status);
 				if (ierror != MPI_SUCCESS)
 					{
@@ -1904,7 +1986,7 @@
 
 	int				i, j, whichNde, n, *tempNums, numAvailTips, intNodeNum, whichElement, 
 					nLongsNeeded;
-	MrBFlt			sum, oldBls[5], newBls[5], tempLength, depth;
+	MrBFlt			sum=0.0, oldBls[5], newBls[5], tempLength, depth;
 	TreeNode		*p, *q, *up, *dn, *a, *b, *c;
 	PolyTree		constraintTree;
 	Tree			*tempTree;
@@ -2337,9 +2419,20 @@
 		CopyToTreeFromTree (tempTree, userTree);
 
 		/* If unrooted tree is needed and usertree is rooted, we unroot it */
-		if (t->isRooted == NO && tempTree->isRooted == YES)
-			{
+		if (t->isRooted == NO && tempTree->isRooted == YES) {
+#if (defined(DEBUG_BUILDSTARTTREE) || defined(DEBUG_DEROOT_USER_TREE)) /* Allen */
+		  if (ShowMCMCTree(tempTree) == ERROR) {
+		    MrBayesPrint("%s   Problem showing rooted tempTree\n",
+				 spacer);
+		  }
+#endif
 			DerootTree (tempTree, outGroupNum);
+#if (defined(DEBUG_BUILDSTARTTREE) || defined(DEBUG_DEROOT_USER_TREE)) /* Allen */
+		  if (ShowMCMCTree(tempTree) == ERROR) {
+		    MrBayesPrint("%s   Problem showing derooted tempTree\n",
+				 spacer);
+		  }
+#endif
 			}
 
 		/* Make sure calculation root is in right place if tree is unrooted */
@@ -2349,9 +2442,12 @@
 			}
 
 		/* Prune tree if some taxa are deleted from this analysis */
-		if (PruneTree (tempTree) == ERROR)
-			{
+		if (PruneTree (tempTree) == ERROR) {
 			MrBayesPrint ("%s   Error in pruning user tree\n", spacer);
+		  if (ShowMCMCTree(tempTree) == ERROR) { /* Allen */
+		    MrBayesPrint("%s   Problem showing tempTree w/error\n",
+				 spacer);
+		  }
 			goto errorExit;
 			}
 			
@@ -2403,12 +2499,14 @@
 				goto errorExit;
 				}
 			for (i=0; i<chainParams.numStartPerts; i++)
-				{
+				{ /* removed _&&_ p->anc->anc == NULL from below -
+				     was not obeying constraints with perturbations - Allen */
 				do
 					{
 					whichElement = (int)(RandomNumber(seed) * t->nIntNodes);
 					p = t->intDownPass[whichElement];
-					} while (p->anc->anc == NULL && p->isLocked == YES);
+					} while ((p->isLocked == YES) || (p->anc == NULL) ||
+						 (p->anc->anc == NULL));
 				q = p->anc;
 				a  = p->left;
 				b  = p->right;
@@ -2420,6 +2518,7 @@
 					{
 					c  = q->left;
 					}
+				if ((t->isRooted != NO) || (userBrlensDef == YES)) /* - Allen */
 				sum = a->length + b->length + c->length + p->length + q->length;
 				if (RandomNumber(seed) < 0.5)
 					{
@@ -2449,8 +2548,13 @@
 					q->right = a;
 					a->anc  = q;
 					}
-				if (sum > 0.0)
-					{
+				if ((t->isRooted == NO) && (userBrlensDef != YES)) { /* - Allen */
+				  a->length = 0.1;
+				  b->length = 0.1;
+				  c->length = 0.1;
+				  p->length = 0.1;
+				  q->length = 0.1;
+				} else if (sum > (BRLENS_MIN*5)) { /* - Allen */
 					/* We also need to change the branch lengths. However, we don't
 					   change the overall tree length this way. */
 					oldBls[0] = a->length / sum;
@@ -2464,8 +2568,19 @@
 					c->length = newBls[2] * sum;
 					p->length = newBls[3] * sum;
 					q->length = newBls[4] * sum;
+					if (t->isRooted == NO) { /* otherwise, tree length a problem... */
+					  if (a->length < BRLENS_MIN) /* roundoff error, usually - Allen */
+					    a->length = BRLENS_MIN;
+					  if (b->length < BRLENS_MIN)
+					    b->length = BRLENS_MIN;
+					  if (c->length < BRLENS_MIN)
+					    c->length = BRLENS_MIN;
+					  if (p->length < BRLENS_MIN)
+					    p->length = BRLENS_MIN;
+					  if (q->length < BRLENS_MIN)
+					    q->length = BRLENS_MIN;
 					}
-				else
+				} else
 					{
 					a->length = BRLENS_MIN;
 					b->length = BRLENS_MIN;
@@ -2698,7 +2813,7 @@
 		posit = m->rateProbStart + (compCharPos[c] - m->compCharStart) * m->numGammaCats;
 				
 		/* calculate the HMM forward probs fi(x) at site x in HMM */
-		if (siteJump[c] <= MAX_SMALL_JUMP)
+		if ((siteJump[c] > 0) && (siteJump[c] <= MAX_SMALL_JUMP)) /* added siteJump[c] > 0 check - Allen */
 			{
 			max = 0.0;
 			for (i=0; i<nRates; i++)
@@ -2711,7 +2826,7 @@
 					max = F[i];
 				}
 			}
-		else if (siteJump[c] < BIG_JUMP)	/* intermediate jump, calculate trans probs */
+		else if ((siteJump[c] > 0) && (siteJump[c] < BIG_JUMP))	/* intermediate jump, calculate trans probs */
 			{
 			MultiplyMatrixNTimes(nRates, markovTi[0], siteJump[c], markovTiN);
 			max = 0.0;
@@ -3053,7 +3168,7 @@
 
 {
 
-	int				a, b, i, j, k, nLongsNeeded;
+	int				a, b, i, j, k, nLongsNeeded, did_flip;
 	long int		*constraintPartition, *mask, *bitsets;
 	TreeNode		*p=NULL;
 	   	
@@ -3061,7 +3176,7 @@
 	nLongsNeeded = (numLocalTaxa / nBitsInALong) + 1;
 
 	bitsets = (long int *) calloc ((t->nNodes + 1)*nLongsNeeded + /*1*/
-                                       1+numLocalTaxa/nBitsInALong+1, 
+                                       1+numLocalTaxa/nBitsInALong+1, /* exactly what is this trying to add?!? - Allen */
                                        sizeof(long int));
 
 	if (!bitsets)
@@ -3102,11 +3217,12 @@
 			}
 		}
 
-	for (a=b=0; a<numDefinedConstraints; a++)
-		{
+	for (a=b=0; a<numDefinedConstraints; a++) {
 		if (t->constraints[a] == NO)
 			continue;
 
+		did_flip = 0;
+
 		/* set bits in partition to add */
 		for (j=0; j<nLongsNeeded; j++)
 			constraintPartition[j] = 0;
@@ -3121,8 +3237,10 @@
 			}
 
 		/* make sure outgroup is outside constrained partition if unrooted tree */
-		if (t->isRooted == NO && IsBitSet(localOutGroup, constraintPartition))
+		if (t->isRooted == NO && IsBitSet(localOutGroup, constraintPartition)) {
 			FlipBits(constraintPartition, nLongsNeeded, mask);
+			did_flip = 1;
+		}
 
 		/* find the locked node */
 		for (i=j=0; i<t->nNodes; i++)
@@ -3136,17 +3254,18 @@
 	
 		if (j != 1)
 			{
-			MrBayesPrint ("%s   Lock number %d is not set or set several places\n", spacer, b);
+			MrBayesPrint ("%s   Lock number %d is not set or set several places (%d)\n", spacer, b, j);
 			goto errorExit;
 			}
 
 		/* check that locked node is correct */
-		for (i=0; i<nLongsNeeded; i++)
-			{
-			if (p->partition[i] != constraintPartition[i]) 
-				{
-				MrBayesPrint ("%s   Lock number %d is set for the wrong node\n", spacer, b);
+		for (i=0; i<nLongsNeeded; i++) {
+		  if (p->partition[i] != constraintPartition[i]) {
+		    MrBayesPrint ("%s   Lock number %d is set for the wrong node (%ld vs %ld; did_flip %d)\n", spacer, b,
+				  p->partition[i], constraintPartition[i], did_flip); /* - Allen */
+#ifdef DEBUG_CONSTRAINTS
 				goto errorExit;
+#endif
 				}
 			}
 
@@ -3573,6 +3692,7 @@
 				{
 				foundIt = YES;
 				p->lockID = b;
+				p->isLocked = YES; /* Allen */
 				if (t->isCalibrated == YES)
 					{
 					p->isDated = YES;
@@ -3604,9 +3724,24 @@
 int CheckTemperature (void)
 
 {
-
+  if ((chainParams.useSA == YES) && (chainParams.numChains > 1)) { /* - Allen */
+    MrBayesPrint ("%s   Cannot use Simulated Annealing with multiple (%d) chains.\n",
+		  spacer, chainParams.numChains);
+    return (ERROR);
+  }
+  if ((chainParams.useAdapt == YES) && (chainParams.numChains > 1)) { /* - Allen */
+    MrBayesPrint ("%s   Cannot use Adapt with multiple (%d) chains.\n",
+		  spacer, chainParams.numChains);
+    return (ERROR);
+  }
+  if ((chainParams.useAdapt == YES) && (chainParams.numRuns > 1)) { /* - Allen */
+    MrBayesPrint ("%s   Cannot use Adapt with multiple (%d) runs.\n",
+		  spacer, chainParams.numRuns);
+    return (ERROR);
+  }
 	if (chainParams.userDefinedTemps == YES)
 	        {
+
 		  if(AreDoublesEqual(chainParams.userTemps[0], 1.0, ETA)==NO)
 			{
 			MrBayesPrint ("%s   The first user-defined temperature must be 1.0.\n", spacer);
@@ -3737,8 +3872,10 @@
 #	if defined DEBUG_COMPRESSDATA
 	if (PrintMatrix() == ERROR)
 		goto errorExit;
+#       if (! defined(MPI_ENABLED))
 	getchar();
 #	endif
+#	endif
 
 	/* set all pointers that will be allocated locally to NULL */
 	isTaken = NULL;
@@ -3961,6 +4098,9 @@
 
 	compMatrixRowSize = newColumn;
 
+	MrBayesPrint("%s   Have %d unique site patterns total\n", spacer, 
+		     numCompressedChars); /* - Allen */
+
 	/* now we know the size, so we can allocate space for the compressed matrix ... */
 	if (memAllocs[ALLOC_COMPMATRIX] == YES)
 		{
@@ -3980,10 +4120,10 @@
 		MrBayesPrint ("%s   numSitesOfPat not free in CompressData\n", spacer);
 		goto errorExit;
 		}
-	numSitesOfPat = (CLFlt *) calloc (numCompressedChars, sizeof(CLFlt));
+	numSitesOfPat = (CLFlt2 *) calloc (numCompressedChars, sizeof(CLFlt2));
 	if (!numSitesOfPat)
 		{
-		MrBayesPrint ("%s   Problem allocating numSitesOfPat (%d)\n", spacer, numCompressedChars * sizeof(MrBFlt));
+		MrBayesPrint ("%s   Problem allocating numSitesOfPat (%d)\n", spacer, numCompressedChars * sizeof(CLFlt2));
 		goto errorExit;
 		}
 	memAllocs[ALLOC_NUMSITESOFPAT] = YES;
@@ -4007,16 +4147,18 @@
 			compMatrix[pos(i,j,compMatrixRowSize)] = tempMatrix[pos(i,j,numLocalChar)];
 
 	for (i=0; i<numCompressedChars; i++)
-		numSitesOfPat[i] = (CLFlt) tempSitesOfPat[i];
+		numSitesOfPat[i] = (CLFlt2) tempSitesOfPat[i];
 
 	for (i=0; i<compMatrixRowSize; i++)
 		origChar[i] = tempChar[i];
 
-#	if defined (DEBUG_COMPRESSDATA)
+#	if (defined(DEBUG_COMPRESSDATA) || defined(DEBUG_COMPRESSDATA2)) /* - Allen */
 	if (PrintCompMatrix() == ERROR)
 		goto errorExit;
+#       if (! defined(MPI_ENABLED))
 	getchar();
 #	endif
+#	endif
 
 	/* free the temporary variables */
 	free (tempSitesOfPat);
@@ -4055,7 +4197,8 @@
 {
 
 	int				c, k;
-	CLFlt			*clL, *clR, *clP, *pL, *pR, *tiPL, *tiPR;
+	CLFlt			*clL, *clR, *clP;
+	CLFlt2                  *pL, *pR, *tiPL, *tiPR;
 	ModelInfo		*m;
 	
 	/* find model settings for this division */
@@ -4109,7 +4252,8 @@
 
 	int				a, b, c, h, i, j, k, shortCut, *lState=NULL, *rState=NULL, catStart,
 					nObsStates, nStates, nStatesSquared, preLikeJump;
-	CLFlt			likeL, likeR, *pL, *pR, *tiPL, *tiPR, *clL, *clR, *clP;
+	CLFlt			likeL, likeR, *clL, *clR, *clP;
+	CLFlt2                  *pL, *pR, *tiPL, *tiPR;
 	ModelInfo		*m;
 	
 	/* find model settings for this division and nStates, nStatesSquared */
@@ -4291,7 +4435,8 @@
 
 {
 	int				c, h, i, j, k, shortCut, *lState=NULL, *rState=NULL;
-	CLFlt			*clL, *clR, *clP, *pL, *pR, *tiPL, *tiPR;
+	CLFlt			*clL, *clR, *clP;
+	CLFlt2                  *pL, *pR, *tiPL, *tiPR;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -4446,7 +4591,7 @@
 
 
 
-
+/* NOTE: CLFlt, CLFlt2 assumptions made in below! - Allen */
 
 #if defined SSE
 /*----------------------------------------------------------------
@@ -4532,7 +4677,8 @@
 {
 
 	int				a, b, c, h, i, j, k, shortCut, *lState=NULL, *rState=NULL, nStates, nStatesSquared;
-	CLFlt			likeL, likeR, *pL, *pR, *tiPL, *tiPR, *clL, *clR, *clP;
+	CLFlt			likeL, likeR, *clL, *clR, *clP;
+	CLFlt2                  *pL, *pR, *tiPL, *tiPR;
 	ModelInfo		*m;
 	
 	/* find model settings for this division and nStates, nStatesSquared */
@@ -4693,7 +4839,8 @@
 {
 
 	int				a, c, h, i, j, k, nStates, nCats;
-	CLFlt			*clL, *clR, *clP, *pL, *pR, *tiPL, *tiPR, likeL, likeR;
+	CLFlt			*clL, *clR, *clP, likeL, likeR;
+	CLFlt2                  *pL, *pR, *tiPL, *tiPR;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -4761,7 +4908,8 @@
 {
 
 	int				c, k;
-	CLFlt			*clL, *clR, *clP, *clA, *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
+	CLFlt			*clL, *clR, *clP, *clA;
+	CLFlt2                  *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
 	ModelInfo		*m;
 
 	/* find model settings for this division */
@@ -4823,8 +4971,8 @@
 
 	int				a, b, c, h, i, j, k, shortCut, *lState=NULL, *rState=NULL, *aState=NULL,
 					catStart, nObsStates, nStates, nStatesSquared, preLikeJump;
-	CLFlt			likeL, likeR, likeA, *clL, *clR, *clP, *clA, *pL, *pR, *pA,
-					*tiPL, *tiPR, *tiPA;
+	CLFlt			likeL, likeR, likeA, *clL, *clR, *clP, *clA;
+	CLFlt2                  *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
 	ModelInfo		*m;
 	
 	/* find model settings for this division and nStates, nStatesSquared */
@@ -5057,7 +5205,8 @@
 
 {
 	int				c, h, i, j, k, shortCut, *lState=NULL, *rState=NULL, *aState=NULL;
-	CLFlt			*clL, *clR, *clP, *clA, *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
+	CLFlt			*clL, *clR, *clP, *clA;
+	CLFlt2                  *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -5270,7 +5419,7 @@
 
 
 
-
+/* NOTE ASSUMPTIONS REGARDING CLFlt, CLFlt2 -Allen */
 
 #if defined SSE
 /*----------------------------------------------------------------
@@ -5371,10 +5520,11 @@
 
 {
 
-	int				a, b, c, h, i, j, k, shortCut, *lState=NULL, *rState=NULL, *aState=NULL,
-					nStates, nStatesSquared;
-	CLFlt			likeL, likeR, likeA, *clL, *clR, *clP, *clA, *pL, *pR, *pA,
-					*tiPL, *tiPR, *tiPA;
+        int			a, b, c, h, i, j, k, shortCut, nStates, nStatesSquared;
+	int                     *lState=NULL, *rState=NULL, *aState=NULL;
+					
+	CLFlt			likeL, likeR, likeA, *clL, *clR, *clP, *clA;
+	CLFlt2                  *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
 	ModelInfo		*m;
 	
 	/* find model settings for this division and nStates, nStatesSquared */
@@ -5576,8 +5726,8 @@
 {
 
 	int				a, c, h, i, j, k, nStates, nCats;
-	CLFlt			*clL, *clR, *clP, *clA, *pL, *pR, *pA, *tiPL, *tiPR, *tiPA,
-					likeL, likeR, likeA;
+	CLFlt			*clL, *clR, *clP, *clA, likeL, likeR, likeA;
+	CLFlt2                  *pL, *pR, *pA, *tiPL, *tiPR, *tiPA;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -5649,7 +5799,8 @@
 {
 
 	int				c, k;
-	CLFlt			*clFA, *clFP, *clDP, *pA, *tiPA, likeUp[2];
+	CLFlt			*clFA, *clFP, *clDP, likeUp[2];
+	CLFlt2                  *pA, *tiPA;
 	ModelInfo		*m;
 	
 	/* find model settings for this division */
@@ -5717,7 +5868,8 @@
 {
 
 	int				a, c, i, j, k, nStates, nStatesSquared, nGammaCats;
-	CLFlt			*clFA, *clFP, *clDP, *pA, *tiPA, *likeUp, sum;
+	CLFlt			*clFA, *clFP, *clDP, *likeUp, sum;
+	CLFlt2                  *pA, *tiPA;
 	ModelInfo		*m;
 	
 	/* find model settings for this division */
@@ -5810,7 +5962,8 @@
 {
 
 	int				c, k, nGammaCats;
-	CLFlt			*clFA, *clFP, *clDP, *pA, *tiPA, likeUp[4];
+	CLFlt			*clFA, *clFP, *clDP, likeUp[4];
+	CLFlt2                  *pA, *tiPA;
 	ModelInfo		*m;
 	
 	/* find model settings for this division */
@@ -5891,7 +6044,8 @@
 {
 
 	int				a, c, i, j, k, nStates, nCats;
-	CLFlt			*clFA, *clFP, *clDP, *pA, *tiPA, likeUp[10], sum;
+	CLFlt			*clFA, *clFP, *clDP, likeUp[10], sum;
+	CLFlt2                  *pA, *tiPA;
 	ModelInfo		*m;
 	
 	/* find model settings for this division */
@@ -6038,12 +6192,11 @@
 |		variation
 |
 -----------------------------------------------------------------*/
-int CondLikeScaler_Gen (TreeNode *p, int division, int chain)
-
+int CondLikeScaler_Gen (TreeNode *p, int division, int chain, int do_full_scaling)
 {
 
-	int				c, i, j, k, n, nStates;
-	CLFlt			scaler, *clP, *scPNew, *scPOld, *lnScaler;
+	int				c, i, j, k, n, nStates, doing_full_scaling;
+	CLFlt			scaler, *clP, *scP, *lnScaler, oldLnScaler;
 	ModelInfo		*m;
 
 	m = &modelSettings[division];
@@ -6052,8 +6205,7 @@
 	/* find conditional likelihood pointer */
 	clP = condLikePtr[chain][p->index] + m->condLikeStart + Bit(division, &p->clSpace[0]) * condLikeRowSize;
 
-	scPNew = nodeScaler[chain] + m->compCharStart + (2 * (p->index - numLocalTaxa) +  Bit(division, &p->clSpace[0]   )) * numCompressedChars;
-	scPOld = nodeScaler[chain] + m->compCharStart + (2 * (p->index - numLocalTaxa) + (Bit(division, &p->clSpace[0])^1)) * numCompressedChars;
+	scP = nodeScaler[chain] + m->compCharStart + ((2*(p->index - numLocalTaxa)) + state[chain]) * numCompressedChars;
 
 	lnScaler = treeScaler[chain] + m->compCharStart + state[chain] * numCompressedChars;
 
@@ -6060,34 +6212,121 @@
 	/* subtract old values; these can also be associated with old scaler nodes */
 	if (Bit(division, &p->scalersSet[0]) == 1)
 		{
-		for (c=0; c<m->numChars; c++)
-			lnScaler[c] -= scPOld[c];	/*subtract old value */
+		  for (c=0; c<m->numChars; c++) {
+		    oldLnScaler = lnScaler[c];
+		    if ((lnScaler[c] > LOG_CL_WARN_THRESH) &&
+			(num_times_scaler_warn3 < MAX_TIMES_SCALER_WARN3)) {
+		      MrBayesPrint("%s   WARNING: for division %d (partition %d) char %d node %d scP is %1.30lf and lnScaler is (before subtracting scP) %1.30lf (numGammaCats is %d, nStates is %d)\n",
+				   spacer, division, division+1, c, p->index, scP[c], lnScaler[c],
+				   m->numGammaCats, nStates);
+		      num_times_scaler_warn3++;
+		      if (num_times_scaler_warn3 >= MAX_TIMES_SCALER_WARN3) {
+			MrBayesPrint("%s   WARNING: not warning again on scaler (3) until after next likelihood print\n",spacer);
+		      }
+		    }
+		    lnScaler[c] -= scP[c];	/*subtract old value */
+		    if (lnScaler[c] > LOG_CL_WARN_THRESH) {
+		      num_times_scaler_adjusted2++;
+		      if (((oldLnScaler > 0.0) ||
+			   (scP[c] > 0.0) ||
+			   (lnScaler[c] > (-0.01*(oldLnScaler+scP[c]))))) {
+			num_times_scaler_adjusted2_big++;
+			if (num_times_scaler_warn4 < MAX_TIMES_SCALER_WARN4) {
+			  MrBayesPrint("%s   WARNING: for division %d (partition %d) char %d node %d scP is %1.30lf and lnScaler is %1.30lf before subtracting scP, %1.30lf after subtracting scP (numGammaCats is %d, nStates is %d); adjusting to 0.0\n",
+				       spacer, division, division+1, c, p->index, scP[c], oldLnScaler,
+				       lnScaler[c], m->numGammaCats, nStates);
+			  num_times_scaler_warn4++;
+			  if (num_times_scaler_warn4 >= MAX_TIMES_SCALER_WARN4) {
+			    MrBayesPrint("%s   WARNING: not warning again on scaler (4) reduction\n",spacer);
+			  }
+			}
+		      }
+		      lnScaler[c] = (CLFlt) 0.0;
+		    }
+		    scP[c] = (CLFlt) 0.0;
+		  }
 		FlipOneBit(division, &p->scalersSet[0]);
 		}
 
-	/* add new values only if currently a scaler node */
+#if defined(ALTSCALING)
+	if ((p->left == NULL) || (! do_full_scaling))
+	  doing_full_scaling = 0;
+	else {
+	  doing_full_scaling = 1;
+	  did_scaling = 1;
+	}
+#else
 	if (p->scalerNode == YES)
-		{
+	  doing_full_scaling = 1;
+	else
+	  doing_full_scaling = 0;
+#endif
+	/* add new values only if currently a scaler node */
+	if (doing_full_scaling) {
 		i = j = 0;
-		for (c=0; c<m->numChars; c++)
-			{
+		for (c=0; c<m->numChars; c++) {
 			scaler = 0.0;
 			for (k=0; k<m->numGammaCats; k++)
 				{
-				for (n=0; n<nStates; n++)
-					{
-					if (clP[i] > scaler)
+				  for (n=0; n<nStates; n++) {
+				    if (clP[i] > scaler) {
+				      if ((clP[i] > CL_WARN_THRESH3) &&
+					  (num_times_cl_warn < MAX_TIMES_CL_WARN) &&
+					  ((lnScaler[c]+log(clP[i])) > SCALER_WARN_THRESH)) { /* Allen */
+					MrBayesPrint("%s   WARNING: for division %d (partition %d) char %d node %d GammaCat %d state %d clP is %1.30lf (numGammaCats is %d, nStates is %d, i is %d, scaler was %1.30lf, lnScaler was %1.30lf)\n",
+						     spacer, division, division+1, c, p->index, k, n, clP[i],
+						     m->numGammaCats, nStates, i, scaler, lnScaler[c]);
+					num_times_cl_warn++;
+					if (num_times_cl_warn >= MAX_TIMES_CL_WARN) {
+					  MrBayesPrint("%s   WARNING: not warning again on clP\n",spacer);
+					}
+				      }
 						scaler = clP[i];
+				    }
 					i++;
 					}
 				}
+#if defined(DEBUG_SCALER_NAN)
+			if ((! (isfinite(scaler))) || (scaler < (CLFlt)(1/(m->numGammaCats*nStates))))
+			  MrBayesPrint("%s   WARNING: for division %d (partition %d) char %d scaler is %1.30lf (numGammaCats is %d, nStates is %d, lnScaler was %1.30lf)\n",
+				       spacer, division, division+1, c, scaler, m->numGammaCats, nStates, lnScaler[c]);
+#endif
 			for (k=0; k<m->numGammaCats; k++)
 				{
 				for (n=0; n<nStates; n++)
 					clP[j++] /= scaler;
 				}
-			scPNew[c] = (CLFlt) log (scaler);	/* store scaler       */
-			lnScaler[c] += scPNew[c];	/* add to site scaler */
+			scP[c] = (CLFlt) log (scaler);	/* store scaler       */
+			if ((lnScaler[c] > LOG_CL_WARN_THRESH) &&
+			    (num_times_scaler_warn2 < MAX_TIMES_SCALER_WARN2)) {
+			  MrBayesPrint("%s   WARNING: for division %d (partition %d) char %d node %d scaler is %1.30lf, scP is %1.30lf, and lnScaler is (before adding scP) %1.30lf (numGammaCats is %d, nStates is %d)\n",
+				       spacer, division, division+1, c, p->index, scaler, scP[c], lnScaler[c],
+				       m->numGammaCats, nStates);
+			  num_times_scaler_warn2++;
+			  if (num_times_scaler_warn2 >= MAX_TIMES_SCALER_WARN2) {
+			    MrBayesPrint("%s   WARNING: not warning again on scaler (2)\n",spacer);
+			  }
+			}
+			lnScaler[c] += scP[c];	/* add to site scaler */
+			if (((scP[c] > 0.0) && (lnScaler[c] > LIKE_EPSILON)) ||
+			    (lnScaler[c] > LOG_CL_WARN_THRESH)) { /* Allen */
+			  num_times_scaler_adjusted++;
+			  if ((lnScaler[c] > LOG_CL_WARN_THRESH) &&
+			      (scaler > CL_WARN_THRESH)) {
+			    num_times_scaler_adjusted_big++;
+			    if (num_times_scaler_warn < MAX_TIMES_SCALER_WARN) {
+			      MrBayesPrint("%s   WARNING: for division %d (partition %d) char %d node %d scaler is %1.30lf, scP is %1.30lf, and lnScaler is %1.30lf (numGammaCats is %d, nStates is %d); reducing lnScaler to %1.6lf\n",
+					   spacer, division, division+1, c, p->index, scaler, scP[c], lnScaler[c],
+					   m->numGammaCats, nStates, (CLFlt) 0.0);
+			      num_times_scaler_warn++;
+			      if (num_times_scaler_warn >= MAX_TIMES_SCALER_WARN) {
+				MrBayesPrint("%s   WARNING: not warning again on scaler reduction\n",spacer);
+			      }
+			    }
+			  }
+			  scP[c] -= lnScaler[c];
+			  lnScaler[c] = (CLFlt) 0.0;
+			}
 			}
 		FlipOneBit(division, &p->scalersSet[0]);
 		}
@@ -6107,10 +6346,10 @@
 |		variation
 |
 -----------------------------------------------------------------*/
-int CondLikeScaler_NUC4 (TreeNode *p, int division, int chain)
+int CondLikeScaler_NUC4 (TreeNode *p, int division, int chain, int do_full_scaling)
 
 {
-	int				c, i, j, k;
+	int				c, i, j, k, doing_full_scaling;
 	CLFlt			scaler, *clP, *scPNew, *scPOld, *lnScaler;
 	ModelInfo		*m;
 	
@@ -6132,8 +6371,22 @@
 		FlipOneBit(division,&p->scalersSet[0]);	/* clear flag marking scalers set */
 		}
 
-	/* multiply by new values only if currently a scaler node */
+#if defined(ALTSCALING)
+	if ((p->left == NULL) || (! do_full_scaling))
+	  doing_full_scaling = 0;
+	else {
+	  doing_full_scaling = 1;
+	  did_scaling = 1;
+	}
+#else
 	if (p->scalerNode == YES)
+	  doing_full_scaling = 1;
+	else
+	  doing_full_scaling = 0;
+#endif
+
+	/* multiply by new values only if currently a scaler node */
+	if (doing_full_scaling)
 		{
 		i = j = 0;
 		for (c=0; c<m->numChars; c++)
@@ -6184,11 +6437,11 @@
 |		variation using fast log approximation
 |
 -----------------------------------------------------------------*/
-int CondLikeScaler_NUC4_fast (TreeNode *p, int division, int chain)
+int CondLikeScaler_NUC4_fast (TreeNode *p, int division, int chain, int do_full_scaling)
 
 {
-	int				c, i, j, k, x, *expPtr, log_2;
-	CLFlt			scaler, *scPNew, *scPOld, *lnScaler, *clP;
+	int				c, i, j, k, x, *expPtr, log_2, doing_full_scaling;
+	CLFlt			scaler, *clP, *scPNew, *scPOld, *lnScaler;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -6209,8 +6462,22 @@
 		FlipOneBit(division,&p->scalersSet[0]);	/* clear flag marking scalers set */
 		}
 
-	/* multiply by new values only if currently a scaler node */
+#if defined(ALTSCALING)
+	if ((p->left == NULL) || (! do_full_scaling))
+	  doing_full_scaling = 0;
+	else {
+	  doing_full_scaling = 1;
+	  did_scaling = 1;
+	}
+#else
 	if (p->scalerNode == YES)
+	  doing_full_scaling = 1;
+	else
+	  doing_full_scaling = 0;
+#endif
+
+	/* multiply by new values only if currently a scaler node */
+	if (doing_full_scaling)
 		{
 		i = j = 0;
 		for (c=0; c<m->numChars; c++)
@@ -6281,11 +6548,11 @@
 |	CondLikeScaler_NY98: codon model with omega variation
 |
 -----------------------------------------------------------------*/
-int CondLikeScaler_NY98 (TreeNode *p, int division, int chain)
+int CondLikeScaler_NY98 (TreeNode *p, int division, int chain, int do_full_scaling)
 
 {
 
-	int				c, i, j, k, n, nStates;
+	int				c, i, j, k, n, nStates, doing_full_scaling;
 	CLFlt			scaler, *clP, *scPNew, *scPOld, *lnScaler;
 	ModelInfo		*m;
 
@@ -6308,8 +6575,22 @@
 		FlipOneBit(division,&p->scalersSet[0]);
 		}
 
-	/* add new values only if currently a scaler node */
+#if defined(ALTSCALING)
+	if ((p->left == NULL) || (! do_full_scaling))
+	  doing_full_scaling = 0;
+	else {
+	  doing_full_scaling = 1;
+	  did_scaling = 1;
+	}
+#else
 	if (p->scalerNode == YES)
+	  doing_full_scaling = 1;
+	else
+	  doing_full_scaling = 0;
+#endif
+
+	/* add new values only if currently a scaler node */
+	if (doing_full_scaling)
 		{
 		i = j = 0;
 		for (c=0; c<m->numChars; c++)
@@ -6350,11 +6631,11 @@
 |		rate variation
 |
 -----------------------------------------------------------------*/
-int CondLikeScaler_Std (TreeNode *p, int division, int chain)
+int CondLikeScaler_Std (TreeNode *p, int division, int chain, int do_full_scaling)
 
 {
 
-	int				a, c, i, j, k, nStates, numReps;
+	int				a, c, i, j, k, nStates, numReps, doing_full_scaling;
 	CLFlt			scaler, *clP, *scPNew, *scPOld, *lnScaler;
 	ModelInfo		*m;
 	
@@ -6376,8 +6657,22 @@
 		FlipOneBit (division, &p->scalersSet[0]);
 		}
 
-	/* add new values only if currently a scaler node */
+#if defined(ALTSCALING)
+	if ((p->left == NULL) || (! do_full_scaling))
+	  doing_full_scaling = 0;
+	else {
+	  doing_full_scaling = 1;
+	  did_scaling = 1;
+	}
+#else
 	if (p->scalerNode == YES)
+	  doing_full_scaling = 1;
+	else
+	  doing_full_scaling = 0;
+#endif
+
+	/* add new values only if currently a scaler node */
+	if (doing_full_scaling)
 		{
 		i = j = 0;
 		for (c=0; c<m->numChars; c++)
@@ -7603,7 +7898,7 @@
 			else if (expecting == Expecting(NUMBER))
 				{
 				sscanf (tkn, "%d", &tempI);
-				runIDSeed = tempI;
+				runIDSeed = tempI; /* should the below be runIDSeed, not swapSeed? XXX - Allen */
 				MrBayesPrint ("%s   Setting run ID [stamp] seed to %ld [for GRID use]\n", spacer, swapSeed);
 				expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
 				}
@@ -7957,6 +8252,12 @@
 					free(tempStr);
 					return (ERROR);
 					}
+				if ((chainParams.useAdapt == YES) && (tempI > 1)) {
+				  MrBayesPrint ("%s   Too many runs (maximum of 1 run with Adapt)\n",
+						spacer);
+				  free(tempStr);
+				  return (ERROR);
+				}
 				chainParams.numRuns = tempI;
 				fileNameChanged = YES;
 				MrBayesPrint ("%s   Setting number of runs to %d\n", spacer, chainParams.numRuns);
@@ -7988,6 +8289,18 @@
 					free(tempStr);
 					return (ERROR);
 					}
+				if ((chainParams.useSA == YES) && (tempI > 1)) {
+				  MrBayesPrint ("%s   Too many chains (maximum of 1 chain with Simulated Annealing)\n",
+						spacer);
+				  free(tempStr);
+				  return (ERROR);
+				}
+				if ((chainParams.useAdapt == YES) && (tempI > 1)) {
+				  MrBayesPrint ("%s   Too many chains (maximum of 1 chain with Adapt)\n",
+						spacer);
+				  free(tempStr);
+				  return (ERROR);
+				}
 				chainParams.numChains = tempI;
 				MrBayesPrint ("%s   Setting number of chains to %d\n", spacer, chainParams.numChains);
 				expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
@@ -7999,9 +8312,18 @@
 				}
 			}
 		/* set Temp (chainTemp) ***************************************************************/
-		else if (!strcmp(parmName, "Temp"))
-			{
-			if (expecting == Expecting(EQUALSIGN))
+		else if (!strcmp(parmName, "Temp")) {
+		  if (chainParams.useSA == YES) {
+		    MrBayesPrint ("%s   Chain temps are incompatible with current Simulated Annealing implementation\n",
+				  spacer);
+		    free(tempStr);
+		    return (ERROR);
+		  } else if (chainParams.useAdapt == YES) {
+		    MrBayesPrint ("%s   Chain temps are incompatible with current Adapt implementation\n",
+				  spacer);
+		    free(tempStr);
+		    return (ERROR);
+		  } else if (expecting == Expecting(EQUALSIGN))
 				{
 				tempIndex = 0;
 				expecting = Expecting(NUMBER) | Expecting(LEFTPAR);
@@ -8047,9 +8369,80 @@
 				free(tempStr);
 				return (ERROR);
 				}
+		/* set SimAnneal (useSA) - Allen */
+		} else if (!strcmp(parmName, "SimAnneal")) {
+		  if (expecting == Expecting(EQUALSIGN))
+		    expecting = Expecting(ALPHA);
+		  else if (expecting == Expecting(ALPHA)) {
+		    if (IsArgValid(tkn, tempStr) == NO_ERROR) {
+		      if (!strcmp(tempStr, "Yes")) {
+			if (chainParams.numChains == 4) { /* default */
+			  MrBayesPrint ("%s   Setting chains to 1 due to SimAnneal\n", spacer);
+			  chainParams.numChains = 1;
+			} else if (chainParams.numChains > 1) {
+			  MrBayesPrint ("%s   %d chains invalid for SimAnneal\n", spacer, chainParams.numChains);
+			  free (tempStr);
+			  return (ERROR);
+			}
+			chainParams.useSA = YES;
+		      } else
+			chainParams.useSA = NO;
+		    } else {
+		      MrBayesPrint ("%s   Invalid argument for SimAnneal\n", spacer);
+		      free (tempStr);
+		      return (ERROR);
+		    }
+		    if (chainParams.useSA == YES)
+		      MrBayesPrint ("%s   Using (single-chain) Simulated Annealing\n", spacer);
+		    else
+		      MrBayesPrint ("%s   Not using (single-chain) Simulated Annealing\n", spacer);
+		    expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
+		  } else {
+		    free(tempStr);
+		    return (ERROR);
+		  }
+		/* set Adapt (useAdapt) - Allen */
+		} else if (!strcmp(parmName, "Adapt")) {
+		  if (expecting == Expecting(EQUALSIGN))
+		    expecting = Expecting(ALPHA);
+		  else if (expecting == Expecting(ALPHA)) {
+		    if (IsArgValid(tkn, tempStr) == NO_ERROR) {
+		      if (!strcmp(tempStr, "Yes")) {
+			if (chainParams.numChains == 4) { /* default */
+			  MrBayesPrint ("%s   Setting chains to 1 due to Adapt\n", spacer);
+			  chainParams.numChains = 1;
+			} else if (chainParams.numChains > 1) {
+			  MrBayesPrint ("%s   %d chains invalid for Adapt\n", spacer, chainParams.numChains);
+			  free (tempStr);
+			  return (ERROR);
+			}
+			if (chainParams.numRuns == 2) { /* default */
+			  MrBayesPrint ("%s   Setting runs to 1 due to Adapt\n", spacer);
+			  chainParams.numRuns = 1;
+			} else if (chainParams.numRuns > 1) {
+			  MrBayesPrint ("%s   %d runs invalid for Adapt\n", spacer, chainParams.numRuns);
+			  free (tempStr);
+			  return (ERROR);
+			}
+			chainParams.useAdapt = YES;
+		      } else
+			chainParams.useAdapt = NO;
+		    } else {
+		      MrBayesPrint ("%s   Invalid argument for Adapt\n", spacer);
+		      free (tempStr);
+		      return (ERROR);
+		    }
+		    if (chainParams.useAdapt == YES)
+		      MrBayesPrint ("%s   Using (single-chain) Adaptation\n", spacer);
+		    else
+		      MrBayesPrint ("%s   Not using (single-chain) Adaptation\n", spacer);
+		    expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
+		  } else {
+		    free(tempStr);
+		    return (ERROR);
 			}
 		/* set Reweight (weightScheme) ********************************************************/
-		else if (!strcmp(parmName, "Reweight"))
+		} else if (!strcmp(parmName, "Reweight"))
 			{
 			if (expecting == Expecting(EQUALSIGN))
 				expecting = Expecting(LEFTPAR);
@@ -8218,14 +8611,19 @@
 					free(tempStr);
 					return (ERROR);
 					}
-				if (tempD > 0.50)
-					{
-					MrBayesPrint ("%s   Burnin fraction too high (> 0.50)\n", spacer);
+				if (tempD > 0.90)
+				  { /* changed from 0.50 - Allen */
+					MrBayesPrint ("%s   Burnin fraction too high (> 0.90)\n", spacer);
 					free(tempStr);
 					return (ERROR);
 					}
 				chainParams.burninFraction = tempD;
-				MrBayesPrint ("%s   Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
+				if (((tempD > 0.50) && (chainParams.useSA != YES)) || (tempD > 0.75))
+				  MrBayesPrint ("%s   Setting burnin fraction to %.3f (this is high!)\n", spacer,
+						chainParams.burninFraction);
+				else
+				  MrBayesPrint ("%s   Setting burnin fraction to %.3f\n", spacer,
+						chainParams.burninFraction);
 				expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
 				}
 			else 
@@ -8287,7 +8685,7 @@
 					return (ERROR);
 					}
 				chainParams.stopVal = tempD;
-				MrBayesPrint ("%s   Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
+				MrBayesPrint ("%s   Setting stopval to %.6f\n", spacer, chainParams.stopVal);
 				expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
 				}
 			else 
@@ -8667,7 +9065,7 @@
 {
 
 	int			b, c, i, j, k, n, chn, index, nOfParams, tempInt;
-	MrBFlt		sum, *bs, *value, *subValue, symDir[10], scaler;
+	MrBFlt		sum, *bs, *value, *subValue, symDir[10], scaler, tempFlt;
 
 	Param		*p;
 	ModelInfo	*m;
@@ -8690,6 +9088,8 @@
 		MrBayesPrint ("%s   paramValues not free in FillNormalParams\n", spacer);
 		return ERROR;
 		}
+	MrBayesPrint("%s   Allocating paramValues for %d params\n", spacer,
+		     nOfParams); /* -Allen */
 	paramValues = (MrBFlt *) calloc (nOfParams, sizeof(MrBFlt));
 	if (!paramValues)
 		{
@@ -9008,16 +9408,12 @@
 
 					/* if user has not set dirichlet parameters, go with default */
 					/* overall variance equals number of states */
-					if (mp->numDirParams == 0)
-						for (i=0; i<mp->nStates; i++)
+					if (mp->numDirParams == 0) {
+						for (i=0; i<mp->nStates; i++) {
 							value[i] = mp->stateFreqsDir[i] = 1.0;
-					else
-						for (i=0; i<mp->nStates; i++)
-							value[i] = mp->stateFreqsDir[i];
-
 					/* now fill in subvalues */
-					for (i=0; i<mp->nStates; i++)
 						subValue[i] =  (1.0 / mp->nStates);
+						}
 #					if defined ASYMMETRY
 					if (mp->dataType == RESTRICTION)
 						{
@@ -9025,9 +9421,16 @@
 							subValue[i] =  0.5;
 						}
 #					endif
+					} else { /* Allen */
+						sum = 0.0;
+						for (i=0; i<mp->nStates; i++) {
+							value[i] = mp->stateFreqsDir[i];
+							sum += mp->stateFreqsDir[i];
 					}
-
-				else if (p->paramId == PI_USER)
+						for (i=0; i<mp->nStates; i++)
+							subValue[i] = mp->stateFreqsDir[i]/sum;
+					}
+				} else if (p->paramId == PI_USER)
 					{
 					for (i=0; i<mp->nStates; i++)
 						subValue[i] =  mp->stateFreqsFix[i];
@@ -9090,8 +9493,10 @@
 					{
 					if (GetEmpiricalFreqs (p->relParts, p->nRelParts) == ERROR)
 						return (ERROR);
-					for (i=0; i<mp->nStates; i++)
+					  for (i=0; i<mp->nStates; i++) {
 						value[i] = empiricalFreqs[i];
+						subValue[i] = empiricalFreqs[i]; /* Allen */
+					  }
 					}
 
 				else if (p->paramId == PI_EQUAL)
@@ -9154,12 +9559,14 @@
 				/* first get hyperprior */
 				if (p->paramId == SHAPE_UNI)
 					{
-					value[0] = 0.5;
-					if (value[0] < mp->shapeUni[0] || value[0] > mp->shapeUni[1])
-						value[0] = mp->shapeUni[0] + (mp->shapeUni[1] - mp->shapeUni[0]) *  0.5;
+					  value[0] = 0.65; /* Grishin distance - Allen */
+					  if (value[0] < mp->shapeUni[0])
+					    value[0] = mp->shapeUni[0] + (mp->shapeUni[1] - mp->shapeUni[0]) *  0.1;
+					  else if (value[0] > mp->shapeUni[1])
+					    value[0] = mp->shapeUni[0] + (mp->shapeUni[1] - mp->shapeUni[0]) *  0.9;
 					}
 				else if (p->paramId == SHAPE_EXP)
-					value[0] = 0.5;
+					value[0] = 0.65;
 				else if (p->paramId == SHAPE_FIX)
 					value[0] = mp->shapeFix;
 				/* now fill in rates */
@@ -9169,9 +9576,8 @@
 			else if (p->paramType == P_PINVAR)
 				{
 				/* Fill in pInvar ***************************************************************************************/
-				if (p->paramId == PINVAR_UNI)
-					value[0] = 0.0;
-
+				if (p->paramId == PINVAR_UNI) /* was 0.0 - Allen */
+				  value[0] = mp->pInvarUni[0] + (mp->pInvarUni[1] - mp->pInvarUni[0]) *  0.1;
 				else if (p->paramId == PINVAR_FIX)
 					value[0] =  mp->pInvarFix;
 				}
@@ -9178,21 +9584,23 @@
 			else if (p->paramType == P_CORREL)
 				{
 				/* Fill in correlation parameter of adgamma model *******************************************************/
-				if (p->paramId == CORREL_UNI)
+				  if (p->paramId == CORREL_UNI) { /* was just 0.0 - Allen */
 					value[0] = 0.0;
-
-				else if (p->paramId == CORREL_FIX)
+				    if (value[0] < mp->corrUni[0])
+				      value[0] = mp->corrUni[0];
+				    else if (value[0] > mp->corrUni[1])
+				      value[0] = mp->corrUni[1];
+				  } else if (p->paramId == CORREL_FIX)
 					value[0] =  mp->corrFix;
 				
 				/* Fill in correlation matrices */
 				AutodGamma (subValue, value[0], mp->numGammaCats);
 				}
-			else if (p->paramType == P_SWITCH)
-				{
+			else if (p->paramType == P_SWITCH) {
 				/* Fill in switchRates for covarion model ***************************************************************/
 				for (j=0; j<2; j++)
 					{
-					if (p->paramId == SWITCH_UNI)
+					if ((p->paramId == SWITCH_UNI) || (p->paramId == SWITCH_PINVUNI))
 						value[j] = RandomNumber(seed) * (mp->covswitchUni[1] - mp->covswitchUni[0]) + mp->covswitchUni[0];
 
 					else if (p->paramId == SWITCH_EXP)
@@ -9201,8 +9609,22 @@
 					else if (p->paramId == SWITCH_FIX)
 						value[j] = mp->covswitchFix[j];
 					}
+				if (p->paramId == SWITCH_PINVUNI) { /* - Allen */
+				  if ((value[0]+value[1]) < LIKE_EPSILON) {
+				    value[0] += 0.5*(mp->covswitchUni[1] - mp->covswitchUni[0]);
+				    value[1] += 0.5*(mp->covswitchUni[1] - mp->covswitchUni[0]);
 				}
-			else if (p->paramType == P_RATEMULT)
+				  tempFlt = value[1]/(value[0]+value[1]);
+				  if (tempFlt > mp->covswitchMaxPinvar) {
+				    value[1] = (mp->covswitchMaxPinvar*value[0])/(1-mp->covswitchMaxPinvar);
+				    if (value[1] < mp->covswitchUni[0]) {
+				      value[1] = mp->covswitchUni[0];
+				      value[0] =
+					(value[1]*(1-mp->covswitchMaxPinvar))/mp->covswitchMaxPinvar;
+				    }
+				  }
+				}
+			} else if (p->paramType == P_RATEMULT)
 				{
 				/* Fill in rateMult for division rates ******************************************************************/
 				for (j=0; j<p->nValues; j++)
@@ -9359,10 +9781,10 @@
 		goto errorExit;
 		}
 	memAllocs[ALLOC_NUMSITESOFPAT] = NO;
-	numSitesOfPat = (CLFlt *) realloc((void *) numSitesOfPat, numCompressedChars * chainParams.numChains * sizeof(MrBFlt));
+	numSitesOfPat = (CLFlt2 *) realloc((void *) numSitesOfPat, numCompressedChars * chainParams.numChains * sizeof(CLFlt2));
 	if (!numSitesOfPat)
 		{
-		MrBayesPrint ("%s   Problem reallocating numSitesOfPat (%d)\n", spacer, numCompressedChars * chainParams.numChains * sizeof(MrBFlt));
+		MrBayesPrint ("%s   Problem reallocating numSitesOfPat (%d)\n", spacer, numCompressedChars * chainParams.numChains * sizeof(CLFlt2));
 		goto errorExit;
 		}
 	memAllocs[ALLOC_NUMSITESOFPAT] = YES;
@@ -10563,7 +10985,7 @@
 	
 	int				c, i, n, division;
 	long			*pL, *pR, *pP, *pA, x;
-	CLFlt			*nSitesOfPat;
+	CLFlt2			*nSitesOfPat;
 	TreeNode        *p;
 	ModelInfo		*m;
 
@@ -10636,7 +11058,7 @@
 	
 	int				c, i, n, division;
 	long			*pL, *pR, *pP, x;
-	CLFlt			*nSitesOfPat; 
+	CLFlt2			*nSitesOfPat; 
 	TreeNode		*p;
 	ModelInfo		*m;
 
@@ -10685,7 +11107,7 @@
 	
 	int				c, i, n, division;
 	long			*pL, *pR, *pP, *pA, x;
-	CLFlt			*nSitesOfPat;
+	CLFlt2			*nSitesOfPat;
 	MrBFlt			length;
 	TreeNode		*p;
 	ModelInfo		*m;
@@ -11144,7 +11566,7 @@
 	int				i, j, c, d;
 	long			*cl, *clL, *clR, *clA, *clU, *clM, x, y, z;
 	CLFlt			temp1, temp2=0.0;
-	CLFlt			*nSitesOfPat;
+	CLFlt2			*nSitesOfPat;
 	TreeNode		*p, *q;
 	ModelInfo		*m;
 
@@ -11417,6 +11839,7 @@
 
 	/* For now, we wonly allow random swaps in the MPI version. Other schemes require
 	   tagging of messages, or a dedicated server node doing message processing.      */
+	/* ? that is the opposite of what the below code does... XXX - Allen */
 	(*swapA) = (int) (RandomNumber(&swapSeed) * chainParams.numChains);
 	(*swapB) = (int) (RandomNumber(&swapSeed) * (chainParams.numChains - 1));
 	if ((*swapB) == (*swapA))
@@ -11424,14 +11847,12 @@
 
 #	else
 
-	if (chainParams.swapAdjacentOnly == NO)
-		{
+	if (chainParams.swapAdjacentOnly == NO) {
 		(*swapA) = (int) (RandomNumber(&swapSeed) * chainParams.numChains);
 		(*swapB) = (int) (RandomNumber(&swapSeed) * (chainParams.numChains - 1));
 		if ((*swapB) == (*swapA))
 			(*swapB) = chainParams.numChains - 1;
-		}
-	else
+	} else
 		{
 		(*swapA) = (int) (RandomNumber(&swapSeed) * (chainParams.numChains - 1));
 		(*swapB) = (*swapA) + 1;
@@ -11523,7 +11944,8 @@
 
 	/* check if conditional likelihoods are needed */
 	if (oneMatSize > 0)
-		MrBayesPrint ("%s   Initializing conditional likelihoods for internal nodes\n", spacer);
+	  MrBayesPrint ("%s   Initializing conditional likelihoods for internal nodes (oneMatSize is %d)\n",
+			spacer, oneMatSize); /* - Allen */
 	else
 		return NO_ERROR;
 
@@ -11553,6 +11975,8 @@
 		MrBayesPrint ("%s   the conditional likelihoods for the chain. Please try\n", spacer);
 		MrBayesPrint ("%s   allocating more memory to the program, or running the\n", spacer);
 		MrBayesPrint ("%s   data on a computer with more memory.\n", spacer);
+		MrBayesPrint ("%s   numLocalChains: %d; oneMatSize: %d; nNodes: %d.\n",
+			      spacer,numLocalChains,oneMatSize,nNodes); /* Allen */
 		return ERROR;
 		}	
 	else
@@ -11578,6 +12002,8 @@
 		MrBayesPrint ("%s   the conditional likelihoods for the chain. Please try\n", spacer);
 		MrBayesPrint ("%s   allocating more memory to the program, or running the\n", spacer);
 		MrBayesPrint ("%s   data on a computer with more memory.\n", spacer);
+                MrBayesPrint ("%s   numLocalChains: %d; nScalerNodes: %d; numCompressedChars: %d.\n",
+			      spacer,numLocalChains,nScalerNodes,numCompressedChars); /* Allen */
 		if (treeScalerSpace) 
 			free (treeScalerSpace);
 		if (treeScaler) 
@@ -11803,15 +12229,15 @@
 		return ERROR;
 		}
 #if defined SSE
-	tiProbSpace = (CLFlt *) ALIGNEDMALLOC (2 * numLocalChains * nNodes * tiProbRowSize * sizeof(CLFlt), 16);
+	tiProbSpace = (CLFlt2 *) ALIGNEDMALLOC (2 * numLocalChains * nNodes * tiProbRowSize * sizeof(CLFlt2), 16);
 #else
-	tiProbSpace = (CLFlt *) SafeMalloc (2 * numLocalChains * nNodes * tiProbRowSize * sizeof(CLFlt));
+	tiProbSpace = (CLFlt2 *) SafeMalloc (2 * numLocalChains * nNodes * tiProbRowSize * sizeof(CLFlt2));
 	
 	/* PAUL: initialization of ti data, otherwise (according to valgrind) uninitialized data is used in some of
 	 * the likelihood functions. This might be a bug! */
-	/*memset(tiProbSpace, 0, 2 * numLocalChains * nNodes * tiProbRowSize * sizeof(CLFlt));*/
+	/*memset(tiProbSpace, 0, 2 * numLocalChains * nNodes * tiProbRowSize * sizeof(CLFlt2));*/
 #endif
-	tiProbs = (CLFlt **) SafeMalloc (numLocalChains * sizeof(CLFlt *));
+	tiProbs = (CLFlt2 **) SafeMalloc (numLocalChains * sizeof(CLFlt2 *));
 	if (!tiProbSpace || !tiProbs)
 		{
 		MrBayesPrint ("%s   The program ran out of memory while trying to allocate\n", spacer);
@@ -11886,7 +12312,10 @@
 			MrBayesPrint ("%s   The program ran out of memory while trying to allocate\n", spacer);
 			MrBayesPrint ("%s   the eigen values and vectors for the chain. Please try\n", spacer);
 			MrBayesPrint ("%s   allocating more memory to the program, or running the\n", spacer);
-			MrBayesPrint ("%s   data on a computer with more memory.\n", spacer);
+			MrBayesPrint ("%s   data on a computer with more memory. (2*%d*%d*%d=%d (%d/chain) bytes)\n",
+				      spacer, numLocalChains,cijkRowSize,sizeof(MrBFlt),
+				      2*numLocalChains*cijkRowSize*sizeof(MrBFlt),
+				      2*cijkRowSize*sizeof(MrBFlt));
 			if (cijkSpace) 
 				free (cijkSpace);
 			if (cijks) 
@@ -12651,9 +13080,10 @@
 
 {
 
-	int			c, d, i, j, k, s, maxRates, numReps, oneMatSize,
+  int			c, d, i, j, k, s, maxRates, numReps, oneMatSize, nObs,
 				corrModel[MAX_NUM_DIVS];
 	long int	*charBits;
+	MrBFlt          *bs, sum;
 	CLFlt		*cL;
 	ModelInfo	*m;
 	ModelParams *mp=NULL;
@@ -12690,7 +13120,7 @@
 
 	/* check if conditional likelihoods are needed */
 	if (oneMatSize > 0)
-		MrBayesPrint ("%s   Initializing conditional likelihoods for terminals\n", spacer);
+		MrBayesPrint ("%s   Initializing conditional likelihoods for terminals (oneMatSize is %d)\n", spacer, oneMatSize);
 	else
 		return NO_ERROR;
 
@@ -12833,7 +13263,7 @@
 
 	/* allocate and set space for terminal conditional likelihoods unless	*/
 	/* we are going to augment data											*/
-	for (d=0; d<numCurrentDivisions; d++)
+	for (d=0; d<numCurrentDivisions; d++) /* umm... mp = &modelParams[d] needed? - Allen */
 		if (!strcmp(mp->augmentData, "Yes"))
 			return NO_ERROR;
 
@@ -12854,7 +13284,9 @@
 		MrBayesPrint ("%s   The program ran out of memory while trying to allocate\n", spacer);
 		MrBayesPrint ("%s   the tip conditional likelihoods for the chain. Please try\n", spacer);
 		MrBayesPrint ("%s   allocating more memory to the program, or running the\n", spacer);
-		MrBayesPrint ("%s   data on a computer with more memory.\n", spacer);
+		MrBayesPrint ("%s   data on a computer with more memory (oneMatSize is %d).\n", spacer,
+			      oneMatSize);
+		/* - Allen */
 		return ERROR;
 		}	
 	
@@ -12863,8 +13295,7 @@
 		termCondLikes[i] = 0.0;
 
 	/* fill in conditional likelihoods */
-	for (d=0; d<numCurrentDivisions; d++)
-		{
+	for (d=0; d<numCurrentDivisions; d++) {
 		m = &modelSettings[d];
 		mp = &modelParams[d];
 		
@@ -12897,26 +13328,38 @@
 			/* TODO: fill in from compressed matrix */
 			
 			}
-		else
-			{
+		else { /* not doing this for standard - not sure as to proper statefreqs et al - Allen */
+		  bs = GetParamSubVals (m->stateFreq, 0, 0); /* as if chain 0 */
 			numReps = m->numModelStates / mp->nStates;
 			numReps *= m->numTiCats;
 
-			for (i=0; i<numLocalTaxa; i++)
-				{
+		  for (i=0; i<numLocalTaxa; i++) {
 				charBits = parsMatrix + m->parsMatrixStart + i * parsMatrixRowSize;
 				cL = termCondLikes + m->condLikeStart + i * condLikeRowSize;
-				for (c=m->compCharStart; c<m->compCharStop; c++)
-					{
-					for (k=0; k<numReps; k++)
-						{
-						for (s=0; s<mp->nStates; s++)
-							{
+		    for (c=m->compCharStart; c<m->compCharStop; c++) {
+		      for (k=0; k<numReps; k++) {
+			nObs = 0;
+			sum = 0.0;
+			for (s=0; s<mp->nStates; s++) {
+			  if (IsBitSet(s, charBits)) {
+			    nObs++;
+			    sum += bs[s];
+			  }
+			}
+			if ((nObs > 1) && (sum > LIKE_EPSILON)) {
+			  for (s=0; s<mp->nStates; s++) {
+			    if (IsBitSet(s, charBits))
+			      (*cL) = bs[s]/sum;
+			    cL++;
+			  }
+			} else {
+			  for (s=0; s<mp->nStates; s++) {
 							if (IsBitSet(s, charBits))
 								(*cL) = 1.0;
 							cL++;
 							}
 						}
+		      }
 					charBits += m->nParsIntsPerSite;
 					}
 				}
@@ -13452,10 +13894,11 @@
 int Likelihood_Gen (TreeNode *p, int division, int chain, MrBFlt *lnL, int whichSitePats)
 
 {
-	int				c, j, k, nStates, hasPInvar;
-	MrBFlt			s01, s10, probOn, probOff, *swr;
+	int				c, j, k, nStates, hasPInvar, toreturn = NO_ERROR;
+	MrBFlt			s01=50.0, s10=50.0, probOn, probOff, *swr;
 	MrBFlt			covBF[40], freq, *bs, *clInvar=NULL, like, likeI, pInvar=0.0, lnLike;
-	CLFlt			*clP, *lnScaler, *nSitesOfPat;
+	CLFlt			*clP, *lnScaler;
+	CLFlt2                  *nSitesOfPat;
 	ModelInfo		*m;
 	
 	/* find model settings and nStates, pInvar, invar cond likes */
@@ -13492,7 +13935,10 @@
 		for (j=0; j<nStates/2; j++)
 			{
 			covBF[j] = bs[j] * probOn;
-			covBF[j+nStates/2] = bs[j] * probOff;
+			/* try to reduce numerical inaccuracy - Allen */
+			covBF[j+(nStates/2)] = bs[j] - covBF[j];
+			if (covBF[j+(nStates/2)] < LIKE_EPSILON)
+			  covBF[j+(nStates/2)] = bs[j] * probOff;
 			}
 
 		/* finally set bs pointer to adjusted values */
@@ -13531,15 +13977,50 @@
 			like *= freq;
 
 			/* check against LIKE_EPSILON (values close to zero are problematic) */
-			if (like < LIKE_EPSILON)
-				{
-				MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
+			if (like < LIKE_EPSILON) {
+#if defined(ALTSCALING)
+			  had_like_epsilon = 1;
+			  if (! did_scaling)
+			    return ERROR;
+#endif
+		          if (m->switchRates != NULL)
+			    MrBayesPrint("%s   WARNING: In LIKE_EPSILON - for division %d (partition %d) char %d has like = %1.30lf (s01 is %2.15lf, s10 is %2.15lf, lnScaler = %lf)\n",
+					 spacer, division, division+1, c, like, s01, s10, lnScaler[c]);
+			  else
+			    MrBayesPrint("%s   WARNING: In LIKE_EPSILON - for division %d (partition %d) char %d has like = %1.30lf (lnScaler = %lf)\n",
+					 spacer, division, division+1, c, like, lnScaler[c]);
+			  toreturn = ERROR;
+			  /* failing to add in would be saying that is very _likely_... - Allen */
+			  (*lnL) += (lnScaler[c] + log(LIKE_EPSILON)) * nSitesOfPat[c];
+			} else {
+			  lnLike = lnScaler[c] + log(like);
+
+#if defined(ALTSCALING)
+			  if ((! did_scaling) && (lnLike > LIKE_EPSILON)) {
+			    if (m->switchRates != NULL)
+			      MrBayesPrint("%s   WARNING: Likelihood_Gen - for division %d (partition %d) char %d has like %1.30lf lnLike %1.30lf (s01 %2.15lf, s10 %2.15lf, lnScaler %lf)\n",
+					   spacer, division, division+1, c, like, lnLike, s01, s10, lnScaler[c]);
+			    else
+			      MrBayesPrint("%s   WARNING: Likelihood_Gen - for division %d (partition %d) char %d has like %1.30lf lnLike %1.30lf (lnScaler %lf)\n",
+					   spacer, division, division+1, c, like, lnScaler[c]);
+			    had_like_epsilon = 1;
 				return ERROR;
 				}
+#endif
+
+#ifdef DEBUG_SCALER_NAN
+			  if (! (isfinite(lnLike))) {
+			    if (m->switchRates != NULL)
+			      MrBayesPrint("%s   WARNING: For division %d (partition %d) char %d has lnLike %1.30lf (like = %1.30lf, s01 is %2.15lf, s10 is %2.15lf, lnScaler = %lf)\n",
+					   spacer, division, division+1, c, lnLike, like, s01, s10, lnScaler[c]);
 			else	
-				{
-				(*lnL) += (lnScaler[c] +  log(like)) * nSitesOfPat[c];
+			      MrBayesPrint("%s   WARNING: For division %d (partition %d) char %d has lnLike %1.30lf (like = %1.30lf, lnScaler = %lf)\n",
+					   spacer, division, division+1, c, lnLike, like, lnScaler[c]);
+			    return ERROR;
 				}
+#endif
+			  (*lnL) += lnLike * nSitesOfPat[c];
+			}
 			}
 		}
 	else
@@ -13573,24 +14054,43 @@
 			else
 				lnLike = log (like + (likeI / exp (lnScaler[c]))) + lnScaler[c];
 
+#ifdef DEBUG_SCALER_NAN
+			  if (! (isfinite(lnLike))) {
+			    	MrBayesPrint("%s   WARNING: For division %d (partition %d) char %d has lnLike %1.30lf (like = %1.30lf, likeI = %1.30lf, lnScaler = %lf)\n",
+					   spacer, division, division+1, c, lnLike, like, likeI, lnScaler[c]);
+				return ERROR;
+			  }
+#endif
 			/* check against LIKE_EPSILON (values close to zero are problematic) */
-			if (like < LIKE_EPSILON)
-				{
+			  if (((lnScaler[c] >= -200) || (likeI <= 1E-70)) && (like < LIKE_EPSILON)) {
+			    /* see "forget about like" above for why additional conditions - Allen */
 #ifdef DEBUG_OUTPUT
 				printf ("lnScaler[%d] = %lf likeI = %lf\n", c, lnScaler[c], likeI);
 #endif
-				MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
+#if defined(ALTSCALING)
+			    had_like_epsilon = 1;
+			    if (! did_scaling)
+			      return ERROR;
+#endif
+			    MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n",
+					  spacer, division, c, like);
 				getchar();
+			    /* leaving "return ERROR" because am not seeing this and don't want
+			       to alter without being able to test - Allen */
 				return ERROR;
-				}
-			else	
-				{
+#if defined(ALTSCALING)
+			  } else if ((! did_scaling) && (lnLike > LIKE_EPSILON)) {
+			    MrBayesPrint ("%s   WARNING: Likelihood_Gen - for division %d char %d has like %1.30lf, likeI %1.30lf, lnScaler %lf, lnLike %1.30lf\n",
+					  spacer, like, likeI, lnScaler[c], lnLike);
+			    had_like_epsilon = 1;
+			    return ERROR;
+#endif
+			  } else
 				(*lnL) += lnLike * nSitesOfPat[c];
 				}
 			}		
-		}
 		
-	return NO_ERROR;
+	return toreturn; /* admittedly, return value is actually ignored in all cases... XXX - Allen */
 	
 }
 
@@ -13610,7 +14110,8 @@
 
 	int				c, i, k, hasPInvar;
 	MrBFlt			freq, likeI, *bs, like, *clInvar=NULL, pInvar=0.0, lnLike;
-	CLFlt			*clP, *lnScaler, *nSitesOfPat;
+	CLFlt			*clP, *lnScaler;
+	CLFlt2                  *nSitesOfPat;
 	ModelInfo		*m;
 	
 	/* find model settings and pInvar, invar cond likes */
@@ -13663,6 +14164,11 @@
 			/* check against LIKE_EPSILON (values close to zero are problematic) */
 			if (like < LIKE_EPSILON)
 				{
+#if defined(ALTSCALING)
+				had_like_epsilon = 1;
+				if (! did_scaling)
+				  return ERROR;
+#endif
 				MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
 				return ERROR;
 				}
@@ -13707,6 +14213,11 @@
 			/* check against LIKE_EPSILON (values close to zero are problematic) */
 			if (like < LIKE_EPSILON)
 				{
+#if defined(ALTSCALING)
+				  had_like_epsilon = 1;
+				  if (! did_scaling)
+				    return ERROR;
+#endif
 				MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
 				return ERROR;
 				}
@@ -13736,7 +14247,8 @@
 
 	int				c, j, k, nStates;
 	MrBFlt			catLike, like, *bs, *omegaCatFreq;
-	CLFlt			*clP, *lnScaler, *nSitesOfPat;
+	CLFlt			*clP, *lnScaler;
+	CLFlt2                  *nSitesOfPat;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -13775,6 +14287,11 @@
 		/* check against LIKE_EPSILON (values close to zero are problematic) */
 		if (like < LIKE_EPSILON)
 			{
+#if defined(ALTSCALING)
+			  had_like_epsilon = 1;
+			  if (! did_scaling)
+			    return ERROR;
+#endif
 			MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
 			return ERROR;
 			}
@@ -13804,7 +14321,8 @@
 
 	int				c, k;
 	MrBFlt			*bs, freq, like, pUnobserved, pObserved;
-	CLFlt			*clP, *lnScaler, *nSitesOfPat;
+	CLFlt			*clP, *lnScaler;
+	CLFlt2                  *nSitesOfPat;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -13845,6 +14363,11 @@
 	pObserved =  1.0 - pUnobserved;
 	if (pObserved < LIKE_EPSILON)
 		{
+#if defined(ALTSCALING)
+		  had_like_epsilon = 1;
+		  if (! did_scaling)
+		    return ERROR;
+#endif
 		MrBayesPrint ("%s   WARNING: p(Observed) < LIKE_EPSILON - for division %d p(Observed) = %1.30lf\n", spacer, division, pObserved);
 		return ERROR;
 		}
@@ -13860,6 +14383,11 @@
 		/* check against LIKE_EPSILON (values close to zero are problematic) */
 		if (like < LIKE_EPSILON)
 			{
+#if defined(ALTSCALING)
+			  had_like_epsilon = 1;
+			  if (! did_scaling)
+			    return ERROR;
+#endif
 			MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
 			return ERROR;
 			}
@@ -13893,7 +14421,8 @@
 	int				b, c, j, k, nBetaCats, nGammaCats, nStates;
 	MrBFlt			catLike, catFreq, gammaFreq, like, *bs, *bsBase,
 					pUnobserved, pObserved;
-	CLFlt			*clP, *lnScaler, *nSitesOfPat;
+	CLFlt			*clP, *lnScaler;
+	CLFlt2                  *nSitesOfPat;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -13957,6 +14486,11 @@
 			/* check against LIKE_EPSILON (values close to zero are problematic) */
 			if (like < LIKE_EPSILON)
 				{
+#if defined(ALTSCALING)
+				  had_like_epsilon = 1;
+				  if (! did_scaling)
+				    return ERROR;
+#endif
 				MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
 				return ERROR;
 				}
@@ -14033,6 +14567,11 @@
 			/* check against LIKE_EPSILON (values close to zero are problematic) */
 			if (like < LIKE_EPSILON)
 				{
+#if defined(ALTSCALING)
+				  had_like_epsilon = 1;
+				  if (! did_scaling)
+				    return ERROR;
+#endif
 				MrBayesPrint ("%s   WARNING: In LIKE_EPSILON - for division %d char %d has like = %1.30lf\n", spacer, division, c, like);
 				return ERROR;
 				}
@@ -14087,7 +14626,8 @@
 	int				c, i, nStates;
 	long			done, *pL, *pR, *pP, *pA, *oldpP, x;
 	CLFlt			nParsChars, treeLength;
-	CLFlt			length, *nSitesOfPat, *newNodeLength, oldNodeLength;
+	CLFlt			length, *newNodeLength, oldNodeLength;
+	CLFlt2                  *nSitesOfPat;
 	Tree			*t;
 	ModelInfo		*m;
 
@@ -14268,7 +14808,7 @@
 	int				c, i, *nStates;
 	long			*pL, *pR, *pP, *pA, *oldpP, x;
 	CLFlt			*treeLength;
-	CLFlt			*nSitesOfPat;
+	CLFlt2			*nSitesOfPat;
 	Tree			*t;
 	ModelInfo		*m;
 
@@ -14389,8 +14930,11 @@
 			{
 			tree = GetTree(m->brlens, chain, state[chain]);
 
-			if (m->parsModelId == NO)
-				{
+			if (m->parsModelId == NO) {
+#if defined(ALTSCALING)
+			  had_like_epsilon = 0;
+			  did_scaling = 0;
+#endif
 				for (i=0; i<tree->nIntNodes; i++)
 					{
 					p = tree->intDownPass[i];
@@ -14432,8 +14976,13 @@
 							m->CondLikeDown (p, d, chain);
 
 #						if !defined (DEBUG_NOSCALING)
+#if defined(ALTSCALING)
+						if (Bit(d,&p->scalersSet[0]) == 1)
+						  m->CondLikeScaler (p, d, chain, 0);
+#else
 						if (p->scalerNode == YES || Bit(d,&p->scalersSet[0]) == 1)
-							m->CondLikeScaler (p, d, chain);
+						  m->CondLikeScaler (p, d, chain, 1);
+#endif
 #						endif
 						}
 					}
@@ -14441,6 +14990,30 @@
 				
 			lnL = 0.0;
 			m->Likelihood (tree->root->left, d, chain, &lnL, (chainId[chain] % chainParams.numChains));
+#if defined(ALTSCALING)
+			if ((m->parsModelId == NO) && had_like_epsilon) {
+			  had_like_epsilon = 0; /* reset */
+			  for (i=0; i<tree->nIntNodes; i++) {
+			    p = tree->intDownPass[i];
+			    if ((p->upDateCl == YES) && (p->left != NULL)) /* suspect latter check not necc... - Allen */
+			      m->CondLikeScaler (p, d, chain, 1);
+			  }
+			  if (! did_scaling) {
+			    MrBayesPrint("%s   WARNING: For division %d (partition %d) chain %d no scaling done despite repeat\n",
+					 spacer, d, d+1, chain);
+			    did_scaling = 1; /* so that Likelihood_Gen will do the right thing */
+			  }
+			  lnL = 0.0;
+			  m->Likelihood (tree->root->left, d, chain, &lnL, (chainId[chain] % chainParams.numChains));
+			  MrBayesPrint("%s   WARNING: For division %d (partition %d) chain %d lnL is %1.30lf after scaling (chainLnLike before adding is %1.30lf)\n",
+				       spacer, d, d+1, chain, lnL, chainLnLike);
+			}
+#endif
+#ifdef DEBUG_SCALER_NAN
+			if (! (isfinite(lnL)))
+			  MrBayesPrint("%s   WARNING: For division %d (partition %d) chain %d lnL is %1.30lf (chainLnLike is %1.30lf)\n",
+				       spacer, d, d+1, chain, lnL, chainLnLike);
+#endif
 			m->lnLike[2*chain + state[chain]] =  lnL;
 			}
 		chainLnLike += m->lnLike[2*chain + state[chain]];
@@ -14568,7 +15141,7 @@
 				{
 				alphaDir = mp->revMatDir;
 				sum = 0.0;
-				for (i=0; i<6; i++)
+				for (i=0; i<6; i++) /* ASSUMES NUC? XXX - Allen */
 					sum += st[i];
 				for (i=0; i<6; i++)
 					newProp[i] = st[i] / sum;
@@ -14645,7 +15218,9 @@
 			lnPrior += log(1.0) - log(mp->corrUni[1] - mp->corrUni[0]);
 			}
 		else if (p->paramType == P_SWITCH)
-			{
+			{ /* Note that, since covswitchUni values are 0-100, the below will usually be subtracting
+			     a positive number, unlike P_PINVAR, which will be taking the log of a number below 1 -
+			     giving a negative value - and subtracting it, which will increase the lnPrior - Allen */
 			/* switching rate parameter of covarion model */
 			mp = &modelParams[p->relParts[0]];
 			if (p->paramId == SWITCH_UNI)
@@ -14656,8 +15231,12 @@
 				{
 				lnPrior += log(mp->covswitchExp) - mp->covswitchExp * st[0];
 				}
-			}
-		else if (p->paramType == P_RATEMULT)
+			else if (p->paramId == SWITCH_PINVUNI) { /* best guess! - Allen */
+			  x = mp->covswitchUni[0]/(mp->covswitchUni[0]+mp->covswitchUni[1]); /* min pinvar */
+			  lnPrior += log(1.0) - log(mp->covswitchUni[1] - mp->covswitchUni[0])
+			    - log(mp->covswitchMaxPinvar - x);
+			} /* Have decided not to try doing SWITCH_PINVEXP due to uncertainty - Allen */
+		} else if (p->paramType == P_RATEMULT)
 			{
 			
 			
@@ -14934,11 +15513,12 @@
 {
 
 	int				i, j, k, nNodes;
-	MrBFlt			*ct, tempD, lastCoalescenceTime, coalescenceTime, intervalLength;
+	MrBFlt			tempD, lastCoalescenceTime, coalescenceTime, intervalLength;
+	CLFlt2                  *ct;
 	TreeNode		*p;
 
 	/* allocate space for the coalescence times */
-	ct = (MrBFlt *)SafeMalloc((size_t) (t->nIntNodes) * sizeof(MrBFlt));
+	ct = (CLFlt2 *)SafeMalloc((size_t) (t->nIntNodes) * sizeof(CLFlt2));
 	if (!ct)
 		{
 		printf ("\n   ERROR: Problem allocating ct\n");
@@ -14950,7 +15530,7 @@
 		{
 		p = t->intDownPass[i];
 		if (p->anc != NULL)
-			ct[j++] = p->nodeDepth;
+			ct[j++] = (CLFlt2) p->nodeDepth;
 		}
 	nNodes = j;
 
@@ -15219,7 +15799,7 @@
 	/* Change correlation parameter (-1, 1) of adgamma model */
 
 	int			i, isValidP;
-	MrBFlt		oldP, newP, window, minP, maxP, ran, *markovTiValues;
+	MrBFlt		oldP, newP, window, minP, maxP, wasNewP, ran, *markovTiValues;
 	ModelParams *mp;
 
 	/* get size of window, centered on current rho */
@@ -15244,13 +15824,26 @@
 	
 	/* check that new value is valid */
 	isValidP = NO;
-	do
-		{
-		if (newP < minP)
+	do {
+	  if (newP < minP) {
+	    wasNewP = newP;
 			newP = 2* minP - newP;
-		else if (newP > maxP)
+	    if (newP > maxP) { /* avoid infinite loops - Allen */
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_Adgamma): newP %2.6lf orig %2.6lf tmp %2.6lf minP %2.6lf maxP %2.6lf window %2.6lf ran %2.6lf\n",
+                           spacer, newP, oldP, wasNewP, minP, maxP, window, ran);
+              return(NO_ERROR);
+	    }
+	  } else if (newP > maxP) {
+	    wasNewP = newP;
 			newP = 2 * maxP - newP;
-		else
+	    if (newP < minP) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_Adgamma): newP %2.6lf orig %2.6lf tmp %2.6lf minP %2.6lf maxP %2.6lf window %2.6lf ran %2.6lf\n",
+                           spacer, newP, oldP, wasNewP, minP, maxP, window, ran);
+              return(NO_ERROR);
+	    }
+	  } else
 			isValidP = YES;
 		} while (isValidP == NO);
 
@@ -16019,7 +16612,7 @@
 
 	/* change one branch length */
 
-	MrBFlt		tuning, maxV, minV, m, newM, brlensPrExp=0.0;
+	MrBFlt		tuning, maxV, minV, m, newM, tmpM, brlensPrExp=0.0;
 	TreeNode	*p;
 	ModelParams *mp;
 	Tree		*t;
@@ -16048,7 +16641,8 @@
 	do
 		{
 		p = t->allDownPass[(int)(RandomNumber(seed)*t->nNodes)];
-		} while (p->anc == NULL || (t->isRooted == YES && p->anc->anc == NULL));
+		} while (p->anc == NULL || (t->isRooted == YES && p->anc->anc == NULL) ||
+			 (p->length <= 0.0)); /* - Allen */
 
 	/* determine new length */
 	m = p->length;
@@ -16055,12 +16649,32 @@
 	newM = m * exp(tuning * (RandomNumber(seed) - 0.5));
 
 	/* reflect new length if necessary */
-	while (newM < minV || newM > maxV)
-		{
-		if (newM < minV)
+	/* Putting in below code so that will not go into an infinite loop - Allen */
+	while (newM < minV || newM > maxV) {
+	  if (newM < minV) {
+	    tmpM = newM;
 			newM = minV * minV / newM;
-		else if (newM > maxV)
+	    if (newM > maxV) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_BrLen): m %1.6lf orig newM %1.6lf newM %1.6lf minV %1.6lf maxV %1.6lf\n",
+			   spacer, m, tmpM, newM, minV, maxV);
+	      return(NO_ERROR);
+	    }
+	  } else if (newM > maxV) {
+	    tmpM = newM;
 			newM = maxV * maxV / newM;
+	    if (newM < minV) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_BrLen): m %1.6lf orig newM %1.6lf newM %1.6lf minV %1.6lf maxV %1.6lf\n",
+			   spacer, m, tmpM, newM, minV, maxV);
+	      return(NO_ERROR);
+	    }
+	  }
+	}
+	/* getting some extremely small, rather useless moves - Allen */
+	if (fabs(newM - p->length) < BRLEN_EPSILON) {
+	  abortMove = YES;
+	  return(NO_ERROR);
 		}
 	p->length = newM;
 
@@ -16697,7 +17311,7 @@
 					b->isLocked = YES;
 					u->isLocked = NO;
 					b->lockID = u->lockID;
-					u->lockID = 0;
+					u->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			else 
@@ -16733,7 +17347,7 @@
 						u->isLocked = YES;
 						b->isLocked = NO;
 						u->lockID = b->lockID;
-						b->lockID = 0;
+						b->lockID = -1; /* to -1 from 0 - Allen */
 						}
 					}
 				}
@@ -17233,7 +17847,7 @@
 					b->isLocked = YES;
 					u->isLocked = NO;
 					b->lockID = u->lockID;
-					u->lockID = 0;
+					u->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			else 
@@ -17269,7 +17883,7 @@
 						u->isLocked = YES;
 						b->isLocked = NO;
 						u->lockID = b->lockID;
-						b->lockID = 0;
+						b->lockID = -1; /* to -1 from 0 - Allen */
 						}
 					}
 				}
@@ -17731,7 +18345,7 @@
 				b->isLocked = YES;
 				u->isLocked = NO;
 				b->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			}
 		else 
@@ -17765,7 +18379,7 @@
 					u->isLocked = YES;
 					a->isLocked = NO;
 					u->lockID = a->lockID;
-					a->lockID = 0;
+					a->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			}
@@ -18572,13 +19186,27 @@
 	/* record c length and adjust with multiplier using reflection */
 	m = c->length;
 	x = c->length * exp(tuning * (RandomNumber(seed) - 0.5));		/* save the modified dangling branch for later use */
-	while (x < minV || x > maxV)
-		{
-		if (x < minV)
+	while (x < minV || x > maxV) {
+	  if (x < minV) {
+	    y = x;
 			x = minV * minV / x;
-		else if (x > maxV)
+	    if (x > maxV) {
+              abortMove = YES;
+              MrBayesPrint("%s   ERROR (Move_ExtTBR; 1): x %2.6lf orig %2.6lf tmp %2.6lf minV %2.6lf maxV %2.6lf tuning %2.6lf\n",
+                           spacer, x, m, y, minV, maxV, tuning);
+              return(NO_ERROR);
+            }
+	  } else if (x > maxV) {
+	    y = x;
 			x = maxV * maxV / x;
+	    if (x < minV) {
+              abortMove = YES;
+              MrBayesPrint("%s   ERROR (Move_ExtTBR; 1): x %2.6lf orig %2.6lf tmp %2.6lf minV %2.6lf maxV %2.6lf tuning %2.6lf\n",
+                           spacer, x, m, y, minV, maxV, tuning);
+              return(NO_ERROR);
+            }
 		}
+	}
 	
 	/* calculate proposal and prior ratio based on length modification */
 	(*lnProposalRatio) = log (x / m);
@@ -18588,13 +19216,27 @@
 	/* record v length and adjust with multiplier using reflection*/
 	m = v->length;
 	v->length *= exp(tuning * (RandomNumber(seed) - 0.5));
-	while (v->length < minV || v->length > maxV)
-		{
-		if (v->length < minV)
+	while (v->length < minV || v->length > maxV) {
+	  if (v->length < minV) {
+	    y = v->length;
 			v->length = minV * minV / v->length;
-		else if (v->length > maxV)
+	    if (v->length > maxV) {
+              abortMove = YES;
+              MrBayesPrint("%s   ERROR (Move_ExtTBR; 2): len %2.6lf orig %2.6lf tmp %2.6lf minV %2.6lf maxV %2.6lf tuning %2.6lf\n",
+                           spacer, v->length, m, y, minV, maxV, tuning);
+              return(NO_ERROR);
+            }
+	  } else if (v->length > maxV) {
+	    y = v->length;
 			v->length = maxV * maxV / v->length;
+	    if (v->length < minV) {
+              abortMove = YES;
+              MrBayesPrint("%s   ERROR (Move_ExtTBR; 2): len %2.6lf orig %2.6lf tmp %2.6lf minV %2.6lf maxV %2.6lf tuning %2.6lf\n",
+                           spacer, v->length, m, y, minV, maxV, tuning);
+              return(NO_ERROR);
 		}
+	  }
+	}
 	v->upDateTi = YES;
 
 	/* adjust proposal and prior ratio based on length modification */
@@ -18648,16 +19290,29 @@
 
 	/* adjust length of branch to be modified */
 	/* if it is not the root branch of a rooted tree */
-	if (t->isRooted == NO || u->anc->anc != NULL) 
-		{
+	if (t->isRooted == NO || u->anc->anc != NULL) {
 		m = u->length;
 		u->length *= exp(tuning * (RandomNumber(seed) - 0.5));
-		while (u->length < minV || u->length > maxV)
-			{
-			if (u->length < minV)
+	  while (u->length < minV || u->length > maxV) {
+	    if (u->length < minV) {
+	      y = u->length;
 				u->length = minV * minV / u->length;
-			else if (u->length > maxV)
+	      if (u->length > maxV) {
+		abortMove = YES;
+		MrBayesPrint("%s   ERROR (Move_ExtTBR; 3): len %2.6lf orig %2.6lf tmp %2.6lf minV %2.6lf maxV %2.6lf tuning %2.6lf\n",
+			     spacer, u->length, m, y, minV, maxV, tuning);
+		return(NO_ERROR);
+	      }
+	    } else if (u->length > maxV) {
+	      y = u->length;
 				u->length = maxV * maxV / u->length;
+	      if (u->length < minV) {
+		abortMove = YES;
+		MrBayesPrint("%s   ERROR (Move_ExtTBR; 3): len %2.6lf orig %2.6lf tmp %2.6lf minV %2.6lf maxV %2.6lf tuning %2.6lf\n",
+			     spacer, u->length, m, y, minV, maxV, tuning);
+		return(NO_ERROR);
+	      }
+	    }
 			}
 
 		/* adjust proposal and prior ratio based on length modification */
@@ -18708,7 +19363,7 @@
 				b->isLocked = YES;
 				u->isLocked = NO;
 				b->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			}
 		else 
@@ -18744,7 +19399,7 @@
 					u->isLocked = YES;
 					b->isLocked = NO;
 					u->lockID = b->lockID;
-					b->lockID = 0;
+					b->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			}
@@ -19138,7 +19793,7 @@
 				b->isLocked = YES;
 				u->isLocked = NO;
 				b->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			}
 		else 
@@ -19174,7 +19829,7 @@
 					u->isLocked = YES;
 					b->isLocked = NO;
 					u->lockID = b->lockID;
-					b->lockID = 0;
+					b->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			}
@@ -19662,7 +20317,7 @@
 				b->isLocked = YES;
 				u->isLocked = NO;
 				b->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			}
 		else 
@@ -19698,7 +20353,7 @@
 					u->isLocked = YES;
 					b->isLocked = NO;
 					u->lockID = b->lockID;
-					b->lockID = 0;
+					b->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			}
@@ -20236,7 +20891,7 @@
 				b->isLocked = YES;
 				u->isLocked = NO;
 				b->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			}
 		else 
@@ -20272,7 +20927,7 @@
 					u->isLocked = YES;
 					b->isLocked = NO;
 					u->lockID = b->lockID;
-					b->lockID = 0;
+					b->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			}
@@ -20854,7 +21509,7 @@
 				b->isLocked = YES;
 				u->isLocked = NO;
 				b->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			}
 		else 
@@ -20890,7 +21545,7 @@
 					u->isLocked = YES;
 					b->isLocked = NO;
 					u->lockID = b->lockID;
-					b->lockID = 0;
+					b->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			}
@@ -21164,7 +21819,6 @@
 #endif
 
 	return (NO_ERROR);
-	
 }
 
 
@@ -21178,10 +21829,9 @@
 	/* change gamma shape parameter using multiplier */
 	
 	int			i, isAPriorExp, isValidA;
-	MrBFlt		oldA, newA, minA, maxA, alphaExp=0.0, ran, factor, tuning, *gammaRates;
+	MrBFlt		oldA, newA, wasNewA, minA, maxA, alphaExp=0.0, ran, factor, tuning, *gammaRates;
 	ModelParams *mp;
 
-	/* get tuning parameter */
 	tuning = mvp[0];
 
 	/* get model params */
@@ -21216,13 +21866,27 @@
 
 	/* check validity */
 	isValidA = NO;
-	do
-		{
-		if (newA < minA)
+	/* Putting in below code so that will not go into an infinite loop (have seen this) - Allen */
+	do {
+	  if (newA < minA) {
+	    wasNewA = newA;
 			newA = minA * minA / newA;
-		else if (newA > maxA)
+	    if (newA > maxA) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_GammaShape_M): newA %2.6lf orig %2.6lf tmp %2.6lf minA %2.6lf maxA %2.6lf tuning %2.6lf factor %2.6lf ran %2.6lf\n",
+                           spacer, newA, oldA, wasNewA, minA, maxA, tuning, factor, ran);
+              return(NO_ERROR);
+	    }
+	  } else if (newA > maxA) {
+	    wasNewA = newA;
 			newA = maxA * maxA / newA;
-		else
+	    if (newA < minA) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_GammaShape_M): newA %2.6lf orig %2.6lf tmp %2.6lf minA %2.6lf maxA %2.6lf tuning %2.6lf factor %2.6lf ran %2.6lf\n",
+                           spacer, newA, oldA, wasNewA, minA, maxA, tuning, factor, ran);
+              return(NO_ERROR);
+	    }
+	  } else
 			isValidA = YES;
 		} while (isValidA == NO);
 
@@ -21583,7 +22247,7 @@
 				v->isLocked = YES;
 				u->isLocked = NO;
 				v->lockID = u->lockID;
-				u->lockID = 0;
+				u->lockID = -1; /* to -1 from 0 - Allen */
 				}
 			if (v->left == d)
 				v->right = u;
@@ -22468,7 +23132,7 @@
 int Move_NodeSlider (Param *param, int chain, long *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp)
 
 {
-	MrBFlt		tuning, maxV, minV, oldM, newM, brlensPrExp=0.0, newMin, newMax, oldMin, oldMax;
+	MrBFlt		tuning, maxV, minV, oldM, newM, brlensPrExp=0.0, newMin, newMax, oldMin, oldMax, tmpM;
 	TreeNode	*p, *q;
 	ModelParams *mp;
 	Tree		*t;
@@ -22508,12 +23172,29 @@
 	/* determine new length */
 	oldM = (q->length + p->length);
 	newM = oldM * exp(tuning * (RandomNumber(seed) - 0.5));
-	while (newM < 2.0 * minV || newM > 2.0 * maxV)
-		{
-		if (newM < 2.0 * minV)
-			newM = 4.0 * minV * minV / newM;
-		else if (newM > 2.0 * maxV)
-			newM = 4.0 * maxV * maxV / newM;
+	/* Putting in below code so that will not go into an infinite loop - Allen */
+	while ((newM < (2.0*minV)) || (newM > (2.0*maxV))) {
+	  if (newM < (2.0*minV)) {
+	    tmpM = newM;
+	    newM = (4.0 * minV * minV)/newM;
+	    if (newM > (2.0*maxV)) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_NodeSlider): oldM %1.6lf (%1.6lf %1.6lf) orig newM %1.6lf newM %1.6lf 2*minV %1.6lf 2*maxV %1.6lf\n",
+			   spacer, oldM, q->length, p->length, tmpM, newM, (2.0*minV),
+			   (2.0*maxV));
+	      return(NO_ERROR);
+	    }
+	  } else if (newM > (2.0*maxV)) {
+	    tmpM = newM;
+	    newM = (4.0 * maxV * maxV)/newM;
+	    if (newM < (2.0*minV)) {
+	      abortMove = YES;
+	      MrBayesPrint("%s   ERROR (Move_NodeSlider): oldM %1.6lf (%1.6lf %1.6lf) orig newM %1.6lf newM %1.6lf 2*minV %1.6lf 2*maxV %1.6lf\n",
+			   spacer, oldM, q->length, p->length, tmpM, newM, (2.0*minV),
+			   (2.0*maxV));
+	      return(NO_ERROR);
+	    }
+	  }
 		}
 
 	/* determine new lengths of p and q */
@@ -22549,7 +23230,148 @@
 }
 
 
+/*-----------------------------------------------------------------------------------
+|
+|	Move_NodeSlider2: move the position of one node without changing topology
+|
+-------------------------------------------------------------------------------------*/
+
+int Move_NodeSlider2 (Param *param, int chain, long *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp) {
+  MrBFlt      tuning, maxV, minV, oldM, newM, brlensPrExp=0.0, newMin, newMax, oldMin, oldMax, tmpM;
+  MrBFlt      oldProp, newProp, minProp, maxProp, window, tmpProp;
+  int         isValidP;
+  TreeNode    *p, *q;
+  ModelParams *mp;
+  Tree	      *t;
+
+  tuning = mvp[0]; /* Larget & Simon's tuning parameter lambda */
+  window = mvp[1]; /* window size for proportion - Allen */
+  
+  mp = &modelParams[param->relParts[0]];
+  
+  /* max and min brlen */
+  if (param->paramId == BRLENS_UNI) {
+    minV = mp->brlensUni[0] > BRLENS_MIN ? mp->brlensUni[0] : BRLENS_MIN;
+    maxV = mp->brlensUni[1];
+  } else {
+    minV = BRLENS_MIN;
+    maxV = BRLENS_MAX;
+    brlensPrExp = mp->brlensExp;
+  }
+
+  /* get tree */
+  t = GetTree (param, chain, state[chain]);
+  
+  /* pick an interior branch */
+  do {
+    p = t->intDownPass[(int)(RandomNumber(seed)*t->nIntNodes)];
+  } while (p->anc == NULL || (t->isRooted == YES && p->anc->anc == NULL));
+
+  /* pick one descendant branch */
+  if (RandomNumber(seed) < 0.5)
+    q = p->left;
+  else
+    q = p->right;
+  
+  /* determine new length */
+  oldM = (q->length + p->length);
+  oldProp = q->length/oldM;
+  newM = oldM * exp(tuning * (RandomNumber(seed) - 0.5));
+  /* Putting in below code so that will not go into an infinite loop - Allen */
+  while ((newM < (2.0*minV)) || (newM > (2.0*maxV))) {
+    if (newM < (2.0*minV)) {
+      tmpM = newM;
+      newM = (4.0 * minV * minV)/newM;
+      if (newM > (2.0*maxV)) {
+	abortMove = YES;
+	MrBayesPrint("%s   ERROR (Move_NodeSlider2): oldM %1.6lf (%1.6lf %1.6lf) orig newM %1.6lf newM %1.6lf 2*minV %1.6lf 2*maxV %1.6lf\n",
+		     spacer, oldM, q->length, p->length, tmpM, newM, (2.0*minV),
+		     (2.0*maxV));
+	return(NO_ERROR);
+      }
+    } else if (newM > (2.0*maxV)) {
+      tmpM = newM;
+      newM = (4.0 * maxV * maxV)/newM;
+      if (newM < (2.0*minV)) {
+	abortMove = YES;
+	MrBayesPrint("%s   ERROR (Move_NodeSlider2): oldM %1.6lf (%1.6lf %1.6lf) orig newM %1.6lf newM %1.6lf 2*minV %1.6lf 2*maxV %1.6lf\n",
+		     spacer, oldM, q->length, p->length, tmpM, newM, (2.0*minV),
+		     (2.0*maxV));
+	return(NO_ERROR);
+      }
+    }
+  }
+
+  /* determine new lengths of p and q */
+  newMin = minV > newM - maxV ? minV : newM - maxV;
+  newMax = maxV < newM - minV ? maxV : newM - minV;
+  oldMin = minV > oldM - maxV ? minV : oldM - maxV;
+  oldMax = maxV > oldM - minV ? maxV : oldM - minV;
+  minProp = newMin/newM;
+  maxProp = newMax/newM;
+  
+  if ((newMax - newMin) < minV)
+    q->length = newMin + RandomNumber(seed) * (newMax - newMin);
+  else if (minProp < maxProp) {
+    /* do a sliding window - Allen */
+    newProp = oldProp + window * (RandomNumber(seed) - 0.5);
+    /* check validity */
+    isValidP = NO;
+    do {
+      if (newProp < minProp) {
+	tmpProp = newProp;
+	newProp = 2*minProp - newProp;
+	if (newProp > maxProp) {
+	  abortMove = YES;
+	  MrBayesPrint("%s   ERROR (Move_NodeSlider2): oldM %1.6lf (%1.6lf %1.6lf) newM %1.6lf oldProp %1.6lf orig newProp %1.6lf newProp %1.6lf minProp %1.6lf maxProp %1.6lf\n",
+		       spacer, oldM, q->length, p->length, newM, oldProp, tmpProp, newProp, minProp, maxProp);
+	  return(NO_ERROR);
+	}
+      } else if (newProp > maxProp) {
+	tmpProp = newProp;
+	newProp = 2*maxProp - newProp;
+	if (newProp < minProp) {
+	  abortMove = YES;
+	  MrBayesPrint("%s   ERROR (Move_NodeSlider2): oldM %1.6lf (%1.6lf %1.6lf) newM %1.6lf oldProp %1.6lf orig newProp %1.6lf newProp %1.6lf minProp %1.6lf maxProp %1.6lf\n",
+		       spacer, oldM, q->length, p->length, newM, oldProp, tmpProp, newProp, minProp, maxProp);
+	  return(NO_ERROR);
+	}
+      } else
+	isValidP = YES;
+    } while (isValidP == NO);
+
+    q->length = newM*newProp;
+    if (q->length < newMin)
+      q->length = newMin;
+    else if (q->length > newMax)
+      q->length = newMax;
+  } else /* just in case is due to roundoff error */
+    q->length = (newMin+newMax)/2;
+
+  p->length = newM - q->length;
+
+  /* set flags for update of transition probabilities at p and q */
+  p->upDateTi = YES;
+  q->upDateTi = YES;
+  p->upDateCl = YES;
+  
+  /* set flags for update of cond likes from p->anc and down to root */
+  while (p->anc != NULL)
+    {
+      p = p->anc;
+      p->upDateCl = YES;
+    }
+  
+  /* calculate proposal ratio */
+  (*lnProposalRatio) = log(newM / oldM) + log ((newMax - newMin) / (oldMax - oldMin));
+  
+  /* update prior if exponential prior on branch lengths */
+  if (param->paramId == BRLENS_EXP)
+    (*lnPriorRatio) = brlensPrExp * (oldM - newM);
+  
+  return (NO_ERROR);
 
+}
 
 
 /*-----------------------------------------------------------------------------------
@@ -24173,7 +24995,7 @@
 					b->isLocked = YES;
 					u->isLocked = NO;
 					b->lockID = u->lockID;
-					u->lockID = 0;
+					u->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			else 
@@ -24220,7 +25042,7 @@
 						u->isLocked = YES;
 						b->isLocked = NO;
 						u->lockID = b->lockID;
-						b->lockID = 0;
+						b->lockID = -1; /* to -1 from 0 - Allen */
 						}
 					}
 				}
@@ -24753,7 +25575,7 @@
 					b->isLocked = YES;
 					u->isLocked = NO;
 					b->lockID = u->lockID;
-					u->lockID = 0;
+					u->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			else 
@@ -24800,7 +25622,7 @@
 						u->isLocked = YES;
 						b->isLocked = NO;
 						u->lockID = b->lockID;
-						b->lockID = 0;
+						b->lockID = -1; /* to -1 from 0 - Allen */
 						}
 					}
 				}
@@ -25362,7 +26184,7 @@
 					b->isLocked = YES;
 					u->isLocked = NO;
 					b->lockID = u->lockID;
-					u->lockID = 0;
+					u->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			else 
@@ -25409,7 +26231,7 @@
 						u->isLocked = YES;
 						b->isLocked = NO;
 						u->lockID = b->lockID;
-						b->lockID = 0;
+						b->lockID = -1; /* to -1 from 0 - Allen */
 						}
 					}
 				}
@@ -26043,7 +26865,7 @@
 					b->isLocked = YES;
 					u->isLocked = NO;
 					b->lockID = u->lockID;
-					u->lockID = 0;
+					u->lockID = -1; /* to -1 from 0 - Allen */
 					}
 				}
 			else 
@@ -26090,7 +26912,7 @@
 						u->isLocked = YES;
 						b->isLocked = NO;
 						u->lockID = b->lockID;
-						b->lockID = 0;
+						b->lockID = -1; /* to -1 from 0 - Allen */
 						}
 					}
 				}
@@ -27480,16 +28302,16 @@
 
 	/* change pi */
 	int			i, nStates;
-	MrBFlt		dirichletParameters[64], *newPi, *oldPi, symDirAlphai, sum, alphaPi, x, y;
+	MrBFlt		dirichletParameters[64], *newPi, *oldPi, *origDirAlpha, sum, alphaPi, x, y;
 
 	/* get the values we need */
 	nStates = param->nSubValues;
-	symDirAlphai = *GetParamVals(param, chain, state[chain]);
+	origDirAlpha = GetParamVals(param, chain, state[chain]);
 	newPi = GetParamSubVals (param, chain, state[chain]);
 	oldPi = GetParamSubVals (param, chain, state[chain] ^ 1);
 	
 #	if defined ASYMMETRY
-	nStates /= 2;
+	nStates /= 2; /* Shouldn't this be only for the RESTRICTION datatype? - Allen */
 	if (RandomNumber (seed) < 0.5)
 		{
 		newPi += 2;
@@ -27499,9 +28321,18 @@
 	
 	alphaPi = mvp[0];
 
-	/* multiply old values with some large number to get new values close to the old ones */
+	sum = 0.0;
 	for (i=0; i<nStates; i++)
+		sum += origDirAlpha[i];
+	sum /= nStates; /* so for default of equal to # of states will be 1 - Allen */
+
+	/* multiply old values with some large number to get new values close to the old ones */
+	for (i=0; i<nStates; i++) {
+		if ((sum != 1) && (AreDoublesEqual(oldPi[i],(origDirAlpha[i]/(sum*nStates)),0.0001) != NO))
+			dirichletParameters[i] = oldPi[i] * alphaPi * sum;
+		else
 		dirichletParameters[i] = oldPi[i] * alphaPi;
+	}
 
 	DirichletRandomVariable (dirichletParameters, newPi, nStates, seed);
 
@@ -27522,24 +28353,46 @@
 	x = LnGamma(sum);
 	for (i=0; i<nStates; i++)
 		x -= LnGamma(newPi[i]*alphaPi);
-	for (i=0; i<nStates; i++)
+	for (i=0; i<nStates; i++) {
+		if (oldPi[i] >= 0.0001) /* Allen */
 		x += (newPi[i]*alphaPi-1.0)*log(oldPi[i]);
+		else
+			x += (newPi[i]*alphaPi-1.0)*log(0.0001);
+	}
 	sum = 0.0;
-	for (i=0; i<nStates; i++)
+	for (i=0; i<nStates; i++) {
+		if (oldPi[i] >= 0.0001)
 		sum += oldPi[i]*alphaPi;
+		else
+			sum += 0.0001*alphaPi;
+	}
 	y = LnGamma(sum);
-	for (i=0; i<nStates; i++)
+	for (i=0; i<nStates; i++) {
+		if (oldPi[i] >= 0.0001) {
 		y -= LnGamma(oldPi[i]*alphaPi);
-	for (i=0; i<nStates; i++)
 		y += (oldPi[i]*alphaPi-1.0)*log(newPi[i]);
+		} else {
+			y -= LnGamma(0.0001*alphaPi);
+			y += (0.0001*alphaPi-1.0)*log(newPi[i]);
+		}
+	}
 	(*lnProposalRatio) = x - y;
 
-	/* get prior ratio */
-	y = x = 0.0;					/* the Gamma part of the prior is the same */
+	/* get prior ratio */ /* The below is my best guess of how it was supposed to be done; it appears
+			         that the old version was copied from Move_StatefreqsSymDirMultistate
+				 and not changed properly; should really try to figure out how to
+				 get it to matter more depending on origDirAlpha! XXX - Allen */
+	y = x = sum = 0.0;				/* the Gamma part of the prior is the same */
 	for (i=0; i<nStates; i++)
-		x += (symDirAlphai-1.0)*log(newPi[i]);
+		sum += origDirAlpha[i];
 	for (i=0; i<nStates; i++)
-		y += (symDirAlphai-1.0)*log(oldPi[i]);
+		x += ((origDirAlpha[i]/sum)-1)*log(newPi[i]);
+	for (i=0; i<nStates; i++) {
+		if (oldPi[i] >= 0.0001)
+			y += ((origDirAlpha[i]/sum)-1)*log(oldPi[i]);
+		else
+			y += ((origDirAlpha[i]/sum)-1)*log(0.0001);
+	}
 	(*lnPriorRatio) = x - y;
 		
 	for (i=0; i<param->nRelParts; i++)
@@ -27660,7 +28513,7 @@
 	/* change switch rate of covarion model using sliding window */
 	
 	int			i, isSPriorExp, isValidS, whichRate;
-	MrBFlt		oldS, newS, window, minS, maxS, sExp=0.0, ran, *value;
+	MrBFlt		oldS, newS, window, minS, maxS, sExp=0.0, ran, *value, tmpS;
 	ModelParams *mp;
 
 	/* decide which switching rate to change */
@@ -27675,15 +28528,65 @@
 	/* get model params */
 	mp = &modelParams[param->relParts[0]];
 	
+	/* get old values of switching rates */
+	value = GetParamVals(param, chain, state[chain]);
+
+
 	/* get minimum and maximum values for switching rate */
-	if (param->paramId == SWITCH_UNI)
-		{
+	if (param->paramId == SWITCH_UNI) {
 		minS = mp->covswitchUni[0];
 		maxS = mp->covswitchUni[1];
 		isSPriorExp = NO;
+        } else if (param->paramId == SWITCH_PINVUNI) { /* - Allen */
+	  minS = mp->covswitchUni[0];
+	  maxS = mp->covswitchUni[1];
+	  isSPriorExp = NO;
+	  if (whichRate == 1) {
+	    /* maximum for s10 that will give mp->covswitchMaxPinvar stationary ProbOff */
+	    tmpS = (value[0]*mp->covswitchMaxPinvar)/(1-mp->covswitchMaxPinvar);
+	    if (maxS > tmpS) {
+	      if (tmpS < minS) { /* Shouldn't happen! */
+	        MrBayesPrint ("%s   ERROR: tmpS %1.2lf minS %1.2lf s01 %1.2lf\n", spacer, tmpS, minS,
+			      value[0]);
+		return (ERROR);
+	      } else if (tmpS == minS) {
+		whichRate = 0;
+		tmpS = (value[1]*(1-mp->covswitchMaxPinvar))/mp->covswitchMaxPinvar;
+		if (minS < tmpS) {
+		  if (tmpS >= maxS) {
+		    MrBayesPrint ("%s   ERROR: tmpS %1.2lf maxS %1.2lf s01 %1.2lf s10 %1.2lf\n", spacer, tmpS, maxS,
+				  value[0],value[1]);
+		    return (ERROR);
+		  } else
+		    minS = tmpS;
 		}
-	else
-		{
+	      } else
+		maxS = tmpS;
+            }
+	  } else {
+	    /* minimum for s01 that will give mp->covswitchMaxPinvar stationary ProbOff */
+	    tmpS = (value[1]*(1-mp->covswitchMaxPinvar))/mp->covswitchMaxPinvar;
+	    if (minS < tmpS) {
+	      if (tmpS > maxS) {
+		MrBayesPrint ("%s   ERROR: tmpS %1.2lf maxS %1.2lf s10 %1.2lf\n", spacer, tmpS, maxS,
+			      value[1]);
+		return (ERROR);
+	      } else if (tmpS == maxS) {
+		whichRate = 1;
+		tmpS = (value[0]*mp->covswitchMaxPinvar)/(1-mp->covswitchMaxPinvar);
+		if (maxS > tmpS) {
+		  if (tmpS <= minS) { /* Shouldn't happen! */
+		    MrBayesPrint ("%s   ERROR: tmpS %1.2lf minS %1.2lf s01 %1.2lf s10 %1.2lf\n", spacer, tmpS, minS,
+				  value[0],value[1]);
+		    return (ERROR);
+		  } else
+		    maxS = tmpS;
+		}
+	      } else
+	        minS = tmpS;
+            }
+          }
+	} else {
 		minS = 0.01;
 		maxS = KAPPA_MAX;
 		sExp = mp->covswitchExp;
@@ -27691,7 +28594,6 @@
 		}
 
 	/* get old value of switching rate */
-	value = GetParamVals(param, chain, state[chain]);
 	newS = oldS = value[whichRate];
 
 	/* change value for switching rate */
@@ -29277,7 +30179,8 @@
 {
 
 	int				c, i, k, nStates;
-	MrBFlt			*bs, freq;
+	MrBFlt			*bs, freq, covBF[40];
+	MrBFlt                  s01=50.0, s10=50.0, probOn, probOff, *swr;
 	CLFlt			*clFP, *cL, sum;
 	char			*tempStr;
 	int             tempStrSize = 200;
@@ -29299,6 +30202,29 @@
 	/* find number of states */
 	nStates = m->numModelStates;
 	
+	/* if covarion model, adjust base frequencies */
+        if (m->switchRates != NULL) { /* code from Likelihood_Gen - Allen */
+	  /* find the stationary frequencies */
+	  swr = GetParamVals(m->switchRates, chain, state[chain]);
+	  s01 = swr[0];
+	  s10 = swr[1];
+	  probOn = s01 / (s01 + s10);
+	  probOff =  1.0 - probOn;
+	  
+	  /* now adjust the base frequencies; on-state stored first in cond likes */
+	  for (i=0; i<(nStates/2); i++) {
+	    covBF[i] = bs[i] * probOn;
+	    /* try to reduce numerical inaccuracy - Allen */
+	    covBF[i+(nStates/2)] = bs[i] - covBF[i];
+	    if (covBF[i+(nStates/2)] < LIKE_EPSILON)
+	      covBF[i+(nStates/2)] = bs[i] * probOff;
+	  }
+	  
+	  /* finally set bs pointer to adjusted values */
+	  bs = covBF;
+	}
+
+
 	/* find frequencies of rate categories */
 	freq =  1.0 /  m->numGammaCats;
 	
@@ -29334,18 +30260,23 @@
 		}
 
 	/* print the resulting conditional likelihoods cycling over uncompressed chars */
-	for (c=0; c<numChar; c++)
-		{
+	for (c=0; c<numChar; c++) {
 		if (charInfo[c].isExcluded == YES || charInfo[c].partitionId[partitionNum-1] != division+1)
 			continue;
 		i = compCharPos[c] - m->compCharStart;
 		cL = ancStateCondLikes + (i*nStates);
-		for (i=0; i<nStates; i++)
-			{
+		if (m->switchRates != NULL) {
+		  for (i=0; i<nStates/2; i++) { /* may later work on indicating covarion on/off... XXX - Allen */
+		    SaveSprintf (&tempStr, &tempStrSize, "%f\t", (cL[i]+cL[i+(nStates/2)]));
+		    if (AddToPrintString (tempStr) == ERROR) return (ERROR);
+		  }
+		} else {
+		  for (i=0; i<nStates; i++) {
 			SaveSprintf (&tempStr, &tempStrSize, "%f\t", cL[i]);
 			if (AddToPrintString (tempStr) == ERROR) return (ERROR);
 			}
 		}
+	}
 	free (tempStr);
 	return NO_ERROR;	
 
@@ -29625,6 +30556,11 @@
 	Tree			*t;
 	TreeNode		*p;
 
+#	if defined (MPI_ENABLED)
+	if (proc_id != 0)
+		return (NO_ERROR);
+#	endif
+
 	if (!chainCondLikes || !termCondLikes)
 		return ERROR;
 
@@ -29686,10 +30622,11 @@
 				MrBayesPrint ("\n");
 				}
 			MrBayesPrint("\n");
+#	if (!defined (MPI_ENABLED))
 			if ((c-m->condLikeStart)% 8 == 0)
 				getchar();
+#endif
 			}
-				
 		}	/* next division */
 
 	return NO_ERROR;
@@ -29722,6 +30659,11 @@
 	extern char		WhichStand (int x);
 
 
+#	if defined (MPI_ENABLED)
+	if (proc_id != 0)
+		return (NO_ERROR);
+#	endif
+
 	if (!compMatrix)
 		return ERROR;
 
@@ -29842,8 +30784,10 @@
 				}
 			MrBayesPrint ("\n\n");
 			}
+#if (! defined(MPI_ENABLED))
 		MrBayesPrint ("Press return to continue\n");
 		getchar();
+#endif
 		}	/* next division */
 
 	return NO_ERROR;
@@ -29872,6 +30816,10 @@
 	extern char		WhichRes (int x);
 	extern char		WhichStand (int x);
 
+#	if defined (MPI_ENABLED)
+	if (proc_id != 0)
+		return (NO_ERROR);
+#	endif
 
 	if (!matrix)
 		return ERROR;
@@ -30065,8 +31013,10 @@
 			{
 			for (i=0; i<chainParams.numChains; i++)
 				{
-				for (j=i+1; j<chainParams.numChains; j++)
-					{
+				  for (j=i+1; j<chainParams.numChains; j++) {
+				    if (swapInfo[n][j][i] < 1)
+				      MrBayesPrintf (fpMcmc, "\tN/A");
+				    else
 					MrBayesPrintf (fpMcmc, "\t%.6f", (MrBFlt) swapInfo[n][i][j] / (MrBFlt) swapInfo[n][j][i]);
 					}
 				}
@@ -30241,7 +31191,7 @@
 	MrBFlt			freq, siteLike, invLike, catLike, pInvar=0.0, *bs,
 					*clInvar=NULL, *catRate, baseRate;
 	MrBFlt			s01, s10, probOn, probOff, *swr, covBF[40];
-	CLFlt			*lnScaler, *clP, *siteRates;
+	CLFlt			*clP, *siteRates, *lnScaler;
 	char			*tempStr;
 	int             tempStrSize = 200;
 	ModelInfo		*m;
@@ -31438,8 +32388,8 @@
 					maxNumExchanges = swapInfo[n][i][j];
 		SaveSprintf(&tempStr, &tempStrSize, "%d", maxNumExchanges);
 		maxLen = (int) strlen(tempStr);
-		if (maxLen < 4)
-			maxLen = 4;
+		if (maxLen < 8) /* from 4 - Allen */
+			maxLen = 8;
 			
 		reweightingChars = NO;
 		if ((chainParams.weightScheme[0] + chainParams.weightScheme[1]) > 0.00001)
@@ -31480,10 +32430,10 @@
 						MrBayesPrint ("%*c%s ", maxLen-3, ' ', " N/A");
 						}
 					else
-						{
-						SaveSprintf(&tempStr, &tempStrSize, "%1.2lf", (MrBFlt)swapInfo[n][i][j]/swapInfo[n][j][i]);
+					  { /* from %1.2lf - Allen */
+						SaveSprintf(&tempStr, &tempStrSize, "%1.6lf", (MrBFlt)swapInfo[n][i][j]/swapInfo[n][j][i]);
 						len = (int) strlen(tempStr);
-						MrBayesPrint ("%*c%1.2lf ", maxLen-len+1, ' ', (MrBFlt)swapInfo[n][i][j]/swapInfo[n][j][i]);
+						MrBayesPrint ("%*c%1.6lf ", maxLen-len+1, ' ', (MrBFlt)swapInfo[n][i][j]/swapInfo[n][j][i]);
 						}
 					}
 				else if (i == j)
@@ -31629,7 +32579,7 @@
 ---------------------------------------------------*/
 void PrintTiProbs (MrBFlt *tP, MrBFlt *bs, int nStates)
 
-{
+{ /* Why is it MrBFlt not CLFlt? - Allen */
 	int		i, j;
 	MrBFlt	*tiP, sum;
 
@@ -31776,12 +32726,22 @@
 }
 
 
-void PrintToScreen (int curGen, time_t endingT, time_t startingT)
+void PrintToScreen (int curGen, int movesAccepted, time_t endingT, time_t startingT)
 {
 
 	int			i, chn, nHours, nMins, nSecs;
 	MrBFlt		timePerGen;
-
+#if defined(UNIX_VERSION)
+	time_t endtime_secs; /* UNIX_VERSION b/c uncertain on time handling otherwise - Allen */
+	struct tm *tmbuf;
+	static char *dayname[] = {"Sun", "Mon", "Tue", "Wed", "Thu", "Fri",
+				  "Sat"};
+	static char *monname[] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
+				  "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
+	/* My apologies about the fixed English-language names,
+	   but asctime's inbuilt return was problematic (I use Perl for
+	   string-handling programs...) - Allen */
+#endif
 #	if defined (MPI_ENABLED)
 	int			numLocalColdChains, numFirstAndLastCold;
 	
@@ -31790,7 +32750,7 @@
 		MrBayesPrint ("\n");
 		MrBayesPrint ("%s   Chain results:\n\n", spacer);
 		}
-	MrBayesPrint ("%s   %4d -- ", spacer, curGen);
+	MrBayesPrint ("%s   %5d -- ", spacer, curGen);
 	numLocalColdChains = numFirstAndLastCold = 0;
 	for (chn=0; chn<numLocalChains; chn++)
 		{
@@ -31841,8 +32801,7 @@
 	else
 		MrBayesPrint ("[...%d remote chains...] ", (chainParams.numChains*chainParams.numRuns) - numLocalChains);
 
-	if (curGen > 1)
-		{
+	if (curGen > 1) {
 		timePerGen = (MrBFlt) ((MrBFlt)(endingT-startingT)/(MrBFlt)curGen);
 		nSecs = (int)((chainParams.numGen - curGen) * timePerGen);
 		nHours = (int)nSecs / 3600;
@@ -31850,6 +32809,20 @@
 		nMins = nSecs / 60; 
 		nSecs -= nMins * 60;
 		MrBayesPrint ("-- %d:%0.2d:%0.2d", nHours, nMins, nSecs);
+#if defined(UNIX_VERSION)
+		if ((nHours > 0) || (nMins > 0)) {
+		  endtime_secs = time(0) + (time_t)((chainParams.numGen - curGen) * timePerGen);
+		  tmbuf = localtime(&endtime_secs); /* Code from Perl - Allen */
+		  MrBayesPrint(" -- Estimated Endtime: %s %s %2d %02d:%02d:%02d %d",
+			       dayname[tmbuf->tm_wday],
+			       monname[tmbuf->tm_mon],
+			       tmbuf->tm_mday,
+			       tmbuf->tm_hour,
+			       tmbuf->tm_min,
+			       tmbuf->tm_sec,
+			       tmbuf->tm_year + 1900);
+		}
+#endif
 		}
 	MrBayesPrint ("\n");
 	fflush (stdout);
@@ -31861,7 +32834,8 @@
 		MrBayesPrint ("\n");
 		MrBayesPrint ("%s   Chain results:\n\n", spacer);
 		}
-	MrBayesPrint ("%s   %5d -- ", spacer, curGen);
+	MrBayesPrint ("%s   %5d (%5d; %6.4lf %%) -- ", spacer, curGen, movesAccepted,
+		      (double)(((double)movesAccepted*100.0)/(double)curGen));
 	if (numLocalChains == 1)
 		MrBayesPrint ("%1.3lf ", curLnL[0]);
 	else
@@ -31900,8 +32874,7 @@
 			}
 		}
 		
-	if (curGen > 1)
-		{
+	if (curGen > 1) {
 		timePerGen = (MrBFlt) ((MrBFlt)(endingT-startingT)/(MrBFlt)curGen);
 		nSecs = (int)((chainParams.numGen - curGen) * timePerGen);
 		nHours = (int)nSecs / 3600;
@@ -31909,6 +32882,20 @@
 		nMins = nSecs / 60; 
 		nSecs -= nMins * 60;
 		MrBayesPrint ("-- %d:%0.2d:%0.2d", nHours, nMins, nSecs);
+#if defined(UNIX_VERSION)
+		if ((nHours > 0) || (nMins > 0)) {
+		  endtime_secs = time(0) + (time_t)((chainParams.numGen - curGen) * timePerGen);
+		  tmbuf = localtime(&endtime_secs); /* Code from Perl - Allen */
+		  MrBayesPrint(" -- Estimated Endtime: %s %s %2d %02d:%02d:%02d %d",
+			       dayname[tmbuf->tm_wday],
+			       monname[tmbuf->tm_mon],
+			       tmbuf->tm_mday,
+			       tmbuf->tm_hour,
+			       tmbuf->tm_min,
+			       tmbuf->tm_sec,
+			       tmbuf->tm_year + 1900);
+		}
+#endif
 		}
 	MrBayesPrint ("\n");
 	
@@ -32960,7 +33947,7 @@
 
 
 
-
+#if !defined(ALTSCALING)
 /*-------------------------------------------------------------------
 |
 |	ResetScalers: reset scaler nodes of all trees of all chains
@@ -32983,11 +33970,12 @@
 			t = GetTreeFromIndex (i, chn, state[chn]);
 		
 			/* set the node depth value of terminal nodes to zero */
-			for (n=0; n<t->nNodes; n++)
-				{
+			for (n=0; n<t->nNodes; n++) {
 				p = t->allDownPass[n];
-				if (p->left == NULL)
+			  if (p->left == NULL) {
 					p->x = 0;
+			    p->scalerNode = NO; /* just in case - Allen */
+			  }
 				}
 
 			/* loop over interior nodes */
@@ -33010,8 +33998,7 @@
 
 	return NO_ERROR;
 }
-
-
+#endif
 
 
 
@@ -33084,8 +34071,15 @@
 {
 	
 	int		i, j, k, n, chn, swapA=0, swapB=0, whichMove, acceptMove, lastDiagnostics, stopChain,
-			nErrors;
+			nErrors, burninGen=250000, partBurninGen1, partBurninGen2, partBurninGen=125000, tempZero, tempUse,
+	                tempUse1, tempUse2, adaptNarrow=NO, adaptWiden=NO, timesAdaptNarrow=0,timesAdaptWiden=0,
+	                movesAccepted=0, timesAdaptNarrowAccepted=0, timesAdaptNarrowRefused=0, adaptOnlyIfBack=NO,
+	                timesAdaptWidenAccepted=0, timesAdaptWidenRefused=0, adaptDidNarrow=NO, adaptDidWiden=NO,
+	                didReportPartBurnin=NO, didReportBurnin=NO;
 	MrBFlt		r=0.0, lnLikelihoodRatio, lnPriorRatio, lnProposalRatio, lnLike=0.0, lnPrior=0.0, f, CPUTime;
+	MrBFlt          mult_short=138.463, mult_long=69.2314, tempSubZero, effectiveGen, successRatio1, minParam;
+	MrBFlt          maxParam, adaptMult, adaptCFactor, newParam, currParam, minParamMult, maxParamMult, newParamMult;
+	MrBFlt          minAdaptCFactor, successRatio2, successRatioUse, origAdaptCFactor;
 	MCMCMove	*theMove, *mv;
 	time_t		startingT, endingT, stoppingT1, stoppingT2;
 	clock_t		previousCPUTime, currentCPUTime;
@@ -33368,12 +34362,77 @@
 #		endif
 		}
 
+	/* Initialize generation cutoffs, mults for Simulated Annealing, Adapt - Allen */
+	if (chainParams.useAdapt == YES) { /* setting these here to eliminate complaints from GCC */
+	  /* "higher" temps usable */
+	  tempZero = (int) ((1.0/(1.0+0.6))*1000.0);
+	  tempSubZero = 1.0 - (1.0/(1.0+0.6));
+	} else {
+	  /* 1/(1+0.2) as starting temperature (0.2 is default Temp for 1st of heated multiple chains) */
+	  tempZero = (int) ((1.0/(1.0+0.2))*1000.0);
+	  tempSubZero = 1.0 - (1.0/(1.0+0.2));
+	}
+	if ((chainParams.useSA == YES) ||
+	    (chainParams.useAdapt == YES)) { /* note that C casts from floating to integer truncate */
+	  if (chainParams.relativeBurnin == YES) {
+	    burninGen = (int) (chainParams.numGen*chainParams.burninFraction);
+	    if (chainParams.burninFraction < 0.25) {
+	      MrBayesPrint ("%s   burninGen is %d (burninFraction is %1.3lf); cannot do Simulated Annealing or Adapt\n",
+			    spacer, burninGen, chainParams.burninFraction);
+	      return ERROR;
+	    }
+	    partBurninGen1 = (int) (chainParams.numGen*chainParams.burninFraction*0.5);
+	    partBurninGen2 = (int) (chainParams.numGen*(chainParams.burninFraction-0.25));
+	  } else {
+	    burninGen = chainParams.chainBurnIn;
+	    if ((burninGen*4) < chainParams.numGen) {
+	      MrBayesPrint ("%s   burninGen is %d but numGen is %d; cannot do Simulated Annealing or Adapt\n",
+			    spacer, burninGen, chainParams.numGen);
+	      return ERROR;
+	    }
+	    partBurninGen1 = (int) (chainParams.chainBurnIn/2);
+	    partBurninGen2 = (int) (chainParams.chainBurnIn - (chainParams.numGen*0.25));
+	  }
+
+	  partBurninGen = partBurninGen1;
+	  if (partBurninGen2 > partBurninGen1)
+	    partBurninGen = partBurninGen2;
+
+	  if (burninGen <= partBurninGen) {
+	    MrBayesPrint ("%s   partBurninGen would be %d (burninGen is %d); cannot do Simulated Annealing or Adapt\n",
+			  spacer, partBurninGen, burninGen);
+	    return ERROR;
+	  }
+
+	  MrBayesPrint("%s   partBurninGen is %d, burninGen is %d\n",
+		       spacer, partBurninGen, burninGen);
+
+	  if (chainParams.useSA == YES) {
+	    if (partBurninGen <= 2) {
+	      MrBayesPrint ("%s   partBurninGen would be %d (burninGen is %d); cannot do Simulated Annealing\n",
+			    spacer, partBurninGen, burninGen);
+	      return ERROR;
+	    }
+
+	    MrBayesPrint("%s   tempSubZero is %.6lf, tempZero/1000 is %.3lf\n",
+			 spacer, tempSubZero, (double) ((double)tempZero/1000.0));
+
+	    /* such that at partBurninGen, will be 0 (mult 1) w/rounding, but starts at tempZero */
+	    mult_short = (0.001-tempSubZero)/(partBurninGen-1);
+	    /* ditto for burninGen */
+	    mult_long = (0.001-tempSubZero)/(burninGen-1);
+	    MrBayesPrint ("%s   mult_short is %.6lf, mult_long is %.6lf\n",
+			  spacer,mult_short,mult_long);
+	  }
+	}
 	/* get chain IDs */
 	GetChainIds ();
 	
 	/* initialize likelihoods and prior                  */
 	/* touch everything and calculate initial cond likes */
+#if !defined(ALTSCALING)
 	ResetScalers ();
+#endif
 	TouchAllPartitions ();
 	for (chn=0; chn<numLocalChains; chn++)
 		{
@@ -33413,6 +34472,10 @@
 	for (i=0; i<chainParams.numRuns; i++)
 		maxLnL0[i] = -100000000.0;
 
+#if defined(MPI_ENABLED) && defined(DEBUG_MPI) && LAM_MPI
+	MPIL_Trace_on();
+#endif
+
 	startingT=time(0);
 	CPUTime = 0.0;
 	previousCPUTime = clock();
@@ -33421,11 +34484,27 @@
 		{
 
 		currentCPUTime = clock();
-		if (currentCPUTime - previousCPUTime > 10 * CLOCKS_PER_SEC)
+		if (currentCPUTime < previousCPUTime) { /* wraparound - Allen */
+		  if ((currentCPUTime < 0) && (previousCPUTime > 0)) { /* Approximating */
+		    if (previousCPUTime > (-1*currentCPUTime)) {
+		      CPUTime += (2*(previousCPUTime + currentCPUTime)) / (MrBFlt) (CLOCKS_PER_SEC);
+		    } else {
+		      CPUTime += (-2*(previousCPUTime + currentCPUTime)) / (MrBFlt) (CLOCKS_PER_SEC);
+		    }
+		    previousCPUTime = currentCPUTime;
+		  } else {
+		    if ((previousCPUTime - currentCPUTime) >= CLOCKS_PER_SEC) {
+		      CPUTime += (previousCPUTime - currentCPUTime) / (MrBFlt) (CLOCKS_PER_SEC);
+		      previousCPUTime = currentCPUTime;
+		    }
+		  }
+		} else {
+		  if ((currentCPUTime - previousCPUTime) > 10 * CLOCKS_PER_SEC)
 			{
 			CPUTime += (currentCPUTime - previousCPUTime) / (MrBFlt) (CLOCKS_PER_SEC);
 			previousCPUTime = currentCPUTime;
 			}
+		}
 
 		/*! AbortRun is set by the signal handler when it receives a CTRL-C */
 #if 		defined  (MPI_ENABLED)
@@ -33435,11 +34514,14 @@
 		if (confirmAbortRun == YES && RequestAbortRun() == 1)
 			return ABORT;
 #endif
+
+#if !defined(ALTSCALING)
 		/* Refresh scalers every SCALER_REFRESH_FREQ generations.                */
 		/* It is done before copying so we know it will take effect immediately. */
 		/* However, the actual scalers are recalculated only when really needed. */
 		if (n % SCALER_REFRESH_FREQ == 0)
 			ResetScalers();
+#endif
 
 		for (chn=0; chn<numLocalChains; chn++)
 			{
@@ -33474,6 +34558,387 @@
 			for (i=0; i<theMove->parm->nRelParts; i++)
 				modelSettings[theMove->parm->relParts[i]].upDateCl = YES;
 
+			i = chainId[chn];
+
+			adaptNarrow = NO;
+			adaptWiden = NO;
+			adaptDidNarrow = NO;
+			adaptDidWiden = NO;
+			adaptOnlyIfBack = NO;
+			successRatioUse = 0.4; /* silence GCC */
+
+			if ((chainParams.useAdapt == YES) && (theMove->doAdaptParam == YES) && /* - Allen */
+			    (theMove->nTried[0] >= 50) &&
+			    (theMove->nTriedAfterAdapt >= 20)) {
+			  j = theMove->whichAdaptParam;
+			  minParam = theMove->minAdapt;
+			  maxParam = theMove->maxAdapt;
+			  successRatio1 = ((MrBFlt)(theMove->nAccepted[0]))/((MrBFlt)(theMove->nTried[0]));
+			  successRatio2 = ((MrBFlt)(theMove->nAcceptedAfterAdapt))/
+			    ((MrBFlt)(theMove->nTriedAfterAdapt)); /* just recent */
+			  if ((successRatio1 < 0.7) &&
+			      (successRatio2 < 0.7) &&
+			      (theMove->tuningParam[j] > minParam)) {
+			    if ((0.1 - successRatio1) > LIKE_EPSILON) {
+			      if (((0.1 - successRatio2) > LIKE_EPSILON) &&
+				  (theMove->nTriedAfterAdapt >= 50)) {
+				adaptNarrow = YES;
+				successRatioUse = MinimumValue(successRatio1,successRatio2);
+			      } else {
+				adaptNarrow = YES;
+				adaptOnlyIfBack = YES;
+				if (n >= burninGen)
+				  successRatioUse = successRatio1;
+				else
+				  successRatioUse = MinimumValue(successRatio1,successRatio2);
+			      }
+			    } else if (((0.1 - successRatio2) > LIKE_EPSILON) &&
+				       (theMove->nTriedAfterAdapt >= 50)) {
+			      adaptNarrow = YES;
+			      adaptOnlyIfBack = YES;
+			      if (n >= burninGen)
+				successRatioUse = successRatio2;
+			      else
+				successRatioUse = MinimumValue(successRatio1,successRatio2);
+			    }
+			  } else if ((successRatio1 > 0.1) &&
+				     (successRatio2 > 0.1) &&
+				     (theMove->tuningParam[j] < maxParam)) {
+			    if ((successRatio1 - 0.7) > LIKE_EPSILON) {
+			      if (((successRatio2 - 0.7) > LIKE_EPSILON) &&
+				  (theMove->nTriedAfterAdapt >= 50)) {
+				adaptWiden = YES;
+				successRatioUse = MaximumValue(successRatio1,successRatio2);
+			      } else {
+				adaptWiden = YES;
+				adaptOnlyIfBack = YES;
+				if (n >= burninGen)
+				  successRatioUse = successRatio1;
+				else
+				  successRatioUse = MaximumValue(successRatio1,successRatio2);
+			      }
+			    } else if (((successRatio2 - 0.7) > LIKE_EPSILON) &&
+				       (theMove->nTriedAfterAdapt >= 50)) {
+			      adaptWiden = YES;
+			      adaptOnlyIfBack = YES;
+			      if (n >= burninGen)
+				successRatioUse = successRatio2;
+			      else
+				successRatioUse = MaximumValue(successRatio1,successRatio2);
+			    }
+			  }
+
+			  /* do not alter after partBurninGen unless moving back toward original */
+			  if (n >= partBurninGen)
+			    adaptOnlyIfBack = YES;
+									  
+			  if (adaptOnlyIfBack == YES) {
+			    if (adaptNarrow == YES) {
+			      minParam = theMove->oldTuningParam;
+			      if (theMove->tuningParam[j] <= minParam)
+				adaptNarrow = NO;
+			    } else if (adaptWiden == YES) {
+			      maxParam = theMove->oldTuningParam;
+			      if (theMove->tuningParam[j] >= maxParam)
+				adaptWiden = NO;
+			    }
+			  }
+
+			  if ((theMove->adaptMultiplier == YES) &&
+			      ((adaptNarrow == YES) || (adaptWiden == YES))) {
+			    currParam = exp(theMove->tuningParam[j]/2);
+			    minParamMult = theMove->minAdaptMult;
+			    maxParamMult = theMove->maxAdaptMult;
+			    if (adaptOnlyIfBack == YES) {
+			      if (adaptNarrow == YES)
+				minParamMult = theMove->oldTuningParamMult;
+			      else
+				maxParamMult = theMove->oldTuningParamMult;
+			    }
+
+			    if ((adaptNarrow == YES) && (currParam > minParamMult)) {
+			      if ((currParam > theMove->oldTuningParamMult) && (n < burninGen)) {
+				/* would get to min in 2 moves of this level - slows down as get closer to min */
+				adaptCFactor = sqrt(currParam/theMove->minAdaptMult)-1.0;
+				origAdaptCFactor = CBRT(currParam/theMove->minAdaptMult)-1.0;
+			      } else {
+				/* would get to min in 3 moves of this level - slows down as get closer to min */
+				origAdaptCFactor = adaptCFactor = CBRT(currParam/theMove->minAdaptMult)-1.0;
+			      }
+			      if (currParam > theMove->oldTuningParamMult) {
+				adaptCFactor = MaximumValue(adaptCFactor,((currParam/theMove->oldTuningParamMult)-1.0));
+				origAdaptCFactor =
+				  MaximumValue(origAdaptCFactor,((currParam/theMove->oldTuningParamMult)-1.0));
+			      }
+			      minAdaptCFactor = MinimumValue(theMove->minAdaptCFactor,((currParam/minParamMult)-1.0));
+			      if (adaptCFactor < minAdaptCFactor)
+				adaptCFactor = minAdaptCFactor;
+			      if (origAdaptCFactor < minAdaptCFactor)
+				origAdaptCFactor = minAdaptCFactor;
+			      adaptMult = 1.0/(1.0 + adaptCFactor*((0.1 - successRatioUse)/0.1));
+			      newParamMult = adaptMult*currParam;
+			      if (newParamMult < currParam) { /* roundoff error is not improbable... */
+				if ((origAdaptCFactor < adaptCFactor) && (newParamMult < theMove->oldTuningParamMult)) {
+				  adaptMult = 1.0/(1.0 + origAdaptCFactor*((0.1 - successRatioUse)/0.1));
+				  newParamMult = MinimumValue(theMove->oldTuningParamMult,(adaptMult*currParam));
+				}
+				if (newParamMult < minParamMult)
+				  newParamMult = minParamMult;
+				newParam = 2.0*log(newParamMult);
+				if (newParam < minParam)
+				  newParam = minParam;
+				if (newParam < theMove->tuningParam[j]) {
+				  theMove->tuningParam[j] = newParam;
+				  timesAdaptNarrow++;
+				  adaptDidNarrow = YES;
+				  theMove->nTriedAfterAdapt = 0;
+				  theMove->nAcceptedAfterAdapt = 0;
+				} else if (theMove->tuningParam[j] > maxParam) {
+				  theMove->tuningParam[j] = maxParam;
+				  timesAdaptNarrow++;
+				  adaptDidNarrow = YES;
+				  theMove->nTriedAfterAdapt = 0;
+				  theMove->nAcceptedAfterAdapt = 0;
+				} else if (newParam > theMove->maxAdapt) {
+				  MrBayesPrint ("%s   Error (%2.15lf over maxAdapt %2.15lf; Mult; curr %2.15lf (%2.15lf); min %2.15lf; max %2.15lf) in Adapt for %s\n",
+						spacer, newParam, theMove->maxAdapt, theMove->tuningParam[j], 
+						currParam, minParam, maxParam, theMove->name);
+				  newParam = sqrt(minParam*theMove->tuningParam[j]);
+				  if (newParam > minParam)
+				    newParam = minParam;
+				  if (newParam < theMove->tuningParam[j]) {
+				    theMove->tuningParam[j] = newParam;
+				    timesAdaptNarrow++;
+				    adaptDidNarrow = YES;
+				    theMove->nTriedAfterAdapt = 0;
+				    theMove->nAcceptedAfterAdapt = 0;
+				  }
+				}
+			      } else if (theMove->tuningParam[j] > maxParam) {
+				theMove->tuningParam[j] = maxParam;
+				timesAdaptNarrow++;
+				adaptDidNarrow = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      }
+			    } else if ((adaptWiden == YES) && (currParam < maxParamMult)) {
+			      if ((currParam < theMove->oldTuningParamMult) && (n < burninGen)) {
+				adaptCFactor = sqrt(theMove->maxAdaptMult/currParam)-1.0;
+				origAdaptCFactor = CBRT(theMove->maxAdaptMult/currParam)-1.0;
+			      } else
+				origAdaptCFactor = adaptCFactor = CBRT(theMove->maxAdaptMult/currParam)-1.0;
+			      if (currParam < theMove->oldTuningParamMult) {
+				adaptCFactor = MaximumValue(adaptCFactor,((theMove->oldTuningParamMult/currParam)-1.0));
+				origAdaptCFactor =
+				  MaximumValue(origAdaptCFactor,((theMove->oldTuningParamMult/currParam)-1.0));
+			      }
+			      minAdaptCFactor = MinimumValue(theMove->minAdaptCFactor,((maxParamMult/currParam)-1.0));
+			      if (adaptCFactor < minAdaptCFactor)
+				adaptCFactor = minAdaptCFactor;
+			      if (origAdaptCFactor < minAdaptCFactor)
+				origAdaptCFactor = minAdaptCFactor;
+			      adaptMult = 1.0 + adaptCFactor*((successRatioUse - 0.7)/0.3);
+			      newParamMult = adaptMult*currParam;
+			      if (newParamMult > currParam) {
+				if ((origAdaptCFactor < adaptCFactor) && (newParamMult > theMove->oldTuningParamMult)) {
+				  adaptMult = 1.0 + origAdaptCFactor*((successRatioUse - 0.7)/0.3);
+				  newParamMult = MaximumValue(theMove->oldTuningParamMult,(adaptMult*currParam));
+				}
+				if (newParamMult > maxParamMult)
+				  newParamMult = maxParamMult;
+				newParam = 2.0*log(newParamMult);
+				if (newParam > maxParam)
+				  newParam = maxParam;
+				if (newParam > theMove->tuningParam[j]) {
+				  theMove->tuningParam[j] = newParam;
+				  timesAdaptWiden++;
+				  adaptDidWiden = YES;
+				  theMove->nTriedAfterAdapt = 0;
+				  theMove->nAcceptedAfterAdapt = 0;
+				} else if (theMove->tuningParam[j] < (minParam+theMove->minAdapt)) {
+				  newParam = minParam + theMove->minAdapt;
+				  if (newParam > maxParam)
+				    newParam = maxParam;
+				  if (newParam > theMove->tuningParam[j]) {
+				    theMove->tuningParam[j] = newParam;
+				    timesAdaptWiden++;
+				    adaptDidWiden = YES;
+				    theMove->nTriedAfterAdapt = 0;
+				    theMove->nAcceptedAfterAdapt = 0;
+				  }
+				} else if (newParam < theMove->minAdapt) {
+				  MrBayesPrint ("%s   Error (%2.15lf under minAdapt %2.15lf; Mult; curr %2.15lf (%2.15lf); min %2.15lf; max %2.15lf) in Adapt for %s\n",
+						spacer, newParam, theMove->minAdapt, theMove->tuningParam[j], currParam,
+						minParam, maxParam, theMove->name);
+				  newParam = sqrt(maxParam*theMove->tuningParam[j]);
+				  if (newParam > maxParam)
+				    newParam = maxParam;
+				  if (newParam > theMove->tuningParam[j]) {
+				    theMove->tuningParam[j] = newParam;
+				    timesAdaptWiden++;
+				    adaptDidWiden = YES;
+				    theMove->nTriedAfterAdapt = 0;
+				    theMove->nAcceptedAfterAdapt = 0;
+				  }
+				}
+			      } else if (theMove->tuningParam[j] < (minParam+theMove->minAdapt)) {
+				newParam = minParam + theMove->minAdapt;
+				if (newParam > maxParam)
+				  newParam = maxParam;
+				if (newParam > theMove->tuningParam[j]) {
+				  theMove->tuningParam[j] = newParam;
+				  timesAdaptWiden++;
+				  adaptDidWiden = YES;
+				  theMove->nTriedAfterAdapt = 0;
+				  theMove->nAcceptedAfterAdapt = 0;
+				}
+			      }
+			    }
+			  } else if (adaptNarrow == YES) {
+			    if ((theMove->tuningParam[j] > theMove->oldTuningParam) && (n < burninGen)) {
+			      adaptCFactor = sqrt(theMove->tuningParam[j]/theMove->minAdapt)-1.0;
+			      origAdaptCFactor = CBRT(theMove->tuningParam[j]/theMove->minAdapt)-1.0;
+			    } else
+			      origAdaptCFactor = adaptCFactor = CBRT(theMove->tuningParam[j]/theMove->minAdapt)-1.0;
+			    if (theMove->tuningParam[j] > theMove->oldTuningParam) {
+			      adaptCFactor =
+				MaximumValue(adaptCFactor,((theMove->tuningParam[j]/theMove->oldTuningParam)-1.0));
+			      origAdaptCFactor =
+				MaximumValue(origAdaptCFactor,((theMove->tuningParam[j]/theMove->oldTuningParam)-1.0));
+			    }
+			    if (theMove->tuningParam[j] > minParam)
+			      minAdaptCFactor =
+				MinimumValue(theMove->minAdaptCFactor,((theMove->tuningParam[j]/minParam)-1.0));
+			    else if (theMove->tuningParam[j] > theMove->minAdapt) /* should be! */
+			      minAdaptCFactor =
+				MinimumValue(theMove->minAdaptCFactor,((theMove->tuningParam[j]/theMove->minAdapt)-1.0));
+			    else
+			      minAdaptCFactor = theMove->minAdaptCFactor;
+			    if (adaptCFactor < minAdaptCFactor)
+			      adaptCFactor = minAdaptCFactor;
+			    if (origAdaptCFactor < minAdaptCFactor)
+			      origAdaptCFactor = minAdaptCFactor;
+			    adaptMult = 1.0/(1.0 + adaptCFactor*((0.1 - successRatioUse)/0.1));
+			    newParam = adaptMult*theMove->tuningParam[j];
+			    if (newParam < minParam)
+			      newParam = minParam;
+			    if (newParam < theMove->tuningParam[j]) {
+			      if ((origAdaptCFactor < adaptCFactor) && (newParam < theMove->oldTuningParam)) {
+				adaptMult = 1.0/(1.0 + origAdaptCFactor*((0.1 - successRatioUse)/0.1));
+				newParam = MinimumValue(theMove->oldTuningParam,(adaptMult*theMove->tuningParam[j]));
+			      }
+			      if (newParam < theMove->tuningParam[j]) {
+				theMove->tuningParam[j] = newParam;
+				timesAdaptNarrow++;
+				adaptDidNarrow = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      } else if (theMove->tuningParam[j] > maxParam) {
+				theMove->tuningParam[j] = maxParam;
+				timesAdaptNarrow++;
+				adaptDidNarrow = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      }
+			    } else if (theMove->tuningParam[j] > maxParam) {
+			      theMove->tuningParam[j] = maxParam;
+			      timesAdaptNarrow++;
+			      adaptDidNarrow = YES;
+			      theMove->nTriedAfterAdapt = 0;
+			      theMove->nAcceptedAfterAdapt = 0;
+			    } else if (newParam > theMove->maxAdapt) {
+			      MrBayesPrint ("%s   Error (%2.15lf over maxAdapt %2.15lf; curr %2.15lf) in Adapt for %s\n",
+					    spacer, newParam, theMove->maxAdapt, theMove->tuningParam[j], theMove->name);
+			      newParam = (minParam+theMove->tuningParam[j])/2.0;
+			      if (newParam > minParam)
+				newParam = minParam;
+			      if (newParam < theMove->tuningParam[j]) {
+				theMove->tuningParam[j] = newParam;
+				timesAdaptNarrow++;
+				adaptDidNarrow = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      }
+			    }
+			  } else if (adaptWiden == YES) {
+			    if ((theMove->tuningParam[j] < theMove->oldTuningParam) && (n < burninGen)) {
+			      adaptCFactor = sqrt(theMove->maxAdapt/theMove->tuningParam[j])-1.0;
+			      origAdaptCFactor = CBRT(theMove->maxAdapt/theMove->tuningParam[j])-1.0;
+			    } else
+			      origAdaptCFactor = adaptCFactor = CBRT(theMove->maxAdapt/theMove->tuningParam[j])-1.0;
+			    if (theMove->tuningParam[j] < theMove->oldTuningParam) {
+			      adaptCFactor =
+				MaximumValue(adaptCFactor,((theMove->oldTuningParam/theMove->tuningParam[j])-1.0));
+			      origAdaptCFactor =
+				MaximumValue(origAdaptCFactor,((theMove->oldTuningParam/theMove->tuningParam[j])-1.0));
+			    }
+			    if (maxParam > theMove->tuningParam[j])
+			      minAdaptCFactor =
+				MinimumValue(theMove->minAdaptCFactor,((maxParam/theMove->tuningParam[j])-1.0));
+			    else if (theMove->maxAdapt > theMove->tuningParam[j])
+			      minAdaptCFactor =
+				MinimumValue(theMove->minAdaptCFactor,((theMove->maxAdapt/theMove->tuningParam[j])-1.0));
+			    else
+			      minAdaptCFactor = theMove->minAdaptCFactor;
+			    if (adaptCFactor < minAdaptCFactor)
+			      adaptCFactor = minAdaptCFactor;
+			    if (origAdaptCFactor < minAdaptCFactor)
+			      origAdaptCFactor = minAdaptCFactor;
+			    adaptMult = 1.0 + adaptCFactor*((successRatioUse - 0.7)/0.3);
+			    newParam = adaptMult*theMove->tuningParam[j];
+			    if (newParam > maxParam)
+			      newParam = maxParam;
+			    if (newParam > theMove->tuningParam[j]) {
+			      if ((origAdaptCFactor < adaptCFactor) && (newParam > theMove->oldTuningParam)) {
+				adaptMult = 1.0 + origAdaptCFactor*((successRatioUse - 0.7)/0.3);
+				newParam = MaximumValue(theMove->oldTuningParam,(adaptMult*theMove->tuningParam[j]));
+			      }
+			      if (newParam > theMove->tuningParam[j]) {
+				theMove->tuningParam[j] = newParam;
+				timesAdaptWiden++;
+				adaptDidWiden = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      } else if (theMove->tuningParam[j] < (minParam+theMove->minAdapt)) {
+				newParam = minParam + theMove->minAdapt;
+				if (newParam > maxParam)
+				  newParam = maxParam;
+				if (newParam > theMove->tuningParam[j]) {
+				  theMove->tuningParam[j] = newParam;
+				  timesAdaptWiden++;
+				  adaptDidWiden = YES;
+				  theMove->nTriedAfterAdapt = 0;
+				  theMove->nAcceptedAfterAdapt = 0;
+				}
+			      }
+			    } else if (theMove->tuningParam[j] < (minParam+theMove->minAdapt)) {
+			      newParam = minParam + theMove->minAdapt;
+			      if (newParam > maxParam)
+				newParam = maxParam;
+			      if (newParam > theMove->tuningParam[j]) {
+				theMove->tuningParam[j] = newParam;
+				timesAdaptWiden++;
+				adaptDidWiden = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      }
+			    } else if (newParam < theMove->minAdapt) {
+			      MrBayesPrint ("%s   Error (%2.15lf under minAdapt %2.15lf; curr %2.15lf) in Adapt for %s\n",
+					    spacer, newParam, theMove->minAdapt, theMove->tuningParam[j], theMove->name);
+			      newParam = (maxParam+theMove->tuningParam[j])/2.0;
+			      if (newParam > maxParam)
+				newParam = maxParam;
+			      if (newParam > theMove->tuningParam[j]) {
+				theMove->tuningParam[j] = newParam;
+				timesAdaptWiden++;
+				adaptDidWiden = YES;
+				theMove->nTriedAfterAdapt = 0;
+				theMove->nAcceptedAfterAdapt = 0;
+			      }
+			    }
+			  }
+			}
+
 			/*TouchAllPartitions();*/    /* for debugging copying shortcuts */
 
 			/* make move */
@@ -33488,16 +34953,104 @@
 				}
 
 			/* calculate likelihood ratio */
-			if (abortMove == NO)
-				{
+			if (abortMove == NO) {
 				  /*TouchAllTrees(chn);*/  /* for debugging copying shortcuts */
 				lnLike = LogLike(chn);
 				lnLikelihoodRatio = lnLike - curLnL[chn];
 				lnPrior = curLnPr[chn] + lnPriorRatio;
 
+				if (chainParams.useSA == YES) {
+				  if (n >= burninGen) {
+				    /* no temperature adjustment */
+				  } else if ((n*mult_short) <= 1.0) {
+				    lnLikelihoodRatio *= (tempZero/1000.0);
+				    lnPriorRatio      *= (tempZero/1000.0);
+				  } else if ((lnPriorRatio+lnLikelihoodRatio) > 0.0) {
+				    /* we would like at least 30% refusals */
+				    effectiveGen = ((n-movesAccepted)/0.3);
+				    if (effectiveGen <= 1.0)
+				      tempUse1 = tempZero;
+				    else if (effectiveGen >= partBurninGen)
+				      tempUse1 = 1000;
+				    else {
+				      tempUse1 = (int) (tempZero - (MrBFlt)(mult_short*effectiveGen*1000.0));
+				      if (tempUse1 > 999)
+					tempUse1 = 999;
+				    }
+
+				    tempUse = tempUse1;
+
+				    if (tempUse < 999) {
+				      if (n <= 1)
+					tempUse2 = tempZero;
+				      else
+					tempUse2 = (int) (tempZero - (MrBFlt)(mult_long*(MrBFlt)(n)*1000.0));
+				      if (tempUse2 > 999)
+					tempUse2 = 999;
+				      if (tempUse2 > tempUse)
+					tempUse = tempUse2;
+				    }
+
+				    if (tempUse < 1000) {
+				      if (tempUse < tempZero) {
+					lnLikelihoodRatio *= (tempZero/1000.0);
+					lnPriorRatio      *= (tempZero/1000.0);
+				      } else {
+					lnLikelihoodRatio *= (tempUse/1000.0);
+					lnPriorRatio      *= (tempUse/1000.0);
+				      }
+				    }
+				  } else if ((lnPriorRatio+lnLikelihoodRatio) < 0.0) {
+				    /* admittedly, some prior work (on Adaptive Simulated Annealing - L. Ingber)
+				       indicates that would be best to only decrease effective tempeature
+				       when have had acceptance... may work on later, although admittedly
+				       this somewhat contradicts "annealing" - if no acceptances, is
+				       a sign the system is at "equilibrium"... XXX - Allen */
+				    /* we would like at least 10% acceptances - thus 1/0.1 = 10; + 1.0 is
+				       because n includes this generation */
+				    if ((chainParams.useAdapt == YES) &&
+					(theMove->doAdaptParam == YES))
+				      effectiveGen = (movesAccepted*10.0) + 1.0;
+				    else
+				      effectiveGen = ((theMove->nAccepted[i]*10.0)/theMove->relProposalProb) + 1.0;
+				    if ((effectiveGen*mult_short) <= 1.0)
+				      tempUse1 = tempZero;
+				    else if (effectiveGen >= partBurninGen)
+				      tempUse1 = 1000;
+				    else {
+				      tempUse1 = (int) (tempZero - (MrBFlt)(mult_short*effectiveGen*1000.0));
+				      if (tempUse1 > 999)
+					tempUse1 = 999;
+				    }
+
+				    tempUse = tempUse1;
+
+				    if (tempUse < 999) {
+				      if ((n*mult_long) <= 1.0)
+					tempUse2 = tempZero;
+				      else
+					tempUse2 = (int) (tempZero - (MrBFlt)(mult_long*(MrBFlt)(n)*1000.0));
+				      if (tempUse2 > 999)
+					tempUse2 = 999;
+				      if (tempUse2 > tempUse)
+					tempUse = tempUse2;
+				    }
+
+				    if (tempUse < 1000) {
+				      if (tempUse < tempZero) {
+					lnLikelihoodRatio *= (tempZero/1000.0);
+					lnPriorRatio      *= (tempZero/1000.0);
+				      } else {
+					lnLikelihoodRatio *= (tempUse/1000.0);
+					lnPriorRatio      *= (tempUse/1000.0);
+				      }
+				    }
+				  }
+				} else {
 				/* heat */
-				lnLikelihoodRatio *= Temperature (chainId[chn]);
-				lnPriorRatio      *= Temperature (chainId[chn]);
+				  lnLikelihoodRatio *= Temperature (i);
+				  lnPriorRatio      *= Temperature (i);
+				}
 			
 				/* calculate the acceptance probability */
 				if (lnLikelihoodRatio + lnPriorRatio + lnProposalRatio < -100.0)
@@ -33510,12 +35063,47 @@
 
 			/* decide to accept or reject the move */
 			acceptMove = NO;
-			i = chainId[chn];
+			if (abortMove == NO) {
 			theMove->nTried[i]++;
-			if (abortMove == NO && RandomNumber(seed) < r)
-				{
+			  theMove->nTriedAfterAdapt++;
+			  if ((chainParams.useAdapt == YES) && (theMove->doAdaptParam == YES) &&
+			      ((adaptDidNarrow == YES) || (adaptDidWiden == YES))) {
+			    j = theMove->whichAdaptParam;
+			    if (theMove->minParamTried > theMove->tuningParam[j])
+			      theMove->minParamTried = theMove->tuningParam[j];
+			    if (theMove->maxParamTried < theMove->tuningParam[j])
+			      theMove->maxParamTried = theMove->tuningParam[j];
+			  }
+			  if (RandomNumber(seed) < r) {
 				acceptMove = YES;
 				theMove->nAccepted[i]++;
+			    movesAccepted++;
+			    theMove->nAcceptedAfterAdapt++;
+			    if ((chainParams.useAdapt == YES) && (theMove->doAdaptParam == YES) &&
+				((adaptDidNarrow == YES) || (adaptDidWiden == YES))) {
+			      j = theMove->whichAdaptParam;
+			      if (theMove->minParamAccepted > theMove->tuningParam[j])
+				theMove->minParamAccepted = theMove->tuningParam[j];
+			      /* yes, it can be both above the max and below the min so far (starting vals) */
+			      if (theMove->maxParamAccepted < theMove->tuningParam[j])
+				theMove->maxParamAccepted = theMove->tuningParam[j];
+			      if (adaptDidNarrow)
+				timesAdaptNarrowAccepted++;
+			      else
+				timesAdaptWidenAccepted++;
+			    }
+			  } else if ((chainParams.useAdapt == YES) && (theMove->doAdaptParam == YES) &&
+				     ((adaptDidNarrow == YES) || (adaptDidWiden == YES))) {
+			    j = theMove->whichAdaptParam;
+			    if (theMove->minParamRefused > theMove->tuningParam[j])
+			      theMove->minParamRefused = theMove->tuningParam[j];
+			    if (theMove->maxParamRefused < theMove->tuningParam[j])
+			      theMove->maxParamRefused = theMove->tuningParam[j];
+			    if (adaptDidNarrow)
+			      timesAdaptNarrowRefused++;
+			    else
+			      timesAdaptWidenRefused++;
+			  }
 				}
 
 			/* update the chain */
@@ -33542,13 +35130,20 @@
 		/* attempt swap(s) */
 		if (chainParams.numChains > 1 && n % chainParams.swapFreq == 0)
 			{
-			for (i = 0; i<chainParams.numRuns; i++)
-				{
-				for (j = 0; j<chainParams.numSwaps; j++)
+			  for (i = 0; i<chainParams.numRuns; i++) {
+			    /* the below had previously lacked a bracket... XXX - Allen */
+			    for (j = 0; j<chainParams.numSwaps; j++) {
+#ifdef DEBUG_MPI
+			      MrBayesPrint ("%s   Getting swappers (run %d iteration %d)\n",spacer,i,j);
+#endif
 				GetSwappers (&swapA, &swapB, i);
-				if (AttemptSwap (swapA, swapB, seed) == ERROR)
-					{
+#ifdef DEBUG_MPI
+			      MrBayesPrint ("%s   Swappers are %d and %d\n",spacer,swapA,swapB);
+#endif
+			      if (AttemptSwap (swapA, swapB, seed) == ERROR) {
 #					if defined (MPI_ENABLED)
+				if (proc_id == 0)
+				  MrBayesPrint ("%s   Unsuccessful swap of states\n", spacer);
 					nErrors++;
 #					else
 					MrBayesPrint ("%s   Unsuccessful swap of states\n", spacer);
@@ -33555,16 +35150,92 @@
 					return ERROR;
 #					endif
 					}
+#ifdef DEBUG_MPI
+			      MrBayesPrint ("%s   Tried swapping states\n", spacer);
+#endif
+			    }
 				}
 			}
 
 		/* print information to screen */
-		if (n == 1 || n % chainParams.printFreq == 0)
-			PrintToScreen(n, time(0), startingT);
+		if (n == 1 || n % chainParams.printFreq == 0) {
+		  PrintToScreen(n, movesAccepted, time(0), startingT);
+		  if ((timesAdaptNarrow > 0) || (timesAdaptWiden > 0)) {
+		    MrBayesPrint("%s  Adapt narrowed windows %d times (%d accepted, %d refused), widened %d times (%d accepted, %d refused)\n",
+				 spacer,timesAdaptNarrow,timesAdaptNarrowAccepted,timesAdaptNarrowRefused,
+				 timesAdaptWiden,timesAdaptWidenAccepted,timesAdaptWidenRefused);
+		    timesAdaptNarrow = timesAdaptNarrowAccepted = timesAdaptNarrowRefused = 0;
+		    timesAdaptWiden = timesAdaptWidenAccepted = timesAdaptWidenRefused = 0;
+		  }
+		  if (num_times_scaler_adjusted > 0) { /* - Allen */
+#		    if defined (MPI_ENABLED)
+		    if (proc_id == 0)
+#endif
+		      MrBayesPrint ("%s  WARNING: Scaler adjusted %d times (%d big)\n",spacer,
+				    num_times_scaler_adjusted, num_times_scaler_adjusted_big);
+		    num_times_scaler_adjusted = 0;
+		    num_times_scaler_adjusted_big = 0;
+		  }
+		  if (num_times_scaler_adjusted2 > 0) { /* - Allen */
+#		    if defined (MPI_ENABLED)
+		    if (proc_id == 0)
+#endif
+		      MrBayesPrint ("%s  WARNING: Scaler adjusted (2) %d times (%d big)\n",spacer,
+				    num_times_scaler_adjusted2, num_times_scaler_adjusted2_big);
+		    num_times_scaler_adjusted2 = 0;
+		    num_times_scaler_adjusted2_big = 0;
+		  }
+		  num_times_scaler_warn3 = 0;
+		  if (((chainParams.useSA == YES) ||
+		       (chainParams.useAdapt == YES)) &&
+		      (chainParams.numChains == 1) && (chainParams.numRuns == 1) &&
+		      (((n >= burninGen) && (didReportBurnin == NO)) ||
+		       ((n >= partBurninGen) && (didReportPartBurnin == NO)))) {
+		    MrBayesPrint ("%s   Current acceptance rates for the moves:\n", spacer);
+		    MrBayesPrint ("%s      With prob.  Chain accepted changes to\n", spacer);
+		  
+		    for (i=0; i<numMoves; i++) {
+		      mv = &moves[i];
+		      if (mv->nTried[0] == 0)
+			MrBayesPrint ("%s      not tried   param. %d %s\n", spacer, mv->parm->index+1, mv->name);
+		      else {
+			MrBayesPrint ("%s       %6.2f %%   param. %d %s\n", spacer,
+				      100.0*(MrBFlt)(mv->nAccepted[0])/(MrBFlt)(mv->nTried[0]), mv->parm->index+1,
+				      mv->name);
+			if ((chainParams.useAdapt == YES) && (mv->doAdaptParam == YES)) {
+			  if (mv->nTriedAfterAdapt > 0)
+			    MrBayesPrint("%s       %6.2f %%   param. %d %s after adapt\n", spacer,
+					 100.0*(MrBFlt)(mv->nAcceptedAfterAdapt)/(MrBFlt)(mv->nTriedAfterAdapt),
+					 mv->parm->index+1, mv->name);
+			  else
+			    MrBayesPrint("%s       not tried  param. %d %s after adapt\n", spacer,
+					 mv->parm->index+1, mv->name);
+			  k = mv->whichAdaptParam;
+			  MrBayesPrint("%s     Adapt:             %6.6f to       %6.6f\n",spacer,
+				       mv->oldTuningParam,mv->tuningParam[k]);
+			  MrBayesPrint("%s     Adapt tried:       %6.6f to       %6.6f\n",spacer,
+				       mv->minParamTried, mv->maxParamTried);
+			  MrBayesPrint("%s     Adapt accepted:    %6.6f to       %6.6f\n",spacer,
+				     mv->minParamAccepted, mv->maxParamAccepted);
+			  MrBayesPrint("%s     Adapt refused:     %6.6f to       %6.6f\n",spacer,
+				       mv->minParamRefused, mv->maxParamRefused);
+			}
+		      }
+		    }
+		    if (n >= burninGen) {
+		      didReportPartBurnin = YES;
+		      didReportBurnin = YES;
+		    } else
+		      didReportPartBurnin = YES;
+		  }
+		}
 
 		/* print information to files */
 		if (n == 1 || n == chainParams.numGen || n % chainParams.sampleFreq == 0)
 			{
+#if defined(DEBUG_MPI)
+			  MrBayesPrint ("%s  WARNING: Attempting printing of states\n",spacer);
+#endif
 			if (PrintStatesToFiles (n) == ERROR)
 				{
 #				if defined (MPI_ENABLED)
@@ -33581,6 +35252,9 @@
 			if (chainParams.numRuns > 1 && ((n > 1 && chainParams.relativeBurnin == YES)
 				|| (n >= chainParams.chainBurnIn * chainParams.sampleFreq && chainParams.relativeBurnin == NO)))
 				{
+#if defined(DEBUG_MPI)
+				  MrBayesPrint ("%s  WARNING: Checking tree samples\n",spacer);
+#endif
 				/* we need some space for coming output */
 				MrBayesPrint ("\n");
 				/* first reassemble files if mpi version */				
@@ -33651,6 +35325,9 @@
 				if (chainParams.allComps == YES)
 					PrintTopConvInfo ();
 #				if defined (MPI_ENABLED)
+#                               if defined(DEBUG_MPI)
+				MrBayesPrint ("%s  WARNING: Attempting stopChain broadcast\n",spacer);
+#endif
 					}
 				ierror = MPI_Bcast (&stopChain, 1, MPI_INT, 0, MPI_COMM_WORLD);
 				if (ierror != MPI_SUCCESS)
@@ -33661,6 +35338,9 @@
 #				endif
 				}
 			/* part of the following function needs to be performed by all MPI processors */
+#                       if defined(DEBUG_MPI)
+			MrBayesPrint ("%s  WARNING: Attempting printing of MCMC data\n",spacer);
+#                       endif
 			if (PrintMCMCDiagnosticsToFile (n) == ERROR)
 				nErrors++;
 #			if defined (MPI_ENABLED)
@@ -33688,7 +35368,19 @@
 			{
 			  stoppingT1 = time(0); 
 			  currentCPUTime = clock();
+			  if (currentCPUTime < previousCPUTime) { /* wraparound - Allen */
+			    if ((currentCPUTime < 0) && (previousCPUTime > 0)) { /* Approximating */
+			      if (previousCPUTime > (-1*currentCPUTime)) {
+				CPUTime += (2*(previousCPUTime + currentCPUTime)) / (MrBFlt) (CLOCKS_PER_SEC);
+			      } else {
+				CPUTime += (-2*(previousCPUTime + currentCPUTime)) / (MrBFlt) (CLOCKS_PER_SEC);
+			      }
+			    } else {
+			      CPUTime += (previousCPUTime - currentCPUTime) / (MrBFlt) (CLOCKS_PER_SEC);
+			    }
+			  } else {
 			  CPUTime += (currentCPUTime - previousCPUTime) / (MrBFlt) CLOCKS_PER_SEC;
+			  }
 			  previousCPUTime = currentCPUTime;
 			  chainParams.numGen += ExtendChainQuery ();
 			  stoppingT2 = time(0);
@@ -33713,7 +35405,22 @@
 		} /* end run chain */
 	endingT = time(0);
 	currentCPUTime = clock();
+#if defined(MPI_ENABLED) && defined(DEBUG_MPI) && LAM_MPI
+	MPIL_Trace_off();
+#endif
+	if (currentCPUTime < previousCPUTime) { /* wraparound - Allen */
+	  if ((currentCPUTime < 0) && (previousCPUTime > 0)) { /* Approximating */
+	    if (previousCPUTime > (-1*currentCPUTime)) {
+	      CPUTime += (2*(previousCPUTime + currentCPUTime)) / (MrBFlt) (CLOCKS_PER_SEC);
+	    } else {
+	      CPUTime += (-2*(previousCPUTime + currentCPUTime)) / (MrBFlt) (CLOCKS_PER_SEC);
+	    }
+	  } else {
+	    CPUTime += (previousCPUTime - currentCPUTime) / (MrBFlt) (CLOCKS_PER_SEC);
+	  }
+	} else {
 	CPUTime += (currentCPUTime - previousCPUTime) / (MrBFlt) CLOCKS_PER_SEC;
+	}
 
 	CloseMBPrintFiles (); /* redundant because files closed in FreeChainMemory but kept here as a safeguard in case of future changes */
 
@@ -33805,8 +35512,7 @@
 #	endif
 
 	/* print acceptance rates for the moves */
-	for (j=0; j<chainParams.numChains*chainParams.numRuns; j++)
-		{
+	for (j=0; j<chainParams.numChains*chainParams.numRuns; j++) {
 		if (chainParams.numChains == 1)
 			{
 			if (chainParams.numRuns == 1)
@@ -33829,20 +35535,39 @@
 				MrBayesPrint ("%s   Acceptance rates for the moves in chain %d of run %d (heated):\n", spacer, j%chainParams.numChains+1, j/chainParams.numChains+1);
 			}
 
-		if (j % chainParams.numChains == 0 || chainParams.allChains == YES)
-			{
+		if (j % chainParams.numChains == 0 || chainParams.allChains == YES) {
 			MrBayesPrint ("%s      With prob.  Chain accepted changes to\n", spacer);
 
-			for (i=0; i<numMoves; i++)
-				{
+		  for (i=0; i<numMoves; i++) {
 				mv = &moves[i];
-				if (mv->nTried[0] == 0)
+		    if (mv->nTried[j] == 0)
 					MrBayesPrint ("%s      not tried   param. %d %s\n", spacer, mv->parm->index+1, mv->name);
+		    else {
+		      MrBayesPrint ("%s       %6.2f %%   param. %d %s\n", spacer,
+				    100.0*(MrBFlt)(mv->nAccepted[j])/(MrBFlt)(mv->nTried[j]), mv->parm->index+1,
+				    mv->name);
+		      if ((chainParams.useAdapt == YES) && (mv->doAdaptParam == YES)) {
+			if (mv->nTriedAfterAdapt > 0)
+			  MrBayesPrint("%s       %6.2f %%   param. %d %s after adapt\n", spacer,
+				       100.0*(MrBFlt)(mv->nAcceptedAfterAdapt)/(MrBFlt)(mv->nTriedAfterAdapt),
+				       mv->parm->index+1, mv->name);
 				else
-					MrBayesPrint ("%s       %6.2f %%   param. %d %s\n", spacer, 100.0*(MrBFlt)(mv->nAccepted[j])/(MrBFlt)(mv->nTried[j]), mv->parm->index+1, mv->name);
+			  MrBayesPrint("%s       not tried  param. %d %s after adapt\n", spacer,
+				       mv->parm->index+1, mv->name);
+			k = mv->whichAdaptParam;
+			MrBayesPrint("%s     Adapt:             %6.6f to       %6.6f\n",spacer,
+				     mv->oldTuningParam,mv->tuningParam[k]);
+			MrBayesPrint("%s     Adapt tried:       %6.6f to       %6.6f\n",spacer,
+				     mv->minParamTried, mv->maxParamTried);
+			MrBayesPrint("%s     Adapt accepted:    %6.6f to       %6.6f\n",spacer,
+				     mv->minParamAccepted, mv->maxParamAccepted);
+			MrBayesPrint("%s     Adapt refused:     %6.6f to       %6.6f\n",spacer,
+				     mv->minParamRefused, mv->maxParamRefused);
 				}
 			}
 		}
+		}
+	}
 
 	/* output information on the success of the chain state swap proposals */
 	if (PrintSwapInfo () == ERROR)
@@ -33899,16 +35624,52 @@
 int SaveSprintf(char **target, int *targetLen, char *fmt, ...) {
   va_list    argp;
   int        len,retval;
+#if !defined(VISUAL)
+  char  *tempStr;
+#endif
 
   va_start(argp, fmt);
 #ifdef VISUAL
   len = _vsnprintf(NULL, 0, fmt, argp);
 #else
+# ifdef _SGI_SOURCE
+   len = _xpg5_vsnprintf(NULL, 0, fmt, argp);
+   if (len < 0) {
+     va_end(argp);
+     va_start(argp, fmt);
+     tempStr = (char *) SafeMalloc((size_t) (*targetLen * sizeof(char)));
+     if (tempStr) {
+      len = _xpg5_vsnprintf(tempStr, *targetLen, fmt, argp);
+      free(tempStr);
+     }
+   }
+# else
   len = vsnprintf(NULL, 0, fmt, argp);
+   if (len < 0) {
+     va_end(argp);
+     va_start(argp, fmt);
+     tempStr = (char *) SafeMalloc((size_t) (*targetLen * sizeof(char)));
+     if (tempStr) {
+      len = vsnprintf(tempStr, *targetLen, fmt, argp);
+      free(tempStr);
+     }
+   }
+# endif
+#endif
+
+  if (len < 0) {
+    fprintf(stderr, "Negative length for printing to /%s/ in SaveSprintf\n",*target);
+#ifdef BUFSIZ
+    len = *targetLen + BUFSIZ;
+#else
+    len = *targetLen + 4096;
 #endif
+  }
 
   va_end(argp);
 
+  len++; /* len from vsnprintf does not include trailing null - Allen */
+
   if(len>*targetLen) 
         {
 /* 	  fprintf(stderr, "* increasing buffer from %d to %d bytes\n", *targetLen, len+TARGETLENDELTA); */
@@ -33917,7 +35678,11 @@
         }
 
   va_start(argp,fmt);
+#ifdef VISUAL
   retval=vsprintf(*target, fmt, argp);
+#else
+  retval=vsnprintf(*target, *targetLen, fmt, argp);
+#endif
   va_end(argp);
 
 /*   fprintf(stderr, "* savesprintf /%s/\n",*target); */
@@ -35504,11 +37270,19 @@
 					}
 
 			/* find the parameter x prior type */
+			if (mp->dataType == PROTEIN) { /* Allen */
+				if ((!strcmp(mp->aaModelPr, "Fixed")) &&
+				    (!strcmp(mp->aaModel, "Gtr")) &&
+				    (!strcmp(mp->aaRevMatPr,"Dirichlet")))
+					p->paramId = REVMAT_DIR;
+				else
+					p->paramId = REVMAT_FIX;
+			} else {
 			if (!strcmp(mp->revMatPr,"Dirichlet"))
 				p->paramId = REVMAT_DIR;
 			else
 				p->paramId = REVMAT_FIX;
-
+			}
 			if (p->paramId != REVMAT_FIX)
 				p->printParam = YES;
 			if (FillRelPartsString (p, partString) == YES)
@@ -35548,7 +37322,6 @@
 								{
 								SaveSprintf(&tempStr, &tempStrSize, "\tr(%c<->%c)", StateCode_AA(n1), StateCode_AA(n2));
 								strcat (p->paramName, tempStr);
-								strcat (p->paramName, partString);
 								}
 							}
 						}
@@ -36157,7 +37930,12 @@
 			/* find the parameter x prior type */
 			if (!strcmp(mp->covSwitchPr,"Uniform"))
 				p->paramId = SWITCH_UNI;
-			else if (!strcmp(mp->covSwitchPr,"Exponential"))
+			else if (!strcmp(mp->covSwitchPr,"PinvarUniform")) { /* - Allen */
+			  if (mp->covswitchMaxPinvar < 1)
+			    p->paramId = SWITCH_PINVUNI;
+			  else
+			    p->paramId = SWITCH_UNI;
+			} else if (!strcmp(mp->covSwitchPr,"Exponential"))
 				p->paramId = SWITCH_EXP;
 			else
 				p->paramId = SWITCH_FIX;
@@ -37036,11 +38814,11 @@
 		/* check if we should calculate ancestral states */
 		if (!strcmp(mp->inferAncStates,"Yes"))
 			{
-			if (mp->dataType == PROTEIN && !strcmp(mp->covarionModel, "No"))
+			if (mp->dataType == PROTEIN) /* am implementing covarion anc states - Allen */
 				m->printAncStates = YES;
 			else if (mp->dataType == DNA || mp->dataType == RNA)
 				{
-				if (!strcmp(mp->nucModel,"4by4") && !strcmp(mp->covarionModel, "No"))
+				if (!strcmp(mp->nucModel,"4by4"))
 					m->printAncStates = YES;
 				}
 			else if (mp->dataType == STANDARD || mp->dataType == RESTRICTION)
@@ -37087,7 +38865,7 @@
 
 {
 	
-	int			i, j, j1, k, n;
+	int			i, j, j1, k, n, has_covarion;
 	MrBFlt		sum;
 
 	Param		*p;
@@ -37109,9 +38887,22 @@
 		for (i=0; i<NUM_MOVE_TYPES; i++)
 			for (j=0; j<moveTypes[i].nApplicable; j++)
 				if (moveTypes[i].applicableTo[j] == p->paramId)
-					if (p->relProposalProb * moveTypes[i].relProposalProb > 0.00001)
+			    if (p->relProposalProb * moveTypes[i].relProposalProb > 0.00001) {
+			      if (moveTypes[i].covarion_wanted != 0) {
+				has_covarion = -1;
+				/* gamma has diff acceptances with, w/out covarion - Allen */
+				for (j1=0; j1<p->nRelParts; j1++) {
+				  if (!strcmp(modelParams[p->relParts[j1]].covarionModel, "Yes")) {
+				    has_covarion = 1;
+				    break;
+				  }
+				}
+				if (has_covarion != moveTypes[i].covarion_wanted)
+				  continue;
+			      }
 						numMoves++;
 		}
+		}
 
 	/* then allocate space for moves */
 	if (memAllocs[ALLOC_MOVES] == YES)
@@ -37153,6 +38944,18 @@
 					{
 					if (p->relProposalProb * moveTypes[i].relProposalProb < 0.00001)
 						continue;
+					if (moveTypes[i].covarion_wanted != 0) {
+					  has_covarion = -1;
+					  /* gamma has diff acceptances with, w/out covarion - Allen */
+					  for (j1=0; j1<p->nRelParts; j1++) {
+					    if (!strcmp(modelParams[p->relParts[j1]].covarionModel, "Yes")) {
+					      has_covarion = 1;
+					      break;
+					    }
+					  }
+					  if (has_covarion != moveTypes[i].covarion_wanted)
+					    continue;
+					}
 					moves[n].name = moveTypes[i].name;
 					moves[n].shortName = moveTypes[i].shortName;
 					moves[n].relProposalProb = p->relProposalProb * moveTypes[i].relProposalProb;
@@ -37160,6 +38963,45 @@
 					moves[n].cumProposalProb = sum;
 					moves[n].parm = p;
 					moves[n].moveFxn = moveTypes[i].moveFxn;
+					moves[n].doAdaptParam = moveTypes[i].doAdaptParam;
+					moves[n].whichAdaptParam = moveTypes[i].whichAdaptParam;
+					if ((moves[n].doAdaptParam == YES) &&
+					    (moveTypes[i].minimum[moves[n].whichAdaptParam] < LIKE_EPSILON)) {
+					  MrBayesPrint ("%s   WARNING: %s cannot be adaptive due to too-low minimum %lf\n",
+							spacer, moves[n].name,
+							moveTypes[i].minimum[moves[n].whichAdaptParam]);
+					  moves[n].doAdaptParam = NO;
+					}
+					moves[n].adaptMultiplier = moveTypes[i].adaptMultiplier;
+					if (moves[n].doAdaptParam == YES) {
+					  moves[n].oldTuningParam = moveTypes[i].tuningParam[moves[n].whichAdaptParam];
+					  moves[n].minParamTried = moves[n].oldTuningParam;
+					  moves[n].maxParamTried = moves[n].oldTuningParam;
+					  moves[n].minAdapt = moveTypes[i].minimum[moves[n].whichAdaptParam];
+					  if (p->paramId == SHAPE_UNI)
+					    moves[n].maxAdapt =
+					      MinimumValue(moveTypes[i].maximum[moves[n].whichAdaptParam],
+							   (log(modelParams[p->relParts[0]].shapeUni[1]/
+								modelParams[p->relParts[0]].shapeUni[0])*2.0));
+					  else if ((p->paramId == SWITCH_UNI) || (p->paramId == SWITCH_PINVUNI))
+					    moves[n].maxAdapt =
+					      MinimumValue(moveTypes[i].maximum[moves[n].whichAdaptParam],
+							   ((modelParams[p->relParts[0]].covswitchUni[1] -
+							     modelParams[p->relParts[0]].covswitchUni[0])*2.0));
+					  else
+					    moves[n].maxAdapt = moveTypes[i].maximum[moves[n].whichAdaptParam];
+					  moves[n].minParamRefused = moves[n].maxAdapt;
+					  moves[n].maxParamRefused = moves[n].minAdapt;
+					  moves[n].minParamAccepted = moves[n].maxAdapt;
+					  moves[n].maxParamAccepted = moves[n].minAdapt;
+					  if (moves[n].adaptMultiplier == YES) { /* precalculate */
+					    moves[n].minAdaptMult = exp(moves[n].minAdapt/2.0);
+					    moves[n].maxAdaptMult = exp(moves[n].maxAdapt/2.0);
+					  }
+					}
+					moves[n].minAdaptCFactor = moveTypes[i].minAdaptCFactor;
+					moves[n].nTriedAfterAdapt = 0;
+					moves[n].nAcceptedAfterAdapt = 0;
 					for (j1=0; j1<chainParams.numRuns; j1++)
 						moves[n].nAccepted[j1] = moves[n].nTried[j1] = 0;
 					moves[n].tuningParam[0] = moveTypes[i].tuningParam[0];
@@ -38316,6 +40158,12 @@
 		mt->name = mt->shortName = mt->nameTuning[0] = mt->nameTuning[1] = "";
 		mt->relProposalProb = 0.0;
 		mt->parsimonyBased = NO;
+		mt->doAdaptParam = NO;
+		mt->whichAdaptParam = 0;
+		mt->adaptMultiplier = NO;
+		/* if had acceptance rate of 1.0, would mult by 1/0.7 */
+		mt->minAdaptCFactor = (1.0/0.7)-1.0;
+		mt->covarion_wanted = 0;
 		}
 
 	/* All moves are in alphabetic order of move function name. This list should match
@@ -38343,9 +40191,9 @@
 	mt->applicableTo[0] = CORREL_UNI;
 	mt->nApplicable = 1;
 	mt->moveFxn = &Move_Adgamma;
-	mt->relProposalProb = 1.0;
+	mt->relProposalProb = 3.0; /* changed from 1.0 - Allen */
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 0.5;  /* window size */
+	mt->tuningParam[0] = 0.25;  /* window size - changed to 0.25 due to usage of 0<->1 range - Allen */
 	mt->minimum[0] = 0.001;
 	mt->maximum[0] = 1.999;
 	mt->parsimonyBased = NO;
@@ -38400,12 +40248,20 @@
 	mt->applicableTo[1] = BRLENS_EXP;
 	mt->nApplicable = 2;
 	mt->moveFxn = &Move_BrLen;
+#ifdef NOTOPOLOGY /* - Allen */
+	mt->relProposalProb = 10.0;
+#else
 	mt->relProposalProb = 0.0;
+#endif
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 2.0 * log (2.0);  /* lambda */
-	mt->minimum[0] = 0.00001;
-	mt->maximum[0] = 10000000.0;
+	mt->tuningParam[0] = 2.0 * log (20.0);  /* lambda; was log(2.0) - Allen */
+	mt->minimum[0] = 2.0 * log((BRLEN_EPSILON+BRLENS_MIN)/BRLENS_MIN);
+	mt->maximum[0] = log(BRLENS_MAX/BRLENS_MIN)*2.0; /* original would cause problems with exp - Allen */
 	mt->parsimonyBased = NO;
+	mt->doAdaptParam = YES;
+	mt->whichAdaptParam = 0;
+	mt->adaptMultiplier = YES;
+	/* leaving minAdaptCFactor at default for now */
 	mt->level = STANDARD_USER;
 
 	/* Move_Extinction */
@@ -38534,14 +40390,18 @@
 	mt->applicableTo[1] = TOPOLOGY_NCL_CONSTRAINED_HOMO;
 	mt->nApplicable = 2;
 	mt->moveFxn = &Move_ExtTBR;
-	mt->relProposalProb = 15.0;
+#ifdef NOTOPOLOGY /* -Allen */
+	mt->relProposalProb = 0.0;
+#else
+	mt->relProposalProb = 22.0; /* changed from 15.0 - Allen */
+#endif
 	mt->numTuningParams = 2;
-	mt->tuningParam[0] = 0.8;  /* extension probability */
-	mt->tuningParam[1] = 2.0 * log (1.6);  /* lambda */
+	mt->tuningParam[0] = 0.54;  /* extension probability */ /* changed from 0.8 - Allen */
+	mt->tuningParam[1] = 2.0 * log (1.048);  /* lambda */ /* changed from 2.0 * log(1.6) - Allen */
 	mt->minimum[0] = 0.00001;
 	mt->maximum[0] = 0.99;
 	mt->minimum[1] = 0.00001;
-	mt->maximum[1] = 10000000.0;
+	mt->maximum[1] = 200.0; /* default changed due to potential problems with exp - Allen */
 	mt->parsimonyBased = NO;
 	mt->level = STANDARD_USER;
 
@@ -38631,7 +40491,7 @@
 
 	/* Move_GammaShape_M */
 	mt = &moveTypes[i++];
-	mt->name = "(gamma shape) with multiplier";
+	mt->name = "(gamma shape) with multiplier (non-covarion)";
 	mt->nameTuning[0] = "multiplier tuning parameter (lambda)";
 	mt->shortName = "mult";
 	mt->applicableTo[0] = SHAPE_UNI;
@@ -38640,10 +40500,37 @@
 	mt->moveFxn = &Move_GammaShape_M;
 	mt->relProposalProb = 1.0;
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 2.0 * log (1.5);  /* lambda */
+	mt->tuningParam[0] = 2.0 * log(1.5); /* lambda for non-covarion - Allen */
+	mt->minimum[0] = 2.0 * log(0.65/(0.65-MIN_SHAPE_PARAM)); /* 0.65 is equiv to Grishin distance - Allen */
+	/* orig 10000000.0 as max caused problems with exp - Allen */
+	mt->maximum[0] = log(MAX_SHAPE_PARAM/MIN_SHAPE_PARAM)*2.0;
+	mt->parsimonyBased = NO;
+	mt->doAdaptParam = YES;
+	mt->whichAdaptParam = 0;
+	mt->adaptMultiplier = YES;
+	mt->covarion_wanted = -1;
+	/* leaving minAdaptCFactor as per default */
+	mt->level = STANDARD_USER;
+
+	mt = &moveTypes[i++];
+	mt->name = "(gamma shape) with multiplier (covarion)";
+	mt->nameTuning[0] = "multiplier tuning parameter (lambda)";
+	mt->shortName = "mult_with_covarion";
+	mt->applicableTo[0] = SHAPE_UNI;
+	mt->applicableTo[1] = SHAPE_EXP;
+	mt->nApplicable = 2;
+	mt->moveFxn = &Move_GammaShape_M; /* same function as above! - Allen */
+	mt->relProposalProb = 1.5; /* was 1.0 - Allen */
+	mt->numTuningParams = 1;
+	mt->tuningParam[0] = 2.0 * log(1.00003); /* lambda for covarion - Allen */
 	mt->minimum[0] = 0.00001;
-	mt->maximum[0] = 10000000.0;
+	mt->maximum[0] = log(MAX_SHAPE_PARAM/MIN_SHAPE_PARAM)*2.0; /* see above */
 	mt->parsimonyBased = NO;
+	mt->doAdaptParam = YES;
+	mt->whichAdaptParam = 0;
+	mt->adaptMultiplier = YES;
+	mt->covarion_wanted = 1;
+	/* leaving minAdaptCFactor as per default */
 	mt->level = STANDARD_USER;
 
 	/* Move_Growth */
@@ -38673,11 +40560,15 @@
 	mt->applicableTo[1] = TOPOLOGY_NCL_CONSTRAINED_HOMO;
 	mt->nApplicable = 2;
 	mt->moveFxn = &Move_Local;
-	mt->relProposalProb = 5.0;
+#ifdef NOTOPOLOGY /* - Allen */
+	mt->relProposalProb = 0.0;	
+#else
+	mt->relProposalProb = 14.0; /* changed from 5.0 - Allen */
+#endif
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 2.0 * log (1.1);  /* lambda */
+	mt->tuningParam[0] = 2.0 * log (1.00105);  /* lambda */ /* changed from log(1.1) - Allen */
 	mt->minimum[0] = 0.00001;
-	mt->maximum[0] = 10000000.0;
+	mt->maximum[0] = 200.0; /* changed from 10000000.0 due to problems with exp - Allen */
 	mt->parsimonyBased = NO;
 	mt->level = STANDARD_USER;
 
@@ -38738,10 +40629,43 @@
 	mt->moveFxn = &Move_NodeSlider;
 	mt->relProposalProb = 0.0;
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 2.0 * log (1.1);  /* lambda */
-	mt->minimum[0] = 0.00001;
-	mt->maximum[0] = 10000000.0;
+	mt->tuningParam[0] = 2.0 * log(1.02);  /* lambda; was log(1.1) - Allen */
+	mt->minimum[0] = 2.0 * log((BRLEN_EPSILON+BRLENS_MIN)/BRLENS_MIN);
+	mt->maximum[0] = log(BRLENS_MAX/BRLENS_MIN)*2.0; /* see Move_BrLen */
+	mt->parsimonyBased = NO;
+	mt->doAdaptParam = NO; /* turning off b/c is too insensitive to tuning param; use SimAnneal instead */
+	mt->whichAdaptParam = 0;
+	mt->adaptMultiplier = YES;
+	mt->minAdaptCFactor = 2; /* guess! move is relatively insensitive to tuning parameter... - Allen */
+	mt->level = STANDARD_USER;
+
+	/* Move_NodeSlider2 */
+	mt = &moveTypes[i++];
+	mt->name = "(branch lengths) with nodeslider2";
+	mt->nameTuning[0] = "multiplier tuning parameter (lambda)";
+	mt->nameTuning[1] = "sliding window size for proportion";
+	mt->shortName = "nslider2";
+	mt->applicableTo[0] = BRLENS_UNI;
+	mt->applicableTo[1] = BRLENS_EXP;
+	mt->nApplicable = 2;
+	mt->moveFxn = &Move_NodeSlider2;
+#ifdef NOTOPOLOGY /* - Allen */
+	mt->relProposalProb = 27.0;
+#else
+	mt->relProposalProb = 0.0;
+#endif
+	mt->numTuningParams = 2;
+	mt->tuningParam[0] = 2.0 * log(1.011);  /* lambda; was log(1.1) - Allen */
+	mt->minimum[0] = 2.0 * log((BRLEN_EPSILON+BRLENS_MIN)/BRLENS_MIN);
+	mt->maximum[0] = log(BRLENS_MAX/BRLENS_MIN)*2.0; /* see Move_BrLen */
+	mt->tuningParam[1] = 0.05;
+	mt->minimum[1] = 0.002;
+	mt->maximum[1] = 0.999;
 	mt->parsimonyBased = NO;
+	mt->doAdaptParam = YES;
+	mt->whichAdaptParam = 1;
+	mt->adaptMultiplier = NO;
+	/* leaving minAdaptCFactor at default */
 	mt->level = STANDARD_USER;
 
 	/* Move_NodeSliderClock */
@@ -38976,12 +40900,15 @@
 	mt->applicableTo[0] = PINVAR_UNI;
 	mt->nApplicable = 1;
 	mt->moveFxn = &Move_Pinvar;
-	mt->relProposalProb = 1.0;
+	mt->relProposalProb = 0.5; /* was 1.0 - Allen */
 	mt->numTuningParams = 1;
 	mt->tuningParam[0] = 0.1;  /* window size */
-	mt->minimum[0] = 0.001;
+	mt->minimum[0] = 0.002; /* is mult by max of 0.5 - Allen */
 	mt->maximum[0] = 0.999;
 	mt->parsimonyBased = NO;
+	mt->doAdaptParam = YES; /* - Allen */
+	mt->whichAdaptParam = 0;
+	/* leaving minAdaptCFactor at default */
 	mt->level = STANDARD_USER;
 
 	/* Move_RanSPR1 */
@@ -39076,9 +41003,13 @@
 	mt->applicableTo[0] = RATEMULT_DIR;
 	mt->nApplicable = 1;
 	mt->moveFxn = &Move_RateMult_Dir;
-	mt->relProposalProb = 1.0;
+#ifdef NOTOPOLOGY
+	mt->relProposalProb = 5.0;
+#else
+	mt->relProposalProb = 10.0; /* was 1.0; should really be proportional to num of rate partitions... XXX - Allen */
+#endif
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 500.0; /* alphaPi */
+	mt->tuningParam[0] = 5000.0; /* alphaPi; was 500.0 - Allen */
 	mt->minimum[0] = 0.00000001;
 	mt->maximum[0] = 10000000.0;
 	mt->parsimonyBased = NO;
@@ -39092,7 +41023,7 @@
 	mt->applicableTo[0] = REVMAT_DIR;
 	mt->nApplicable = 1;
 	mt->moveFxn = &Move_Revmat_Dir;
-	mt->relProposalProb = 1.0;
+	mt->relProposalProb = 3.0; /* was 1.0 - Allen */
 	mt->numTuningParams = 1;
 	mt->tuningParam[0] = 500.0;  /* alphaPi */
 	mt->minimum[0] = 0.0001;
@@ -39164,7 +41095,7 @@
 	mt->applicableTo[0] = PI_DIR;
 	mt->nApplicable = 1;
 	mt->moveFxn = &Move_Statefreqs;
-	mt->relProposalProb = 1.0;
+	mt->relProposalProb = 0.5; /* changed from 1.0 - Allen */
 	mt->numTuningParams = 1;
 	mt->tuningParam[0] = 300.0; /* so-called "alphaPi" */
 	mt->minimum[0] = 0.001;
@@ -39196,15 +41127,19 @@
 	mt->nameTuning[0] = "sliding window size";
 	mt->shortName = "slider";
 	mt->applicableTo[0] = SWITCH_UNI;
-	mt->applicableTo[1] = SWITCH_EXP;
-	mt->nApplicable = 2;
+	mt->applicableTo[1] = SWITCH_PINVUNI;
+	mt->applicableTo[2] = SWITCH_EXP;
+	mt->nApplicable = 3;
 	mt->moveFxn = &Move_SwitchRate;
-	mt->relProposalProb = 1.0;
+	mt->relProposalProb = 3.0; /* was 1.0 - Allen */
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 1.0;  /* window size */
+	mt->tuningParam[0] = 0.0006;  /* window size */ /* was 1.0 - Allen */
 	mt->minimum[0] = 0.00001;
-	mt->maximum[0] = 1000.0;
+	mt->maximum[0] = 200.0; /* window size more sensible max */
 	mt->parsimonyBased = NO;
+	mt->doAdaptParam = YES; /* XXX - Allen */
+	mt->whichAdaptParam = 0;
+	/* leaving minAdaptCFactor at default */
 	mt->level = STANDARD_USER;
 
 	/* Move_SwitchRate_M */
@@ -39218,7 +41153,7 @@
 	mt->moveFxn = &Move_SwitchRate_M;
 	mt->relProposalProb = 0.0;
 	mt->numTuningParams = 1;
-	mt->tuningParam[0] = 2.0 * log (1.5);  /* lambda */
+	mt->tuningParam[0] = 2.0 * log (1.6);  /* lambda */
 	mt->minimum[0] = 0.00001;
 	mt->maximum[0] = 10000000.0;
 	mt->parsimonyBased = NO;
@@ -39930,7 +41865,7 @@
 	int			i, j, k, index;
 	MrBFlt		t, u, x, z, beta, bigPi_j[4], pij, bigPij,
 				*catRate, baseRate, theRate, *pis;
-	CLFlt		*tiP;
+	CLFlt2		*tiP;
 	ModelInfo	*m;
 #	if defined SSE
 	float		temp[16];
@@ -40055,7 +41990,7 @@
 	register int	i, j, k, n, s, index;
 	MrBFlt			t, *catRate, baseRate, *eigenValues, *cijk, 
 					EigValexp[64], sum, *ptr, theRate, correctionFactor;
-	CLFlt			*tiP;
+	CLFlt2			*tiP;
 	ModelInfo		*m;
 #if defined SSE
 	float			temp[16];
@@ -40199,8 +42134,8 @@
 	
 	register int	i, j, k, n, s, index;
 	int				sizeOfSingleCijk;
-	MrBFlt			t, *eigenValues, *cijk, EigValexp[64], sum, *ptr, correctionFactor;
-	CLFlt			*tiP;
+	MrBFlt			t, *eigenValues, *cijk, EigValexp[64], sum, *ptr, correctionFactor, baseRate;
+	CLFlt2			*tiP;
 	ModelInfo		*m;
 	
 	m = &modelSettings[division];
@@ -40235,10 +42170,13 @@
 	   cijk information for these models) */
 	sizeOfSingleCijk = m->nCijk / m->nCijkParts;
 
+	/* get base rate - Allen */
+	baseRate =  GetRate (division, chain);
+
 	/* fill in values */
 	for (k=index=0; k<m->nCijkParts; k++)
 		{
-		t =  p->length * correctionFactor;
+		t =  p->length * baseRate * correctionFactor;
 		for (s=0; s<n; s++)
 			EigValexp[s] =  exp(eigenValues[s] * t);
 
@@ -40294,7 +42232,7 @@
 	int			i, j, k, index;
 	MrBFlt		t, kap, u, w, x, y, z, beta, bigPi_j[4], pij, bigPij, *pis,
 				*catRate, baseRate, theRate;
-	CLFlt		*tiP;
+	CLFlt2		*tiP;
 	ModelInfo	*m;
 #if defined SSE
 	float		temp[16];
@@ -40422,7 +42360,8 @@
 	
 	int			i, j, k, index;
 	MrBFlt		*catRate, baseRate;
-	CLFlt		pNoChange, pChange, *tiP;
+	CLFlt		pNoChange, pChange;
+	CLFlt2          *tiP;
 	ModelInfo	*m;
 #if defined SSE
 	float		temp[16];
@@ -40517,7 +42456,7 @@
 	int			k, index;
 	MrBFlt		baseRate, eV, mu, theRate, v,
 				*bs, *catRate;
-	CLFlt		*tiP;
+	CLFlt2		*tiP;
 	ModelInfo	*m;
 	
 	/* find model settings for the division */
@@ -40574,7 +42513,8 @@
 	MrBFlt		v, eV1, eV2, eV3, eV4, eV5, *catRate,
 				baseRate, theRate, pi, f1, f2, f3, f4, f5, f6, f7, root, EigValexp[10],
 				*eigenValues, *cijk, sum, *bs, mu;
-	CLFlt		pNoChange, pChange, *tiP;
+	CLFlt		pNoChange, pChange;
+	CLFlt2          *tiP;
 	ModelInfo	*m;
 #	if defined (DEBUG_TIPROBS_STD)
 	int			index3;
diff -NaurdbB -xCVS -xconfigure$ -xconfig.h.in -xMakefile.in mrbayes-3.1.2.old/model.c mrbayes-3.1.2.new/model.c
--- mrbayes-3.1.2.old/model.c	Fri Dec 23 11:26:29 EST 2005
+++ mrbayes-3.1.2.new/model.c	Mon Aug 28 19:36:50 EDT 2006
@@ -364,7 +364,11 @@
 							MrBayesPrint ("%s                     Switching rates have independent uniform dist-\n", spacer);
 							MrBayesPrint ("%s                     ributions on the interval (%1.2lf,%1.2lf).\n", spacer, modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1]);
 							}
-						else if (!strcmp(modelParams[i].covSwitchPr,"Exponential"))
+						else if (!strcmp(modelParams[i].covSwitchPr,"PinvarUniform")) { /* - Allen */
+						  MrBayesPrint ("%s                     Switching rates have independent uniform distributions\n", spacer);
+						  MrBayesPrint ("%s                     on the interval (%1.2lf,%1.2lf) with a pinvar limit of %1.2lf.\n", spacer,
+								modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1],modelParams[i].covswitchMaxPinvar);
+						} else if (!strcmp(modelParams[i].covSwitchPr,"Exponential"))
 							{
 							MrBayesPrint ("%s                     Switching rates have independent exponential\n", spacer);
 							MrBayesPrint ("%s                     distributions with parameters (%1.2lf).\n", spacer, modelParams[i].covswitchExp);
@@ -574,7 +578,7 @@
 								if (!strcmp(modelParams[i].shapePr,"Uniform"))
 									{
 									MrBayesPrint ("%s                     Gamma shape parameter is uniformly dist-\n", spacer);
-									MrBayesPrint ("%s                     ributed on the interval (%1.2lf,%1.2lf).\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
+									MrBayesPrint ("%s                     ributed on the interval (%1.4lf,%1.2lf).\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
 									}
 								else if (!strcmp(modelParams[i].shapePr,"Exponential"))
 									{
@@ -583,7 +587,7 @@
 									}
 								else
 									{
-									MrBayesPrint ("%s                     Gamma shape parameter is fixed to %1.2lf.\n", spacer, modelParams[i].shapeFix);
+									MrBayesPrint ("%s                     Gamma shape parameter is fixed to %1.4lf.\n", spacer, modelParams[i].shapeFix);
 									}
 								}
 							if ((!strcmp(modelParams[i].ratesModel, "Propinv") || !strcmp(modelParams[i].ratesModel, "Invgamma")) && !strcmp(modelParams[i].covarionModel, "No"))
@@ -787,9 +791,9 @@
 					
 					if (j == P_TRATIO)
 						MrBayesPrint (" Parameter  = Tratio\n", spacer);
-					else if (j == P_REVMAT)
+					else if (j == P_REVMAT) {
 						MrBayesPrint (" Parameter  = Revmat\n", spacer);
-					else if (j == P_OMEGA)
+					} else if (j == P_OMEGA)
 						MrBayesPrint (" Parameter  = Omega\n", spacer);
 					else if (j == P_PI)
 						MrBayesPrint (" Parameter  = Statefreq\n", spacer);
@@ -834,8 +838,8 @@
 						else
 							MrBayesPrint ("%s            Prior      = Fixed(%1.2lf)\n", spacer, modelParams[i].tRatioFix);
 						}
-					else if (j == P_REVMAT)
-						{
+					else if (j == P_REVMAT) { /* Allen */
+						if (modelParams[i].dataType == DNA || modelParams[i].dataType == RNA) {
 						if (!strcmp(modelParams[i].revMatPr,"Dirichlet"))
 							MrBayesPrint ("%s            Prior      = Dirichlet(%1.2lf,%1.2lf,%1.2lf,%1.2lf,%1.2lf,%1.2lf)\n", spacer, 
 							modelParams[i].revMatDir[0], modelParams[i].revMatDir[1], modelParams[i].revMatDir[2],
@@ -844,8 +848,19 @@
 							MrBayesPrint ("%s            Prior      = Fixed(%1.2lf,%1.2lf,%1.2lf,%1.2lf,%1.2lf,%1.2lf)\n", spacer, 
 							modelParams[i].revMatFix[0], modelParams[i].revMatFix[1], modelParams[i].revMatFix[2],
 							modelParams[i].revMatFix[3], modelParams[i].revMatFix[4], modelParams[i].revMatFix[5]);
-						}
-					else if (j == P_OMEGA)
+						} else if (modelParams[i].dataType == PROTEIN) {
+							if ((!strcmp(modelParams[i].aaModelPr, "Fixed")) &&
+				    			    (!strcmp(modelParams[i].aaModel, "Gtr")) &&
+				    			    (!strcmp(modelParams[i].aaRevMatPr,"Dirichlet")))
+								MrBayesPrint ("%s            Prior      = Dirichlet(...)\n", spacer);
+							else if (!strcmp(modelParams[i].aaModelPr, "Mixed"))
+								MrBayesPrint ("%s            Prior      = Mixed (see below)\n", spacer);
+							else
+								MrBayesPrint ("%s            Prior      = Fixed(%s)\n",spacer,
+									     modelParams[i].aaModel);
+						} else
+							MrBayesPrint("%s            Prior      = %s\n",spacer,modelParams[i].revMatPr);  
+					} else if (j == P_OMEGA)
 						{
 						if (!strcmp(modelParams[i].omegaVar,"Equal"))
 							{
@@ -987,11 +1002,11 @@
 					else if (j == P_SHAPE)
 						{
 						if (!strcmp(modelParams[i].shapePr,"Uniform"))
-							MrBayesPrint ("%s            Prior      = Uniform(%1.2lf,%1.2lf)\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
+							MrBayesPrint ("%s            Prior      = Uniform(%1.4lf,%1.2lf)\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
 						else if (!strcmp(modelParams[i].shapePr,"Exponential"))
 							MrBayesPrint ("%s            Prior      = Exponential(%1.2lf)\n", spacer, modelParams[i].shapeExp);
 						else
-							MrBayesPrint ("%s            Prior      = Fixed(%1.2lf)\n", spacer, modelParams[i].shapeFix);
+							MrBayesPrint ("%s            Prior      = Fixed(%1.4lf)\n", spacer, modelParams[i].shapeFix);
 						}
 					else if (j == P_PINVAR)
 						{
@@ -1011,6 +1026,9 @@
 						{
 						if (!strcmp(modelParams[i].covSwitchPr,"Uniform"))
 							MrBayesPrint ("%s            Prior      = Uniform(%1.2lf,%1.2lf)\n", spacer, modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1]);
+						else if (!strcmp(modelParams[i].covSwitchPr,"PinvarUniform"))
+							MrBayesPrint ("%s            Prior      = PinvarUniform(%1.2lf,%1.2lf,%1.2lf)\n",
+								      spacer, modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1], modelParams[i].covswitchMaxPinvar);
 						else if (!strcmp(modelParams[i].covSwitchPr,"Exponential"))
 							MrBayesPrint ("%s            Prior      = Exponential(%1.2lf)\n", spacer, modelParams[i].covswitchExp);
 						else
@@ -3231,13 +3249,13 @@
 									}
 								if (modelParams[i].shapeUni[0] < MIN_SHAPE_PARAM)
 									{
-									MrBayesPrint ("%s   Lower value for uniform cannot be less than %1.2lf\n", spacer, MIN_SHAPE_PARAM);
+									MrBayesPrint ("%s   Lower value for uniform cannot be less than %1.4lf\n", spacer, MIN_SHAPE_PARAM);
 									return (ERROR);
 									}
 								if (nApplied == 0 && numCurrentDivisions == 1)
-									MrBayesPrint ("%s   Setting Shapepr to Uniform(%1.2lf,%1.2lf)\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
+									MrBayesPrint ("%s   Setting Shapepr to Uniform(%1.4lf,%1.2lf)\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1]);
 								else
-									MrBayesPrint ("%s   Setting Shapepr to Uniform(%1.2lf,%1.2lf) for partition %d\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1], i+1);
+									MrBayesPrint ("%s   Setting Shapepr to Uniform(%1.4lf,%1.2lf) for partition %d\n", spacer, modelParams[i].shapeUni[0], modelParams[i].shapeUni[1], i+1);
 								expecting  = Expecting(RIGHTPAR);
 								}
 							}
@@ -3262,13 +3280,13 @@
 								}
 							if (modelParams[i].shapeFix < MIN_SHAPE_PARAM)
 								{
-								MrBayesPrint ("%s   Shape parameter cannot be less than %1.2lf\n", spacer, MIN_SHAPE_PARAM);
+								MrBayesPrint ("%s   Shape parameter cannot be less than %1.4lf\n", spacer, MIN_SHAPE_PARAM);
 								return (ERROR);
 								}
 							if (nApplied == 0 && numCurrentDivisions == 1)
-								MrBayesPrint ("%s   Setting Shapepr to Fixed(%1.2lf)\n", spacer, modelParams[i].shapeFix);
+								MrBayesPrint ("%s   Setting Shapepr to Fixed(%1.4lf)\n", spacer, modelParams[i].shapeFix);
 							else
-								MrBayesPrint ("%s   Setting Shapepr to Fixed(%1.2lf) for partition %d\n", spacer, modelParams[i].shapeFix, i+1);
+								MrBayesPrint ("%s   Setting Shapepr to Fixed(%1.4lf) for partition %d\n", spacer, modelParams[i].shapeFix, i+1);
 							expecting  = Expecting(RIGHTPAR);
 							}
 						}
@@ -3296,8 +3314,8 @@
 					{
 					nApplied = NumActiveParts ();
 					for (i=0; i<numCurrentDivisions; i++)
-						{
-						if ((activeParts[i] == YES || nApplied == 0) && (modelParams[i].dataType == DNA || modelParams[i].dataType == RNA))
+						{ /* Had lacked PROTEIN as possible; not adding RESTRICTION or STANDARD since am not familiar with details of those - Allen */
+						if ((activeParts[i] == YES || nApplied == 0) && (modelParams[i].dataType == DNA || modelParams[i].dataType == RNA || modelParams[i].dataType == PROTEIN))
 							strcpy(modelParams[i].pInvarPr, tempStr);
 						}
 					}
@@ -3318,8 +3336,8 @@
 				{
 				nApplied = NumActiveParts ();
 				for (i=0; i<numCurrentDivisions; i++)
-					{
-					if ((activeParts[i] == YES || nApplied == 0) && (modelParams[i].dataType == DNA || modelParams[i].dataType == RNA))
+					{ /* Had lacked PROTEIN as possible; not adding RESTRICTION or STANDARD since am not familiar with details of those - Allen */
+					if ((activeParts[i] == YES || nApplied == 0) && (modelParams[i].dataType == DNA || modelParams[i].dataType == RNA || modelParams[i].dataType == PROTEIN))
 						{
 						if (!strcmp(modelParams[i].pInvarPr,"Uniform"))
 							{
@@ -3762,7 +3780,7 @@
 								{
 								if (modelParams[i].covswitchUni[0] >= modelParams[i].covswitchUni[1])
 									{
-									MrBayesPrint ("%s   Lower value for uniform should be greater than upper value\n", spacer);
+									MrBayesPrint ("%s   Lower value for uniform should be less than upper value\n", spacer);
 									return (ERROR);
 									}
 								if (nApplied == 0 && numCurrentDivisions == 1)
@@ -3772,7 +3790,40 @@
 								expecting  = Expecting(RIGHTPAR);
 								}
 							}
-						else if (!strcmp(modelParams[i].covSwitchPr,"Exponential"))
+						else if (!strcmp(modelParams[i].covSwitchPr,"PinvarUniform")) { /* - Allen */
+						  sscanf (tkn, "%lf", &tempD);
+						  if (numVars[i] > 1)
+						    modelParams[i].covswitchMaxPinvar = tempD;
+                                                  else
+						    modelParams[i].covswitchUni[numVars[i]] = tempD;
+                                                  numVars[i]++;
+						  if (numVars[i] < 3)
+						    expecting  = Expecting(COMMA);
+						  else {
+						    if (modelParams[i].covswitchUni[0] >= modelParams[i].covswitchUni[1]) {
+						      MrBayesPrint ("%s   Lower value for uniform should be less than upper value\n", spacer);
+						      return (ERROR);
+						    } else if (modelParams[i].covswitchMaxPinvar > 1) {
+						      MrBayesPrint ("%s   Covswitchpr max pinvar must be 1 or less\n", spacer);
+						      return (ERROR);
+						    } else if (modelParams[i].covswitchMaxPinvar <= 0) {
+						      MrBayesPrint ("%s   Covswitchpr max pinvar must be greater than 0\n", spacer);
+						      return (ERROR);
+						    } else if (modelParams[i].covswitchMaxPinvar <= (modelParams[i].covswitchUni[0]/(modelParams[i].covswitchUni[0]+modelParams[i].covswitchUni[1]))) {
+						      MrBayesPrint ("%s   Covswitchpr max pinvar must be above %1.4lf with uniform %1.2lf to %1.2lf\n",
+								    spacer, (modelParams[i].covswitchUni[0]/(modelParams[i].covswitchUni[0]+modelParams[i].covswitchUni[1])), modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1]);
+						      return (ERROR);
+						    } else if (modelParams[i].covswitchMaxPinvar == 1)
+						      MrBayesPrint ("%s   Covswitchpr max pinvar of 1 has no effect\n", spacer);
+						    if (nApplied == 0 && numCurrentDivisions == 1)
+						      MrBayesPrint ("%s   Setting Covswitchpr to PinvarUniform(%1.2lf,%1.2lf,%1.4lf)\n",
+								    spacer, modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1], modelParams[i].covswitchMaxPinvar);
+						    else
+						      MrBayesPrint ("%s   Setting Covswitchpr to PinvarUniform(%1.2lf,%1.2lf,%1.4lf) for partition %d\n", spacer,
+								    modelParams[i].covswitchUni[0], modelParams[i].covswitchUni[1], modelParams[i].covswitchMaxPinvar, i+1);
+						    expecting  = Expecting(RIGHTPAR);
+						  }
+						} else if (!strcmp(modelParams[i].covSwitchPr,"Exponential"))
 							{
 							sscanf (tkn, "%lf", &tempD);
 							modelParams[i].covswitchExp = tempD;
@@ -4100,7 +4151,10 @@
 										sum += tempStateFreqs[j];
 									if (AreDoublesEqual (sum, (MrBFlt) 1.0, (MrBFlt) 0.001) == NO)
 										{
-										MrBayesPrint ("%s   State frequencies do not sum to 1.0\n", spacer);
+										  if (nApplied == 0 && numCurrentDivisions == 1)
+										    MrBayesPrint ("%s   State frequencies do not sum to 1.0 (%f)\n", spacer, sum);
+										  else
+										    MrBayesPrint ("%s   State frequencies do not sum to 1.0 for partition %d (%f)\n", spacer, i+1, sum);
 										return (ERROR);
 										}
 									strcpy(modelParams[i].stateFreqsFixType, "User");
@@ -6138,9 +6192,9 @@
 			*isApplic2 = NO;
 
 		/* check that data type is the same for both partitions */
-		if (isFirstNucleotide == YES && isSecondNucleotide == NO)
+		if (isFirstNucleotide != isSecondNucleotide)
 			isSame = NO;
-		if (isFirstProtein == YES && isSecondProtein == NO)
+		if (isFirstProtein != isSecondProtein)
 			isSame = NO;
 
 		/* GTR applies to both part1 and part2. We now need to check if the prior is the same for both. */
@@ -6764,8 +6818,15 @@
 				isSame = NO;
 			if (AreDoublesEqual (modelParams[part1].covswitchUni[1], modelParams[part2].covswitchUni[1], (MrBFlt) 0.00001) == NO)
 				isSame = NO;
-			}
-		else if (!strcmp(modelParams[part1].covSwitchPr,"Exponential") && !strcmp(modelParams[part2].covSwitchPr,"Exponential"))
+		} else if (!strcmp(modelParams[part1].covSwitchPr,"PinvarUniform") && !strcmp(modelParams[part2].covSwitchPr,"PinvarUniform")) {
+		  /* We assume the person has not set covswitchMaxPinvar to 1, which makes PinvarUniform equivalent to Uniform - Allen */
+		  if (AreDoublesEqual (modelParams[part1].covswitchUni[0], modelParams[part2].covswitchUni[0], (MrBFlt) 0.00001) == NO)
+		    isSame = NO;
+		  else if (AreDoublesEqual (modelParams[part1].covswitchUni[1], modelParams[part2].covswitchUni[1], (MrBFlt) 0.00001) == NO)
+		    isSame = NO;
+		  else if (AreDoublesEqual (modelParams[part1].covswitchMaxPinvar, modelParams[part2].covswitchMaxPinvar, (MrBFlt) 0.00001) == NO)
+		    isSame = NO;
+		} else if (!strcmp(modelParams[part1].covSwitchPr,"Exponential") && !strcmp(modelParams[part2].covSwitchPr,"Exponential"))
 			{
 			if (AreDoublesEqual (modelParams[part1].covswitchExp, modelParams[part2].covswitchExp, (MrBFlt) 0.00001) == NO)
 				isSame = NO;
