public class RunConsense {
	
	private TextOutputStream stream;
	private Consensus consensus;
	
//	private String criterion;
	private String consensusType;
	private InformationCriterion criterion;
	private int numModels;
	private Model[] model;
	private int[] order;
	private double confidenceInterval;
	private Vector<Model> confidenceModels;

	private double[] w;
	private double[] cumw;
	private boolean[] isInInterval;

	public static es.uvigo.darwin.jmodeltest.threads.SwingWorker workerConsense;

	// constructor
	public RunConsense(InformationCriterion tcriterion, String tconsensusType,
			double minterval) {

		criterion = tcriterion;
		consensusType = tconsensusType;
		numModels = tcriterion.getNumModels();
		model = ModelTest.getCandidateModels();
		order = new int[numModels];
		confidenceInterval = minterval;
		confidenceModels = new Vector<Model>();

		w = new double[numModels];
		cumw = new double[numModels];
		isInInterval = new boolean[numModels];

		stream = ModelTest.getMainConsole();

		/**
		 * This action listener, called by the "Start" button, effectively forks
		 * the thread that does the work.
		 * 
		 * /* Invoking start() on the SwingWorker causes a new Thread to be
		 * creaconsense that will call construct(), and then finished(). Note
		 * that finished() is called even if the worker is interrupconsense
		 * because we catch the InterrupconsenseException in doConsense().
		 */
		if (ModelTest.buildGUI)
			System.out.println("\nComputing model averaged phylogeny ...");

		buildConfidenceInterval();
		consensus = doConsense();
		printConsensus();

		if (ModelTest.buildGUI) {
			XManager.getInstance().getPane().setCaretPosition(
					XManager.getInstance().getPane().getDocument()
					.getLength());
			System.out.println("OK");
		}

	}

	/**
	 * This method represents the application code that we'd like to run on a
	 * separate thread.
	 */
	// Object doConsense()
	private Consensus doConsense() {

		List<WeightedTree> treeList = new ArrayList<WeightedTree>();

		for (Model m : confidenceModels) {

			try {
				// parse tree
				String tree = m.getTreeString();
				StringReader sr = new StringReader(tree);
				Tree t = new ReadTree(new PushbackReader(sr));

				double weight;
				// set criterion
				switch (criterion.getType()) {
				case InformationCriterion.IC_AIC:
					weight = m.getAICw();
					break;
				case InformationCriterion.IC_AICc:
					weight = m.getAICcw();
					break;
				case InformationCriterion.IC_BIC:
					weight = m.getBICw();
					break;
				case InformationCriterion.IC_DT:
					weight = m.getDTw();
					break;
				default:
					weight = 0.0d;
				}

				treeList.add(new WeightedTree(t, weight));

			} catch (TreeParseException e1) {
				// TODO Auto-generated catch block
				e1.printStackTrace();
			}
		}

		double consensusThreshold;
		if (consensusType.equals("strict")) {
			consensusThreshold = 0.99999999d;
		} else {
			consensusThreshold = 0.5d;
		}
		;
		Consensus consensus = new Consensus(treeList, consensusThreshold,
				Consensus.BRANCH_LENGTHS_MEDIAN);

		return consensus;
	} // doConsense

	/**************
	 * buildConfidenceInterval ************************
	 * 
	 * Builds the confidence interval of selected models and their cumulative
	 * weight
	 * 
	 * The model that just passed the confidence will be or not in the interval
	 * by chance (see below)
	 ****************************************************************/

	private void buildConfidenceInterval() {
		int i;
		Model tmodel;
		tmodel = model[0];

		order = criterion.order;

		// set alias
		for (i = 0; i < numModels; i++) {
			tmodel = model[order[i]];

			w[i] = criterion.getWeight(tmodel);
			cumw[i] = criterion.getCumWeight(tmodel);
			
		}

		// construct the confidence interval for models
		if (confidenceInterval == 1.0) {
			for (i = 0; i < numModels; i++) {
				tmodel = model[order[i]];
				isInInterval[i] = true;
				confidenceModels.add(tmodel);
			}
		} else {
			for (i = 0; i < numModels; i++) {
				tmodel = model[order[i]];

				// System.out.print("name=" + tmodel.name + " w=" + w[i] +
				// " cumw=" + cumw[i]);

				if (cumw[i] <= confidenceInterval) {
					isInInterval[i] = true;
					confidenceModels.add(tmodel);
				} else
					break;
			}

			// lets decide whether the model that just passed the confidence
			// interval should be included (suggested by John Huelsenbeck)
			double probOut = (tmodel.getCumAICw() - confidenceInterval)
					/ tmodel.getAICw();
			double probIn = 1.0 - probOut;
			Random generator = new Random();
			double randomNumber = generator.nextDouble();
			if (randomNumber <= probIn) {
				isInInterval[i] = true;
				confidenceModels.add(tmodel);
			} else
				isInInterval[i] = false;
		}
	}

	public double getConfidenceInterval() {
		return confidenceInterval;
	}
	
	public String getConsensusType() {
		return consensusType;
	}
	
	public List<Model> getConfidenceModels() {
		return confidenceModels;
	}
	
	public Tree getConsensus() {
		return consensus.getConsensusTree();
	}
	
	private void printConsensus() {
		
		double consensusThreshold = consensusType.equals("50% majority rule")?0.5:1.0;
		
		// print results for best AIC model
		stream.println(" ");stream.println(" ");stream.println(" ");
		stream.println("---------------------------------------------------------------");
		stream.println("*                                                             *");
		stream.println("*                    MODEL AVERAGED PHYLOGENY                 *");
		stream.println("*                                                             *");
		stream.println("---------------------------------------------------------------");

		if (criterion.getType() == InformationCriterion.IC_DT)
			Utilities
					.printRed("\nWarning: The DT weights used for this model averaged phylogeny are very gross"
							+ " and should be used with caution. See the program documentation.\n");
        stream.println(" ");		
        stream.println("Selection criterion: . . . . " + criterion);
        stream.print("Confidence interval: . . . . "); 
        stream.printf("%4.2f\n", confidenceInterval);
        stream.println("Consensus type:. . . . . . . " + consensusType);
        stream.println(" ");

     // print confidence set
		stream.println(" ");
		stream.print("Using " + confidenceModels.size() + " models in the ");
		stream.printf("%4.2f ", confidenceInterval);
		stream.print("confidence interval = ");
		for (Model m : confidenceModels) {
			stream.print(m.getName() + " ");
		}
		stream.println(" ");
		
        Set<FixedBitSet> keySet = consensus.getCladeSupport().keySet();
        List<FixedBitSet> splitsInConsensus = new ArrayList<FixedBitSet>();
        List<FixedBitSet> splitsOutFromConsensus = new ArrayList<FixedBitSet>();

        for (FixedBitSet fbs : keySet) {
            if (fbs.cardinality() > 1) {
                double psupport = (1.0 * consensus.getCladeSupport().get(fbs)) / 1.0;
                if (psupport < consensusThreshold) {
                    splitsOutFromConsensus.add(fbs);
                } else {
                    splitsInConsensus.add(fbs);
                }
            }
        }

        stream.println(" ");
        stream.println("Species in order:");
        stream.println(" ");

        int numTaxa = consensus.getIdGroup().getIdCount();
        for (int i = 0; i < numTaxa; i++) {
            Identifier id = consensus.getIdGroup().getIdentifier(i);
            stream.println("    " + (i + 1) + ". " + id.getName());
        }
 
        stream.println(" ");
        stream.println("Bipartitions included in the consensus tree");
        stream.println(" ");

        stream.println(consensus.getSetsIncluded());
        stream.println(" ");
        
        Tree consensusTree = consensus.getConsensusTree();
        StringWriter sw = new StringWriter();
        PrintWriter pw = new PrintWriter(sw);
        TreeUtilities.printASCII(consensusTree, pw);
        pw.flush();
        stream.println(sw);
        stream.println(" ");
        String newickTree = TreeUtilities.toNewick(consensusTree, true, true, true);
        stream.println(newickTree);
        stream.println(" ");
        stream.println("Note: this tree is unrooted. Branch lengths are the expected number of "
        		+ "substitutions per site. Labels next to parentheses represent phylogenetic "
        		+ "uncertainty due to model selection (see documentation)");
	}

} // class RunConsense public abstract class RunPhyml extends Observable implements Observer {

	// Set of variables for tuning the guided search algorithm
	private static final boolean filterFrequencies = true;
	private static final boolean filterRateMatrix = true;
	private static final boolean filterRateVariation = true;

	protected ApplicationOptions options;
	protected Model[] models;
	protected Model gtrModel = null;

	public static final String[] COMPATIBLE_VERSIONS = {
		"20130103", "20131022", 
		"20141009", "20141029", 
		"20150501", "20151222"};

	public static String PHYML_VERSION = "3.0";

	public static String PHYML_TREE_SUFFIX = "_phyml_tree_";
	public static String PHYML_STATS_SUFFIX = "_phyml_stats_";

	public static File phymlBinary;
	public static String phymlBinaryStr;
	private static String CURRENT_DIRECTORY = ModelTestConfiguration.PATH;
	private static boolean PHYML_GLOBAL = ModelTestConfiguration.isGlobalPhymlBinary();
	public static String PHYML_PATH = CURRENT_DIRECTORY + "exe/phyml/";
	private static boolean compatiblePhyml = false;
	
	protected Observer progress;

	public static boolean isCompatible() {
		return compatiblePhyml;
	}
	
	public static boolean checkBinary() {
		boolean canExecute = false;
    if (!ModelTestConfiguration.isGlobalPhymlBinary()) {
      if (!RunPhyml.phymlBinary.exists()) {
        if (ModelTest.MPJ_ME == 0)
          Utilities.printRed("ERROR: PhyML binary cannot be found: "
              + RunPhyml.phymlBinary.getAbsolutePath() + "\n");
      } else if (!RunPhyml.phymlBinary.canExecute()) {
        if (ModelTest.MPJ_ME == 0)
          Utilities.printRed(
              "ERROR: PhyML binary exists, but it cannot be executed: "
                  + RunPhyml.phymlBinary.getAbsolutePath() + "\n");
      } else if (!RunPhyml.isCompatible()) {
        if (ModelTest.MPJ_ME == 0) {
          Utilities.printRed(
              "WARNING: PhyML binary is not in the list of compatibility: \n");
          Utilities.printRed(RunPhyml.phymlBinary.getAbsolutePath() + " v"
              + RunPhyml.PHYML_VERSION + "\n");
          Utilities.printRed("Compatible versions: ");
          for (int i = 0; i < RunPhyml.COMPATIBLE_VERSIONS.length; i++)
            Utilities.printBlue(RunPhyml.COMPATIBLE_VERSIONS[i] + " ");
          Utilities.printBlue("\n");
          Utilities.printRed(
              "jModelTest will try to continue execution anyway, but it might fail.\n");
        }
        canExecute = true;
      } else {
        if (ModelTest.MPJ_ME == 0)
          Utilities.printBlue(
              "PhyML binary: " + RunPhyml.phymlBinary.getAbsolutePath() + " v"
                  + RunPhyml.PHYML_VERSION + "\n");
        canExecute = true;
      }
    }
		return canExecute;
	}
	private static boolean checkPhymlCompatibility(String binary) {
		boolean binaryFound = false;
		if (Utilities.findCurrentOS() != Utilities.OS_LINUX)
			return false;
		// get process and execute command line
		String cmd[] = {binary, "--version"};
		Runtime rt = Runtime.getRuntime();
		Process proc;
		try {
			proc = rt.exec(cmd, null, null);
			InputStream stdin = proc.getInputStream();
			InputStreamReader isr = new InputStreamReader(stdin);
			BufferedReader br = new BufferedReader(isr);

			String line = null;
			while ( (line = br.readLine()) != null) {
				if (line.toLowerCase().contains("phyml version")) {
					String[] linesplit = line.trim().replace(".", "").split(" ");
					String version = linesplit[linesplit.length - 1];
					PHYML_VERSION = version;
					if (Arrays.asList(COMPATIBLE_VERSIONS).contains(version)) {
						compatiblePhyml = true;
					} else {
						compatiblePhyml = false;
					}
					binaryFound = true;
				}
			}
			br.close();
			isr.close();
			stdin.close();
			proc.destroy();
		} catch (IOException e) {
			return false;
		}
		return binaryFound;
	}
	
	public static boolean isPhymlGlobal() {
		return PHYML_GLOBAL;
	}
	
	static {
		if (PHYML_GLOBAL) {
			PHYML_PATH = "";
			phymlBinaryStr = "phyml";
		} else {
			/* check the local paths */
			String path = ModelTestConfiguration.getExeDir();
			if (!path.startsWith(File.separator)) {
				PHYML_PATH = CURRENT_DIRECTORY + File.separator + path;
			} else {
				PHYML_PATH = path;
			}
			if (!PHYML_PATH.endsWith(File.separator)) {
				PHYML_PATH += File.separator;
			}
			phymlBinary = new File(PHYML_PATH + "phyml");
			if (phymlBinary.exists() && phymlBinary.canExecute()) {
				phymlBinaryStr = phymlBinary.getAbsolutePath();
				compatiblePhyml = true;
			} else {
				phymlBinaryStr = PHYML_PATH + Utilities.getBinaryVersion();
			}
			/* Check if binary exists and is compatible */
			phymlBinary = new File(phymlBinaryStr);
			if (!checkPhymlCompatibility(phymlBinaryStr)) {
				/* Check for system wide PhyML */
				if (!(phymlBinary.exists() && phymlBinary.canExecute()))
				{
					if(checkPhymlCompatibility("/usr/bin/phyml")) {
						PHYML_PATH = "/usr/bin";
						phymlBinary = new File(PHYML_PATH + "/phyml");
						phymlBinaryStr = phymlBinary.getAbsolutePath();						
					}
				}
			}
		}
	}
	public RunPhyml(Observer progress, ApplicationOptions options,
			Model[] models) {
		if (models != null)
			this.models = new Model[models.length];
		else
			this.models = new Model[0];
		for (int i = 0; i < this.models.length; i++)
			this.models[i] = models[i];
		this.options = options;
		this.progress = progress;
		this.addObserver(progress);
		Arrays.sort(this.models, new ModelComparator());
	}

	public void execute() {
		// remove stuff from exe directories before starting
		deleteFiles();
		printSettings(ModelTest.getMainConsole());

		// locate GTR model
		String searchFor;
		int gtrParams;
		if (options.doI && options.doG) {
			searchFor="GTR+I+G";
			gtrParams = 10;
		} else if(options.doI) {
			searchFor="GTR+I";
			gtrParams = 9;
		} else if(options.doG) {
			searchFor="GTR+G";
			gtrParams = 9;
		} else {
			searchFor="GTR";
			gtrParams = 8;
		}
		for (int i = (models.length - 1); i >= 0; i--) {
			if (models[i].getName().startsWith(searchFor)) {
				gtrModel = models[i];
				break;
			}
		}
		if (gtrModel == null) {
			gtrModel = new Model(0, searchFor, "012345", gtrParams, false, false, false, true, options.doI, options.doG, 2, 4);
		}
		// estimate a NJ-JC tree if needed
		if (options.fixedTopology) {
			Model jcModel = null;
			for (Model model : models) {
				if (model.getName().equals("JC")) {
					jcModel = model;
					break;
				}
			}

			if (jcModel != null) {
				notifyObservers(ProgressInfo.BASE_TREE_INIT, 0, jcModel, null);

				PhymlSingleModel jcModelPhyml = new PhymlSingleModel(jcModel,
						0, true, false, options);
				jcModelPhyml.addObserver(this);
				jcModelPhyml.run();

				// create JCtree file
				TextOutputStream JCtreeFile = new TextOutputStream(options
						.getTreeFile().getAbsolutePath(), false);
				JCtreeFile.print(jcModel.getTreeString() + "\n");
				JCtreeFile.close();

				notifyObservers(ProgressInfo.BASE_TREE_COMPUTED, 0, jcModel,
						null);
			}

		}

		if (options.isGuidedSearch()) {

			if (gtrModel != null) {
				// compute GTR model
				notifyObservers(ProgressInfo.GTR_OPTIMIZATION_INIT, models.length, gtrModel, null);
				PhymlSingleModel gtrPhymlModel = new PhymlSingleModel(
						gtrModel, 0, false, false, options);
				gtrPhymlModel.run();
				notifyObservers(ProgressInfo.GTR_OPTIMIZATION_COMPLETED, models.length, gtrModel, null);
				
				GuidedSearchManager gsm = new GuidedSearchManager(
						options.getGuidedSearchThreshold(), gtrModel,
						filterFrequencies, filterRateMatrix,
						filterRateVariation);
	
				models = gsm.filterModels(models);
				ModelTest.setCandidateModels(models);
			} else {
				notifyObservers(ProgressInfo.GTR_NOT_FOUND, models.length, models[0], null);
			}
		}

		// compute likelihood scores for all models

		notifyObservers(ProgressInfo.OPTIMIZATION_INIT, 0, models[0], null);

		doPhyml();

	}

	public void executeIgnoreGaps(Model[] models) {
		notifyObservers(ProgressInfo.REOPTIMIZATION_INIT, models.length, models[0], null);
		parallelExecute(models, true);
		notifyObservers(ProgressInfo.REOPTIMIZATION_COMPLETED, models.length, null, null);
	}
	
	protected boolean parallelExecute(Model models[], boolean ignoreGaps) {
		
		ExecutorService threadPool = Executors.newFixedThreadPool(options
				.getNumberOfThreads());
		Collection<Callable<Object>> c = new ArrayList<Callable<Object>>();
		int current = 0;
		for (Model model : models) {
			if (model != null) {
				PhymlSingleModel psm = new PhymlSingleModel(model, current, false,
						ignoreGaps, options);
				psm.addObserver(this);
				c.add(Executors.callable(psm));
	
				current++;
			}
		}

		Collection<Future<Object>> futures = null;
		try {
			futures = threadPool.invokeAll(c);
		} catch (InterruptedException e) {
			notifyObservers(ProgressInfo.INTERRUPTED, 0, null, null);
		}

		if (futures != null) {
			for (Future<Object> f : futures) {
				try {
					f.get();
				} catch (InterruptedException ex) {
					notifyObservers(ProgressInfo.INTERRUPTED, 0, null, null);
					ex.printStackTrace();
					return false;
				} catch (ExecutionException ex) {
					// Internal exception while computing model.
					// Let's continue with errors
					ex.printStackTrace();
					return false;
				}
			}
		}
		return true;
	}
	
	/***************************
	 * printSettings ***************************** * Prints the settings for the
	 * likelihood calculation * * *
	 ***********************************************************************/

	protected void printSettings(TextOutputStream stream) {

		stream.println(" ");
		stream.println(" ");
		stream.println("---------------------------------------------------------------");
		stream.println("*                                                             *");
		stream.println("*        COMPUTATION OF LIKELIHOOD SCORES WITH PHYML          *");
		stream.println("*                                                             *");
		stream.println("---------------------------------------------------------------");
		stream.println(" ");
		stream.println("::Settings::");
		stream.println(" ");
		stream.println(" Phyml version = " + PHYML_VERSION);
		stream.println(" Phyml binary = " + phymlBinary.getName());
		stream.println(" Phyml path = " + phymlBinary.getAbsolutePath()
				.substring(0, phymlBinary.getAbsolutePath().lastIndexOf(File.separator)) + File.separator);
		stream.println(" Candidate models = " + models.length);
		stream.print("   number of substitution schemes = ");

		if (options.getSubstTypeCode() == 0)
			stream.println("3");
		else if (options.getSubstTypeCode() == 1)
			stream.println("5");
		else if (options.getSubstTypeCode() == 2)
			stream.println("7");
		else
			stream.println("11");

		if (options.doF)
			stream.println("   including models with equal/unequal base frequencies (+F)");
		else
			stream.println("   including only models with equal base frequencies");

		if (options.doI)
			stream.println("   including models with/without a proportion of invariable sites (+I)");
		else
			stream.println("   including only models without a proportion of invariable sites");

		if (options.doG)
			stream.println("   including models with/without rate variation among sites (+G)"
					+ " (nCat = " + options.numGammaCat + ")");
		else
			stream.println("   including only models without rate variation among sites");

		stream.print(" Optimized free parameters (K) =");
		stream.print(" substitution parameters");
		if (options.countBLasParameters)
			stream.print(" + " + options.getNumBranches() + " branch lengths");
		if (options.optimizeMLTopology)
			stream.print(" + topology");
		stream.println(" ");

		stream.print(" Base tree for likelihood calculations = ");
		if (options.userTopologyExists) {
			stream.println("fixed user tree topology.");
			stream.println(" ");
			stream.print("User tree " + "("
					+ options.getInputTreeFile().getName() + ") = ");
			stream.println(options.getUserTree());
			stream.println(" ");
		} else if (options.fixedTopology) {
			stream.println("fixed BIONJ-JC tree topology");
		} else if (options.optimizeMLTopology) {
			stream.println("ML tree");
		} else {
			stream.println("BIONJ tree");
		}

		if (options.optimizeMLTopology) {
			stream.print(" Tree topology search operation = ");
			switch (options.treeSearchOperations) {
			case NNI:
				stream.println("NNI");
				break;
			case SPR:
				stream.println("SPR");
				break;
			case BEST:
				stream.println("BEST");
				break;
			}
		}
		if (options.isClusteringSearch()) {
			stream.println(" Using hill-climbing hierarchical clustering");
		}
		
		if (options.isGuidedSearch()) {
			stream.println(" Using heuristic model filtering ");
		}
		stream.println(" ");
	}

	protected abstract Object doPhyml();

	/***********************************************************************
	 * interruptThread
	 * 
	 * Interrupts the calculation of likelihood scores
	 ***********************************************************************/
	public void interruptThread() {
		notifyObservers(ProgressInfo.OPTIMIZATION_COMPLETED_INTERRUPTED, 0,
				null, null);
		// workerPhyml.interrupt();
		// IS_INTERRUPTED = true;
	}

	private void deleteFiles() {
		/* phymlFolder */
		//if (options.getLogFile() != null)
			options.getLogFile().delete();
	}

	protected void notifyObservers(int type, int value, Model model,
			String message) {
		setChanged();
		notifyObservers(new ProgressInfo(type, value, model, message));
	}

	@Override
	public void update(Observable o, Object arg) {
		setChanged();
		notifyObservers(arg);
	}

} // class RunPhyml public class RunPhymlMPJ extends RunPhyml {


	private List<Model> myModels;

	// Synchronization package variables.
	// Thread safe under current operation. Keep in mind.
	volatile Model rootModel = null;
	volatile boolean rootModelRequest = false;

	public RunPhymlMPJ(Observer progress, ApplicationOptions options,
			Model[] models) {
		super(progress, options, models);
		// this.deleteObserver(progress);
		myModels = new ArrayList<Model>();

	}

	public void distribute() {

		List<Model> modelList = Arrays.asList(models);
		Distributor distributor = new Distributor(modelList, this,
				ModelTest.MPJ_ME, ModelTest.MPJ_SIZE);
		distributor.addObserver(progress);
		Thread distributorThread = new Thread(distributor);
		distributorThread.start();
		request();

		for (Model model : ModelTest.getCandidateModels()) {
			model.update(modelList.get(modelList.indexOf(model)));
		}

		notifyObservers(ProgressInfo.OPTIMIZATION_COMPLETED_OK, models.length,
				null, null);
	}

	public void request() {

		List<PhymlSingleModel> phymlEstimatorList = new ArrayList<PhymlSingleModel>();

		Model[] lastComputedModel = new Model[1];
		while (true) {
			// send request to root
			Model[] modelToReceive = null;
			Model model = null;
			if (ModelTest.MPJ_ME > 0) {
				Request modelRequest = MPI.COMM_WORLD.Isend(lastComputedModel,
						0, 1, MPI.OBJECT, 0, Distributor.TAG_SEND_REQUEST);
				// prepare reception
				modelToReceive = new Model[1];
				// wait for request
				modelRequest.Wait();
				// receive model
				Request modelReceive = MPI.COMM_WORLD.Irecv(modelToReceive, 0,
						1, MPI.OBJECT, 0, Distributor.TAG_SEND_MODEL);
				modelReceive.Wait();
				model = modelToReceive[0];
			} else {
				// This strategy is an easy way to avoid the problem of
				// thread-safety in MPJ-Express. It works correctly, but
				// it also causes to introduce coupling between this class
				// and Distributor having to define two volatile attributes:
				// rootModelRequest and rootModel.
				rootModelRequest = true;
				while (rootModelRequest) {
					try {
						Thread.sleep(200);
					} catch (InterruptedException e) {
						throw new InternalException(
								"Thread interrupted");
					}
				}
				model = rootModel;
			}
			if (model == null)
				break;
			else {
				// compute
				myModels.add(model);
				PhymlSingleModel runenv = new PhymlSingleModel(model, 0, false,
						false, options);
				runenv.addObserver(this);

				if (!runenv.compute())
					throw new InternalException("Optimization error");

				phymlEstimatorList.add(runenv);
				lastComputedModel[0] = runenv.getModel();
			}
		}

		// endTime = System.currentTimeMillis();
	}

	protected Object doPhyml() {
		return null;
	} // doPhyml

	public void execute() {

		if (ModelTest.MPJ_ME == 0) {
			printSettings(ModelTest.getMainConsole());

			// TODO: Send topology to each processor
			// estimate a NJ-JC tree if needed
			if (options.fixedTopology) {
				notifyObservers(ProgressInfo.BASE_TREE_INIT, 0, models[0], null);

				PhymlSingleModel jcModel = new PhymlSingleModel(models[0], 0,
						true, false, options);
				jcModel.run();

				// create JCtree file
				TextOutputStream JCtreeFile = new TextOutputStream(options
						.getTreeFile().getAbsolutePath(), false);
				JCtreeFile.print(models[0].getTreeString() + "\n");
				JCtreeFile.close();

				options.setUserTree(models[0].getTreeString());

				notifyObservers(ProgressInfo.BASE_TREE_COMPUTED, 0, models[0],
						null);

			}

			// compute likelihood scores for all models
			System.out.println("computing likelihood scores for " + models.length
					+ " models with Phyml " + PHYML_VERSION);
		}

		// sincronize ApplicationOptions from root
		ApplicationOptions[] optionsBCast = new ApplicationOptions[1];
		optionsBCast[0] = options;
		MPI.COMM_WORLD.Bcast(optionsBCast, 0, 1, MPI.OBJECT, 0);
		this.options = optionsBCast[0];
		ApplicationOptions.setInstance(this.options);

		if (ModelTest.MPJ_ME == 0) {
			notifyObservers(ProgressInfo.OPTIMIZATION_INIT, 0,
					models[0], null);
			distribute();
		} else {
			try {
				this.options.buildWorkFiles();
			} catch (IOException e) {
				e.printStackTrace();
			}

			request();
		}
	}

	@Override
	public void update(Observable o, Object arg) {
		// Ignore runtime messages
		// setChanged();
		// notifyObservers(arg);
	}
} // class RunPhyml public class XManager {

	private static XManager instance;
	
	public static final Font FONT_CONSOLE, FONT_LABEL, 
		FONT_MENU, FONT_STATUS, FONT_SLIDER, FONT_LABEL_BIG,
		FONT_TABULAR;
	public static final Color LABEL_BLUE_COLOR = new Color(102,102,153);
	public static final Color LABEL_FAIL_COLOR = new Color(153, 0, 0);
	public static final Color LABEL_GREEN_COLOR = new Color(0, 153, 0);
	public static final Color PANEL_BORDER_COLOR = new Color(153, 153, 153);
	public static final Color MENU_COLOR = new Color(130,130,150);//(199, 199, 220);
	public static final Color MENU_FG_COLOR = new Color(0,0,0);
	public static final Color PANE_BACK_COLOR = Color.WHITE;
	public static final Color INNER_BORDER_COLOR = new Color(182,182,182); //182
	public static final Color OUTER_BORDER_COLOR = new Color(89,89,89); //89
	public static final Color DARK_GRAY_COLOR = new Color(89,89,89); //89
	public static final Color STATUS_BACK_COLOR = new Color(220,220,220); //220
	
	public static final Point MAIN_LOCATION = new Point(281, 80);
	
	public static SimpleAttributeSet blueText;
	public static SimpleAttributeSet redText;
	public static SimpleAttributeSet blackText;

	public FrameMain frame;
	public FrameResults resultsFrame;
	private FramePreferences preferencesFrame;
	public static boolean resultsFrameBuilt = false;
	private JTextPane PANE;

	static {

		if (Utilities.isWindows() == false) {
			FONT_CONSOLE = new Font(Font.MONOSPACED, 0, 12);
		} else {
			FONT_CONSOLE = new Font("Lucida Console", 0, 12);
		}
		FONT_LABEL = new Font("Application", 1, 10);
		FONT_LABEL_BIG = new Font("Application", 0, 12);
		FONT_SLIDER = new Font("Application", 1, 9);
		FONT_MENU = new Font("Dialog", 0, 9);
		FONT_STATUS = new Font("Dialog", 0, 9);
		FONT_TABULAR = new java.awt.Font("Verdana", 0, 12);

		blackText = new SimpleAttributeSet();
		redText = new SimpleAttributeSet();
		blueText = new SimpleAttributeSet();

		StyleConstants.setForeground(blueText, Color.blue);
		StyleConstants.setForeground(redText, Color.red);
		StyleConstants.setForeground(blackText, Color.black);
	}

	private XManager() {

		try {
			try {
				UIManager.setLookAndFeel(UIManager
						.getSystemLookAndFeelClassName());
			} catch (Exception e) {
				e.printStackTrace();
			}

			frame = new FrameMain();
			frame.initComponents();
			frame.setVisible(true);
			PANE = frame.getMainEditorPane();

			ModelTest.setMainConsole(new TextOutputStream(new PrintStream(
					new DocumentOutputStream(PANE.getDocument()))));

			ModelTest.setCurrentOutStream(ModelTest.getMainConsole());
			ModelTest.printHeader(ModelTest.getMainConsole());
			ModelTest.printNotice(ModelTest.getMainConsole());
			ModelTest.printCitation(ModelTest.getMainConsole());

			// check expiration date
			// CheckExpiration (frame);
		} catch (Exception e) {
			e.printStackTrace();
		}

	}

	/****************************
	 * CheckExpiration*************************** * Checks for expiration and
	 * quits if needed * * *
	 ***********************************************************************/
	public static void CheckExpiration(JFrame theframe) {
		java.util.Calendar now = java.util.Calendar.getInstance();

		if ((now.get(Calendar.MONTH) != Calendar.JUNE && now.get(Calendar.MONTH) != Calendar.DECEMBER)
				|| now.get(Calendar.YEAR) != 2007) {
			JOptionPane.showMessageDialog(theframe,
					"Program has expired! \n    Bye...", "jModelTest warning",
					JOptionPane.WARNING_MESSAGE);
			theframe.dispose();
			System.exit(0);
		}
	}

	public void setPane(JTextPane pane) {
		PANE = pane;
	}

	public JTextPane getPane() {
		return PANE;
	}

	/****************************
	 * buildFrameResults **************************** * Builds the frame that
	 * displays the model selection results * *
	 ************************************************************************/

	public void buildFrameResults() {
		try {
			if (resultsFrameBuilt) {
				resultsFrame.dispose();
			}
			resultsFrame = new FrameResults();
			resultsFrame.initComponents();
			resultsFrame.setVisible(false);
			resultsFrameBuilt = true;
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public static XManager getInstance() {
		if (instance == null) {
			instance = new XManager();
		}
		return instance;
	}
	
	public void loadFramePreferences() {
		if (preferencesFrame == null) {
			preferencesFrame = new FramePreferences();			
		}
		preferencesFrame.setVisible(true);
	}
	
	public void setLikeLabelText(String text) {
		frame.setLikeLabelText(text);
	}
	
	public void setDataLabelText(String text) {
		frame.setDataLabelText(text);
	}
	
	public void setLikeLabelColor(Color color) {
		frame.setLikeLabelColor(color);
	}
	
	public void enableMenuAIC(boolean enabled) {
		frame.enableMenuAIC(enabled);
	}
	
	public void enableMenuBIC(boolean enabled) {
		frame.enableMenuBIC(enabled);
	}
	
	public void enableMenuDT(boolean enabled) {
		frame.enableMenuDT(enabled);
	}
	
	public void enableMenuhLRT(boolean enabled) {
		frame.enableMenuhLRT(enabled);
	}
	
	public void enableMenuAveraging(boolean enabled) {
		frame.enableMenuAveraging(enabled);
	}
	
	public void enableMenuShowModelTable(boolean enabled) {
		frame.enableMenuShowModelTable(enabled);
	}
	
	public void enableMenuHtmlOutput(boolean enabled) {
		frame.enableMenuHtmlOutput(enabled);
	}
	
	public void selectedMenuResultsBLasParameters(boolean selected) {
		frame.selectedMenuResultsBLasParameters(selected);
	}
	
	public static ImageIcon makeIcon(String imageName, String altText) {
		//Look for the image.
	    String imgLocation = "icons/" + imageName
	        + ".gif";
	    URL imageURL = XManager.class.getResource(imgLocation);
	    ImageIcon icon = null;
	    if (imageURL != null) {
	    	icon = new ImageIcon(imageURL, altText);
	    }
	    return icon;
	}
	
	public static JButton makeIconButton(String imageName,
		      String toolTipText, String altText) {

		    ImageIcon icon = makeIcon(imageName, altText);
		    
		    //Create and initialize the button.
		    JButton button = new JButton();
		    button.setToolTipText(toolTipText);
		    button.setLayout(null);
		    button.setOpaque(false);
		    button.setSize(20,20);
		    button.setVisible(true);

		    if (icon != null) { //image found
		      button.setIcon(icon);
		    } else { //no image found
		      button.setText(altText);
		    }

		    return button;
	  }
} public abstract class AlignmentReader {
	private static ApplicationOptions options = ApplicationOptions
			.getInstance();;

	/****************************
	 * readDataFile **************************** * Reads the input file and gets
	 * the number of taxa and alignment * length * *
	 ***********************************************************************/

	static public void getHeader(String infilenameComplete) {
		// needs the complete path to the file
		try {
			TextInputStream in = new TextInputStream(infilenameComplete);
			String line = in.readLine();
			in.close();
			StringTokenizer reader = new StringTokenizer(line);
			options.setNumTaxa(Integer.parseInt(reader.nextToken()));
			options.setNumSites(Integer.parseInt(reader.nextToken()));
		} catch (FileNotFoundException e) {
			JOptionPane.showMessageDialog(new JFrame(),
					"Could not read the input alignment", "jModelTest error",
					JOptionPane.ERROR_MESSAGE);
		}

		if (options.getNumTaxa() <= 4)
			JOptionPane.showMessageDialog(new JFrame(),
					"The number of taxa does not seem to be correct: "
							+ options.getNumTaxa(), "jModelTest error",
					JOptionPane.ERROR_MESSAGE);

		if (options.getNumSites() <= 1)
			JOptionPane.showMessageDialog(new JFrame(),
					"The number of sites does not seem to be correct: "
							+ options.getNumTaxa(), "jModelTest error",
					JOptionPane.ERROR_MESSAGE);
	}

	/**
	 * Creates an alignment instance.
	 * 
	 * @param out
	 *            the out
	 * @param pr
	 *            the pushback reader which contains the alignment
	 * @param debug
	 *            the debug state
	 * 
	 * @return the alignment
	 * 
	 * @throws AlignmentParseException
	 *             the alignment parse exception
	 * @throws IOException
	 *             Signals that an I/O exception has occured.
	 */
	public static Alignment createAlignment(PrintWriter out, PushbackReader pr,
			boolean debug) throws AlignmentParseException, IOException {

		if (debug) {
			out.println("");
			out.println("**********************************************************");
			out.println("  Reading alignment...");
		}
		Alignment alignment = null;
		try {
			alignment = new ReadAlignment(pr);
		} catch (pal.alignment.AlignmentParseException e) {
			throw new AlignmentParseException(e.getMessage());
		}

		List<String> seqNames = new ArrayList<String>(
				alignment.getSequenceCount());
		for (int i = 0; i < alignment.getSequenceCount(); i++) {
			seqNames.add(alignment.getIdentifier(i).getName());
		}
		
		if (debug) {
			for (int i = 0; i < alignment.getSequenceCount(); i++) {
				out.println("    Sequence #" + (i + 1) + ": "
						+ alignment.getIdentifier(i).getName());
			}
			out.println("   Alignment contains " + alignment.getSequenceCount()
					+ " sequences of length " + alignment.getSiteCount());
			out.println("");
			out.println("**********************************************************");
			out.println("");
		}

		return alignment;
	}
} // class AlignmentReader public abstract class HtmlReporter {

	private static String[] TEMPLATE_DIRS = { "resources" };
	private static String[] TEMPLATE_FILES = {
			"resources" + File.separator + "style.css",
			"resources" + File.separator + "homeIcon.gif",
			"resources" + File.separator + "topIcon.gif",
			"resources" + File.separator + "logo0.png" };
	private static Map<String, Object> datamodel;
	private static File LOG_DIR;
	private static File IMAGES_DIR;
	private static String RESOURCES_DIR;
	
	static {
		String strLogDir = ModelTestConfiguration.getLogDir();
		if (!strLogDir.startsWith(File.separator)) {
			strLogDir = ModelTestConfiguration.PATH + strLogDir;
		}
		LOG_DIR = new File(strLogDir);
		IMAGES_DIR = new File(strLogDir + File.separator + "images");
		RESOURCES_DIR = ModelTestConfiguration.PATH + "resources" + File.separator + "template";
	}
	
	public static void buildReport(ApplicationOptions options, Model models[],
			File mOutputFile) {
		buildReport(options, models, mOutputFile, null);
	}
	
	public static void buildReport(ApplicationOptions options, Model models[],
			File mOutputFile, TreeSummary summary) {

		File outputFile;
		if (mOutputFile != null) {
			if (!(mOutputFile.getName().endsWith(".htm") || mOutputFile.getName()
					.endsWith(".html"))) {
				outputFile = new File(mOutputFile.getAbsolutePath() + ".html");
			} else {
				outputFile = mOutputFile;
			}
		} else {
			outputFile = new File(LOG_DIR.getPath() + File.separator
					+ options.getInputFile().getName()
					+ ".jmodeltest." + options.getExecutionName()
					+ ".html");
		}
		
		// Add the values in the datamodel
		datamodel = new HashMap<String, Object>();
		java.util.Date current_time = new java.util.Date();
		datamodel.put("date", current_time.toString());
		datamodel.put("system", System.getProperty("os.name") + " "
				+ System.getProperty("os.version") + ", arch: "
				+ System.getProperty("os.arch") + ", bits: "
				+ System.getProperty("sun.arch.data.model") + ", numcores: "
				+ Runtime.getRuntime().availableProcessors());
		
		fillInWithOptions(options);
		fillInWithSortedModels(models);
		datamodel.put("isTopologiesSummary", summary!=null ? new Integer(1) : new Integer(0));
		if (summary != null) {
			fillInWithTopologies(summary, options);
		}

		if (options.doAIC) {
			Collection<Map<String, String>> aicModels = new ArrayList<Map<String, String>>();
			Map<String, String> bestAicModel = new HashMap<String, String>();
			fillInWIthInformationCriterion(ModelTest.getMyAIC(), aicModels,
					bestAicModel);
			datamodel.put("aicModels", aicModels);
			datamodel.put("bestAicModel", bestAicModel);
			datamodel.put("aicConfidenceCount", ModelTest.getMyAIC()
					.getConfidenceModels().size());
			StringBuffer aicConfModels = new StringBuffer();
			for (Model model : ModelTest.getMyAIC().getConfidenceModels())
				aicConfModels.append(model.getName() + " ");
			datamodel.put("aicConfidenceList", aicConfModels.toString());
			if (options.writePAUPblock) {
				ByteArrayOutputStream baos = new ByteArrayOutputStream();
				PrintStream ps = new PrintStream(baos);
				TextOutputStream strOutput = new TextOutputStream(ps);
			    ModelTest.WritePaupBlock(strOutput, "AIC", ModelTest.getMyAIC().getMinModel());
			    try {
			    	String pblock = baos.toString("UTF8");
			    	pblock = pblock.replaceAll("\n", "<br/>");
					datamodel.put("aicPaup", pblock);
				} catch (UnsupportedEncodingException e) {
				}
				
			}
			buildChart(outputFile, ModelTest.getMyAIC());
			datamodel.put("aicEuImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_eu_AIC.png");
			datamodel.put("aicRfImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_rf_AIC.png");
		}

		if (options.doAICc) {
			Collection<Map<String, String>> aiccModels = new ArrayList<Map<String, String>>();
			Map<String, String> bestAiccModel = new HashMap<String, String>();
			fillInWIthInformationCriterion(ModelTest.getMyAICc(), aiccModels,
					bestAiccModel);
			datamodel.put("aiccModels", aiccModels);
			datamodel.put("bestAiccModel", bestAiccModel);
			datamodel.put("aiccConfidenceCount", ModelTest.getMyAICc()
					.getConfidenceModels().size());
			StringBuffer aiccConfModels = new StringBuffer();
			for (Model model : ModelTest.getMyAICc().getConfidenceModels())
				aiccConfModels.append(model.getName() + " ");
			datamodel.put("aiccConfidenceList", aiccConfModels.toString());
			if (options.writePAUPblock) {
				ByteArrayOutputStream baos = new ByteArrayOutputStream();
				PrintStream ps = new PrintStream(baos);
				TextOutputStream strOutput = new TextOutputStream(ps);
			    ModelTest.WritePaupBlock(strOutput, "AICc", ModelTest.getMyAICc().getMinModel());
			    try {
			    	String pblock = baos.toString("UTF8");
			    	pblock = pblock.replaceAll("\n", "<br/>");
					datamodel.put("aiccPaup", pblock);
				} catch (UnsupportedEncodingException e) {
				}
				
			}
			buildChart(outputFile, ModelTest.getMyAICc());
			datamodel.put("aiccEuImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_eu_AICc.png");
			datamodel.put("aiccRfImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_rf_AICc.png");
		}

		if (options.doBIC) {
			Collection<Map<String, String>> bicModels = new ArrayList<Map<String, String>>();
			Map<String, String> bestBicModel = new HashMap<String, String>();
			fillInWIthInformationCriterion(ModelTest.getMyBIC(), bicModels,
					bestBicModel);
			datamodel.put("bicModels", bicModels);
			datamodel.put("bestBicModel", bestBicModel);
			datamodel.put("bicConfidenceCount", ModelTest.getMyBIC()
					.getConfidenceModels().size());
			StringBuffer bicConfModels = new StringBuffer();
			for (Model model : ModelTest.getMyBIC().getConfidenceModels())
				bicConfModels.append(model.getName() + " ");
			datamodel.put("bicConfidenceList", bicConfModels.toString());
			if (options.writePAUPblock) {
				ByteArrayOutputStream baos = new ByteArrayOutputStream();
				PrintStream ps = new PrintStream(baos);
				TextOutputStream strOutput = new TextOutputStream(ps);
			    ModelTest.WritePaupBlock(strOutput, "BIC", ModelTest.getMyBIC().getMinModel());
			    try {
			    	String pblock = baos.toString("UTF8");
			    	pblock = pblock.replaceAll("\n", "<br/>");
					datamodel.put("bicPaup", pblock);
				} catch (UnsupportedEncodingException e) {
				}
				
			}
			buildChart(outputFile, ModelTest.getMyBIC());
			datamodel.put("bicEuImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_eu_BIC.png");
			datamodel.put("bicRfImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_rf_BIC.png");
		}

		if (options.doDT) {
			Collection<Map<String, String>> dtModels = new ArrayList<Map<String, String>>();
			Map<String, String> bestDtModel = new HashMap<String, String>();
			fillInWIthInformationCriterion(ModelTest.getMyDT(), dtModels,
					bestDtModel);
			datamodel.put("dtModels", dtModels);
			datamodel.put("bestDtModel", bestDtModel);
			datamodel.put("dtConfidenceCount", ModelTest.getMyDT()
					.getConfidenceModels().size());
			StringBuffer dtConfModels = new StringBuffer();
			for (Model model : ModelTest.getMyDT().getConfidenceModels())
				dtConfModels.append(model.getName() + " ");
			datamodel.put("dtConfidenceList", dtConfModels.toString());
			if (options.writePAUPblock) {
				ByteArrayOutputStream baos = new ByteArrayOutputStream();
				PrintStream ps = new PrintStream(baos);
				TextOutputStream strOutput = new TextOutputStream(ps);
			    ModelTest.WritePaupBlock(strOutput, "DT", ModelTest.getMyDT().getMinModel());
			    try {
			    	String pblock = baos.toString("UTF8");
			    	pblock = pblock.replaceAll("\n", "<br/>");
					datamodel.put("dtPaup", pblock);
				} catch (UnsupportedEncodingException e) {
				}
				
			}
			buildChart(outputFile, ModelTest.getMyDT());
			datamodel.put("dtEuImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_eu_DT.png");
			datamodel.put("dtRfImagePath",IMAGES_DIR.getName() + File.separator + outputFile.getName() + "_rf_DT.png");
		}

		datamodel.put("doAICAveragedPhylogeny",
				ModelTest.getConsensusAIC() != null ? new Integer(1)
						: new Integer(0));
		if (ModelTest.getConsensusAIC() != null) {
			datamodel.put("aicConsensusTree", TreeUtilities.toNewick(ModelTest
					.getConsensusAIC().getConsensus(), true, true, true));
			datamodel.put("consensusType", ModelTest.getConsensusAIC()
					.getConsensusType());
		}
		datamodel.put("doAICcAveragedPhylogeny",
				ModelTest.getConsensusAICc() != null ? new Integer(1)
						: new Integer(0));
		if (ModelTest.getConsensusAICc() != null) {
			datamodel.put("aiccConsensusTree", TreeUtilities.toNewick(ModelTest
					.getConsensusAICc().getConsensus(), true, true, true));
			datamodel.put("consensusType", ModelTest.getConsensusAICc()
					.getConsensusType());
		}
		datamodel.put("doBICAveragedPhylogeny",
				ModelTest.getConsensusBIC() != null ? new Integer(1)
						: new Integer(0));
		if (ModelTest.getConsensusBIC() != null) {
			datamodel.put("bicConsensusTree", TreeUtilities.toNewick(ModelTest
					.getConsensusBIC().getConsensus(), true, true, true));
			datamodel.put("consensusType", ModelTest.getConsensusBIC()
					.getConsensusType());
		}
		datamodel.put("doDTAveragedPhylogeny",
				ModelTest.getConsensusDT() != null ? new Integer(1)
						: new Integer(0));
		if (ModelTest.getConsensusDT() != null) {
			datamodel.put("dtConsensusTree", TreeUtilities.toNewick(ModelTest
					.getConsensusDT().getConsensus(), true, true, true));
			datamodel.put("consensusType", ModelTest.getConsensusDT()
					.getConsensusType());
		}

		// Process the template using FreeMarker
		try {
			freemarkerDo(datamodel, "index.html", outputFile);
		} catch (Exception e) {
			System.out
					.println("There was a problem building the html log files: "
							+ e.getLocalizedMessage());
		}
		
	}

	// Process a template using FreeMarker and print the results
	static void freemarkerDo(Map<String, Object> datamodel, String template,
			File outputFile) throws Exception {
		if (!LOG_DIR.exists() || !LOG_DIR.isDirectory()) {
			LOG_DIR.delete();
			LOG_DIR.mkdir();
		}
		if (!IMAGES_DIR.exists() || !IMAGES_DIR.isDirectory()) {
			IMAGES_DIR.delete();
			IMAGES_DIR.mkdir();
		}
		
		// Check auxiliary files
		for (String file : TEMPLATE_DIRS) {
			File auxDir = new File(LOG_DIR.getPath() + File.separator + file);
			if (!auxDir.exists()) {
				auxDir.mkdirs();
			}
		}

		for (String file : TEMPLATE_FILES) {
			File auxFile = new File(LOG_DIR.getPath() + File.separator + file);
			if (!auxFile.exists()) {
				File inFile = new File(RESOURCES_DIR + File.separator
						+ file);
				if (inFile.exists()) {
					copyFile(inFile, auxFile);
				}
			}
		}
		
		Configuration cfg = new Configuration();

		cfg.setDirectoryForTemplateLoading(new File(RESOURCES_DIR));

		Template tpl = cfg.getTemplate(template);
		OutputStreamWriter output = new FileWriter(outputFile);

		tpl.process(datamodel, output);
	}

	private static void fillInWithOptions(ApplicationOptions options) {

		StringBuffer arguments = new StringBuffer();
		for (String argument : ModelTest.arguments)
			arguments.append(argument + " ");
		datamodel.put("arguments", arguments);
		datamodel.put("alignName", options.getInputFile());
		datamodel.put("numTaxa", options.getNumTaxa());
		datamodel.put("seqLength", options.getNumSites());
		datamodel.put("phymlVersion", RunPhyml.PHYML_VERSION);
		datamodel.put("phymlBinary", Utilities.getBinaryVersion());
		datamodel.put("candidateModels", ModelTest.getCandidateModels().length);
		if (options.getSubstTypeCode() == 0)
			datamodel.put("substSchemes", "3");
		else if (options.getSubstTypeCode() == 1)
			datamodel.put("substSchemes", "5");
		else if (options.getSubstTypeCode() == 2)
			datamodel.put("substSchemes", "7");
		else
			datamodel.put("substSchemes", "11");

		datamodel
				.put("includeF", options.doF ? new Integer(1) : new Integer(0));
		datamodel
				.put("includeG", options.doG ? new Integer(1) : new Integer(0));
		datamodel
				.put("includeI", options.doI ? new Integer(1) : new Integer(0));
		datamodel.put("isAIC", options.doAIC ? new Integer(1) : new Integer(0));
		datamodel.put("isAICc", options.doAICc ? new Integer(1)
				: new Integer(0));
		datamodel.put("isBIC", options.doBIC ? new Integer(1) : new Integer(0));
		datamodel.put("isDT", options.doDT ? new Integer(1) : new Integer(0));
		datamodel.put("numCat", options.numGammaCat);
		datamodel.put("isPAUP", options.writePAUPblock ? new Integer(1) : new Integer(0));

		StringBuffer optimizedParameters = new StringBuffer(
				"Substitution parameters ");
		if (options.countBLasParameters)
			optimizedParameters.append("+ " + options.getNumBranches()
					+ " branch lengths ");
		if (options.optimizeMLTopology)
			optimizedParameters.append("+ topology");
		datamodel.put("freeParameters", optimizedParameters.toString());

		datamodel.put("userTreeDef",
				options.userTopologyExists ? new Integer(1) : new Integer(0));
		if (options.fixedTopology)
			datamodel.put("baseTree", "Fixed BioNJ");
		else if (options.optimizeMLTopology)
			datamodel.put("baseTree", "Maximum Likelihood");
		else if (options.userTopologyExists) {
			datamodel.put("baseTree", "Fixed user tree topology");
			datamodel.put("userTreeFilename", options.getInputTreeFile()
					.getName());
			datamodel.put("userTree", options.getUserTree());
		} else
			datamodel.put("baseTree", "BioNJ");

		switch (options.treeSearchOperations) {
		case NNI:
			datamodel.put("searchAlgorithm", "NNI");
			break;
		case SPR:
			datamodel.put("searchAlgorithm", "SPR");
			break;
		case BEST:
			datamodel.put("searchAlgorithm", "Best of {NNI, SPR}");
			break;
		}

		datamodel.put("confidenceInterval", String.format(Locale.ENGLISH,
				"%5.2f", options.confidenceInterval * 100));
	}

	private static void fillInWithSortedModels(Model[] models) {

		Collection<Map<String, String>> sortedModels = new ArrayList<Map<String, String>>();
		int index = 1;
		for (Model model : models) {
			Map<String, String> modelMap = new HashMap<String, String>();
			modelMap.put("index", String.valueOf(index++));
			modelMap.put("name", model.getName());
			modelMap.put("partition", model.getPartition());
			modelMap.put("lnl",
					String.format(Locale.ENGLISH, "%5.4f", model.getLnL()));
			modelMap.put("k", String.valueOf(model.getK()));
			modelMap.put(
					"fA",
					model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getfA()) : "-");
			modelMap.put(
					"fC",
					model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getfC()) : "-");
			modelMap.put(
					"fG",
					model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getfG()) : "-");
			modelMap.put(
					"fT",
					model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getfT()) : "-");
			modelMap.put(
					"titv",
					model.ispT() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getTitv()) : "-");
			modelMap.put(
					"rA",
					model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getRa()) : "-");
			modelMap.put(
					"rB",
					model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getRb()) : "-");
			modelMap.put(
					"rC",
					model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getRc()) : "-");
			modelMap.put(
					"rD",
					model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getRd()) : "-");
			modelMap.put(
					"rE",
					model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getRe()) : "-");
			modelMap.put(
					"rF",
					model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getRf()) : "-");
			modelMap.put(
					"pInv",
					model.ispI() ? String.format(Locale.ENGLISH, "%5.4f",
							model.getPinv()) : "-");
			modelMap.put(
					"shape",
					model.ispG() ? String.format(Locale.ENGLISH, "%6.4f",
							model.getShape()) : "-");
			modelMap.put("tree", model.getTreeString());
			sortedModels.add(modelMap);
		}

		datamodel.put("sortedModels", sortedModels);
	}

	private static void fillInWithTopologies(TreeSummary summary, ApplicationOptions options) {
		datamodel.put("numberOfTopologies", summary.getNumberOfTopologies());
		Collection<Map<String, String>> sortedTopologies = new ArrayList<Map<String, String>>();
		for (int index=0; index<summary.getNumberOfTopologies(); index++) {
			Map<String, String> topologyMap = new HashMap<String, String>();
			topologyMap.put("index", String.valueOf(index));
			Tree topology = summary.getTopology(index);
			topologyMap.put("tree", TreeUtilities.toNewick(topology, false, false, false));
			List<Model> models = summary.getModelsByTopology(index);
			StringBuilder sb = new StringBuilder();
			for (Model model : models) {
				sb.append(model.getName() + " ");
			}
			topologyMap.put("models", sb.toString());
			topologyMap.put("models", sb.toString());
			
			if (options.doAIC) {
				topologyMap.put("aicRank", String.valueOf(summary.aicIndexOf(topology)));
				topologyMap.put("aicRF", String.valueOf(summary.aicRfOf(topology)));
				topologyMap.put("aicAvgDistance", Utilities.format(summary.aicAvgDistance(topology), 6, 4, true));
				topologyMap.put("aicVarDistance", Utilities.format(summary.aicVarDistance(topology), 6, 4, true));
				topologyMap.put("aicWeight", Utilities.format(summary.aicWeight(topology), 6, 4, false));
			}
			
			if (options.doAICc) {
				topologyMap.put("aiccRank", String.valueOf(summary.aiccIndexOf(topology)));
				topologyMap.put("aiccRF", String.valueOf(summary.aiccRfOf(topology)));
				topologyMap.put("aiccAvgDistance", Utilities.format(summary.aiccAvgDistance(topology), 6, 4, true));
				topologyMap.put("aiccVarDistance", Utilities.format(summary.aiccVarDistance(topology), 6, 4, true));
				topologyMap.put("aiccWeight", Utilities.format(summary.aiccWeight(topology), 6, 4, false));
			}
			
			if (options.doBIC) {
				topologyMap.put("bicRank", String.valueOf(summary.bicIndexOf(topology)));
				topologyMap.put("bicRF", String.valueOf(summary.bicRfOf(topology)));
				topologyMap.put("bicAvgDistance", Utilities.format(summary.bicAvgDistance(topology), 6, 4, true));
				topologyMap.put("bicVarDistance", Utilities.format(summary.bicVarDistance(topology), 6, 4, true));
				topologyMap.put("bicWeight", Utilities.format(summary.bicWeight(topology), 6, 4, false));
			}
			
			if (options.doDT) {
				topologyMap.put("dtRank", String.valueOf(summary.dtIndexOf(topology)));
				topologyMap.put("dtRF", String.valueOf(summary.dtRfOf(topology)));
				topologyMap.put("dtAvgDistance", Utilities.format(summary.dtAvgDistance(topology), 6, 4, true));
				topologyMap.put("dtVarDistance", Utilities.format(summary.dtVarDistance(topology), 6, 4, true));
				topologyMap.put("dtWeight", Utilities.format(summary.dtWeight(topology), 6, 4, false));
			}
			
			sortedTopologies.add(topologyMap);
		}
		
		datamodel.put("sortedTopologies", sortedTopologies);
	}
	
	private static void fillInWIthInformationCriterion(InformationCriterion ic,
			Collection<Map<String, String>> sortedModels,
			Map<String, String> bestModel) {

		Model model = ic.getModel(0);
		bestModel.put("index", String.valueOf(1));
		bestModel.put("name", model.getName());
		bestModel.put("partition", model.getPartition());
		bestModel.put("lnl",
				String.format(Locale.ENGLISH, "%5.4f", model.getLnL()));
		bestModel.put("k", String.valueOf(model.getK()));
		bestModel.put(
				"fA",
				model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getfA()) : "-");
		bestModel.put(
				"fC",
				model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getfC()) : "-");
		bestModel.put(
				"fG",
				model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getfG()) : "-");
		bestModel.put(
				"fT",
				model.ispF() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getfT()) : "-");
		bestModel.put(
				"titv",
				model.ispT() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getTitv()) : "-");
		bestModel.put(
				"rA",
				model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getRa()) : "-");
		bestModel.put(
				"rB",
				model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getRb()) : "-");
		bestModel.put(
				"rC",
				model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getRc()) : "-");
		bestModel.put(
				"rD",
				model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getRd()) : "-");
		bestModel.put(
				"rE",
				model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getRe()) : "-");
		bestModel.put(
				"rF",
				model.ispR() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getRf()) : "-");
		bestModel.put(
				"pInv",
				model.ispI() ? String.format(Locale.ENGLISH, "%5.4f",
						model.getPinv()) : "-");
		bestModel.put(
				"shape",
				model.ispG() ? String.format(Locale.ENGLISH, "%6.4f",
						model.getShape()) : "-");
		bestModel.put("value",
				String.format(Locale.ENGLISH, "%5.4f", ic.getValue(model)));
		bestModel.put("delta",
				String.format(Locale.ENGLISH, "%5.4f", ic.getDelta(model)));
		bestModel.put("weight",
				String.format(Locale.ENGLISH, "%5.4f", ic.getWeight(model)));
		bestModel.put("tree", model.getTreeString());
		bestModel.put("cumWeight",
				String.format(Locale.ENGLISH, "%5.4f", ic.getCumWeight(model)));
		sortedModels.add(bestModel);
		for (int i = 1; i < ic.getNumModels(); i++) {
			model = ic.getModel(i);
			Map<String, String> modelMap = new HashMap<String, String>();
			modelMap.put("index", String.valueOf(i + 1));
			modelMap.put("name", model.getName());
			modelMap.put("partition", model.getPartition());
			modelMap.put("lnl",
					String.format(Locale.ENGLISH, "%5.4f", model.getLnL()));
			modelMap.put("k", String.valueOf(model.getK()));
			modelMap.put("value",
					String.format(Locale.ENGLISH, "%5.4f", ic.getValue(model)));
			modelMap.put("delta",
					String.format(Locale.ENGLISH, "%5.4f", ic.getDelta(model)));
			modelMap.put("weight",
					String.format(Locale.ENGLISH, "%5.4f", ic.getWeight(model)));
			modelMap.put(
					"cumWeight",
					String.format(Locale.ENGLISH, "%5.4f",
							ic.getCumWeight(model)));
			modelMap.put("tree", model.getTreeString());
			sortedModels.add(modelMap);
		}
	}

	@SuppressWarnings("resource")
	public static void copyFile(File in, File out) throws IOException {
		FileChannel inChannel = new FileInputStream(in).getChannel();
		FileChannel outChannel = new FileOutputStream(out).getChannel();
		try {
			inChannel.transferTo(0, inChannel.size(), outChannel);
		} catch (IOException e) {
			throw e;
		} finally {
			if (inChannel != null)
				inChannel.close();
			if (outChannel != null)
				outChannel.close();
		}
	}

	private static void buildChart(File mOutputFile, InformationCriterion ic) {
		
		/* This line prevents Swing exceptions on headless environments */
		if (GraphicsEnvironment.isHeadless()) { return; }
		
		int width = 500;
		int height = 300;
		try {
			if (!IMAGES_DIR.exists()) {
				IMAGES_DIR.mkdir();
			}
			ChartUtilities.saveChartAsPNG(new File(IMAGES_DIR.getPath() + File.separator + mOutputFile.getName() + "_rf_" + ic + ".png"), 
					RFHistogram.buildRFHistogram(ic),
					width, height);
			ChartUtilities.saveChartAsPNG(new File(IMAGES_DIR.getPath() + File.separator + mOutputFile.getName() + "_eu_" + ic + ".png"), 
					RFHistogram.buildEuclideanHistogram(ic),
					width, height);
		} catch (IOException e) {
//			e.printStackTrace();
		}
	}
} public class TextOutputStream extends OutputStream
{  

   public TextOutputStream(String s, boolean append)
   {
          try
          {
			out = new PrintWriter(new FileOutputStream(s,append));
          }
          catch(IOException e)
          {
                  out = null;
          }
   }

   public TextOutputStream(PrintStream s)
   {  out = new PrintWriter(s);
   }

   /**
    * Constructs a text output stream that sends its output to a file 
    * @param s the file name
    */

   public TextOutputStream(String s)
   {  try
      {  out = new PrintWriter(new FileOutputStream(s));
      }
      catch(IOException e)
      {  out = null;
      }
   }


   /**
    * Tests if the output stream is no longer valid
    * @return <tt>true</tt> if the stream has failed
    */

   public boolean fail()
   {  return out == null;
   }

   public boolean checkError()
   {
          return out.checkError();
   }

   public void flush()
   {
          out.flush();
   }

   /**
    * Closes the stream 
    */

   public void close()
   {  out.close();
   }

   /**
    * Writes one byte to the stream
    * @param b the byte to write
    */

   public void write(int b)
   {  out.write(b);
   }

   /**
    * Prints a string to a stream
    * @param s the string to print
    */

   public void print(String s) 
   { 
          out.print(s); 
          flush();
   }

   /**
    * Prints a character to a stream
    * @param c the character to print
    */

   public void print(char c) { out.print(c); flush(); }

   /**
    * Prints an object to a stream. The <tt>toString</tt> method of the object is invoked
    * @param o the object to print
    */

   public void print(Object o) { out.print(o);  flush();}

   /**
    * Prints a long integer to a stream
    * @param l the number to print
    */

   public void print(long l) { out.print(l); flush(); }

   /**
    * Prints a double precision floating point number to a stream
    * @param l the number to print
    */

   public void print(double x) { out.print(x);  flush();}

   /**
    * Prints a string to a stream, followed by a newline
    * @param s the string to print
    */

   public void println(String s) 
   { 
      out.println(s);
          flush();
   }

   /**
    * Prints a character to a stream, followed by a newline
    * @param c the character to print
    */

   public void println(char c) { out.println(c); flush(); }

   /**
    * Prints an object to a stream, followed by a newline. The <tt>toString</tt> method of the object is invoked
    * @param o the object to print
    */

   public void println(Object o) { out.println(o); flush(); }

   /**
    * Prints a long integer to a stream, followed by a newline
    * @param l the number to print
    */

   public void println(long l) { out.println(l); flush(); }

   /**
    * Prints a double precision floating point number to a stream, followed by a newline
    * @param l the number to print
    */

   public void println(double x) { out.println(x); flush(); }

   /**
    * Prints a newline
    */

   public void println() { out.println(); flush(); }

   private void parseFormat(String s)
   {  width = 0;
      precision = -1;
      pre = "";
      post = "";
      leading_zeroes = false;
      show_plus = false;
      alternate = false;
      show_space = false;
      left_align = false;
      fmt = ' ';

      int length = s.length();
      int parse_state = 0;
      // 0 = prefix, 1 = flags, 2 = width, 3 = precision,
      // 4 = format, 5 = end
      int i = 0;

      while (parse_state == 0)
      {  if (i >= length) parse_state = 5;
         else if (s.charAt(i) == '%')
         {  if (i < length - 1)
            {  if (s.charAt(i + 1) == '%')
               {  pre = pre + '%';
                  i++;
               }
               else
                  parse_state = 1;
            }
            else throw new java.lang.IllegalArgumentException();
         }
         else
            pre = pre + s.charAt(i);
         i++;
      }
      while (parse_state == 1)
      {  if (i >= length) parse_state = 5;
         else if (s.charAt(i) == ' ') show_space = true;
         else if (s.charAt(i) == '-') left_align = true;
         else if (s.charAt(i) == '+') show_plus = true;
         else if (s.charAt(i) == '0') leading_zeroes = true;
         else if (s.charAt(i) == '#') alternate = true;
         else { parse_state = 2; i--; }
         i++;
      }
      while (parse_state == 2)
      {  if (i >= length) parse_state = 5;
         else if ('0' <= s.charAt(i) && s.charAt(i) <= '9')
         {  width = width * 10 + s.charAt(i) - '0';
            i++;
         }
         else if (s.charAt(i) == '.')
         {  parse_state = 3;
            precision = 0;
            i++;
         }
         else
            parse_state = 4;
      }
      while (parse_state == 3)
      {  if (i >= length) parse_state = 5;
         else if ('0' <= s.charAt(i) && s.charAt(i) <= '9')
         {  precision = precision * 10 + s.charAt(i) - '0';
            i++;
         }
         else
            parse_state = 4;
      }
      if (parse_state == 4)
      {  if (i >= length) parse_state = 5;
         else fmt = s.charAt(i);
         i++;
      }
      if (i < length)
         post = s.substring(i, length);
   
   
   }

   /**
    * Prints a double precision floating point number to a stream, with formatting instructions
    * @param fmt the formatting instruction string. 
    * The string has a prefix, a format code and a suffix. The prefix and suffix
    * become part of the formatted output. The format code directs the
    * formatting of the parameter to be formatted. The code has the
    * following structure
    * <ul>
    * <li> a % (required)
    * <li> a modifier (optional)
    * <dl>
    * <dt> + <dd> forces display of + for positive numbers
    * <dt> 0 <dd> show leading zeroes
    * <dt> - <dd> align left in the field
    * <dt> space <dd> prepend a space in front of positive numbers
    * <dt> # <dd> use "alternate" format: Don't suppress trailing zeroes in general floating point format.
    * </dl>
    * <li> an integer denoting field width (optional)
    * <li> a period followed by an integer denoting precision (optional)
    * <li> a format descriptor (required)
    * <dl>
    * <dt>f <dd> floating point number in fixed format
    * <dt>e, E <dd> floating point number in exponential notation (scientific format). The E format results in an uppercase E for the exponent (1.14130E+003), the e format in a lowercase e.
    * <dt>g, G <dd> floating point number in general format (fixed format for small numbers, exponential format for large numbers). Trailing zeroes are suppressed. The G format results in an uppercase E for the exponent (if any), the g format in a lowercase e.
    * </dl>
    * </ul>
    * @param l the number to print
    */

   public void printf(String fmt, double x)
   {
 
      parseFormat(fmt);
      print(form(x));
 
   }

   /**
    * Prints a long integer to a stream, with formatting instructions
    * @param fmt the formatting instruction string. 
    * The string has a prefix, a format code and a suffix. The prefix and suffix
    * become part of the formatted output. The format code directs the
    * formatting of the parameter to be formatted. The code has the
    * following structure
    * <ul>
    * <li> a % (required)
    * <li> a modifier (optional)
    * <dl>
    * <dt> + <dd> forces display of + for positive numbers
    * <dt> 0 <dd> show leading zeroes
    * <dt> - <dd> align left in the field
    * <dt> space <dd> prepend a space in front of positive numbers
    * <dt> # <dd> use "alternate" format: Add 0 or 0x for octal or hexadecimal numbers. 
    * </dl>
    * <li> an integer denoting field width (optional)
    * <li> a format descriptor (required)
    * <dl>
    * <dt>d, i <dd> integer in decimal
    * <dt>x <dd> integer in hexadecimal
    * <dt>o <dd> integer in octal
    * </dl>
    * </ul>
    * @param x the number to print
    */

   public void printf(String fmt, long x)
   {  parseFormat(fmt);
      print(form(x));
   }

   /**
    * Prints a character to a stream, with formatting instructions
    * @param fmt the formatting instruction string. 
    * The string has a prefix, a format code and a suffix. The prefix and suffix
    * become part of the formatted output. The format code directs the
    * formatting of the parameter to be formatted. The code has the
    * following structure
    * <ul>
    * <li> a % (required)
    * <li> a modifier (optional)
    * <dl>
    * <dt> + <dd> forces display of + for positive numbers
    * <dt> 0 <dd> show leading zeroes
    * <dt> - <dd> align left in the field
    * <dt> space <dd> prepend a space in front of positive numbers
    * <dt> # <dd> use "alternate" format: Add 0 or 0x for octal or hexadecimal numbers. 
    * </dl>
    * <li> an integer denoting field width (optional)
    * <li> a format descriptor (required)
    * <dl>
    * <dt>d, i <dd> integer in decimal
    * <dt>x <dd> integer in hexadecimal
    * <dt>o <dd> integer in octal
    * <dt>c <dd> character
    * </dl>
    * </ul>
    * @param x the character to print
    */

   public void printf(String fmt, char x)
   {  parseFormat(fmt);
      print(form(x));
   }

   /**
    * Prints a string to a stream, with formatting instructions
    * @param fmt the formatting instruction string. 
    * The string has a prefix, a format code and a suffix. The prefix and suffix
    * become part of the formatted output. The format code directs the
    * formatting of the parameter to be formatted. The code has the
    * following structure
    * <ul>
    * <li> a % (required)
    * <li> a modifier (optional)
    * <dl>
    * <dt> - <dd> align left in the field
    * </dl>
    * <li> an integer denoting field width (optional)
    * <li> a format descriptor (required)
    * <dl>
    * <dt>s <dd> string
    * </dl>
    * </ul>
    * @param x the character to print
    */

   public void printf(String fmt, String x)
   {  parseFormat(fmt);
      print(form(x));
   }

   private static String repeat(char c, int n)
   {  if (n <= 0) return "";
      StringBuffer s = new StringBuffer(n);
      for (int i = 0; i < n; i++) s.append(c);
      return s.toString();
   }

   private static String convert(long x, int n, int m, String d)
   {  if (x == 0) return "0";
      String r = "";
      while (x != 0)
      {  r = d.charAt((int)(x & m)) + r;
         x = x >>> n;
      }
      return r;
   }

   private String pad(String r)
   {  String p = repeat(' ', width - r.length());
      if (left_align) return pre + r + p + post;
      else return pre + p + r + post;
   }

   private String sign(int s, String r)
   {  String p = "";
      if (s < 0) p = "-";
      else if (s > 0)
      {  if (show_plus) p = "+";
         else if (show_space) p = " ";
      }
      else
      {  if (fmt == 'o' && alternate && r.length() > 0 && r.charAt(0) != '0') p = "0";
         else if (fmt == 'x' && alternate) p = "0x";
         else if (fmt == 'X' && alternate) p = "0X";
      }
      int w = 0;
      if (leading_zeroes)
         w = width;
      else if ((fmt == 'd' || fmt == 'i' || fmt == 'x' || fmt == 'X' || fmt == 'o')
         && precision > 0) w = precision;

      return p + repeat('0', w - p.length() - r.length()) + r;
   }


   private String fixed_format(double d)
   {  String f = "";
      if (d > 0x7FFFFFFFFFFFFFFFL) return exp_format(d);

      long l = (long)(precision == 0 ? d + 0.5 : d);
      f = f + l;

      double fr = d - l; // fractional part
      if (fr >= 1 || fr < 0) return exp_format(d);


	/* DP change 210205 */
	String fraction =  frac_part(fr);
	if (fraction == "0")
	{
	int myInteger = (int) d + 1;
	String trailing_zeroes = "";
    for (int i = 1; i <= precision; i++)
        trailing_zeroes = trailing_zeroes + "0";
	// return "1." + trailing_zeroes; /* DP 150506 error here, it would print 12.99999 as 1.0!*/
		return myInteger + "." + trailing_zeroes;
	  }
	else	
    	return f + fraction;

// return f + frac_part(fr);

	/* DP */
     	
   }


/* DP 210205 : there was a problem when the precision should round to 1.0. Instead it would print 0
i.e.,	printf("%10.4f",0.99998) would print 0.0000
*/

   private String frac_part(double fr)
   // precondition: 0 <= fr < 1
   {  String z = "";
      if (precision > 0)
      {  double factor = 1;
         String leading_zeroes = "";
         for (int i = 1; i <= precision && factor <= 0x7FFFFFFFFFFFFFFFL; i++)
         {  factor *= 10;
            leading_zeroes = leading_zeroes + "0";
         }
  
        long l = (long) (factor * fr + 0.5);
    
	/* DP 210205 added */
		if (l == factor)
			return "0";

  	     z = leading_zeroes + l;
         z = z.substring(z.length() - precision, z.length());

      }

      if (precision > 0 || alternate) z = "." + z;
      
      if ((fmt == 'G' || fmt == 'g') && !alternate)
      // remove trailing zeroes and decimal point
      {  int t = z.length() - 1;
         while (t >= 0 && z.charAt(t) == '0') t--;
         if (t >= 0 && z.charAt(t) == '.') t--;
         z = z.substring(0, t + 1);
      }
      return z;
   }

   private String exp_format(double d)
   {  String f = "";
      int e = 0;
      double dd = d;
      double factor = 1;
      if (d != 0)
      {  while (dd > 10) { e++; factor /= 10; dd = dd / 10; }
         while (dd < 1) { e--; factor *= 10; dd = dd * 10; }
      }
      if ((fmt == 'g' || fmt == 'G') && e >= -4 && e < precision)
         return fixed_format(d);

      d = d * factor;
      f = f + fixed_format(d);

      if (fmt == 'e' || fmt == 'g')
         f = f + "e";
      else
         f = f + "E";

      String p = "000";
      if (e >= 0)
      {  f = f + "+";
         p = p + e;
      }
      else
      {  f = f + "-";
         p = p + (-e);
      }

      return f + p.substring(p.length() - 3, p.length());
   }

   private String form(double x)
   {  String r;
      if (precision < 0) precision = 6;
      int s = 1;
      if (x < 0) { x = -x; s = -1; }
      if (fmt == 'f')
         r = fixed_format(x);
      else if (fmt == 'e' || fmt == 'E' || fmt == 'g' || fmt == 'G')
         r = exp_format(x);
      else throw new java.lang.IllegalArgumentException();

      return pad(sign(s, r));
   }

   private String form(long x)
   {  String r;
      int s = 0;
      if (fmt == 'd' || fmt == 'i')
      {  if (x < 0) 
         {  r = ("" + x).substring(1);
            s = -1; 
         }
         else 
         {  r = "" + x; 
            s = 1;
         }
      }
      else if (fmt == 'o')
         r = convert(x, 3, 7, "01234567");
      else if (fmt == 'x')
         r = convert(x, 4, 15, "0123456789abcdef");
      else if (fmt == 'X')
         r = convert(x, 4, 15, "0123456789ABCDEF");
      else throw new java.lang.IllegalArgumentException();

      return pad(sign(s, r));
   }

   private String form(char c)
   {  if (fmt != 'c')
         throw new java.lang.IllegalArgumentException();

      String r = "" + c;
      return pad(r);
   }

   private String form(String s)
   {  if (fmt != 's')
         throw new java.lang.IllegalArgumentException();
      if (precision >= 0) s = s.substring(0, precision);
      return pad(s);
   }

   private PrintWriter out;

   private int width;
   private int precision;
   private String pre;
   private String post;
   private boolean leading_zeroes;
   private boolean show_plus;
   private boolean alternate;
   private boolean show_space;
   private boolean left_align;
   private char fmt; // one of cdeEfgGiosxXos
} public class Model implements Serializable {

	private static final long serialVersionUID = 4173257281316945987L;

	private String name, partition, treeString;
	private Tree tree;

	private int id;
	private int K;
	private int numTi;
	private int numTv;
	private int numGammaCat;

	private boolean pF, pT, pR, pI, pG, pV;

	private boolean isInAICinterval, isInAICcinterval, isInBICinterval,
			isInDTinterval;

	private double lnL;
	/* Unconstrained LK */
	private double unconstrainedLnL;
	private double lnLIgnoringGaps;

	/* statistical */
	private double AIC, AICd, AICw, cumAICw, uAICd;
	private double AICc, AICcd, AICcw, cumAICcw, uAICcd;
	private double BIC, BICd, BICw, cumBICw, uBICd;
	private double DT, DTd, DTw, cumDTw;

	private double fA, fG, fC, fT, titv, kappa;
	private double Ra, Rb, Rc, Rd, Re, Rf, pinv, shape;

	private long computationTime;

	public Model(int mid, String mname, String mpartition, int mparameters,
			boolean mpF, boolean mpT, boolean mpV, boolean mpR, boolean mpI,
			boolean mpG, int mTi, int mTv) {

		id = mid;
		name = mname;
		partition = mpartition; // model substitution specification (e.g. GTR =
								// 012345)
		K = mparameters;
		pF = mpF; // includes unequal base frequencies
		pT = mpT; // includes ti/tv
		pV = mpV; // only for phyml: includes a Ti/Tv, apart from R(x)
		pR = mpR; // includes Rates R(x)
		pI = mpI; // includes p-inv
		pG = mpG; // includes gamma
		numTi = mTi; // number of transition rates
		numTv = mTv; // number of transversion rates

		lnL = 0.0;
		pinv = 0.0;
		shape = ModelTest.INFINITY;
	}

	public Model(int mid, String mname, String mpartition, int mparameters) {

		id = mid;
		name = mname;
		partition = mpartition; // model substitution specification (e.g. GTR =
								// 012345)
		K = mparameters;
		pF = false; // includes unequal base frequencies
		pT = false; // includes ti/tv
		pV = false; // only for phyml: includes a Ti/Tv, apart from R(x)
		pR = false; // includes Rates R(x)
		pI = false; // includes p-inv
		pG = false; // includes gamma
		numTi = 0; // number of transition rates
		numTv = 0; // number of transversion rates

		lnL = 0.0;
		pinv = 0.0;
		shape = ModelTest.INFINITY;
	}

	/****************************
	 * print ************************************ * Print model components * *
	 ************************************************************************/

	public void print(TextOutputStream stream) {

		stream.println("   Model = " + getName());
		stream.println("   partition = " + getPartition());

		if (getLnL() == 0) {
			stream.println(" OPTIMIZATION FAILED!");
		} else {
			stream.printf("   -lnL = %.4f", getLnL());
			stream.printf("\n   K = %d", getK());
			/*
			 * stream.printf ("\n     optimized substitution parameters = %d",
			 * Model.freeParameters[id]); if (ModelTest.countBLasParameters ==
			 * true) stream.printf ("\n     optimized branch lengths = %d",
			 * ModelTest.numBranches); if (ModelTest.optimizeMLTopology)
			 * stream.printf ("\n     optimized topology = %d", 1);
			 */
			if (ispF()) {
				stream.printf("\n   freqA = %5.4f ", getfA());
				stream.printf("\n   freqC = %5.4f ", getfC());
				stream.printf("\n   freqG = %5.4f ", getfG());
				stream.printf("\n   freqT = %5.4f ", getfT());
			}
			if (ispT()) {
				stream.printf("\n   kappa = %.4f", getKappa());
				stream.printf(" (ti/tv = %.4f)", getTitv());
			}
			if (ispR()) {
				stream.printf("\n   R(a) [AC] =  %.4f", getRa());
				stream.printf("\n   R(b) [AG] =  %.4f", getRb());
				stream.printf("\n   R(c) [AT] =  %.4f", getRc());
				stream.printf("\n   R(d) [CG] =  %.4f", getRd());
				stream.printf("\n   R(e) [CT] =  %.4f", getRe());
				stream.printf("\n   R(f) [GT] =  %.4f", 1.0);
			}
			if (ispI())
				stream.printf("\n   p-inv = %5.4f", getPinv());
			if (ispG()) {
				stream.printf("\n   gamma shape = %6.4f", getShape());
			}
			stream.println(" ");
		}
	}

	public String printForTesting() {
		StringBuilder str = new StringBuilder();
		str.append(getName() + " ");
		str.append(ispF() ? (getfA() + " " + getfC() + " " + getfG() + " "
				+ getfT() + " ") : ("NA NA NA NA "));
		str.append(ispT() ? (getKappa() + " " + getTitv() + " ") : ("NA NA "));

		str.append(ispR() ? (getRa() + " " + getRb() + " " + getRc() + " "
				+ getRd() + " " + getRe() + " " + getRf() + " ")
				: ("NA NA NA NA NA NA "));
		str.append(ispI() ? (getPinv() + " ") : ("NA "));
		str.append(ispG() ? (getShape() + " ") : ("NA "));
		str.append(getPartition());
		return str.toString();
	}

	public void setCumBICw(double cumBICw) {
		this.cumBICw = cumBICw;
	}

	public double getCumBICw() {
		return cumBICw;
	}

	public void setLnL(double lnL) {
		this.lnL = lnL;
	}

	public double getLnL() {
		return lnL;
	}

	public void setLnLIgnoringGaps(double lnLIgnoringGaps) {
		this.lnLIgnoringGaps = lnLIgnoringGaps;
	}

	public double getLnLIgnoringGaps() {
		return lnLIgnoringGaps;
	}

	public double getUnconstrainedLnL() {
		return unconstrainedLnL;
	}

	public void setUnconstrainedLnL(double unconstrainedLnL) {
		this.unconstrainedLnL = unconstrainedLnL;
	}

	public void setAIC(double aIC) {
		AIC = aIC;
	}

	public double getAIC() {
		return AIC;
	}

	public void setAICd(double aICd) {
		AICd = aICd;
	}

	public double getAICd() {
		return AICd;
	}

	public void setAICw(double aICw) {
		AICw = aICw;
	}

	public double getAICw() {
		return AICw;
	}

	public void setUAICd(double uAICd) {
		this.uAICd = uAICd;
	}

	public double getUAICd() {
		return uAICd;
	}

	public void setCumAICw(double cumAICw) {
		this.cumAICw = cumAICw;
	}

	public double getCumAICw() {
		return cumAICw;
	}

	public void setAICc(double aICc) {
		AICc = aICc;
	}

	public double getAICc() {
		return AICc;
	}

	public void setAICcd(double aICcd) {
		AICcd = aICcd;
	}

	public double getAICcd() {
		return AICcd;
	}

	public void setAICcw(double aICcw) {
		AICcw = aICcw;
	}

	public double getAICcw() {
		return AICcw;
	}

	public void setCumAICcw(double cumAICcw) {
		this.cumAICcw = cumAICcw;
	}

	public void setUAICcd(double uAICcd) {
		this.uAICcd = uAICcd;
	}

	public double getUAICcd() {
		return uAICcd;
	}

	public double getCumAICcw() {
		return cumAICcw;
	}

	public void setBICd(double bICd) {
		BICd = bICd;
	}

	public double getBICd() {
		return BICd;
	}

	public void setBIC(double bIC) {
		BIC = bIC;
	}

	public double getBIC() {
		return BIC;
	}

	public void setBICw(double bICw) {
		BICw = bICw;
	}

	public double getBICw() {
		return BICw;
	}

	public void setUBICd(double uBICd) {
		this.uBICd = uBICd;
	}

	public double getUBICd() {
		return uBICd;
	}

	public void setDT(double dT) {
		DT = dT;
	}

	public double getDT() {
		return DT;
	}

	public void setCumDTw(double cumDTw) {
		this.cumDTw = cumDTw;
	}

	public double getCumDTw() {
		return cumDTw;
	}

	public void setDTw(double dTw) {
		DTw = dTw;
	}

	public double getDTw() {
		return DTw;
	}

	public void setDTd(double dTd) {
		DTd = dTd;
	}

	public double getDTd() {
		return DTd;
	}

	public String getName() {
		return name;
	}

	public void setPartition(String partition) {
		this.partition = partition;
	}

	public String getPartition() {
		return partition;
	}

	public Tree getTree() {
		return tree;
	}

	public void setTreeString(String tree) throws TreeParseException {
		this.treeString = tree;

		if (tree == null)
			throw new TreeParseException("Attempting to set a null tree");
		StringReader sr = new StringReader(tree);
		this.tree = new ReadTree(new PushbackReader(sr));
	}

	public String getTreeString() {
		return treeString;
	}

	public void setfG(double fG) {
		this.fG = fG;
	}

	public double getfG() {
		return fG;
	}

	public void setfA(double fA) {
		this.fA = fA;
	}

	public double getfA() {
		return fA;
	}

	public void setfC(double fC) {
		this.fC = fC;
	}

	public double getfC() {
		return fC;
	}

	public void setfT(double fT) {
		this.fT = fT;
	}

	public double getfT() {
		return fT;
	}

	public void setTitv(double titv) {
		this.titv = titv;
	}

	public double getTitv() {
		return titv;
	}

	public void setKappa(double kappa) {
		this.kappa = kappa;
	}

	public double getKappa() {
		return kappa;
	}

	public void setRa(double ra) {
		Ra = ra;
	}

	public double getRa() {
		return Ra;
	}

	public void setShape(double shape) {
		this.shape = shape;
	}

	public double getShape() {
		return shape;
	}

	public void setPinv(double pinv) {
		this.pinv = pinv;
	}

	public double getPinv() {
		return pinv;
	}

	public void setRf(double rf) {
		Rf = rf;
	}

	public double getRf() {
		return Rf;
	}

	public void setRe(double re) {
		Re = re;
	}

	public double getRe() {
		return Re;
	}

	public void setRd(double rd) {
		Rd = rd;
	}

	public double getRd() {
		return Rd;
	}

	public void setRc(double rc) {
		Rc = rc;
	}

	public double getRc() {
		return Rc;
	}

	public void setRb(double rb) {
		Rb = rb;
	}

	public double getRb() {
		return Rb;
	}

	public int getId() {
		return id;
	}

	public void setNumGammaCat(int numGammaCat) {
		this.numGammaCat = numGammaCat;
	}

	public int getNumGammaCat() {
		return numGammaCat;
	}

	public int getK() {
		return K;
	}

	public int getNumTi() {
		return numTi;
	}

	public boolean ispV() {
		return pV;
	}

	public boolean ispG() {
		return pG;
	}

	public boolean ispI() {
		return pI;
	}

	public boolean ispR() {
		return pR;
	}

	public boolean ispT() {
		return pT;
	}

	public boolean ispF() {
		return pF;
	}

	public void setInDTinterval(boolean isInDTinterval) {
		this.isInDTinterval = isInDTinterval;
	}

	public boolean isInDTinterval() {
		return isInDTinterval;
	}

	public void setInBICinterval(boolean isInBICinterval) {
		this.isInBICinterval = isInBICinterval;
	}

	public boolean isInBICinterval() {
		return isInBICinterval;
	}

	public void setInAICcinterval(boolean isInAICcinterval) {
		this.isInAICcinterval = isInAICcinterval;
	}

	public boolean isInAICcinterval() {
		return isInAICcinterval;
	}

	public void setInAICinterval(boolean isInAICinterval) {
		this.isInAICinterval = isInAICinterval;
	}

	public boolean isInAICinterval() {
		return isInAICinterval;
	}

	public int getNumTv() {
		return numTv;
	}

	@Override
	public int hashCode() {
		final int prime = 31;
		int result = 1;
		result = prime * result + (pF ? 1231 : 1237);
		result = prime * result + (pG ? 1231 : 1237);
		result = prime * result + (pI ? 1231 : 1237);
		result = prime * result
				+ ((partition == null) ? 0 : partition.hashCode());
		return result;
	}

	@Override
	public boolean equals(Object obj) {
		if (this == obj)
			return true;
		if (obj == null)
			return false;
		if (getClass() != obj.getClass())
			return false;
		Model other = (Model) obj;
		if (pF != other.pF)
			return false;
		if (pG != other.pG)
			return false;
		if (pI != other.pI)
			return false;
		if (partition == null) {
			if (other.partition != null)
				return false;
		} else if (!partition.equals(other.partition))
			return false;
		return true;
	}

	public long getComputationTime() {
		return computationTime;
	}

	public void setComputationTime(long computationTime) {
		this.computationTime = computationTime;
	}

	public void update(Model model) {
		if (this.equals(model)) {
			lnL = model.lnL;
			numGammaCat = model.numGammaCat;
			shape = model.shape;
			pinv = model.pinv;
			kappa = model.kappa;
			treeString = model.treeString;
			tree = model.tree;
			computationTime = model.computationTime;
		}
	}
} // class Model public class ConsoleProgressObserver implements Observer {

	private TextOutputStream stream;
	/** Timer for calculate the elapsed time **/
	private long startTime;
	private int totalModels;
	private int completedModels = 0;
	private boolean threadScheduling;
	private ApplicationOptions options;

	public ConsoleProgressObserver(ApplicationOptions options) {
		this.startTime = System.currentTimeMillis();
		this.stream = ModelTest.getMainConsole();
		this.threadScheduling = options.threadScheduling;
		this.options=options;
	}

	@Override
	@SuppressWarnings("fallthrough")
	public synchronized void update(Observable o, Object arg) {

		if (arg != null) {
			ProgressInfo info = (ProgressInfo) arg;

			switch (info.getType()) {

			case ProgressInfo.BASE_TREE_INIT:
				stream.println(" ");
				stream.println("Estimating a BIONJ-JC tree ... ");
				break;

			case ProgressInfo.BASE_TREE_COMPUTED:
				stream.println("OK");
				stream.println(info.getModel().getName() + " tree: "
						+ info.getModel().getTreeString());
				break;

			case ProgressInfo.OPTIMIZATION_INIT:
				this.totalModels = options.getNumModels();
				stream.println(" ");stream.println(" ");
				stream.println("::Progress::");
				stream.println(" ");
				stream.println("Model \t\t Exec. Time \t Total Time \t -lnL");
				stream.println("-------------------------------------------------------------------------");
				break;
			
			case ProgressInfo.GTR_OPTIMIZATION_INIT:
				stream.println("[Heuristic search] Optimizing " + info.getModel().getName() + " model");
				break;
				
			case ProgressInfo.GTR_OPTIMIZATION_COMPLETED:
				stream.println("[Heuristic search] Starting model filtering");
				break;
				
			case ProgressInfo.SINGLE_OPTIMIZATION_INIT:
				break;

			case ProgressInfo.SINGLE_OPTIMIZATION_COMPLETED:
				completedModels++;
				stream.print(info.getModel().getName() + "\t");
				double modelLnL;
				if (info.getValue() == ProgressInfo.VALUE_REGULAR_OPTIMIZATION) {
					modelLnL = info.getModel().getLnL(); 
				} else {
					modelLnL = info.getModel().getLnLIgnoringGaps();
				}
				if (info.getModel().getName().length()<8)
					stream.print("\t");
				stream.print(info.getMessage() + "\t" 
						+ Utilities.calculateRuntime(startTime, System.currentTimeMillis()) + "\t" 
						+ String.format(Locale.ENGLISH, "%13.4f", modelLnL));
				if (ModelTest.MPJ_RUN && threadScheduling) {
					stream.println(" ");
				} else {
					if (info.isHeuristicSearch() && info.getValue() == ProgressInfo.VALUE_REGULAR_OPTIMIZATION) {
						stream.println("\t ["+info.getHeuristicStage()+"/6] (" + completedModels + "/" + info.getNumModelsInStage() + ")");
						if (completedModels == info.getNumModelsInStage()) {
							completedModels = 0;
						}
					} else {
						stream.println("\t (" + completedModels + "/" + totalModels + ")");
					}
				}
				
				if (ModelTestConfiguration.isCkpEnabled()) {
					try {
						OutputStream file = new FileOutputStream(
								options.getCkpFile());
						OutputStream buffer = new BufferedOutputStream(file);
						ObjectOutput output = new ObjectOutputStream(buffer);
						output.writeObject(ModelTest.getCandidateModels());
						output.close();
					} catch (IOException ex) {
						System.err.println("Cannot perform output.");
					}
				}
				
				break;
				
			case ProgressInfo.REOPTIMIZATION_INIT:
				this.totalModels = info.getValue();
				this.completedModels = 0;
				this.startTime = System.currentTimeMillis();
				stream.println(" ");stream.println(" ");
				stream.println("Some models should be reoptimized for checking lnL against the unconstrained likelihood");
				stream.println(" ");
				stream.println("Model \t\t Exec. Time \t Total Time\t-lnL w/o gaps");
				stream.println("-------------------------------------------------------------------------");
				break;
				
			case ProgressInfo.REOPTIMIZATION_COMPLETED:

				stream.println(" ");

				stream.println("  Unconstrained -lnL       = " + options.getUnconstrainedLnL());
				stream.println("  Number of patterns found = " + options.getNumPatterns());

				stream.println(" ");
				break;
				
			case ProgressInfo.INTERRUPTED:
				stream.println(" ");
				stream.println("Computation of likelihood scores discontinued ...");
				break;

			case ProgressInfo.ERROR:
				stream.println(info.getMessage());
				stream.println(" ");
				stream.println("Computation of likelihood scores discontinued ...");
				System.exit(-1);
				break;
			case ProgressInfo.ERROR_BINARY_NOEXISTS:
				stream.println("");
				stream.println("ERROR: PhyML binary does not exists: " + info.getMessage());
				stream.println("");
				ModelTest.finalize(ProgressInfo.ERROR_BINARY_NOEXISTS);
				break;
			case ProgressInfo.ERROR_BINARY_NOEXECUTE:
				stream.println("");
				stream.println("ERROR: PhyML binary does not have execution permission: " + info.getMessage());
				stream.println("");
				ModelTest.finalize(ProgressInfo.ERROR_BINARY_NOEXECUTE);
				break;
			case ProgressInfo.OPTIMIZATION_COMPLETED_OK:

				stream.println(" ");
				for (Model model : ModelTest.getCandidateModels()) {
					model.print(ModelTest.getMainConsole());
					ModelTest.getMainConsole().println(" ");
				}

				if (options.isAmbiguous()) {
				stream.println("  Best-fit models should be reoptimized for comparison with unconstrained likelihood");
				} else {
					stream.println("  Unconstrained -lnL       = " + options.getUnconstrainedLnL());
					stream.println("  Number of patterns found = " + options.getNumPatterns());
				}
				stream.println(" ");
				stream.println("Computation of likelihood scores completed. It took "
						+ Utilities.calculateRuntime(startTime,
								System.currentTimeMillis()) + ".");
				stream.println(" ");
				// continue

			case ProgressInfo.OPTIMIZATION_COMPLETED_INTERRUPTED:
				// dispose
				break;
			}
		}

	}

} public class AICc extends InformationCriterion {

	// constructor
	public AICc(boolean mwritePAUPblock, boolean mdoImportances,
			boolean mdoModelAveraging, double minterval) {
		super(mwritePAUPblock, mdoImportances, mdoModelAveraging, minterval);
	}

	/****************************
	 * computeAICc ******************************* * Computes the corrected AIC
	 * (AICc) for small sample sizes (relative * to the number of parameters)
	 * and finds out the model with the * minimum AICc * *
	 ************************************************************************/

	public void compute() {

		boolean sorted;
		int temp2, pass;
		double[] tempw = new double[options.getNumModels()];
		double min, sumExp, cum, temp1;

		// Calculate AICc and minAICc
		min = computeAicc(models[0], options);

		if (doCheckAgainstULK) {
			unconstrainedModel.setAICc(computeSingle(unconstrainedModel));
		}
		minModel = models[0];
		for (Model model : models) {
			model.setAICc(computeAicc(model, options));

			if (model.getAICc() < min) {
				min = model.getAICc();
				minModel = model;
			}
			
			if (doCheckAgainstULK) {
				model.setUAICcd(model.getAICc() - unconstrainedModel.getAICc());
			}
		}

		validResults = min > 0;
		
		// Calculate Akaike differences
		sumExp = 0;
		for (int i = 0; i < numModels; i++) {
			models[i].setAICcd(models[i].getAICc() - minModel.getAICc());
			sumExp += Math.exp(-0.5 * models[i].getAICcd());
		}

		// Calculate Akaike weights
		for (int i = 0; i < numModels; i++) {
			if (models[i].getAICcd() > 1000)
				models[i].setAICcw(0.0);
			else
				models[i].setAICcw(Math.exp(-0.5 * models[i].getAICcd()) / sumExp);
			tempw[i] = models[i].getAICc();
			order[i] = i;
		}

		// Get order for Akaike weights and calculate cumWeigths
		sorted = false;
		pass = 1;
		while (!sorted) {
			sorted = true;
			for (int i = 0; i < (numModels - pass); i++)
				if (tempw[i] > tempw[i + 1]) {
					temp1 = tempw[i + 1];
					tempw[i + 1] = tempw[i];
					tempw[i] = temp1;

					temp2 = order[i + 1];
					order[i + 1] = order[i];
					order[i] = temp2;

					sorted = false;
				}
			pass++;
		}

		cum = 0;
		for (int i = 0; i < numModels; i++) {
			cum += models[order[i]].getAICcw();
			models[order[i]].setCumAICcw(cum);
		}

		// confidence interval
		buildConfidenceInterval();

		// parameter importances
		if (doImportances || doModelAveraging)
			parameterImportance();

		// model averaging
		if (doModelAveraging)
			averageModels();

	}

	public double computeSingle(Model model) {
		return computeAicc(model, options);
	}
	
	public static double computeAicc(Model model, ApplicationOptions options) {
		int K;
		if (options.countBLasParameters)
			K = model.getK();
		else
			K = model.getK() - options.getNumBranches();

		double aicc;
		if ((options.getSampleSize() - K - 1) == 0) {
			aicc = 0;
		} else {
			aicc = 2 * (model.getLnL() + K)
				+ ((double) 2 * K * (K + 1) / (options.getSampleSize() - K - 1));
		}
		return aicc;
	}
	
	public static double computeAicc(double lnL, int k, ApplicationOptions options) {
		int K;
		if (options.countBLasParameters)
			K = k;
		else
			K = k - options.getNumBranches();

		double aicc;
		if ((options.getSampleSize() - K - 1) == 0) {
			aicc = 0;
		} else {
			aicc = 2 * (lnL + K)
				+ ((double) 2 * K * (K + 1) / (options.getSampleSize() - K - 1));
		}
		System.out.println("AICC 2 IS " + aicc);
		return aicc;
	}
	
	protected void printHeader(TextOutputStream stream) {
		stream.println("\n\n\n---------------------------------------------------------------");
		stream.println("*                                                             *");
		stream.println("*        CORRECTED AKAIKE INFORMATION CRITERION (AICc)        *");
		stream.println("*                                                             *");
		stream.println("---------------------------------------------------------------");
		stream.println(" ");
		stream.println(" Sample size: " + options.getSampleSize());
	}
	
	protected void printFooter(TextOutputStream stream) {
		stream.println("-lnL:\tnegative log likelihod");
		stream.println(" K:\tnumber of estimated parameters");
		stream.println(" AICc:\tCorrected Akaike Information Criterion");
		stream.println(" delta:\tAICc difference");
		stream.println(" weight:\tAICc weight");
		stream.println(" cumWeight:\tcumulative AICc weight");
	}
	

	/**************
	 * buildConfidenceInterval ************************
	 * 
	 * Builts the confidence interval of selected models and their cumulative
	 * weight
	 * 
	 * The model that just passed the confidence will be or not in the interval
	 * by chance (see below)
	 ****************************************************************/

	public void buildConfidenceInterval() {
		int i;
		Model tmodel = models[0];
		cumWeight = 0;

		// first construct the confidence interval for models
		if (confidenceInterval >= 1.0d) {
			for (i = 0; i < numModels; i++) {
				tmodel = models[order[i]];
				tmodel.setInAICcinterval(true);
				confidenceModels.add(tmodel);
			}
			cumWeight = 1.0;
		} else {
			for (i = 0; i < numModels; i++) {
				tmodel = models[order[i]];
				if (tmodel.getCumAICcw() <= confidenceInterval) {
					tmodel.setInAICcinterval(true);
					confidenceModels.add(tmodel);
					cumWeight += tmodel.getAICcw();
				} else
					break;
			}

			// lets decide whether the model that just passed the confidence
			// interval should be included (suggested by John Huelsenbeck)
			double probOut = (tmodel.getCumAICcw() - confidenceInterval)
					/ tmodel.getAICcw();
			double probIn = 1.0 - probOut;
			Random generator = new Random();
			double randomNumber = generator.nextDouble();
			if (randomNumber <= probIn) {
				tmodel.setInAICcinterval(true);
				if (!confidenceModels.contains(tmodel)) confidenceModels.add(tmodel);
				cumWeight += tmodel.getAICcw();
			} else
				tmodel.setInAICcinterval(false);

			/*
			 * System.out.print("name=" + tmodel.name + " w=" + tmodel.AICcw +
			 * " cumw=" + tmodel.cumAICw); System.out.print(" in=" + probIn +
			 * " out=" + probOut); System.out.print(" r=" + randomNumber +
			 * " isIn=" + tmodel.isInAICcinterval + " ");
			 */
		}

		/*
		 * System.out.print(confidenceInterval + " confidence interval (" +
		 * numModels + " models) = ["); for (Enumeration
		 * e=confidenceModels.elements(); e.hasMoreElements();) { Model m =
		 * (Model)e.nextElement(); System.out.print(m.name + " "); }
		 * System.out.print("]");
		 */
	}

	public double getMinModelValue() {
		return minModel.getAICc();
	}

	public double getMinModelWeight() {
		return minModel.getAICcw();
	}
	
	@Override
	public double getValue(Model m) {
		return m.getAICc();
	}
	
	@Override
	public double getWeight(Model m) {
		return m.getAICcw();
	}
	
	@Override
	public double getDelta(Model m) {
		return m.getAICcd();
	}
	
	@Override
	public double getUDelta(Model m) {
		return m.getUAICcd();
	}
	
	@Override
	public double setUDelta(Model m) {
		m.setUAICcd(computeAicc(m.getLnLIgnoringGaps(), m.getK(), options)
				- computeAicc(unconstrainedModel.getLnL(), 
						unconstrainedModel.getK(), options));
		return m.getUAICcd();
	}
	
	@Override
	public double getCumWeight(Model m) {
		return m.getCumAICcw();
	}
	
	@Override
	public int getType() {
		return IC_AICc;
	}

} // class AICc public class TreeSummary {

	private int AIC_INDEX = 0;
	private int AICC_INDEX = 1;
	private int BIC_INDEX = 2;
	private int DT_INDEX = 3;
	// number of Information Criteria
	private int IC_COUNT = 4;

	private Tree bestTree[] = new Tree[IC_COUNT];
	private TreeDistancesCache rfDistances = TreeRFDistancesCache.getInstance();
	private TreeDistancesCache euclideanDistances = TreeEuclideanDistancesCache
			.getInstance();
	private Hashtable<Tree, SummaryRow> summary;
	// sorted topologies by criterion
	private List<Tree> topologiesAIC = null, topologiesAICc = null,
			topologiesBIC = null, topologiesDT = null, sortedTopologies = null;

	public TreeSummary(Tree bestAIC, Tree bestAICc, Tree bestBIC, Tree bestDT,
			Model[] models) {
		if (bestAIC != null) {
			this.bestTree[AIC_INDEX] = bestAIC;
			this.topologiesAIC = new ArrayList<Tree>();
		}
		if (bestAICc != null) {
			this.bestTree[AICC_INDEX] = bestAICc;
			this.topologiesAICc = new ArrayList<Tree>();
		}
		if (bestBIC != null) {
			this.bestTree[BIC_INDEX] = bestBIC;
			this.topologiesBIC = new ArrayList<Tree>();
		}
		if (bestDT != null) {
			this.bestTree[DT_INDEX] = bestDT;
			this.topologiesDT = new ArrayList<Tree>();
		}
		this.summary = new Hashtable<Tree, SummaryRow>();

		for (Model model : models) {
			boolean done = false;
			for (Tree tree : summary.keySet()) {
				if (sameTopology(tree, model.getTree())) {
					summary.get(tree).addModel(model);
					done = true;
					break;
				}
			}
			if (!done) {
				// new topology
				SummaryRow newRow = new SummaryRow(model.getTree());
				newRow.addModel(model);
				summary.put(model.getTree(), newRow);
				if (topologiesAIC != null)
					topologiesAIC.add(model.getTree());
				if (topologiesAICc != null)
					topologiesAICc.add(model.getTree());
				if (topologiesBIC != null)
					topologiesBIC.add(model.getTree());
				if (topologiesDT != null)
					topologiesDT.add(model.getTree());
			}
		}

		for (Tree tree : summary.keySet()) {
			summary.get(tree).computeEuclideanDistances();
		}
		// sort by criteria
		if (topologiesAIC != null)
			Collections.sort(topologiesAIC, new AicComparator());
		if (topologiesAICc != null)
			Collections.sort(topologiesAICc, new AiccComparator());
		if (topologiesBIC != null)
			Collections.sort(topologiesBIC, new BicComparator());
		if (topologiesDT != null)
			Collections.sort(topologiesDT, new DtComparator());
		
		sortedTopologies = topologiesAIC != null ? topologiesAIC
				: topologiesBIC != null ? topologiesBIC
						: topologiesAICc != null ? topologiesAICc
								: topologiesDT != null ? topologiesDT
									: new ArrayList<Tree>(summary.keySet());
	}

	public int getNumberOfTopologies() {
		return summary.size();
	}

	public Tree getTopology (int index) {
		return sortedTopologies.get(index);
	}
	
	public List<Model> getModelsByTopology (int index) {
		return summary.get(getTopology(index)).models;
	}
	
	public List<Model> getAICModels(int index) {
		Tree key = topologiesAIC.get(index);
		return summary.get(key).models;
	}

	public List<Model> getBICModels(int index) {
		Tree key = topologiesBIC.get(index);
		return summary.get(key).models;
	}

	public List<Model> getAICcModels(int index) {
		Tree key = topologiesAICc.get(index);
		return summary.get(key).models;
	}

	public List<Model> getDTModels(int index) {
		Tree key = topologiesDT.get(index);
		return summary.get(key).models;
	}

	public int aicIndexOf(Tree tree) {
		return topologiesAIC.indexOf(tree);
	}

	public long aicRfOf(Tree tree) {
		return summary.get(tree).rfDistance[AIC_INDEX];
	}
	
	public double aicAvgDistance(Tree tree) {
		return summary.get(tree).avgEuclideanDistance[AIC_INDEX];
	}
	
	public double aicVarDistance(Tree tree) {
		return summary.get(tree).varEuclideanDistance[AIC_INDEX];
	}
	
	public double aiccWeight(Tree tree) {
		return summary.get(tree).support[AICC_INDEX];
	}
	
	public int aiccIndexOf(Tree tree) {
		return topologiesAICc.indexOf(tree);
	}

	public long aiccRfOf(Tree tree) {
		return summary.get(tree).rfDistance[AICC_INDEX];
	}
	
	public double aiccAvgDistance(Tree tree) {
		return summary.get(tree).avgEuclideanDistance[AICC_INDEX];
	}
	
	public double aiccVarDistance(Tree tree) {
		return summary.get(tree).varEuclideanDistance[AICC_INDEX];
	}
	
	public double aicWeight(Tree tree) {
		return summary.get(tree).support[AIC_INDEX];
	}
	
	public int bicIndexOf(Tree tree) {
		return topologiesBIC.indexOf(tree);
	}

	public long bicRfOf(Tree tree) {
		return summary.get(tree).rfDistance[BIC_INDEX];
	}
	
	public double bicAvgDistance(Tree tree) {
		return summary.get(tree).avgEuclideanDistance[BIC_INDEX];
	}
	
	public double bicVarDistance(Tree tree) {
		return summary.get(tree).varEuclideanDistance[BIC_INDEX];
	}
	
	public double bicWeight(Tree tree) {
		return summary.get(tree).support[BIC_INDEX];
	}
	
	public int dtIndexOf(Tree tree) {
		return topologiesDT.indexOf(tree);
	}

	public long dtRfOf(Tree tree) {
		return summary.get(tree).rfDistance[DT_INDEX];
	}
	
	public double dtAvgDistance(Tree tree) {
		return summary.get(tree).avgEuclideanDistance[DT_INDEX];
	}
	
	public double dtVarDistance(Tree tree) {
		return summary.get(tree).varEuclideanDistance[DT_INDEX];
	}
	
	public double dtWeight(Tree tree) {
		return summary.get(tree).support[DT_INDEX];
	}
	
	private boolean sameTopology(Tree t1, Tree t2) {
		return (rfDistances.getDistance(t1, t2) == 0);
	}

	private long rfDistance(Tree t1, Tree t2) {
		if (t1 == null || t2 == null) {
			return -1;
		} else {
			return Math.round(rfDistances.getDistance(t1, t2));
		}
	}

	class SummaryRow {
		List<Model> models;
		Tree commonTopology;
		long rfDistance[] = new long[] { 0l, 0l, 0l, 0l };
		double avgEuclideanDistance[] = new double[] { 0.0d, 0.0d, 0.0d, 0.0d };
		double varEuclideanDistance[] = new double[] { 0.0d, 0.0d, 0.0d, 0.0d };
		double support[] = new double[] { 0.0d, 0.0d, 0.0d, 0.0d };

		SummaryRow(Tree commonTopology) {
			this.commonTopology = commonTopology;
			for (int i = 0; i < IC_COUNT; i++) {
				this.rfDistance[i] = rfDistance(commonTopology, bestTree[i]);
			}
			this.models = new ArrayList<Model>();
		}

		void addModel(Model model) {
			// check same topology
			if (checkTopology(model)) {
				models.add(model);
				support[AIC_INDEX] += model.getAICw();
				support[AICC_INDEX] += model.getAICcw();
				support[BIC_INDEX] += model.getBICw();
				support[DT_INDEX] += model.getDTw();
			}
		}

		void computeEuclideanDistances() {
			// set variables to zero
			for (int i = 0; i < IC_COUNT; i++) {
				if (bestTree[i] != null) {
					avgEuclideanDistance[i] = 0.0d;
					varEuclideanDistance[i] = 0.0d;
					for (Model model : models) {
						double distance = euclideanDistances.getDistance(
								model.getTree(), bestTree[i]);
						avgEuclideanDistance[i] += distance;
						varEuclideanDistance[i] += distance * distance;
					}
					avgEuclideanDistance[i] /= models.size();
					varEuclideanDistance[i] /= models.size();
					varEuclideanDistance[i] -= avgEuclideanDistance[i]
							* avgEuclideanDistance[i];
				}
			}
		}

		boolean checkTopology(Model model) {
			return (rfDistances.getDistance(commonTopology, model.getTree()) == 0);
		}
	}

	class AicComparator implements Comparator<Tree> {

		@Override
		public int compare(Tree o1, Tree o2) {
			long d1 = rfDistance(o1, bestTree[AIC_INDEX]);
			long d2 = rfDistance(o2, bestTree[AIC_INDEX]);
			return (int) (d1 - d2);
		}

	}

	class AiccComparator implements Comparator<Tree> {

		@Override
		public int compare(Tree o1, Tree o2) {
			long d1 = rfDistance(o1, bestTree[AICC_INDEX]);
			long d2 = rfDistance(o2, bestTree[AICC_INDEX]);
			return (int) (d1 - d2);
		}

	}

	class BicComparator implements Comparator<Tree> {

		@Override
		public int compare(Tree o1, Tree o2) {
			long d1 = rfDistance(o1, bestTree[BIC_INDEX]);
			long d2 = rfDistance(o2, bestTree[BIC_INDEX]);
			return (int) (d1 - d2);
		}

	}

	class DtComparator implements Comparator<Tree> {

		@Override
		public int compare(Tree o1, Tree o2) {
			long d1 = rfDistance(o1, bestTree[DT_INDEX]);
			long d2 = rfDistance(o2, bestTree[DT_INDEX]);
			return (int) (d1 - d2);
		}

	}

	public void print(TextOutputStream stream) {
		int MAX_LEN = 80;
		if (topologiesAIC != null || topologiesBIC != null
				|| topologiesAICc != null || topologiesDT != null) {
			List<Tree> sortedTopologies = topologiesAIC != null ? topologiesAIC
					: topologiesBIC != null ? topologiesBIC
							: topologiesAICc != null ? topologiesAICc
									: topologiesDT;
			stream.println("::Optimized Topologies Summary::");
			stream.println("");
			stream.println("There are " + getNumberOfTopologies()
					+ " different topologies.");
			stream.println("");
			for (int i = 0; i < getNumberOfTopologies(); i++) {
				int index = i + 1;
				Tree currentTree = sortedTopologies.get(i);

				SummaryRow summaryRow = summary.get(currentTree);

				int aicIndex = topologiesAIC != null ? aicIndexOf(currentTree) + 1
						: -1;
				int aiccIndex = topologiesAICc != null ? aiccIndexOf(currentTree) + 1
						: -1;
				int bicIndex = topologiesBIC != null ? bicIndexOf(currentTree) + 1
						: -1;
				int dtIndex = topologiesDT != null ? dtIndexOf(currentTree) + 1
						: -1;
				stream.println("Topology Id: " + index);
				stream.println("\tRank\tWeight\t\t RF\tAvgEucl\t\tVarEucl");
				if (topologiesAIC != null) {
						stream.println(getIcRow("AIC", aicIndex, summaryRow, AIC_INDEX));
				}
				if (topologiesBIC != null) {
					stream.println(getIcRow("BIC", bicIndex, summaryRow, BIC_INDEX));
				}
				if (topologiesAICc != null) {
					stream.println(getIcRow("AICc", aiccIndex, summaryRow, AICC_INDEX));
				}
				if (topologiesDT != null) {
					stream.println(getIcRow("DT", dtIndex, summaryRow, DT_INDEX));
				}
				stream.println("Models supporting:   "
						+ summaryRow.models.size());

				stream.print("                     ");
				int chars = 0;
				for (Model model : summaryRow.models) {
					if (chars >= MAX_LEN) {
						stream.println("");
						stream.print("                     ");
						chars = 0;
					}
					stream.print(model.getName() + " ");
					chars += model.getName().length() + 1;
				}
				stream.println("");
				stream.println("");
			}
		}
	}
	
	private String getIcRow(String name, int rankIndex, SummaryRow summaryRow, int criterionIndex) {
		return name + "\t"
				+ Utilities.format(rankIndex,3,0,false)
				+ "\t"
				+ Utilities
						.asPercent(summaryRow.support[criterionIndex] * 100)
				+ "\t\t" + Utilities.format(summaryRow.rfDistance[criterionIndex],3,0,false)
				+ "\t" + Utilities.format(summaryRow.avgEuclideanDistance[criterionIndex],8,2,true)
				+ "\t" + Utilities.format(summaryRow.varEuclideanDistance[criterionIndex],8,2,true);
	}
} public class TreeUtilities {

    public static final int DEFAULT_COLUMN_WIDTH = 70;
    public static final String TREE_WEIGHT_ATTRIBUTE = "weight";
    public static final String TREE_CLADE_SUPPORT_ATTRIBUTE = "support";
    public static final String TREE_NAME_ATTRIBUTE = "treeName";

	public TreeUtilities() {
	}

	/****************************
	 * readTree ********************************** * Reads a tree in Newick
	 * format, rooted or unroored, binary or with * polytomies, and with or
	 * without branch lengths. Returns the root * node *
	 * 
	 * @throws IOException
	 * @throws TreeParseException
	 *             *
	 ************************************************************************/

	public static Tree readTree(String treefilename) throws IOException,
			TreeParseException {
		Tree tree;
		try {
			tree = new ReadTree(treefilename);
		} catch (TreeParseException e) {
			throw e;
		} catch (IOException e) {
			throw e;
		}
		return tree;
	}

	/**
     * Make sure subtree below node has consistent heights, i.e. node height is higher than it's descendants
     * 
     * @param tree the tree
     * @param node the node
     * 
     * @return height of node
     */
    public static double insureConsistency(Tree tree, Node node) {
        double height = TreeUtilities.safeNodeHeight(tree, node);
        if (node.isLeaf()) {
            return height;
        } else {
            for (int i = 0; i < node.getChildCount(); i++) {
                Node n = node.getChild(i);
                final double childHeight = insureConsistency(tree, n);
                height = Math.max(height, childHeight);
            }
        }

        node.setNodeHeight(height);
        return height;
    }
    
    /**
     * Calculates the number of branches from node to most remote tip.
     * 
     * @param node the starting node
     * 
     * @return the node distance
     */
    public static int nodeDistance(final Node node) {
        if (node.isLeaf()) {
            return 0;
        }

        int d = 0;
        for (int i = 0; i < node.getChildCount(); i++) {
            Node n = node.getChild(i);
            d = Math.max(d, nodeDistance(n));
        }
        return d + 1;
    }
    /**
     * Calculates the safe node height.
     * 
     * @param tree the tree
     * @param node the node
     * 
     * @return the height of the node
     */
    public static double safeNodeHeight(final Tree tree, final Node node) {
        if (node.getNodeHeight() > 0.0) {
            return node.getNodeHeight();
        }
        return TreeUtilities.nodeDistance(node);
    }
    
	private static void putCharAtLevel(PrintWriter out, int level, char c,
            int[] position) {
        int n = position[level] - 1;
        for (int i = 0; i < n; i++) {
            out.print(' ');
        }
        out.print(c);
    }
	
	private static void printlnNodeWithNumberAndLabel(PrintWriter out, Node node, int level,
            int numExternalNodes, boolean[] umbrella, int[] position) {
        for (int i = 0; i < level - 1; i++) {
            if (umbrella[i]) {
                putCharAtLevel(out, i, '|', position);
            } else {
                putCharAtLevel(out, i, ' ', position);
            }
        }

        putCharAtLevel(out, level - 1, '+', position);

        int branchNumber;
        if (node.isLeaf()) {
            branchNumber = node.getNumber() + 1;
        } else {
            branchNumber = node.getNumber() + 1 + numExternalNodes;
        }

        String numberAsString = Integer.toString(branchNumber);

        int numDashs = position[level] - numberAsString.length();

        for (int i = 0; i < numDashs ; i++) {
            out.print('-');
        }
        out.print(numberAsString);

        if (node.isLeaf()) {
            out.println(" " + node.getIdentifier());
        } else {
            if (!node.getIdentifier().equals(Identifier.ANONYMOUS)) {
                out.print("(" + node.getIdentifier() + ")");
            }
            out.println();
        }
    }
	private static void printNodeInASCII(PrintWriter out, Node node, int level, int m, int maxm,
            int numExternalNodes, boolean[] umbrella, int[] position, double proportion, int minLength) {
        position[level] = (int) (node.getBranchLength() * proportion);

        if (position[level] < minLength) {
            position[level] = minLength;
        }

        if (node.isLeaf()) // external branch
        {
            if (m == maxm - 1) {
                umbrella[level - 1] = true;
            }

            printlnNodeWithNumberAndLabel(out, node, level, numExternalNodes, umbrella, position);

            if (m == 0) {
                umbrella[level - 1] = false;
            }
        } else // internal branch
        {
            for (int n = node.getChildCount() - 1; n > -1; n--) {
                printNodeInASCII(out, node.getChild(n), level + 1, n, node.getChildCount(),
                        numExternalNodes, umbrella, position, proportion, minLength);

                if (m == maxm - 1 && n == node.getChildCount() / 2) {
                    umbrella[level - 1] = true;
                }

                if (n != 0) {
                    if (n == node.getChildCount() / 2) {
                        printlnNodeWithNumberAndLabel(out, node, level, numExternalNodes, umbrella, position);
                    } else {
                        for (int i = 0; i < level + 1; i++) {
                            if (umbrella[i]) {
                                putCharAtLevel(out, i, '|', position);
                            } else {
                                putCharAtLevel(out, i, ' ', position);
                            }
                        }
                        out.println();
                    }
                }

                if (m == 0 && n == node.getChildCount() / 2) {
                    umbrella[level - 1] = false;
                }
            }
        }
    }
	
	// Print picture of current tree in ASCII
    public static void printASCII(Tree tree, PrintWriter out) {
        tree.createNodeList();

        int numExternalNodes = tree.getExternalNodeCount();
        int numInternalNodes = tree.getInternalNodeCount();
        int numBranches = numInternalNodes + numExternalNodes - 1;

        boolean[] umbrella = new boolean[numExternalNodes];
        int[] position = new int[numExternalNodes];

        int minLength = (Integer.toString(numBranches)).length() + 1;

        int MAXCOLUMN = 40;
        Node root = tree.getRoot();
        if (root.getNodeHeight() == 0.0) {
            NodeUtils.lengths2Heights(root);
        }
        double proportion = (double) MAXCOLUMN / root.getNodeHeight();

        for (int n = 0; n < numExternalNodes; n++) {
            umbrella[n] = false;
        }

        position[0] = 1;
        for (int i = root.getChildCount() - 1; i > -1; i--) {
            printNodeInASCII(out, root.getChild(i), 1, i, root.getChildCount(),
                    numExternalNodes, umbrella, position, proportion, minLength);
            if (i != 0) {
                putCharAtLevel(out, 0, '|', position);
                out.println();
            }
        }
    }
	public static void printNH(PrintWriter out, Tree tree, Node node,
			boolean printLengths, boolean printInternalLabels,
			boolean printCladeSupport) {

		if (!node.isLeaf()) {
			out.print("(");

			for (int i = 0; i < node.getChildCount(); i++) {
				if (i != 0) {
					out.print(",");
				}

				printNH(out, tree, node.getChild(i), printLengths,
						printInternalLabels, printCladeSupport);
			}

			out.print(")");
		}

		if (!node.isRoot()) {
			if (node.isLeaf() || printInternalLabels) {

				String id = node.getIdentifier().toString();
				out.print(id);
			}

			if (printCladeSupport) {
				if (tree.getAttribute(node, TREE_CLADE_SUPPORT_ATTRIBUTE) != null) {
					double support = (Double) tree.getAttribute(node,
							TREE_CLADE_SUPPORT_ATTRIBUTE);
					out.printf(":" + FormattedOutput.getInstance().getDecimalString(support, 4));
				}
			}

			if (printLengths) {
				out.printf(":" + FormattedOutput.getInstance().getDecimalString(node.getBranchLength(), 10));
			}
		}

	}

	public static String toNewick(Tree tree, boolean printLengths,
			boolean printInternalLabels, boolean printCladeSupport) {

		StringWriter sw = new StringWriter();
		PrintWriter mp = new PrintWriter(sw);

		printNH(mp, tree, tree.getRoot(), printLengths, printInternalLabels,
				printCladeSupport);

		sw.append(';');
		return sw.toString();
	}

	public static double getEuclideanTreeDistance(Tree t1, Tree t2) {
		double sum = 0.0;
        int numberOfInternalNodes = t1.getInternalNodeCount();
        if (numberOfInternalNodes != t2.getInternalNodeCount()) {
            throw new RuntimeException("Different number of internal nodes: " +
                    t1.getInternalNodeCount() + " vs " + t2.getInternalNodeCount());
        }
        int numberOfExternalNodes = t1.getExternalNodeCount();
        if (numberOfExternalNodes != t2.getExternalNodeCount()) {
            throw new RuntimeException("Different number of external nodes: " +
                    t1.getInternalNodeCount() + " vs " + t2.getInternalNodeCount());
        }
        for (int i = 0; i < numberOfInternalNodes; i++) {
            double bl1 = t1.getInternalNode(i).getBranchLength();
            double bl2 = t2.getInternalNode(i).getBranchLength();
            sum += (bl1 - bl2) * (bl1 - bl2);
        }
        for (int i = 0; i < numberOfExternalNodes; i++) {
            double bl1 = t1.getExternalNode(i).getBranchLength();
            double bl2 = t2.getExternalNode(i).getBranchLength();
            sum += (bl1 - bl2) * (bl1 - bl2);
        }
        return Math.sqrt(sum);
	}
	
	public static double getRobinsonFouldsTreeDistance(Tree t1, Tree t2) {
		SplitSystem s1 = SplitUtils.getSplits(t1);
		IdGroup idGroup = s1.getIdGroup();
		SplitSystem s2 = SplitUtils.getSplits(idGroup, t2);

		if (s1.getLabelCount() != s2.getLabelCount())
			throw new IllegalArgumentException("Number of labels must be the same!");

		int ns1 = s1.getSplitCount();
		int ns2 = s2.getSplitCount();

		// number of splits in t1 missing in t2
		int fn = 0;
		for (int i = 0; i < ns1; i++)
		{
			if (!s2.hasSplit(s1.getSplit(i))) fn++;
		}

		// number of splits in t2 missing in t1
		int fp = 0;
		for (int i = 0; i < ns2; i++)
		{
			if (!s1.hasSplit(s2.getSplit(i))) fp++;
		}


		return ((double) fp + (double) fn);
	}
} // end of class public class Simulation {
	private ApplicationOptions options;
	
	private static final double nil = -99999;

	private static final String AIC = "AIC";
	private static final String AICc = "AICc";
	private static final String BIC = "BIC";
	private static final String DT = "DT";
	private static final String[] IC_TYPES = { AIC, AICc, BIC, DT };

	public Simulation(ApplicationOptions options) {
		this.options = options;
	}

	/*************************
	 * runSimulations ****************************** * Organizes all the tasks
	 * that the program needs to carry out * * *
	 ***********************************************************************/
	public void run() {
		int i;
		boolean append;
		/*
		 * .trees: base trees for every model in standard orden
		 * 
		 * .AIC_tre: best tree for the AIC model .AIC_mmi_mtre: model averaged
		 * AIC tree (50% majority rule) .AIC_mmi_atre: model averaged AIC tree
		 * (strict) ... same for AICc, BIC and DT
		 * 
		 * .parms: parameter estimates for every model in standard orden
		 * .AIC_parms: parameter estimates for the AIC model .AIC_imp_parms:
		 * parameter importances obtained with the AIC .AIC_mmi_parms: parameter
		 * importances obtained with the AIC ... same for AICc, BIC and DT
		 */

		// make simulation directories
		String simdir = "sims/";
		new File(simdir).mkdir();

		String outdir = simdir + "out/";
		new File(outdir).mkdir();

		String treedir = simdir + "trees/";
		new File(treedir).mkdir();

		String parmdir = simdir + "parameters/";
		new File(parmdir).mkdir();

		String simsName = options.simulationsName;

		if (options.getInputFile().getName().endsWith("001")
				|| options.getInputFile().getName().endsWith("001.phy")
				|| options.getInputFile().getName().endsWith("001.dat"))
			append = false;
		else
			append = true;

		// outfiles
		TextOutputStream outConsole = ModelTest.setMainConsole(new TextOutputStream(
				outdir + simsName + ".out", append));
		TextOutputStream treeConsole = new TextOutputStream(treedir + simsName
				+ ".trees", append);
		TextOutputStream parameterConsole = new TextOutputStream(parmdir
				+ simsName + ".parms", append);

		// print header information
		ModelTest.printHeader(outConsole);

		// print notice information
		// printNotice(outConsole);

		// check expiration date
		// CheckExpiration (outConsole);

		// print the command line
		outConsole.print("\nArguments =");
		for (i = 0; i < ModelTest.arguments.length; i++)
			outConsole.print(" " + ModelTest.arguments[i]);

		// open data file
		File file = options.getInputFile();
		outConsole.print("\n\nReading data file \"" + file.getName()
				+ "\"...");

		if (file.exists()) {
			try {

//				File outputFile = File
//						.createTempFile("jmodeltest", "input.aln");
				String alnStr = ModelTestService.readAlignment(file, options.getAlignmentFile());

				PushbackReader pr = new PushbackReader(
						new StringReader(alnStr));
				options.setAlignment(AlignmentReader.createAlignment(
						new PrintWriter(System.err),
						pr, 
						true)); // file

				outConsole.println(" OK.");
				outConsole.println("  number of sequences: "
						+ options.getNumTaxa());
				outConsole.println("  number of sites: " + options.getNumSites());
			} catch (Exception e)// file cannot be read correctly
			{
				System.err.println("\nThe specified file \""
						+ file.getName()
						+ "\" cannot be read as an alignment");
				outConsole.println(" failed.\n");
				System.exit(0);
			}
		} else // file does not exist
		{
			System.err.println("\nThe specified file \""
					+ file.getName() + "\" cannot be found");
			outConsole.println(" failed.\n");
			System.exit(0);
		}

		// open tree file if necessary
		if (options.userTopologyExists) {
			File treefile = options.getTreeFile();
			options.setInputTreeFile(file);
			outConsole.print("Reading tree file \"" + treefile.getName()
					+ "\"...");

			// read the tree in
			Tree tree = null;
			try {
				tree = TreeUtilities.readTree(treefile.getAbsolutePath());
			} catch (IOException e) {
				System.err.println("\nThe specified tree file \""
						+ treefile.getName() + "\" cannot be found");
				outConsole.println(" failed.\n");
				System.exit(0);
			} catch (TreeParseException e) {
				System.err.println("\nThe specified file \""
						+ treefile.getName()
						+ "\" cannot be read as valid Newick tree");
				outConsole.println(" failed.\n");
				System.exit(0);
			}
			if (tree != null) {
				options.setUserTree(TreeUtilities.toNewick(tree, true,
						false, false));
				TextOutputStream out = new TextOutputStream(
						options.getTreeFile().getAbsolutePath());
				out.print(TreeUtilities.toNewick(tree, true, false, false));
				out.close();
				outConsole.println(" OK.");
			}

		}

		// print some progress
		System.err.print("Doing \"" + file.getName() + "\" ... ");

		// calculate number of models
		if (options.getSubstTypeCode() == 0)
			options.setNumModels(3);
		else if (options.getSubstTypeCode() == 1)
			options.setNumModels(5);
		else if (options.getSubstTypeCode() == 2)
			options.setNumModels(7);
		else
			options.setNumModels(11);

		if (options.doF)
			options.setNumModels(options.getNumModels() * 2);

		if (options.doI && options.doG)
			options.setNumModels(options.getNumModels() * 4);
		else if (options.doI || options.doG)
			options.setNumModels(options.getNumModels() * 2);

		// build set of models
		options.setCandidateModels();

		// calculate likelihoods with phyml in the command line
		RunPhyml phymlrun = new RunPhymlThread(new ConsoleProgressObserver(options), options, ModelTest.getCandidateModels());
		phymlrun.execute();

		// print all model trees to tree file and parameters to parms file
		if (!append)
			parameterConsole
					.println("data\tname\tln\tK\tfA\tfC\tfG\tfT\tkappa\ttitv\trAC\tAG\trAT\trCG\trCT\trGT\tpinvI\tshapeG\tpinvIG\tshapeIG");
		for (i = 0; i < options.getNumModels(); i++) {
			treeConsole.println(options.getInputFile().getName() + "\t"
					+ ModelTest.getCandidateModels()[i].getName() + "\t"
					+ ModelTest.getCandidateModels()[i].getTreeString());
			printModelLine(
					ModelTest.getCandidateModels()[i],
					parameterConsole,
					options.getInputFile().getName() + "\t"
							+ ModelTest.getCandidateModels()[i].getName(), nil, nil);
		}

		// do AIC if selected
		if (options.doAIC) {
			// TODO: oquhuh
			printSelection("AIC", treedir, simsName, parmdir, outConsole,
					append);
		}

		// do AICc if selected
		if (options.doAICc) {
			printSelection("AICc", treedir, simsName, parmdir, outConsole,
					append);
		}

		// do BIC if selected
		if (options.doBIC) {
			printSelection("BIC", treedir, simsName, parmdir, outConsole,
					append);
		}

		// do DT if selected
		if (options.doDT) {
			printSelection("DT", treedir, simsName, parmdir, outConsole, append);
		}

		// do hLRT if selected
		if (options.doHLRT) {
			HLRT myHLRT = new HLRT(options);
			myHLRT.compute(!options.backwardHLRTSelection,
					options.confidenceLevelHLRT, options.writePAUPblock);
			treeConsole.println(ModelTest.averagedTreeString);
		}

		// do dLRT if selected
		if (options.doDLRT) {
			HLRT myHLRT = new HLRT(options);
			myHLRT.computeDynamical(!options.backwardHLRTSelection,
					options.confidenceLevelHLRT, options.writePAUPblock);
		}

		outConsole.println("\n=> This run has finished.\n");
		
		treeConsole.close();
		
		System.err.println("OK.");
	}

	/****************************
	 * printModelLine *************************** * Print model components in a
	 * tabulated format * *
	 ************************************************************************/

	public void printModelLine(Model model, TextOutputStream stream,
			String whichName, double score, double weight) {

		if (whichName == null)
			stream.printf("%-10s", model.getName());
		else
			stream.printf("%-10s", whichName);

		if (model.getLnL() == 0) {
			stream.println("\tOPTIMIZATION FAILED!");
			System.exit(0);
		} else {
			stream.printf("\t%.4f", model.getLnL());
			stream.printf("\t%d", model.getK());

			if (score != nil)
				stream.printf("\t%.4f ", score);

			if (weight != nil)
				stream.printf("\t%.4f ", weight);

			if (model.ispF()) {
				stream.printf("\t%.4f ", model.getfA());
				stream.printf("\t%.4f ", model.getfC());
				stream.printf("\t%.4f ", model.getfG());
				stream.printf("\t%.4f ", model.getfT());
			} else
				stream.print("\t-\t-\t-\t-");

			if (model.ispT()) {
				stream.printf("\t%.4f", model.getKappa());
				stream.printf("\t%.4f", model.getTitv());
			} else
				stream.print("\t-\t-");

			if (model.ispR()) {
				stream.printf("\t%.4f", model.getRa());
				stream.printf("\t%.4f", model.getRb());
				stream.printf("\t%.4f", model.getRc());
				stream.printf("\t%.4f", model.getRd());
				stream.printf("\t%.4f", model.getRe());
				stream.printf("\t%.4f", 1.0);
			} else
				stream.print("\t-\t-\t-\t-\t-\t-");

			if (model.ispI() && model.ispG()) {
				stream.printf("\t-\t-\t%.4f", model.getPinv());
				stream.printf("\t%.4f", model.getShape());
			} else if (model.ispI())
				stream.printf("\t%.4f\t-\t-\t-", model.getPinv());
			else if (model.ispG())
				stream.printf("\t-\t%.4f\t-\t-", model.getShape());
			else
				stream.print("\t-\t-\t-\t-");

			stream.println(" ");
		}
	}

	/*************************
	 * CheckNAsim ********************************* * Returns formatted value or
	 * NA symbol depending on variable status * * *
	 ***********************************************************************/

	public static String CheckNAsim(double value) {
		if (value == Utilities.NA)
			return "-";
		else {
			String s = String.format(Locale.ENGLISH, "%.4f", value);
			return s;
		}
	}

	public void printSelection(String icType, String treedir, String simsName,
			String parmdir, TextOutputStream outConsole, boolean append) {

		// Validate
		if (!Arrays.asList(IC_TYPES).contains(icType)) {
			// TODO: EXCEPTION!
		}

		TextOutputStream ICtre = new TextOutputStream(treedir + simsName + "."
				+ icType + "_tre", append);
		TextOutputStream IC_mmi_mtre = new TextOutputStream(treedir + simsName
				+ "." + icType + "_mmi_mtre", append);
		TextOutputStream IC_mmi_atre = new TextOutputStream(treedir + simsName
				+ "." + icType + "_mmi_atre", append);

		TextOutputStream IC_parms = new TextOutputStream(parmdir + simsName
				+ "." + icType + "_parms", append);
		TextOutputStream IC_imp_parms = new TextOutputStream(parmdir + simsName
				+ "." + icType + "_imp_parms", append);
		TextOutputStream IC_mmi_parms = new TextOutputStream(parmdir + simsName
				+ "." + icType + "_mmi_parms", append);

		InformationCriterion ic;
		if (icType.equals(AIC)) {
			ic = new AIC(options.writePAUPblock, options.doImportances,
					options.doModelAveraging, options.confidenceInterval);
		} else if (icType.equals("AICc")) {
			ic = new AICc(options.writePAUPblock, options.doImportances,
					options.doModelAveraging, options.confidenceInterval);
		} else if (icType.equals("BIC")) {
			ic = new BIC(options.writePAUPblock, options.doImportances,
					options.doModelAveraging, options.confidenceInterval);
		} else if (icType.equals("DT")) {
			ic = new DT(options.writePAUPblock, options.doImportances,
					options.doModelAveraging, options.confidenceInterval);
		} else {
			ic = null;
			// TODO: EXCEPTION!!!
		}

		ic.compute();
		ic.print(outConsole);
		if (icType.equals(AIC)) {
			ModelTest.setMyAIC((AIC) ic);
		} else if (icType.equals("AICc")) {
			ModelTest.setMyAICc((AICc) ic);
		} else if (icType.equals("BIC")) {
			ModelTest.setMyBIC((BIC) ic);
		} else if (icType.equals("DT")) {
			ModelTest.setMyDT((DT) ic);
		}

		ICtre.println(ic.getMinModel().getTreeString());
		if (!append)
			IC_parms.print("data\tname\tln\tK\tscore\tweigth\tfA\tfC\tfG\tfT\tkappa\ttitv\trAC\tAG\trAT\trCG\trCT\trGT\tpinvI\tshapeG\tpinvIG\tshapeIG\n");
		printModelLine(ic.getMinModel(), IC_parms, options.getInputFile().getName()
				+ "\t" + ic.getMinModel().getName(), ic.getMinModelValue(),
				ic.getMinModelWeight());

		if (options.doAveragedPhylogeny) {
			new RunConsense(ic,	"50% majority rule", 
					options.confidenceInterval);
			IC_mmi_mtre.println(ModelTest.averagedTreeString);
			new RunConsense(ic, "strict",
					options.confidenceInterval);
			IC_mmi_atre.println(ModelTest.averagedTreeString);
		}
		if (options.doImportances) {
			if (!append) {
				IC_imp_parms
						.print("data\tfA\tfC\tfG\tfT\tkappa\ttitv\trAC\tAG\trAT\trCG\trCT\trGT\tpinvI\tshapeG\tpinvIG\tshapeIG\n");
			}
			IC_imp_parms.print(options.getInputFile().getName());
			IC_imp_parms.printf("\t%.4f", ic.getIfA());
			IC_imp_parms.printf("\t%.4f", ic.getIfC());
			IC_imp_parms.printf("\t%.4f", ic.getIfG());
			IC_imp_parms.printf("\t%.4f", ic.getIfT());
			IC_imp_parms.printf("\t%.4f", ic.getIkappa());
			IC_imp_parms.printf("\t%.4f", ic.getItitv());
			IC_imp_parms.printf("\t%.4f", ic.getiRa());
			IC_imp_parms.printf("\t%.4f", ic.getiRb());
			IC_imp_parms.printf("\t%.4f", ic.getiRc());
			IC_imp_parms.printf("\t%.4f", ic.getiRd());
			IC_imp_parms.printf("\t%.4f", ic.getiRe());
			IC_imp_parms.printf("\t%.4f", ic.getiRf());
			IC_imp_parms.printf("\t%.4f", ic.getIpinvI());
			IC_imp_parms.printf("\t%.4f", ic.getIshapeG());
			IC_imp_parms.printf("\t%.4f", ic.getIpinvIG());
			IC_imp_parms.printf("\t%.4f", ic.getIshapeIG());
			IC_imp_parms.println(" ");
		}
		if (options.doModelAveraging) {
			if (!append) {
				IC_mmi_parms
						.print("data\tfA\tfC\tfG\tfT\tkappa\ttitv\trAC\tAG\trAT\trCG\trCT\trGT\tpinvI\tshapeG\tpinvIG\tshapeIG\n");
			}
			IC_mmi_parms.print(options.getInputFile().getName());
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAfA()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAfC()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAfG()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAfT()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAkappa()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAtitv()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getaRa()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getaRb()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getaRc()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getaRd()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getaRe()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getaRf()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getApinvI()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAshapeG()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getApinvIG()));
			IC_mmi_parms.printf("\t%s", CheckNAsim(ic.getAshapeIG()));
			IC_mmi_parms.println(" ");
		}
		
		ICtre.close();
		IC_mmi_mtre.close();
		IC_mmi_atre.close();
		IC_parms.close();
		IC_imp_parms.close();
		IC_mmi_parms.close();
	}
} public final class Utilities {

	public static final int NA = Integer.MIN_VALUE;
	public static final int OS_LINUX = 1;
	public static final int OS_OSX = 2;
	public static final int OS_WINDOWS = 3;

	public Utilities() {
	}

	public static String firstNumericToken(String str) {
		StringTokenizer st = new StringTokenizer(str);
		String token = "";
		boolean found = false;
		while (st.hasMoreTokens() && !found) {
			token = st.nextToken();
			found |= isNumber(token);
		}
		if (found)
			return token;
		else
			return null;
	}
	
	public static String lastToken(String str) {
		StringTokenizer st = new StringTokenizer(str);
		String token = "";
		while (st.hasMoreTokens()) {
			token = st.nextToken();
		}
		return token;
	}

	public static boolean isWindows() {
		if (System.getProperty("os.name").startsWith("Window"))
			return true;
		return false;
	}

	public static int findCurrentOS() {
		String os = System.getProperty("os.name");
		if (os.startsWith("Mac"))
			return OS_OSX;
		else if (os.startsWith("Windows"))
			return OS_WINDOWS;
		else
			return OS_LINUX;
	}

	public static String getBinaryVersion() {
		String arch = System.getProperty("os.arch");
		String bit = System.getProperty("sun.arch.data.model");

		String binaryName = null;
		switch (findCurrentOS()) {

		case OS_OSX:
			if (arch.startsWith("ppc")) {
				System.err
						.println("Sorry, PowerPC architecture is no longer supported");
				System.exit(0);
			} else {
				binaryName = "PhyML_3.0_macOS_i386";
			}
			break;
		case OS_LINUX:
			if (bit.startsWith("64"))
				binaryName = "PhyML_3.0_linux64";
			else
				binaryName = "PhyML_3.0_linux32";
			break;
		case OS_WINDOWS:
			binaryName = "PhyML_3.0_win32.exe";
			break;
		default:
			binaryName = null;
		}

		return binaryName;
	}

	public static String calculateRuntimeMinutes(long startTime, long endTime) {
		long seconds = Math.round((endTime - startTime) / 1000.0);
		int hours = (int) (seconds / 3600.0);
		int rest1 = (int) (seconds % (3600.0));
		int minutes = (int) (rest1 / 60.0);
		seconds = (int) (seconds - (hours * 3600 + minutes * 60));
		String h = "" + hours;
		String m = "" + minutes;
		String s = "" + seconds;
		if (minutes < 10)
			m = "0" + m;
		if (seconds < 10)
			s = "0" + s;
		return h + "h:" + m + ":" + s + "";
	}

	public static String getCurrentTime(String format) {
		Calendar cal = new GregorianCalendar();
	    SimpleDateFormat date_format = new SimpleDateFormat(format);
	    return date_format.format(cal.getTime());
	}
	
	public static String displayRuntime(long time) {
		long decimes = Math.round(time / 100.0);
		int hours = (int) (decimes / 36000.0);
		int rest1 = (int) (decimes % 36000.0);
		int minutes = (int) (rest1 / 600.0);
		int rest2 = (int) (rest1 % 600.0);
		int seconds = (int) (rest2 / 10.0);
		decimes = (int) (decimes - (hours * 36000 + minutes * 600 + seconds * 10));
		String h = "" + hours;
		String m = "" + minutes;
		String s = "" + seconds;
		String d = "0" + decimes;
		if (hours < 10)
			h = "0" + h;
		if (minutes < 10)
			m = "0" + m;
		if (seconds < 10)
			s = "0" + s;
		return h + "h:" + m + ":" + s + ":" + d + "";
	}

	public static String calculateRuntime(long startTime, long endTime) {
		return displayRuntime(endTime - startTime);
	}

	public static boolean isNumber(String s) {
		try {
			Double.parseDouble(s);
			return true;
		} catch (NumberFormatException ex) {
			return false;
		}
	}

	public static void toConsoleEnd() {
		if (ModelTest.buildGUI)
			XManager.getInstance()
					.getPane()
					.setCaretPosition(
							XManager.getInstance().getPane().getDocument()
									.getLength());
	}

	/**
	 * copyFile
	 * 
	 * Copies src file to dst file. If the dst file does not exist, it is
	 * created modified from The Java Developers Almanac 1.4:
	 * http://www.exampledepot.com/egs/java.io/CopyFile.html
	 */

	public static void copyFile(String source, String destination)
			throws IOException {
		File src = new File(source);
		File dst = new File(destination);

		InputStream in = new FileInputStream(src);
		OutputStream out = new FileOutputStream(dst);

		// Transfer bytes from in to out
		byte[] buf = new byte[1024];
		int len;
		while ((len = in.read(buf)) > 0) {
			out.write(buf, 0, len);
		}
		in.close();
		out.close();
	}

	/**
	 * CheckNA
	 * 
	 * If value is NA it prints "-"
	 */
	public static String checkNA(double value) {
		if (value == NA)
			return "  -  ";
		else {
			String s = String.format(Locale.ENGLISH, "%8.4f", value);
			return s;
		}
	}

	/**
	 * deleteFile
	 * 
	 * Delete a given file
	 */
	public static void deleteFile(String fileName) {

		File f;

		try {
			f = new File(fileName);
		} catch (Exception e) {
			throw new IllegalArgumentException("DeleteFiles: deletion of "
					+ fileName + " failed");
		}

		// Make sure the file or directory exists and isn't write protected
		if (f.exists()) {
			if (!f.canWrite())
				throw new IllegalArgumentException(
						"DeleteFiles: write protected: " + fileName);

			// if it is a directory, make sure it is empty
			if (f.isDirectory()) {
				String[] files = f.list();
				if (files.length > 0)
					throw new IllegalArgumentException(
							"DeleteFiles: directory not empty: " + fileName);
			}

			// explicit close associated stream for windows
			//
			// attempt to delete it
			boolean success = f.delete();

			if (!success)
				throw new IllegalArgumentException("DeleteFiles: deletion of "
						+ fileName + " failed");
		}
	}

	private static void printColor(String text, SimpleAttributeSet color) {
		if (ModelTest.buildGUI) {
			try {
				Document doc = XManager.getInstance().getPane().getDocument();
				doc.insertString(doc.getLength(), text, color);
			} catch (javax.swing.text.BadLocationException e) {
				System.err.println("Bad Location Exception");
			}
		} else {
			System.out.print(text);
		}
	}
	
	/**
	 * printRed
	 * 
	 * Prints to the main console text in red
	 */
	public static void printRed(String text) {
		printColor(text, XManager.redText);
	}
	
	/**
	 * printBlue
	 * 
	 * Prints to the main console text in blue
	 */

	public static void printBlue(String text) {
		printColor(text, XManager.blueText);
	}

	/******************
	 * RoundDoubleTo **************************
	 * 
	 * Rounds a double to a number of significant digits
	 *********************************************************/

	public static double roundDoubleTo(double decimal, int dplaces) {
		BigDecimal bd = new BigDecimal(decimal);
		bd = bd.setScale(dplaces, BigDecimal.ROUND_UP);
		return bd.doubleValue();
	}

	public static String asPercent(double decimal) {
		if (decimal > 100.0d) {
			decimal = 100.0d;
		}
		return format(decimal, 6, 2, false) + "%";
	}

	public static String format(double number, int totalLength, int decimalPlaces, boolean exp) {
		StringBuffer sb;
		String format;
		if (exp) {
			format = "%"+totalLength+"."+decimalPlaces+"e";
		} else {
			format = "%"+totalLength+"."+decimalPlaces+"f";
		}
		sb = new StringBuffer(String.format(Locale.ENGLISH, format, number));
		// normalize string to size 6
		for (int i = sb.length(); i < totalLength; i++) {
			sb.insert(0, " ");
		}

		return sb.toString();
	}
	
	/******************
	 * putSlashBeforeSpaces **************************
	 * 
	 * Inserts slash before paces in a path (MacOS X)
	 *********************************************************/

	public static String putSlashBeforeSpaces(String originalPath) {
		if (originalPath == null)
			return "";

		StringBuffer s = new StringBuffer(originalPath);

		for (int i = 0; i < s.length(); i++) {
			if (s.charAt(i) == ' ')
				s.insert(i++, '\\');
		}
		return s.toString();
	}

	/******************
	 * substituteSpaces **************************
	 * 
	 * substitutes %20 in spaces in a path (MacOS X)
	 *********************************************************/

	public static String substituteSpaces(String originalPath) {
		if (originalPath == null)
			return "";

		StringBuffer s = new StringBuffer(originalPath);

		for (int i = 0; i < s.length(); i++) {
			if (s.charAt(i) == ' ') {
				s.setCharAt(i, '0');
				s.insert(i++, '%');
				s.insert(i++, '2');
			}
		}
		return s.toString();
	}

	/******************
	 * putQuotesAroundSpaces **************************
	 * 
	 * Inserts quotes at correct locations in nominated path string. (Windows)
	 *********************************************************/

	public static String putQuotesAroundSpaces(String originalPath) {
		StringTokenizer st = new StringTokenizer(originalPath, File.separator);
		String newPath = "\"";
		while (st.hasMoreTokens()) {
			String token = st.nextToken();
			if (token.contains(" ")) {
				newPath += File.separator + "\"" + token + "\"";
			} else {
				newPath += File.separator + token;
			}
		}
		return newPath.substring(1, newPath.length());
	}

	/******************
	 * specialConcatStringArrays **************************
	 * 
	 * concatenates two string arrays but removinf the first string from the
	 * second
	 *********************************************************/

	public static String[] specialConcatStringArrays(String[] A, String[] B) {
		String[] C = new String[A.length + B.length - 1];
		System.arraycopy(A, 0, C, 0, A.length);
		System.arraycopy(B, 1, C, A.length, B.length - 1);
		return C;
	}

	/**
	 * Calculate invariable sites.
	 * 
	 * @param alignment the alignment
	 * 
	 * @return the invariable sites in the alignment
	 */
	public static int calculateInvariableSites (Alignment alignment) {

		//use this function to estimate a good starting value for the InvariableSites distribution.
		int numSites    = alignment.getSiteCount();
		int numSeqs     = alignment.getSequenceCount();
		int inv         = 0;
		int tmp;
		boolean tmpInv = true;
		for(int i=0; i < numSites; i++) {
			tmp     = indexOfChar(alignment.getData(0,i));
			tmpInv  = true;
			for(int j=0; j < numSeqs; j++) {
				if(indexOfChar(alignment.getData(j,i)) != tmp) { //if at least one difference in column i:
					j = numSeqs;    //we exit this for.
					tmpInv = false; //i is not an invariable site.
				}
			}
			if(tmpInv)
				inv++;
		}
		return inv;
	}
	
	private static int indexOfChar (char c) {
		char[] charSet = {'A', 'C', 'G', 'T', 'a', 'c', 'g', 't'};
		for (int charIndex = 0; charIndex < charSet.length; charIndex++)
			if (charSet[charIndex] == c)
				return charIndex % 4;

		return -1;
	}
	
	public static double calculateShannonSampleSize (Alignment alignment, boolean doJustShannonEntropy) {

		int numSites = alignment.getSiteCount();
		int numTaxa = alignment.getSequenceCount();
		
		//int    pattern[][]    = new int   [numSites][numSeqs];
		double freqs  [][]    = new double[numSites][4];
		byte   state  [][]    = new byte  [numSites][4];
		double siteS  []      = new double[numSites];
		int    sequences[]    = new int[numSites];
		double shannonEntropy = 0.0;

		//We simply count bases at positions and store in state[][]
		for(int i=0; i < numSites; i++) {
			for(int j=0; j < numTaxa; j++) {
				//state[i][indexOfchar(alignment.getData(j,i))]++;
				int index = indexOfChar(alignment.getData(j,i));
				if(index >= 0) {
					state[i][index]++;
					sequences[i]++;
				}
				//state[i][sP.pattern[j][i]]++;
			}
		}

		//For each alignment position, we calculate aminoacid frequencies. And also...
		//For each alignment position, we calculate Shannon Entropy based on previous frequencies.
		for(int i=0; i < numSites; i++) {
			for(int j=0; j < 4; j++) {
				//freqs[i][j] = (double)state[i][j]/(double)numSeqs;
				freqs[i][j] = (double)state[i][j]/(double)sequences[i];
				if(freqs[i][j] > 0)
					siteS[i]   += freqs[i][j]*Math.log(freqs[i][j])/Math.log(2);
			}
		}

		//We sum positions entropies over the whole alignment.
		for(int i=0; i < numSites; i++) {
			shannonEntropy += siteS[i];
		}

		if(doJustShannonEntropy) {
			return -1.0*shannonEntropy; //sum of shannon Entropy positions.
		} else { //if Options.SHANNON_NxL
			double meanShannonEntropy;
			double maxShannonEntropy = 0;
			double normalizedShannonEntropy = 0;
			meanShannonEntropy = -1.0*shannonEntropy/(double)numSites; //mean S for sites
			//let's normalize ShannonEntropy from 0 to 1:
			for(int i=0; i<4; i++) {
				maxShannonEntropy += (1.0/(double)4)*Math.log(1.0/(double)4)/Math.log(2);
			}
			//by this moment we normalize by a "rule of three"
			normalizedShannonEntropy = -1.0*meanShannonEntropy/maxShannonEntropy;
			return (double)numSites*(double)numTaxa*normalizedShannonEntropy; //NxL x averaged Shannon entropy
		}
	}
} // end of class
import java.awt.GraphicsEnvironment;
import java.io.BufferedInputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.IOException;
import java.io.InputStream;
import java.io.ObjectInput;
import java.io.ObjectInputStream;
import java.io.OutputStream;
import java.io.PrintStream;
import java.io.PrintWriter;
import java.io.PushbackReader;
import java.io.StringReader;
import java.net.InetAddress;
import java.net.UnknownHostException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Hashtable;
import java.util.List;
import java.util.Properties;
import java.util.Vector;
import mpi.MPI;
import mpi.MPIException;
import pal.tree.Tree;
import pal.tree.TreeParseException;
import es.uvigo.darwin.jmodeltest.exe.RunConsense;
import es.uvigo.darwin.jmodeltest.exe.RunPhyml;
import es.uvigo.darwin.jmodeltest.exe.RunPhymlHybrid;
import es.uvigo.darwin.jmodeltest.exe.RunPhymlMPJ;
import es.uvigo.darwin.jmodeltest.exe.RunPhymlThread;
import es.uvigo.darwin.jmodeltest.gui.XManager;
import es.uvigo.darwin.jmodeltest.io.AlignmentReader;
import es.uvigo.darwin.jmodeltest.io.HtmlReporter;
import es.uvigo.darwin.jmodeltest.io.TextInputStream;
import es.uvigo.darwin.jmodeltest.io.TextOutputStream;
import es.uvigo.darwin.jmodeltest.model.Model;
import es.uvigo.darwin.jmodeltest.observer.ConsoleProgressObserver;
import es.uvigo.darwin.jmodeltest.selection.AIC;
import es.uvigo.darwin.jmodeltest.selection.AICc;
import es.uvigo.darwin.jmodeltest.selection.BIC;
import es.uvigo.darwin.jmodeltest.selection.DT;
import es.uvigo.darwin.jmodeltest.selection.HLRT;
import es.uvigo.darwin.jmodeltest.selection.InformationCriterion;
import es.uvigo.darwin.jmodeltest.tree.TreeSummary;
import es.uvigo.darwin.jmodeltest.tree.TreeUtilities;
import es.uvigo.darwin.jmodeltest.utilities.Simulation;
import es.uvigo.darwin.jmodeltest.utilities.Utilities;
			stream.printf("%s", model.getName());
			stream.print(") selected by ");
			stream.printf("%s", criterion);
			stream.print("\nwith ");
			stream.printf("%s", programName);
			stream.print(" ");
			stream.printf("%s", CURRENT_VERSION);
			stream.print(" on ");
			java.util.Date current_time = new java.util.Date();
			stream.print(current_time.toString());
			stream.println("]");

			stream.print("\nBEGIN PAUP;");
			stream.print("\nLset");

			/* base frequencies */
			stream.print(" base=");
			if (model.ispF()) {
				stream.print("(");
				stream.printf("%.4f ", model.getfA());
				stream.printf("%.4f ", model.getfC());
				stream.printf("%.4f ", model.getfG());
				/* stream.printf("%.4f",model.fT); */
				stream.print(")");
			} else
				stream.print("equal");

			/* substitution rates */
			if (!model.ispT() && !model.ispR())
				stream.print(" nst=1");
			else if (model.ispT()) {
				stream.print(" nst=2 tratio=");
				stream.printf("%.4f", model.getTitv());
			} else if (model.ispR()) {
				stream.print(" nst=6  rmat=(");
				stream.printf("%.4f ", model.getRa());
				stream.printf("%.4f ", model.getRb());
				stream.printf("%.4f ", model.getRc());
				stream.printf("%.4f ", model.getRd());
				stream.printf("%.4f)", model.getRe());
				/* stream.print("1.0000)"); */
			}

			/* site rate variation */
			stream.print(" rates=");
			if (model.ispG()) {
				stream.print("gamma shape=");
				stream.printf("%.4f", model.getShape());
				stream.print(" ncat=");
				stream.printf("%d", model.getNumGammaCat());
			} else
				stream.print("equal");

			/* invariable sites */
			stream.print(" pinvar=");
			if (model.ispI())
				stream.printf("%.4f", model.getPinv());
			else
				stream.print("0");

			stream.print(";\nEND;");
			stream.print("\n--\n");
		}

		catch (Exception e) {
			e.printStackTrace();
		}
	}

	private void checkInputFiles() {
		// open data file
		File inputFile = options.getInputFile();
		MAIN_CONSOLE.print("\n\nReading data file \"" + inputFile.getName()
				+ "\"...");

		if (inputFile.exists()) {

			try {
				
				String alnString = ModelTestService.readAlignment(inputFile,
						options.getAlignmentFile());
				
				PushbackReader pr = new PushbackReader(
						new StringReader(alnString)
						);
				options.setAlignment(AlignmentReader.createAlignment(
						new PrintWriter(System.err), pr, true)); // file
				
				MAIN_CONSOLE.println(" OK.");
				MAIN_CONSOLE.println("  number of sequences: "
						+ options.getNumTaxa());
				MAIN_CONSOLE.println("  number of sites: "
						+ options.getNumSites());
			} catch (Exception e)// file cannot be read correctly
			{
				System.err.println("\nThe specified file \""
						+ inputFile.getAbsolutePath()
						+ "\" cannot be read as an alignment");
				MAIN_CONSOLE.println(" failed.\n" + e.getMessage());
				throw new InvalidArgumentException.InvalidAlignmentFileException(
						inputFile);
			}
		} else // file does not exist
		{
			System.err.println("\nThe specified file \""
					+ inputFile.getAbsolutePath() + "\" cannot be found");
			MAIN_CONSOLE.println(" failed.\n");
			throw new InvalidArgumentException.UnexistentAlignmentFileException(
					inputFile);
		}

		// open tree file if necessary
		if (options.userTopologyExists) {
			File treefile = options.getInputTreeFile();
			MAIN_CONSOLE.print("Reading tree file \"" + treefile.getName()
					+ "\"...");

			// read the tree in
			Tree tree = null;
			try {
				tree = TreeUtilities.readTree(treefile.getAbsolutePath());