use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class InitMigrationModelConnector method customConnector.
public static boolean customConnector(BeautiDoc doc) {
for (BEASTInterface p : doc.getPartitions("Tree")) {
TreeLikelihood treeLikelihood = (TreeLikelihood) p;
StructuredCoalescentMultiTypeTree tree = (StructuredCoalescentMultiTypeTree) treeLikelihood.treeInput.get();
String pID = BeautiDoc.parsePartition(tree.getID());
SCMigrationModel migModel = (SCMigrationModel) doc.pluginmap.get("migModel.t:" + pID);
SCMigrationModel migModelInit = (SCMigrationModel) doc.pluginmap.get("migModelInit.t:" + pID);
String rateMatrixStr = getParameterString((RealParameter) migModel.rateMatrixInput.get());
String popSizesStr = getParameterString((RealParameter) migModel.popSizesInput.get());
// Ensure model has appropriate number of demes
int uniqueTraitCount = uniqueTraitsInData(tree).size();
StringBuilder rateMatrixStrBuilder = new StringBuilder();
StringBuilder popSizesStrBuilder = new StringBuilder();
migModel.getTypeSet().initAndValidate();
if (migModel.popSizesInput.get().getDimension() != migModel.getNTypes()) {
for (int i = 0; i < migModel.getNTypes(); i++) {
popSizesStrBuilder.append(" 1.0");
for (int j = 0; j < migModel.getNTypes(); j++) {
if (j == i)
continue;
rateMatrixStrBuilder.append(" 1.0");
}
}
popSizesStr = popSizesStrBuilder.toString();
rateMatrixStr = rateMatrixStrBuilder.toString();
((RealParameter) migModel.popSizesInput.get()).setDimension(migModel.getNTypes());
((RealParameter) migModel.popSizesInput.get()).valuesInput.setValue(popSizesStr, (RealParameter) migModel.popSizesInput.get());
((RealParameter) migModel.rateMatrixInput.get()).setDimension(migModel.getNTypes() * (migModel.getNTypes() - 1));
((RealParameter) migModel.rateMatrixInput.get()).valuesInput.setValue(rateMatrixStr, (RealParameter) migModel.rateMatrixInput.get());
((RealParameter) migModel.popSizesInput.get()).initAndValidate();
((RealParameter) migModel.rateMatrixInput.get()).initAndValidate();
migModel.initAndValidate();
}
((RealParameter) migModelInit.popSizesInput.get()).setDimension(migModel.getNTypes());
((RealParameter) migModelInit.popSizesInput.get()).valuesInput.setValue(popSizesStr, (RealParameter) migModelInit.popSizesInput.get());
((RealParameter) migModelInit.rateMatrixInput.get()).setDimension(migModel.getNTypes() * (migModel.getNTypes() - 1));
((RealParameter) migModelInit.rateMatrixInput.get()).valuesInput.setValue(rateMatrixStr, (RealParameter) migModelInit.rateMatrixInput.get());
try {
((RealParameter) migModelInit.popSizesInput.get()).initAndValidate();
((RealParameter) migModelInit.rateMatrixInput.get()).initAndValidate();
migModelInit.initAndValidate();
} catch (Exception ex) {
System.err.println("Error configuring initial migration model.");
}
}
return false;
}
use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class StructuredCoalescentTreeDensity method main.
/**
* Test likelihood result. Duplicate of JUnit test for debugging purposes.
*
* @param argv
* @throws java.lang.Exception
*/
public static void main(String[] argv) throws Exception {
// Assemble test MultiTypeTree:
String newickStr = "(((A[state=1]:0.25)[state=0]:0.25,B[state=0]:0.5)[state=0]:1.5," + "(C[state=0]:1.0,D[state=0]:1.0)[state=0]:1.0)[state=0]:0.0;";
MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
mtTree.initByName("newick", newickStr, "typeLabel", "state", "nTypes", 2);
// Assemble migration model:
RealParameter rateMatrix = new RealParameter();
rateMatrix.initByName("value", "2.0 1.0");
RealParameter popSizes = new RealParameter();
popSizes.initByName("value", "5.0 10.0");
SCMigrationModel migrationModel = new SCMigrationModel();
migrationModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes);
// Set up likelihood instance:
StructuredCoalescentTreeDensity likelihood = new StructuredCoalescentTreeDensity();
likelihood.initByName("migrationModel", migrationModel, "multiTypeTree", mtTree);
// Calculated by hand
double expResult = -16.52831;
double result = likelihood.calculateLogP();
System.out.println("Expected result: " + expResult);
System.out.println("Actual result: " + result);
System.out.println("Difference: " + String.valueOf(result - expResult));
}
use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class Ewing_Test method test2.
@Test
public void test2() throws Exception {
System.out.println("Ewing_test 2");
// Fix seed.
Randomizer.setSeed(42);
// Assemble initial MultiTypeTree
String newickStr = "(((1[&deme=1]:0.5)[&deme=0]:0.5,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
mtTree.initByName("value", newickStr, "typeLabel", "deme");
// Assemble migration model:
RealParameter rateMatrix = new RealParameter("0.1 0.1");
RealParameter popSizes = new RealParameter("7.0 7.0");
SCMigrationModel migModel = new SCMigrationModel();
migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
// Assemble distribution:
StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree, "checkValidity", true);
// Set up state:
State state = new State();
state.initByName("stateNode", mtTree);
// Set up operators:
TypedWilsonBaldingEasy twbOperator = new TypedWilsonBaldingEasy();
twbOperator.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
TypedSubtreeExchangeEasy tsxOperator = new TypedSubtreeExchangeEasy();
tsxOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "isNarrow", true);
MultiTypeUniform mtuOperator = new MultiTypeUniform();
mtuOperator.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
MultiTypeTreeScale mttsOperator = new MultiTypeTreeScale();
mttsOperator.initByName("weight", 1.0, "scaleFactor", 0.8, "useOldTreeScaler", true, "migrationModel", migModel, "multiTypeTree", mtTree);
TypePairBirthDeath tpbdOperator = new TypePairBirthDeath();
tpbdOperator.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
TypeMergeSplit tmsOperator = new TypeMergeSplit();
tmsOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "includeRoot", true);
// Set up stat analysis logger:
MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
// Set up MCMC:
MCMC mcmc = new MCMC();
mcmc.initByName("chainLength", "10000000", "state", state, "distribution", distribution, "operator", twbOperator, "operator", tsxOperator, "operator", mtuOperator, "operator", mttsOperator, "operator", tpbdOperator, "operator", tmsOperator, "logger", logger);
// Run MCMC:
mcmc.run();
System.out.format("height mean = %s\n", logger.getHeightMean());
System.out.format("height var = %s\n", logger.getHeightVar());
System.out.format("height ESS = %s\n", logger.getHeightESS());
// Compare analysis results with truth:
boolean withinTol = (logger.getHeightESS() > 1000) && (Math.abs(logger.getHeightMean() - 23) < 0.5) && (Math.abs(logger.getHeightVar() - 300) < 30.0);
Assert.assertTrue(withinTol);
}
use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class Ewing_Test method test1.
@Test
public void test1() throws Exception {
System.out.println("Ewing_test1");
// Fix seed.
Randomizer.setSeed(1);
// Assemble initial MultiTypeTree
String newickStr = "((1[&deme=0]:1,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
mtTree.initByName("value", newickStr, "typeLabel", "deme");
// Assemble migration model:
RealParameter rateMatrix = new RealParameter("0.1 0.1");
RealParameter popSizes = new RealParameter("7.0 7.0");
SCMigrationModel migModel = new SCMigrationModel();
migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
// Assemble distribution:
StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
// Set up state:
State state = new State();
state.initByName("stateNode", mtTree);
// Set up operators:
TypedWilsonBaldingEasy twbOperator = new TypedWilsonBaldingEasy();
twbOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
TypedSubtreeExchangeEasy tsxOperator = new TypedSubtreeExchangeEasy();
tsxOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "isNarrow", true);
MultiTypeUniform mtuOperator = new MultiTypeUniform();
mtuOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
MultiTypeTreeScale mttsOperator = new MultiTypeTreeScale();
mttsOperator.initByName("weight", 1.0, "scaleFactor", 0.8, "useOldTreeScaler", true, "multiTypeTree", mtTree, "migrationModel", migModel);
TypePairBirthDeath tpbdOperator = new TypePairBirthDeath();
tpbdOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
TypeMergeSplit tmsOperator = new TypeMergeSplit();
tmsOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "includeRoot", true);
// Set up stat analysis logger:
MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
// Logger fileLogger = new Logger();
// TreeHeightLogger thLogger = new TreeHeightLogger();
// thLogger.initByName("tree", mtTree);
// fileLogger.initByName(
// "fileName", "test.log",
// "logEvery", 1000,
// "model", distribution,
// "log", thLogger);
// Set up MCMC:
MCMC mcmc = new MCMC();
mcmc.initByName("chainLength", "10000000", "state", state, "distribution", distribution, "operator", twbOperator, "operator", tsxOperator, "operator", mtuOperator, "operator", mttsOperator, "operator", tpbdOperator, "operator", tmsOperator, "logger", logger);
// Run MCMC:
mcmc.run();
System.out.format("height mean = %s\n", logger.getHeightMean());
System.out.format("height var = %s\n", logger.getHeightVar());
System.out.format("height ESS = %s\n", logger.getHeightESS());
// Compare analysis results with truth:
boolean withinTol = (logger.getHeightESS() > 3000) && (Math.abs(logger.getHeightMean() - 19.15) < 0.5) && (Math.abs(logger.getHeightVar() - 310) < 20);
Assert.assertTrue(withinTol);
}
use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class STX_NR_MTU_TS_Test method test.
@Test
public void test() throws Exception {
System.out.println("STX_NR_MTU_TS test");
// Test passing locally, not on Travis. WHY!?
// Fix seed.
Randomizer.setSeed(53);
// Assemble migration model:
RealParameter rateMatrix = new RealParameter("0.1 0.1");
RealParameter popSizes = new RealParameter("7.0 7.0");
SCMigrationModel migModel = new SCMigrationModel();
migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
// Assemble initial MultiTypeTree
MultiTypeTree mtTree = new StructuredCoalescentMultiTypeTree();
mtTree.initByName("typeLabel", "deme", "migrationModel", migModel, "leafTypes", "1 1 0 0");
// Set up state:
State state = new State();
state.initByName("stateNode", mtTree);
// Assemble distribution:
StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
// Set up operators:
Operator operatorSTX = new TypedSubtreeExchange();
operatorSTX.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
Operator operatorNR = new NodeRetype();
operatorNR.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
Operator operatorMTU = new MultiTypeUniform();
operatorMTU.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
Operator operatorMTTS = new MultiTypeTreeScale();
operatorMTTS.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "scaleFactor", 1.5, "useOldTreeScaler", false);
// Set up stat analysis logger:
MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
// Set up MCMC:
MCMC mcmc = new MCMC();
mcmc.initByName("chainLength", "1000000", "state", state, "distribution", distribution, "operator", operatorSTX, "operator", operatorNR, "operator", operatorMTU, "operator", operatorMTTS, "logger", logger);
// Run MCMC:
mcmc.run();
System.out.format("height mean = %s\n", logger.getHeightMean());
System.out.format("height var = %s\n", logger.getHeightVar());
System.out.format("height ESS = %s\n", logger.getHeightESS());
// Direct simulation:
double[] heights = UtilMethods.getSimulatedHeights(migModel, new IntegerParameter("1 1 0 0"));
double simHeightMean = DiscreteStatistics.mean(heights);
double simHeightVar = DiscreteStatistics.variance(heights);
// Compare analysis results with truth:
boolean withinTol = (logger.getHeightESS() > 500) && (Math.abs(logger.getHeightMean() - simHeightMean) < 2.0) && (Math.abs(logger.getHeightVar() - simHeightVar) < 50);
Assert.assertTrue(withinTol);
}
Aggregations