Search in sources :

Example 1 with AdaptableMCMCOperator

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;
}
Also used : TreeModel(dr.evomodel.tree.TreeModel) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) OperatorSchedule(dr.inference.operators.OperatorSchedule) TreeParameterModel(dr.evomodel.tree.TreeParameterModel) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator) TreeParameterModel(dr.evomodel.tree.TreeParameterModel) Model(dr.inference.model.Model) TreeModel(dr.evomodel.tree.TreeModel) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) Parameter(dr.inference.model.Parameter) MCMCOperator(dr.inference.operators.MCMCOperator) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator)

Example 2 with AdaptableMCMCOperator

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;
}
Also used : OperatorSchedule(dr.inference.operators.OperatorSchedule) Likelihood(dr.inference.model.Likelihood) TreeParameterModel(dr.evomodel.tree.TreeParameterModel) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator) TreeModel(dr.evomodel.tree.TreeModel) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) NodeRef(dr.evolution.tree.NodeRef) TreeParameterModel(dr.evomodel.tree.TreeParameterModel) Model(dr.inference.model.Model) TreeModel(dr.evomodel.tree.TreeModel) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) Parameter(dr.inference.model.Parameter) MCMCOperator(dr.inference.operators.MCMCOperator) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator)

Example 3 with AdaptableMCMCOperator

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;
}
Also used : AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator) OperatorSchedule(dr.inference.operators.OperatorSchedule) MCMCCriterion(dr.inference.mcmc.MCMCCriterion) MCMCOperator(dr.inference.operators.MCMCOperator) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator)

Example 4 with AdaptableMCMCOperator

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;
}
Also used : HashMap(java.util.HashMap) TreeModel(dr.evomodel.tree.TreeModel) DiscretizedBranchRates(dr.evomodel.branchratemodel.DiscretizedBranchRates) FileReader(java.io.FileReader) HashSet(java.util.HashSet) OperatorSchedule(dr.inference.operators.OperatorSchedule) TreeParameterModel(dr.evomodel.tree.TreeParameterModel) IOException(java.io.IOException) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator) BufferedReader(java.io.BufferedReader) TreeParameterModel(dr.evomodel.tree.TreeParameterModel) Model(dr.inference.model.Model) TreeModel(dr.evomodel.tree.TreeModel) ParametricDistributionModel(dr.inference.distribution.ParametricDistributionModel) Parameter(dr.inference.model.Parameter) BranchRates(dr.evolution.tree.BranchRates) DiscretizedBranchRates(dr.evomodel.branchratemodel.DiscretizedBranchRates) MCMCOperator(dr.inference.operators.MCMCOperator) AdaptableMCMCOperator(dr.inference.operators.AdaptableMCMCOperator)

Aggregations

AdaptableMCMCOperator (dr.inference.operators.AdaptableMCMCOperator)4 MCMCOperator (dr.inference.operators.MCMCOperator)4 OperatorSchedule (dr.inference.operators.OperatorSchedule)4 TreeModel (dr.evomodel.tree.TreeModel)3 TreeParameterModel (dr.evomodel.tree.TreeParameterModel)3 Model (dr.inference.model.Model)3 Parameter (dr.inference.model.Parameter)3 DefaultTreeModel (dr.evomodel.tree.DefaultTreeModel)2 BranchRates (dr.evolution.tree.BranchRates)1 NodeRef (dr.evolution.tree.NodeRef)1 DiscretizedBranchRates (dr.evomodel.branchratemodel.DiscretizedBranchRates)1 ParametricDistributionModel (dr.inference.distribution.ParametricDistributionModel)1 MCMCCriterion (dr.inference.mcmc.MCMCCriterion)1 Likelihood (dr.inference.model.Likelihood)1 BufferedReader (java.io.BufferedReader)1 FileReader (java.io.FileReader)1 IOException (java.io.IOException)1 HashMap (java.util.HashMap)1 HashSet (java.util.HashSet)1