use of dr.inference.operators.AdaptableMCMCOperator in project beast-mcmc by beast-dev.
the class BeastCheckpointer method readStateFromFile.
protected long readStateFromFile(File file, MarkovChain markovChain, double[] lnL) {
DoubleParser parser = useFullPrecision ? DoubleParser.HEX : DoubleParser.TEXT;
OperatorSchedule operatorSchedule = markovChain.getSchedule();
long state = -1;
ArrayList<TreeParameterModel> traitModels = new ArrayList<TreeParameterModel>();
try {
FileReader fileIn = new FileReader(file);
BufferedReader in = new BufferedReader(fileIn);
int[] rngState = null;
String line = in.readLine();
String[] fields = line.split("\t");
if (fields[0].equals("rng")) {
// if there is a random number generator state present then load it...
try {
rngState = new int[fields.length - 1];
for (int i = 0; i < rngState.length; i++) {
rngState[i] = Integer.parseInt(fields[i + 1]);
}
} catch (NumberFormatException nfe) {
throw new RuntimeException("Unable to read state number from state file");
}
line = in.readLine();
fields = line.split("\t");
}
try {
if (!fields[0].equals("state")) {
throw new RuntimeException("Unable to read state number from state file");
}
state = Long.parseLong(fields[1]);
} catch (NumberFormatException nfe) {
throw new RuntimeException("Unable to read state number from state file");
}
line = in.readLine();
fields = line.split("\t");
try {
if (!fields[0].equals("lnL")) {
throw new RuntimeException("Unable to read lnL from state file");
}
if (lnL != null) {
lnL[0] = parser.parseDouble(fields[1]);
}
} catch (NumberFormatException nfe) {
throw new RuntimeException("Unable to read lnL from state file");
}
for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
if (!parameter.isImmutable()) {
line = in.readLine();
fields = line.split("\t");
// if (!fields[0].equals(parameter.getParameterName())) {
// System.err.println("Unable to match state parameter: " + fields[0] + ", expecting " + parameter.getParameterName());
// }
int dimension = Integer.parseInt(fields[2]);
if (dimension != parameter.getDimension()) {
System.err.println("Unable to match state parameter dimension: " + dimension + ", expecting " + parameter.getDimension() + " for parameter: " + parameter.getParameterName());
System.err.print("Read from file: ");
for (int i = 0; i < fields.length; i++) {
System.err.print(fields[i] + "\t");
}
System.err.println();
}
if (fields[1].equals("branchRates.categories.rootNodeNumber")) {
// System.out.println("eek");
double value = parser.parseDouble(fields[3]);
parameter.setParameterValue(0, value);
if (DEBUG) {
System.out.println("restoring " + fields[1] + " with value " + value);
}
} else {
if (DEBUG) {
System.out.print("restoring " + fields[1] + " with values ");
}
for (int dim = 0; dim < parameter.getDimension(); dim++) {
try {
parameter.setParameterUntransformedValue(dim, parser.parseDouble(fields[dim + 3]));
} catch (RuntimeException rte) {
System.err.println(rte);
continue;
}
if (DEBUG) {
System.out.print(parser.parseDouble(fields[dim + 3]) + " ");
}
}
if (DEBUG) {
System.out.println();
}
}
}
}
for (int i = 0; i < operatorSchedule.getOperatorCount(); i++) {
// TODO we can no longer assume these are in the right order
// TODO best parse all the "operator" lines and store them so we can mix and match within this for loop
// TODO does not only apply to the operators but also to the parameters
// TODO test using additional tip-date sampling compared to previous run
MCMCOperator operator = operatorSchedule.getOperator(i);
line = in.readLine();
fields = line.split("\t");
if (!fields[1].equals(operator.getOperatorName())) {
throw new RuntimeException("Unable to match " + operator.getOperatorName() + " operator: " + fields[1]);
}
if (fields.length < 4) {
throw new RuntimeException("Operator missing values: " + fields[1]);
}
operator.setAcceptCount(Integer.parseInt(fields[2]));
operator.setRejectCount(Integer.parseInt(fields[3]));
if (operator instanceof AdaptableMCMCOperator) {
if (fields.length != 6) {
throw new RuntimeException("Coercable operator missing parameter: " + fields[1]);
}
((AdaptableMCMCOperator) operator).setAdaptableParameter(parser.parseDouble(fields[4]));
((AdaptableMCMCOperator) operator).setAdaptationCount(Long.parseLong(fields[5]));
}
}
// load the tree models last as we get the node heights from the tree (not the parameters which
// which may not be associated with the right node
Set<String> expectedTreeModelNames = new HashSet<String>();
// store list of TreeModels for debugging purposes
ArrayList<TreeModel> treeModelList = new ArrayList<TreeModel>();
for (Model model : Model.CONNECTED_MODEL_SET) {
if (model instanceof TreeModel) {
if (DEBUG) {
System.out.println("model " + model.getModelName());
}
treeModelList.add((TreeModel) model);
expectedTreeModelNames.add(model.getModelName());
if (DEBUG) {
System.out.println("\nexpectedTreeModelNames:");
for (String s : expectedTreeModelNames) {
System.out.println(s);
}
System.out.println();
}
}
// first add all TreeParameterModels to a list
if (model instanceof TreeParameterModel) {
traitModels.add((TreeParameterModel) model);
}
}
// explicitly link TreeModel (using its unique ID) to a list of TreeParameterModels
// this information is currently not yet used
HashMap<String, ArrayList<TreeParameterModel>> linkedModels = new HashMap<String, ArrayList<TreeParameterModel>>();
for (String name : expectedTreeModelNames) {
ArrayList<TreeParameterModel> tpmList = new ArrayList<TreeParameterModel>();
for (TreeParameterModel tpm : traitModels) {
if (tpm.getTreeModel().getId().equals(name)) {
tpmList.add(tpm);
if (DEBUG) {
System.out.println("TreeModel: " + name + " has been assigned TreeParameterModel: " + tpm.toString());
}
}
}
linkedModels.put(name, tpmList);
}
line = in.readLine();
fields = line.split("\t");
// Read in all (possibly more than one) trees
while (fields[0].equals("tree")) {
if (DEBUG) {
System.out.println("\ntree: " + fields[1]);
}
for (Model model : Model.CONNECTED_MODEL_SET) {
if (model instanceof TreeModel && fields[1].equals(model.getModelName())) {
line = in.readLine();
line = in.readLine();
fields = line.split("\t");
// read number of nodes
int nodeCount = Integer.parseInt(fields[0]);
double[] nodeHeights = new double[nodeCount];
String[] taxaNames = new String[(nodeCount + 1) / 2];
for (int i = 0; i < nodeCount; i++) {
line = in.readLine();
fields = line.split("\t");
nodeHeights[i] = parser.parseDouble(fields[1]);
if (i < taxaNames.length) {
taxaNames[i] = fields[2];
}
}
// on to reading edge information
line = in.readLine();
line = in.readLine();
line = in.readLine();
fields = line.split("\t");
int edgeCount = Integer.parseInt(fields[0]);
if (DEBUG) {
System.out.println("edge count = " + edgeCount);
}
// create data matrix of doubles to store information from list of TreeParameterModels
// size of matrix depends on the number of TreeParameterModels assigned to a TreeModel
double[][] traitValues = new double[linkedModels.get(model.getId()).size()][edgeCount];
// create array to store whether a node is left or right child of its parent
// can be important for certain tree transition kernels
int[] childOrder = new int[edgeCount];
for (int i = 0; i < childOrder.length; i++) {
childOrder[i] = -1;
}
int[] parents = new int[edgeCount];
for (int i = 0; i < edgeCount; i++) {
parents[i] = -1;
}
for (int i = 0; i < edgeCount - 1; i++) {
line = in.readLine();
if (line != null) {
if (DEBUG) {
System.out.println("DEBUG: " + line);
}
fields = line.split("\t");
parents[Integer.parseInt(fields[0])] = Integer.parseInt(fields[1]);
// childOrder[i] = Integer.parseInt(fields[2]);
childOrder[Integer.parseInt(fields[0])] = Integer.parseInt(fields[2]);
for (int j = 0; j < linkedModels.get(model.getId()).size(); j++) {
// traitValues[j][i] = parser.parseDouble(fields[3+j]);
traitValues[j][Integer.parseInt(fields[0])] = parser.parseDouble(fields[3 + j]);
}
}
}
// perform magic with the acquired information
if (DEBUG) {
System.out.println("adopting tree structure");
}
// adopt the loaded tree structure;
((TreeModel) model).beginTreeEdit();
((TreeModel) model).adoptTreeStructure(parents, nodeHeights, childOrder, taxaNames);
if (traitModels.size() > 0) {
System.out.println("adopting " + traitModels.size() + " trait models to treeModel " + ((TreeModel) model).getId());
((TreeModel) model).adoptTraitData(parents, traitModels, traitValues, taxaNames);
}
((TreeModel) model).endTreeEdit();
expectedTreeModelNames.remove(model.getModelName());
}
}
line = in.readLine();
if (line != null) {
fields = line.split("\t");
}
}
if (expectedTreeModelNames.size() > 0) {
StringBuilder sb = new StringBuilder();
for (String notFoundName : expectedTreeModelNames) {
sb.append("Expecting, but unable to match state parameter:" + notFoundName + "\n");
}
throw new RuntimeException("\n" + sb.toString());
}
if (DEBUG) {
System.out.println("\nDouble checking:");
for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
if (parameter.getParameterName().equals("branchRates.categories.rootNodeNumber")) {
System.out.println(parameter.getParameterName() + ": " + parameter.getParameterValue(0));
}
}
System.out.println("\nPrinting trees:");
for (TreeModel tm : treeModelList) {
System.out.println(tm.getId() + ": ");
System.out.println(tm.getNewick());
}
}
if (System.getProperty(BeastCheckpointer.CHECKPOINT_SEED) != null) {
MathUtils.setSeed(Long.parseLong(System.getProperty(BeastCheckpointer.CHECKPOINT_SEED)));
} else if (rngState != null) {
MathUtils.setRandomState(rngState);
}
in.close();
fileIn.close();
// This shouldn't be necessary and if it is then it might be hiding a bug...
/*for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
likelihood.makeDirty();
}*/
} catch (IOException ioe) {
throw new RuntimeException("Unable to read file: " + ioe.getMessage());
}
return state;
}
use of dr.inference.operators.AdaptableMCMCOperator in project beast-mcmc by beast-dev.
the class BeastCheckpointer method writeStateToFile.
protected boolean writeStateToFile(File file, long state, double lnL, MarkovChain markovChain) {
OperatorSchedule operatorSchedule = markovChain.getSchedule();
OutputStream fileOut = null;
try {
fileOut = new FileOutputStream(file);
PrintStream out = useFullPrecision ? new CheckpointPrintStream(fileOut) : new PrintStream(fileOut);
ArrayList<TreeParameterModel> traitModels = new ArrayList<TreeParameterModel>();
int[] rngState = MathUtils.getRandomState();
out.print("rng");
for (int i = 0; i < rngState.length; i++) {
out.print("\t");
out.print(rngState[i]);
}
out.println();
out.print("state\t");
out.println(state);
out.print("lnL\t");
out.println(lnL);
for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
if (!parameter.isImmutable()) {
out.print("parameter");
out.print("\t");
out.print(parameter.getParameterName());
out.print("\t");
out.print(parameter.getDimension());
for (int dim = 0; dim < parameter.getDimension(); dim++) {
out.print("\t");
out.print(parameter.getParameterUntransformedValue(dim));
}
out.print("\n");
}
}
for (int i = 0; i < operatorSchedule.getOperatorCount(); i++) {
MCMCOperator operator = operatorSchedule.getOperator(i);
out.print("operator");
out.print("\t");
out.print(operator.getOperatorName());
out.print("\t");
out.print(operator.getAcceptCount());
out.print("\t");
out.print(operator.getRejectCount());
if (operator instanceof AdaptableMCMCOperator) {
out.print("\t");
out.print(((AdaptableMCMCOperator) operator).getAdaptableParameter());
out.print("\t");
out.print(((AdaptableMCMCOperator) operator).getAdaptationCount());
}
out.println();
}
// check up front if there are any TreeParameterModel objects
for (Model model : Model.CONNECTED_MODEL_SET) {
if (model instanceof TreeParameterModel) {
// System.out.println("\nDetected TreeParameterModel: " + ((TreeParameterModel) model).toString());
traitModels.add((TreeParameterModel) model);
}
}
for (Model model : Model.CONNECTED_MODEL_SET) {
if (model instanceof TreeModel) {
out.print("tree");
out.print("\t");
out.println(model.getModelName());
// replace Newick format by printing general graph structure
// out.println(((TreeModel) model).getNewick());
out.println("#node height taxon");
int nodeCount = ((TreeModel) model).getNodeCount();
out.println(nodeCount);
for (int i = 0; i < nodeCount; i++) {
out.print(((TreeModel) model).getNode(i).getNumber());
out.print("\t");
out.print(((TreeModel) model).getNodeHeight(((TreeModel) model).getNode(i)));
if (((TreeModel) model).isExternal(((TreeModel) model).getNode(i))) {
out.print("\t");
out.print(((TreeModel) model).getNodeTaxon(((TreeModel) model).getNode(i)).getId());
}
out.println();
}
out.println("#edges");
out.println("#child-node parent-node L/R-child traits");
out.println(nodeCount);
for (int i = 0; i < nodeCount; i++) {
NodeRef parent = ((TreeModel) model).getParent(((TreeModel) model).getNode(i));
if (parent != null) {
out.print(((TreeModel) model).getNode(i).getNumber());
out.print("\t");
out.print(((TreeModel) model).getParent(((TreeModel) model).getNode(i)).getNumber());
out.print("\t");
if ((((TreeModel) model).getChild(parent, 0) == ((TreeModel) model).getNode(i))) {
// left child
out.print(0);
} else if ((((TreeModel) model).getChild(parent, 1) == ((TreeModel) model).getNode(i))) {
// right child
out.print(1);
} else {
throw new RuntimeException("Operation currently only supported for nodes with 2 children.");
}
// only print the TreeParameterModel that matches the TreeModel currently being written
for (TreeParameterModel tpm : traitModels) {
if (model == tpm.getTreeModel()) {
out.print("\t");
out.print(tpm.getNodeValue((TreeModel) model, ((TreeModel) model).getNode(i)));
}
}
out.println();
} else {
if (DEBUG) {
System.out.println(((TreeModel) model).getNode(i) + " has no parent.");
}
}
}
}
}
out.close();
fileOut.close();
} catch (IOException ioe) {
System.err.println("Unable to write file: " + ioe.getMessage());
return false;
}
if (DEBUG) {
for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
System.err.println(likelihood.getId() + ": " + likelihood.getLogLikelihood());
}
}
return true;
}
use of dr.inference.operators.AdaptableMCMCOperator in project beast-mcmc by beast-dev.
the class MCMCMC method swapChainTemperatures.
private int swapChainTemperatures() {
if (DEBUG) {
System.out.print("Current scores: ");
for (int i = 0; i < chains.length; i++) {
System.out.print("\t");
if (i == coldChain) {
System.out.print("[");
}
System.out.print(chains[i].getCurrentScore());
if (i == coldChain) {
System.out.print("]");
}
}
System.out.println();
}
int newColdChain = coldChain;
int index1 = MathUtils.nextInt(chains.length);
int index2 = MathUtils.nextInt(chains.length);
while (index1 == index2) {
index2 = MathUtils.nextInt(chains.length);
}
double score1 = chains[index1].getCurrentScore();
MCMCCriterion acceptor1 = ((MCMCCriterion) chains[index1].getAcceptor());
double temperature1 = acceptor1.getTemperature();
double score2 = chains[index2].getCurrentScore();
MCMCCriterion acceptor2 = ((MCMCCriterion) chains[index2].getAcceptor());
double temperature2 = acceptor2.getTemperature();
double logRatio = ((score2 - score1) * temperature1) + ((score1 - score2) * temperature2);
boolean swap = (Math.log(MathUtils.nextDouble()) < logRatio);
if (swap) {
if (DEBUG) {
System.out.println("Swapping chain " + index1 + " and chain " + index2);
}
acceptor1.setTemperature(temperature2);
acceptor2.setTemperature(temperature1);
OperatorSchedule schedule1 = schedules[index1];
OperatorSchedule schedule2 = schedules[index2];
for (int i = 0; i < schedule1.getOperatorCount(); i++) {
MCMCOperator operator1 = schedule1.getOperator(i);
MCMCOperator operator2 = schedule2.getOperator(i);
long tmp = operator1.getAcceptCount();
operator1.setAcceptCount(operator2.getAcceptCount());
operator2.setAcceptCount(tmp);
tmp = operator1.getRejectCount();
operator1.setRejectCount(operator2.getRejectCount());
operator2.setRejectCount(tmp);
double tmp2 = operator1.getSumDeviation();
operator1.setSumDeviation(operator2.getSumDeviation());
operator2.setSumDeviation(tmp2);
if (operator1 instanceof AdaptableMCMCOperator) {
tmp2 = ((AdaptableMCMCOperator) operator1).getAdaptableParameter();
((AdaptableMCMCOperator) operator1).setAdaptableParameter(((AdaptableMCMCOperator) operator2).getAdaptableParameter());
((AdaptableMCMCOperator) operator2).setAdaptableParameter(tmp2);
}
}
if (index1 == coldChain) {
newColdChain = index2;
} else if (index2 == coldChain) {
newColdChain = index1;
}
}
return newColdChain;
}
use of dr.inference.operators.AdaptableMCMCOperator in project beast-mcmc by beast-dev.
the class CheckPointModifier method readStateFromFile.
protected long readStateFromFile(File file, MarkovChain markovChain, double[] lnL) {
OperatorSchedule operatorSchedule = markovChain.getSchedule();
long state = -1;
this.traitModels = new ArrayList<TreeParameterModel>();
try {
FileReader fileIn = new FileReader(file);
BufferedReader in = new BufferedReader(fileIn);
int[] rngState = null;
String line = in.readLine();
String[] fields = line.split("\t");
if (fields[0].equals("rng")) {
// if there is a random number generator state present then load it...
try {
rngState = new int[fields.length - 1];
for (int i = 0; i < rngState.length; i++) {
rngState[i] = Integer.parseInt(fields[i + 1]);
}
} catch (NumberFormatException nfe) {
throw new RuntimeException("Unable to read state number from state file");
}
line = in.readLine();
fields = line.split("\t");
}
try {
if (!fields[0].equals("state")) {
throw new RuntimeException("Unable to read state number from state file");
}
state = Long.parseLong(fields[1]);
} catch (NumberFormatException nfe) {
throw new RuntimeException("Unable to read state number from state file");
}
line = in.readLine();
fields = line.split("\t");
try {
if (!fields[0].equals("lnL")) {
throw new RuntimeException("Unable to read lnL from state file");
}
if (lnL != null) {
lnL[0] = Double.parseDouble(fields[1]);
}
} catch (NumberFormatException nfe) {
throw new RuntimeException("Unable to read lnL from state file");
}
line = in.readLine();
// System.out.println(line);
fields = line.split("\t");
// Tree nodes have numbers as parameter ids
for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
// numbers should be positive but can include zero
if (isTreeNode(parameter.getId()) && isTreeNode(fields[1]) || parameter.getId().equals(fields[1])) {
int dimension = Integer.parseInt(fields[2]);
if (dimension != parameter.getDimension() && !fields[1].equals("branchRates.categories")) {
System.err.println("Unable to match state parameter dimension: " + dimension + ", expecting " + parameter.getDimension() + " for parameter: " + parameter.getParameterName());
System.err.print("Read from file: ");
for (int i = 0; i < fields.length; i++) {
System.err.print(fields[i] + "\t");
}
System.err.println();
}
if (fields[1].equals("branchRates.categories.rootNodeNumber")) {
// System.out.println("eek");
double value = Double.parseDouble(fields[3]);
parameter.setParameterValue(0, value);
if (DEBUG) {
System.out.println("restoring " + fields[1] + " with value " + value);
}
} else {
if (DEBUG) {
System.out.print("restoring " + fields[1] + " with values ");
}
if (fields[1].equals("branchRates.categories")) {
for (int dim = 0; dim < (fields.length - 3); dim++) {
// System.out.println("dim " + dim);
parameter.setParameterValue(dim, Double.parseDouble(fields[dim + 3]));
if (DEBUG) {
System.out.print(Double.parseDouble(fields[dim + 3]) + " ");
}
}
} else {
for (int dim = 0; dim < parameter.getDimension(); dim++) {
parameter.setParameterValue(dim, Double.parseDouble(fields[dim + 3]));
if (DEBUG) {
System.out.print(Double.parseDouble(fields[dim + 3]) + " ");
}
}
}
if (DEBUG) {
System.out.println();
}
}
line = in.readLine();
// System.out.println(line);
fields = line.split("\t");
} else {
if (DEBUG) {
System.out.println(" unable to match " + parameter.getId() + " with " + fields[1] + " (moving on ...)");
}
// there will be more parameters in the connected set than there are lines in the checkpoint file
// TODO keep track of these parameters and print a list of those parameters to screen
}
}
// TODO remove else-clause (and boolean) after multiple rounds of testing
if (IN_MEMORY) {
// first read in all the operator lines from the checkpoint file
// store them in a HashMap (or other structure that allows easy look-up)
HashMap<String, String[]> operatorMap = new HashMap<String, String[]>();
while (fields[0].equals("operator")) {
operatorMap.put(fields[1], fields);
line = in.readLine();
fields = line.split("\t");
}
// then iterate over the operator schedule and look into the HashMap for the information
for (int i = 0; i < operatorSchedule.getOperatorCount(); i++) {
MCMCOperator operator = operatorSchedule.getOperator(i);
String[] lookup = operatorMap.get(operator.getOperatorName());
if (lookup == null) {
// could be additional operator so not necessarily a problem that warrants an exception
if (DEBUG) {
System.out.println("No information found in checkpoint file for operator " + operator.getOperatorName());
}
} else {
// entry was found in stored information
if (DEBUG) {
System.out.println("restoring operator " + operator.getOperatorName() + " with settings from " + lookup[1]);
}
if (lookup.length < 4) {
throw new RuntimeException("Operator missing values: " + lookup[1] + ", length=" + lookup.length);
}
operator.setAcceptCount(Integer.parseInt(lookup[2]));
operator.setRejectCount(Integer.parseInt(lookup[3]));
if (operator instanceof AdaptableMCMCOperator) {
if (lookup.length != 6) {
throw new RuntimeException("Coercable operator missing parameter: " + lookup[1]);
}
((AdaptableMCMCOperator) operator).setAdaptableParameter(Double.parseDouble(lookup[4]));
}
}
// don't forget to remove the entry from the hash map
operatorMap.remove(operator.getOperatorName());
}
if (DEBUG) {
System.out.println("Number of entries left in stored operator map: " + operatorMap.size());
}
} else {
// No changes needed for loading in operators
for (int i = 0; i < operatorSchedule.getOperatorCount(); i++) {
MCMCOperator operator = operatorSchedule.getOperator(i);
if (DEBUG) {
System.out.println("restoring operator " + operator.getOperatorName() + " with settings from " + fields[1]);
}
if (!fields[1].equals(operator.getOperatorName())) {
throw new RuntimeException("Unable to match operator: " + fields[1] + " vs. " + operator.getOperatorName());
}
if (fields.length < 4) {
throw new RuntimeException("Operator missing values: " + fields[1] + ", length=" + fields.length);
}
operator.setAcceptCount(Integer.parseInt(fields[2]));
operator.setRejectCount(Integer.parseInt(fields[3]));
if (operator instanceof AdaptableMCMCOperator) {
if (fields.length != 6) {
throw new RuntimeException("Coercable operator missing parameter: " + fields[1]);
}
((AdaptableMCMCOperator) operator).setAdaptableParameter(Double.parseDouble(fields[4]));
}
line = in.readLine();
fields = line.split("\t");
}
}
// load the tree models last as we get the node heights from the tree (not the parameters which
// which may not be associated with the right node
Set<String> expectedTreeModelNames = new HashSet<String>();
for (Model model : Model.CONNECTED_MODEL_SET) {
if (model instanceof TreeModel) {
expectedTreeModelNames.add(model.getModelName());
}
if (model instanceof TreeParameterModel) {
this.traitModels.add((TreeParameterModel) model);
}
if (model instanceof BranchRates) {
this.rateModel = (BranchRates) model;
}
// e.g. MixtureModelBranchRates uses an array of ParametricDistributionModel
if (model instanceof DiscretizedBranchRates) {
parDistMod = ((DiscretizedBranchRates) model).getParametricDistributionModel();
}
}
while (fields[0].equals("tree")) {
for (Model model : Model.CONNECTED_MODEL_SET) {
if (model instanceof TreeModel && fields[1].equals(model.getModelName())) {
// AR: Can we not just add them to a Flexible tree and then make a new TreeModel
// taking that in the constructor?
// internally, we have a tree with all the taxa
// externally, i.e. in the checkpoint file, we have a tree representation comprising
// a subset of the full taxa set
// write method that adjusts the internal representation, i.e. the one in the connected
// set, according to the checkpoint file and a distance-based approach to position
// the additional taxa
// first read in all the data from the checkpoint file
line = in.readLine();
line = in.readLine();
fields = line.split("\t");
// read number of nodes
int nodeCount = Integer.parseInt(fields[0]);
double[] nodeHeights = new double[nodeCount];
String[] taxaNames = new String[(nodeCount + 1) / 2];
for (int i = 0; i < nodeCount; i++) {
line = in.readLine();
fields = line.split("\t");
nodeHeights[i] = Double.parseDouble(fields[1]);
if (i < taxaNames.length) {
taxaNames[i] = fields[2];
}
}
// on to reading edge information
line = in.readLine();
line = in.readLine();
line = in.readLine();
fields = line.split("\t");
int edgeCount = Integer.parseInt(fields[0]);
// create data matrix of doubles to store information from list of TreeParameterModels
double[][] traitValues = new double[traitModels.size()][edgeCount];
// create array to store whether a node is left or right child of its parent
// can be important for certain tree transition kernels
int[] childOrder = new int[edgeCount];
for (int i = 0; i < childOrder.length; i++) {
childOrder[i] = -1;
}
int[] parents = new int[edgeCount];
for (int i = 0; i < edgeCount; i++) {
parents[i] = -1;
}
for (int i = 0; i < edgeCount; i++) {
line = in.readLine();
if (line != null) {
fields = line.split("\t");
parents[Integer.parseInt(fields[0])] = Integer.parseInt(fields[1]);
childOrder[i] = Integer.parseInt(fields[2]);
for (int j = 0; j < traitModels.size(); j++) {
traitValues[j][i] = Double.parseDouble(fields[3 + j]);
}
}
}
// perform magic with the acquired information
// CheckPointTreeModifier modifyTree = new CheckPointTreeModifier((TreeModel) model);
this.modifyTree = new CheckPointTreeModifier((TreeModel) model, parDistMod);
// this.modifyTree = new CheckPointTreeModifier((TreeModel) model);
modifyTree.adoptTreeStructure(parents, nodeHeights, childOrder, taxaNames);
if (traitModels.size() > 0) {
modifyTree.adoptTraitData(parents, this.traitModels, traitValues);
}
// adopt the loaded tree structure; this does not yet copy the traits on the branches
// ((TreeModel) model).beginTreeEdit();
// ((TreeModel) model).adoptTreeStructure(parents, nodeHeights, childOrder);
// ((TreeModel) model).endTreeEdit();
expectedTreeModelNames.remove(model.getModelName());
}
}
line = in.readLine();
if (line != null) {
fields = line.split("\t");
}
}
if (expectedTreeModelNames.size() > 0) {
StringBuilder sb = new StringBuilder();
for (String notFoundName : expectedTreeModelNames) {
sb.append("Expecting, but unable to match state parameter:" + notFoundName + "\n");
}
throw new RuntimeException(sb.toString());
}
in.close();
fileIn.close();
} catch (IOException ioe) {
throw new RuntimeException("Unable to read file: " + ioe.getMessage());
}
return state;
}
Aggregations