public class ApplicationOptions implements Serializable {

	private static final long serialVersionUID = -3961572952922591321L;
	private static final int AMBIGUOUS_DATATYPE_STATE = 4;
	private static final int DEFAULT_SEED = 12345;

	/** Tree topology search algorithms */
	public static enum TreeSearch {
		NNI, SPR, BEST
	};

	private static ApplicationOptions instance;
	private static Vector<String> testOrder;

	private File inputDataFile, inputTreeFile;
	private File alignmentFile, treeFile;
	private File ckpFile;
	private File logFile;
	private Alignment alignment;
	private boolean isAmbiguous;
	private boolean forceCheckULnL = false;
	private double unconstrainedLnL = 0.0d;
	private int numPatterns = 0;

	// Newick user tree
	private String userTree;

	/*
	 * Number of threads for shared memory execution or static thread scheduling
	 */
	private int numberOfThreads = Runtime.getRuntime().availableProcessors();
	/*
	 * Number of threads for dynamic thread scheduling
	 */
	private File machinesFile;

	public boolean threadScheduling = false;

	public boolean doAIC = false;
	public boolean doAICc = false;
	public boolean doBIC = false;
	public boolean doDT = false;
	public boolean doDLRT = false;
	public boolean doHLRT = false;
	// whether to include the parameter base frequencies
	public boolean doF = false;
	public boolean doI = false; // whether to include the parameter pinv
	public boolean doG = false; // whether to include gamma
	public int numGammaCat = 4;

	public boolean backwardHLRTSelection = false;
	public double confidenceLevelHLRT = 0.01;

	public boolean writePAUPblock = false;
	public boolean doImportances = false;
	public boolean doModelAveraging = false;
	public boolean doAveragedPhylogeny = false;
	public double confidenceInterval = 1.0; // by default include all models

	// for specific simulations
	public boolean doingSimulations = false;
	public String simulationsName = "";

	// whether to use the same user fixed topology for all calculations
	public boolean userTopologyExists = false;
	// whether to use the same BIONJ-JC fixed tree for all calculations
	public boolean fixedTopology = false;
	// whether to optimize the BIONJ-model tree by ML
	public boolean optimizeMLTopology = true;

	public TreeSearch treeSearchOperations = TreeSearch.BEST;

	// Threshold for the guided search mode. A QST == 0.0 means no model
	// but the GTR one is optimized. A high QST means the whole set of
	// set of models will be optimized.
	private double guidedSearchThreshold = 0.0d;
	private boolean doClusteringSearch = false;
	private int heuristicInformationCriterion = InformationCriterion.IC_BIC;

	private int numSites;
	private int numTaxa;
	private int numBranches;
	private int numInvariableSites;
	private int numModels;

  private int rngSeed = DEFAULT_SEED;

	private String executionName;
	
	public String consensusType = "50% majority rule"; // consensus type
														// for model
														// averaged
														// phylogeny
	public boolean countBLasParameters = true; // whether to count branch
	// lengths as parameters

	private int substTypeCode = 0; // number of substitution types to
									// consider
	private ApplicationOptions() { }

	public void createCkpFile() {
		if (ckpFile != null && ckpFile.exists())
			return;
		
		try {
			if (executionName == null) {
				executionName = Utilities.getCurrentTime("yyyyMMddHHmmss");
			}
			if (ModelTestConfiguration.isCkpEnabled() && getInputFile() != null) {
				ckpFile = new File(ModelTestConfiguration.getLogDir() + File.separator + getInputFile().getName() 
						+ "." + executionName + ".ckp");
			} else {
				ckpFile = File.createTempFile("jmodeltest", ".ckp");
				ckpFile.deleteOnExit();
			}
		} catch (IOException e) {
			System.err.println("Error creating checkpointing file");
		}
	}
	
	public void createLogFile() {
		try {
			if (executionName == null) {
				executionName = Utilities.getCurrentTime("yyyyMMddHHmmss");
			}
			if (ModelTestConfiguration.isPhymlLogEnabled() && getInputFile() != null) {
				logFile = new File(ModelTestConfiguration.getLogDir() + File.separator + getInputFile().getName() 
						+ ".phyml." + executionName + ".log");
				if (!logFile.canWrite())
				{
					logFile = null;
				}
			} 
			if (logFile == null) {
				logFile = File.createTempFile("jmodeltest-phyml", ".log");
				logFile.deleteOnExit();
			}
		} catch (IOException e) {
			System.err.println("Error creating PhyML log file");
		}
	}
	
	public void deleteLogFile() {
		if (!ModelTestConfiguration.isPhymlLogEnabled()) {
			logFile.delete();
		}
	}
	
	// This method copies the input files into the scratch.
	// It is important to speed up the I/O in distributed memory.
	public void buildWorkFiles() throws IOException {
		File workDataFile = getAlignmentFile();
		if (!workDataFile.exists()) {
			workDataFile.createNewFile();
			workDataFile.deleteOnExit();
		}
		try {
			String alnStr = ModelTestService.readAlignment(inputDataFile, alignmentFile);

			PushbackReader pr = new PushbackReader(
					new StringReader(alnStr));
			setAlignment(
					AlignmentReader.createAlignment(
							new PrintWriter(System.err), pr, true));
		} catch (AlignmentParseException e) {
			e.printStackTrace();
		}

		if (userTree != null) {
			File workTreeFile = getTreeFile();
			if (!workTreeFile.exists()) {
				workTreeFile.createNewFile();
				workTreeFile.deleteOnExit();
			}
			TextOutputStream out = new TextOutputStream(getTreeFile()
					.getAbsolutePath());
			out.print(getUserTree());
			out.close();
		}
	}

	public static ApplicationOptions getInstance() {
		if (instance == null) {
			instance = new ApplicationOptions();
		}
		return instance;
	}

	/****************************
	 * setCandidateModels *********************** * Build the set of candidate
	 * models * *
	 ************************************************************************/

	public void setCandidateModels() {
		boolean includeModel;

		// fill in list of models
		ModelTest.setCandidateModels(new Model[getNumModels()]);

		if (substTypeCode < 4) {
			int j = 0;
			for (int i = 0; i < ModelTest.MAX_NUM_MODELS; i++) {
				includeModel = true;
				if (ModelConstants.substType[i] > substTypeCode)
					includeModel = false;
				if (!doF && !ModelConstants.equalBaseFrequencies[i]) // is F
					includeModel = false;
				if (!doI
						&& (ModelConstants.rateVariation[i] == 1 || ModelConstants.rateVariation[i] == 3))
					includeModel = false;
				if (!doG
						&& (ModelConstants.rateVariation[i] == 2 || ModelConstants.rateVariation[i] == 3))
					includeModel = false;

				if (includeModel) {
					loadModelConstraints(ModelTest.getCandidateModel(j), j, i);
					j++;
				}
			}
		} else {
			int j = 0;
			for (int ratesCount = 1; ratesCount <= 6; ratesCount++) {
				for (String partition : ModelConstants.fullModelSet
						.get(ratesCount)) {
					boolean baseFrequencies;
					int rateVariation;
					for (int k = 0; k < 8; k++) {
						baseFrequencies = (k / 4 == 1);
						rateVariation = (k % 4);
						includeModel = true;
						if (!doF && baseFrequencies)
							includeModel = false;
						if (!doI && (rateVariation == 1 || rateVariation == 3))
							includeModel = false;
						if (!doG && (rateVariation == 2 || rateVariation == 3))
							includeModel = false;
						if (includeModel) {
							loadModelConstraints(
									ModelTest.getCandidateModel(j), j,
									partition, baseFrequencies, rateVariation,
									ratesCount);
							j++;
						}
					}
				}
			}
		}

		if (ModelTest.buildGUI || ModelTest.testingOrder == null) {
			// set set of hypotheses for hLRTs in default order
			testOrder = new Vector<String>();
			// we need to reinitialize the hypotheses list in case it existed
			// before..
			testOrder.removeAllElements();

			if (doF)
				testOrder.add("freq");

			switch (substTypeCode) {
			case 0:
				testOrder.add("titv");
				testOrder.add("2ti4tv");
				break;
			case 1:
				testOrder.add("titv");
				testOrder.add("2ti");
				testOrder.add("2tv");
				break;
			case 2:
			case 3:
			case 4:
				testOrder.add("titv");
				testOrder.add("2ti");
				testOrder.add("2tv");
				testOrder.add("4tv");
			}

			if (doI)
				testOrder.add("pinv");

			if (doG)
				testOrder.add("gamma");

			// ModelTest.testingOrder.removeAllElements();

			ModelTest.testingOrder = testOrder;
		}
	}

	/************************
	 * loadModelConstraints ************************ * Loads initial definitions
	 * of models * * *
	 ***********************************************************************/

	public void loadModelConstraints(Model currentModel, int order, int modelNo) {
		int modelParameters, BL;
		boolean pF, pT, pV, pR, pI, pG;

		// initialize
		pT = pV = pR = pI = pG = false;

		pF = !ModelConstants.equalBaseFrequencies[modelNo];

		if (ModelConstants.numTransitions[modelNo] == 1
				&& ModelConstants.numTransversions[modelNo] == 1) {
			pT = true;
		} else if (ModelConstants.numTransitions[modelNo] > 1
				|| ModelConstants.numTransversions[modelNo] > 1) {
			pR = true;
		}

		/*
		 * find out first wich models in phyml give titv if
		 * (Model.numTransitions[modelNo] == 2 ||
		 * Model.numTransversions[modelNo] == 2) { pV = true; }
		 */

		switch (ModelConstants.rateVariation[modelNo]) {
		case 1:
			pI = true;
			break;
		case 2:
			pG = true;
			break;
		case 3:
			pI = true;
			pG = true;
			break;
		}

		if (countBLasParameters)
			BL = getNumBranches();
		else
			BL = 0;

		if (optimizeMLTopology)
			modelParameters = ModelConstants.freeParameters[modelNo] + BL + 1;
		else
			modelParameters = ModelConstants.freeParameters[modelNo] + BL;

		ModelTest.getCandidateModels()[order] = new Model(order + 1,
				ModelConstants.modelName[modelNo],
				ModelConstants.modelCode[modelNo], modelParameters, pF, pT, pV,
				pR, pI, pG, ModelConstants.numTransitions[modelNo],
				ModelConstants.numTransversions[modelNo]);
		
		if (ModelTest.getLoadedModels() != null && ModelTest.getLoadedModels().length > 0) {
			for (Model model : ModelTest.getLoadedModels()) {
				if (model.getName().equals(ModelTest.getCandidateModels()[order].getName()) 
						&& model.getLnL() > 0.0) {
					/* load from checkpoint */
					Model otherModel = ModelTest.getCandidateModels()[order];
					otherModel.setLnL(model.getLnL());
					otherModel.setLnLIgnoringGaps(model.getLnLIgnoringGaps());
					otherModel.setShape(model.getShape());
					otherModel.setfA(model.getfA());
					otherModel.setfC(model.getfC());
					otherModel.setfG(model.getfG());
					otherModel.setfT(model.getfT());
					otherModel.setRa(model.getRa());
					otherModel.setRb(model.getRb());
					otherModel.setRc(model.getRc());
					otherModel.setRd(model.getRd());
					otherModel.setRe(model.getRe());
					otherModel.setRf(model.getRf());
					otherModel.setPinv(model.getPinv());
					otherModel.setKappa(model.getKappa());
					otherModel.setTitv(model.getTitv());
					try {
						otherModel.setTreeString(model.getTreeString());
					} catch (TreeParseException e) {
						e.printStackTrace();
					}
					model.setLnL(0.0);
					break;
				}
			}
		}
	}

	public void loadModelConstraints(Model currentModel, int order,
			String partition, boolean baseFrequencies, int rateVariation,
			int ratesCount) {
		int modelParameters, BL, nTi, nTv;
		boolean pF, pT, pV, pR, pI, pG;

		pT = pV = pR = pI = pG = false;

		pF = baseFrequencies;

		switch (rateVariation) {
		case 1:
			pI = true;
			break;
		case 2:
			pG = true;
			break;
		case 3:
			pI = true;
			pG = true;
			break;
		}

		if (countBLasParameters)
			BL = getNumBranches();
		else
			BL = 0;

		modelParameters = (ratesCount - 1) + (pF ? 3 : 0) + (pI ? 1 : 0)
				+ (pG ? 1 : 0) + BL + (optimizeMLTopology ? 1 : 0);

		String modelName;

		if (ModelConstants.codeList.contains(partition)) {
			int index = ModelConstants.codeList.indexOf(partition);
			if (pF)
				index += 4;
			if (pI)
				index += 1;
			if (pG)
				index += 2;
			modelName = ModelConstants.modelName[index];
			nTi = ModelConstants.numTransitions[index];
			nTv = ModelConstants.numTransversions[index];
			if (nTi == 1 && nTv == 1) {
				pT = true;
			} else if (nTi > 1 || nTv > 1) {
				pR = true;
			}
		} else {
			modelName = partition + (pI ? "+I" : "") + (pG ? "+G" : "")
					+ (pF ? "+F" : "");
			nTi = 0;// ModelConstants.getNumberOfTransitions(partition);
			nTv = 0;// ModelConstants.getNumberOfTransversions(partition);
			pT = false;
			pR = true;
		}
		ModelTest.getCandidateModels()[order] = new Model(order + 1, modelName,
				partition, modelParameters, pF, pT, pV, pR, pI, pG, 0, 0);
		
		if (ModelTest.getLoadedModels() != null && ModelTest.getLoadedModels().length > 0) {
			for (Model model : ModelTest.getLoadedModels()) {
				if (model.getName().equals(ModelTest.getCandidateModels()[order].getName()) 
						&& model.getLnL() > 0.0) {
					/* load from checkpoint */
					Model otherModel = ModelTest.getCandidateModels()[order];
					otherModel.setLnL(model.getLnL());
					otherModel.setLnLIgnoringGaps(model.getLnLIgnoringGaps());
					otherModel.setShape(model.getShape());
					otherModel.setfA(model.getfA());
					otherModel.setfC(model.getfC());
					otherModel.setfG(model.getfG());
					otherModel.setfT(model.getfT());
					otherModel.setRa(model.getRa());
					otherModel.setRb(model.getRb());
					otherModel.setRc(model.getRc());
					otherModel.setRd(model.getRd());
					otherModel.setRe(model.getRe());
					otherModel.setRf(model.getRf());
					otherModel.setPinv(model.getPinv());
					otherModel.setKappa(model.getKappa());
					otherModel.setTitv(model.getTitv());
					try {
						otherModel.setTreeString(model.getTreeString());
					} catch (TreeParseException e) {
						e.printStackTrace();
					}
					model.setLnL(0.0);
					break;
				}
			}
		}
	}

	public File getInputFile() {
		return inputDataFile;
	}

	public File getInputTreeFile() {
		return inputTreeFile;
	}

	public void setInputFile(File file) {
		this.inputDataFile = file;
	}

	public void setInputTreeFile(File file) {
		this.inputTreeFile = file;
	}

	public File getAlignmentFile() {
		if (alignmentFile == null) {
			try {
				alignmentFile = File.createTempFile("jmodeltest", ".phy");
				alignmentFile.deleteOnExit();
			} catch (IOException e) {
				alignmentFile = inputDataFile;
			}
		}
		return alignmentFile;
	}

	public Alignment getAlignment() {
		return alignment;
	}

	public void setAlignment(Alignment alignment) {
		this.alignment = alignment;
		setNumTaxa(alignment.getSequenceCount());
		setNumSites(alignment.getSiteCount());
		setNumBranches(2 * numTaxa - 3);
		setNumInvariableSites(Utilities.calculateInvariableSites(alignment));

		DataType dt = alignment.getDataType();
		// check ambiguity
		isAmbiguous = false;
		if (alignment.getDataType().isAmbiguous()) {
			isAmbiguous = true;
		} else {
			for (int i = 0; i < numTaxa; i++) {
				String seq = alignment.getAlignedSequenceString(i);
				if (seq.indexOf(Alignment.GAP) >= 0) {
					isAmbiguous = true;
					break;
				}
				for (int j = 0; j < seq.length(); j++) {
					if (dt.getState(seq.charAt(j)) == AMBIGUOUS_DATATYPE_STATE) {
						isAmbiguous = true;
						break;
					}
				}
			}
		}
	}

	public boolean isForceCheckULnL() {
		return forceCheckULnL;
	}

	public void setForceCheckULnL(boolean forceCheckULnL) {
		this.forceCheckULnL = forceCheckULnL;
	}

	public double getUnconstrainedLnL() {
		return unconstrainedLnL;
	}

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

	public int getNumPatterns() {
		return numPatterns;
	}

	public void setNumPatterns(int numPatterns) {
		this.numPatterns = numPatterns;
	}

	public boolean isAmbiguous() {
		return isAmbiguous;
	}

	public String getExecutionName() {
		return executionName;
	}

	public void setExecutionName(String executionName) {
		this.executionName = executionName;
	}

	public File getTreeFile() {
		if (treeFile == null) {
			try {
				treeFile = File.createTempFile("jmodeltest", ".tree");
				treeFile.deleteOnExit();
			} catch (IOException e) {
				treeFile = inputTreeFile;
			}

		}
		return treeFile;
	}

	public File getCkpFile() {
		return ckpFile;
	}
	
	public File getLogFile() {
		return logFile;
	}

	public String getUserTree() {
		return userTree;
	}

	public void setUserTree(String tree) {
		this.userTree = tree;
	}

	/**
	 * @param substTypeCode
	 *            the substTypeCode to set
	 */
	public void setSubstTypeCode(int substTypeCode) {
		this.substTypeCode = substTypeCode;
		this.doClusteringSearch = (substTypeCode == 4);
	}

	/**
	 * @return the substTypeCode
	 */
	public int getSubstTypeCode() {
		return substTypeCode;
	}

	public double getSampleSize() {
		return numSites;
	}

	public static void setInstance(ApplicationOptions newInstance) {
		File alignmentFile = instance.alignmentFile;
		File treeFile = instance.treeFile;
		instance = newInstance;
		instance.alignmentFile = alignmentFile;
		instance.treeFile = treeFile;
	}

	public int getNumberOfThreads() {
		return numberOfThreads;
	}

	public void setNumberOfThreads(int numberOfThreads) {
		this.numberOfThreads = numberOfThreads;
	}

	public boolean isGuidedSearch() {
		return guidedSearchThreshold > 1e-6;
	}

	public double getGuidedSearchThreshold() {
		return guidedSearchThreshold;
	}

	public void setGuidedSearchThreshold(double guidedSearchThreshold) {
		this.guidedSearchThreshold = guidedSearchThreshold;
	}

	public void setMachinesFile(File machinesFile) throws FileNotFoundException {
		if (ModelTest.HOSTS_TABLE.containsKey(ModelTest.getHostname())) {

			setNumberOfThreads(
					ModelTest.HOSTS_TABLE.get(ModelTest.getHostname()));

		} else {
			System.err.println("");
			System.err.println("WARNING: Machines File format is wrong.");
			System.err.println("         This host: " + ModelTest.getHostname()
					+ " does not exist");
			System.err.println("         Using a single thread");
			System.err.println("");
			this.numberOfThreads = 1;
		}
	}

	public File getMachinesFile() {
		return machinesFile;
	}

	public boolean isClusteringSearch() {
		return doClusteringSearch;
	}

	public void setClusteringSearch(boolean doClusteringSearch) {
		this.doClusteringSearch = doClusteringSearch;
	}

	public int getHeuristicInformationCriterion() {
		return heuristicInformationCriterion;
	}

	public void setHeuristicInformationCriterion(
			int heuristicInformationCriterion) {
		this.heuristicInformationCriterion = heuristicInformationCriterion;
	}

	public int getNumSites() {
		return numSites;
	}

	public void setNumSites(int numSites) {
		this.numSites = numSites;
	}

	public int getNumTaxa() {
		return numTaxa;
	}

	public void setNumTaxa(int numTaxa) {
		this.numTaxa = numTaxa;
	}

	public int getNumBranches() {
		return numBranches;
	}

	public void setNumBranches(int numBranches) {
		this.numBranches = numBranches;
	}

	public int getNumInvariableSites() {
		return numInvariableSites;
	}

	public void setNumInvariableSites(int numInvariableSites) {
		this.numInvariableSites = numInvariableSites;
	}

	public int getNumModels() {
		return numModels;
	}

	public void setNumModels(int numModels) {
		this.numModels = numModels;
	}

  public int getRngSeed() {
    return rngSeed;
  }

  public void setRngSeed(int rngSeed) {
    this.rngSeed = rngSeed;
  }
} public class ModelTest {

	private ApplicationOptions options = ApplicationOptions.getInstance();

	/** The MPJ rank of the process. It is only useful if MPJ is running. */
	public static int MPJ_ME;
	/** The MPJ size of the communicator. It is only useful if MPJ is running. */
	public static int MPJ_SIZE;
	/** The MPJ running state. */
	public static boolean MPJ_RUN;

	public static final Properties USERDEF_PROPERTIES;
	
	// application constant definitions
	public static final int PRECISION = 4;
	public static final double INFINITY = 9999;
	public static final int MAX_NUM_MODELS = 88;
	public static final int MAX_NAME = 60;
	public static final String CURRENT_VERSION = "2.1.10 v20160303";
	public static final String programName = ("jModeltest");
	public static final String URL = "https://github.com/ddarriba/jmodeltest2";
	public static final String WIKI = "https://github.com/ddarriba/jmodeltest2";
	public static final String DISCUSSION_GROUP = "http://groups.google.com/group/jmodeltest";
	public static String CONFIG_FILE = "conf/jmodeltest.conf";
	public static final String UNKNOWN_HOSTNAME = "UNKNOWN";

	private static TextOutputStream MAIN_CONSOLE;
	private static TextOutputStream PHYML_CONSOLE = null;
	private static TextOutputStream CURRENT_OUT_STREAM;

	public static String[] arguments;

	private static boolean AICwasCalculated = false;
	private static boolean AICcwasCalculated = false;
	private static boolean BICwasCalculated = false;
	private static boolean DTwasCalculated = false;

	public static Vector<String> testingOrder; // order of the hLRTs
	public static String averagedTreeString; // model-averaged phylogeny in
												// Newick format
	public static enum ExecMode {
		CONSOLE, GUI
	};
	public static ExecMode execMode;
	
	private static AIC myAIC;
	private static AICc myAICc;
	private static BIC myBIC;
	private static DT myDT;
	private static HLRT myHLRT;

	private static RunConsense consensusAIC;
	private static RunConsense consensusAICc;
	private static RunConsense consensusBIC;
	private static RunConsense consensusDT;

	private static Model[] candidateModels;
	private static Model[] loadedModels;
	private static Model minAIC, minAICc, minBIC, minDT, minHLRT, minDLRT;

	private static String hostname;
	public static Hashtable<String, Integer> HOSTS_TABLE;

	// We can work under a GUI or in the command line
	public static boolean buildGUI = true;

	static {
		USERDEF_PROPERTIES = new Properties();
		InetAddress addr;
		try {
			addr = InetAddress.getLocalHost();
			// Get hostname
			hostname = addr.getHostName();
		} catch (UnknownHostException e) {
			hostname = UNKNOWN_HOSTNAME;
			System.err.println("WARNING: This host is unknown");
			// WARN AND DO NOTHING
		}
	}

	// constructor with GUI
	public ModelTest() {
		if (!GraphicsEnvironment.isHeadless()) {
			options.createLogFile();
			execMode = ExecMode.GUI;
			XManager.getInstance();
			// Check binary
	    	RunPhyml.checkBinary();
		} else {
			System.err.println("");
			System.err.println("ERROR: You are trying to run a GUI interface in a headless server.");
			finalize(-1);
		}
	}

	// constructor without GUI
	public ModelTest(String[] arg) {
		try {
			// open mainConsole
			MAIN_CONSOLE = new TextOutputStream(System.out);
			execMode = ExecMode.CONSOLE;
			ParseArguments();
			
			if (MPJ_ME == 0) {
				// print header information
				printHeader(MAIN_CONSOLE);

				// print citation information
				printCitation(MAIN_CONSOLE);

				// print notice information
				printNotice(MAIN_CONSOLE);
			}
			
			// Check binary
			if (!RunPhyml.checkBinary())
			{
				finalize(-1);
			}
        	
			options.createLogFile();
			options.createCkpFile();
			if (options.doingSimulations) {
				Simulation sim = new Simulation(options);
				sim.run();
			} else
				runCommandLine();
		} catch (Exception e) {
			e.printStackTrace();
		}
		finalize(0);
	}

	/****************************
	 * main ************************************ * Starts the application * * *
	 ***********************************************************************/

	public static boolean loadCheckpoint (File ckpFile) {
		try {
			InputStream file = new FileInputStream(ckpFile);
			InputStream buffer = new BufferedInputStream(file);
			ObjectInput input = new ObjectInputStream(buffer);
			loadedModels = (Model[]) input.readObject();
			int numModels = 0;
			for (Model model : loadedModels) {
				if (model.getLnL() > 0.0) numModels++;
			}
			MAIN_CONSOLE.println(" ok!");
			MAIN_CONSOLE.println("Loaded "+ numModels +" models");
			ApplicationOptions.getInstance().setNumModels(loadedModels.length);
			input.close();
		} catch (ClassNotFoundException ex) {
			MAIN_CONSOLE.println(" cannot perform input.");
			return false;
		} catch (IOException ex) {
			MAIN_CONSOLE.println(" cannot perform input.");
			return false;
		}
		return true;
	}
	
	public static void main(String[] args) {
		// initializing MPJ environment (if available)
		System.err.println("[MPI] Testing MPI environment... (" + hostname
				+ ")");
		try {
			arguments = MPI.Init(args);
			MPJ_ME = MPI.COMM_WORLD.Rank();
			MPJ_SIZE = MPI.COMM_WORLD.Size();
			System.err.println("[MPI] ... OK! [" + hostname + " (" + MPJ_ME
          + ")]");
			MPJ_RUN = true;
		} catch (MPIException e) {
			System.err.println("[MPI] Proceed without MPI");
			MPJ_ME = 0;
			MPJ_SIZE = 1;
			MPJ_RUN = false;
			arguments = args;
		} catch (Exception e) {
			System.err.println("[MPI] Proceed without MPI");
			MPJ_ME = 0;
			MPJ_SIZE = 1;
			MPJ_RUN = false;
			arguments = args;
		} catch (ExceptionInInitializerError e) {
			System.err.println("[MPI] Initializer error!");
			System.err.println(e.getMessage());
			System.exit(-1);
			MPJ_ME = 0;
			MPJ_SIZE = 1;
			MPJ_RUN = false;
			arguments = args;
		} catch (NoClassDefFoundError e) {
			System.err.println("[MPI] Proceed without MPI");
			MPJ_ME = 0;
			MPJ_SIZE = 1;
			MPJ_RUN = false;
			arguments = args;
		}

		if (arguments.length < 1) {
			buildGUI = true;
			new ModelTest();
		} else {
			buildGUI = false;
			new ModelTest(arguments);
		}

		ApplicationOptions.getInstance().getLogFile().delete();
	}

	/****************************
	 * runCommandLine ************************** * Organizes all the tasks that
	 * the program needs to carry out * * *
	 ***********************************************************************/
	public void runCommandLine() {

		if (MPJ_ME == 0) {
			// print the command line
			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE.print("Arguments =");
			for (int i = 0; i < arguments.length; i++)
				MAIN_CONSOLE.print(" " + arguments[i]);

			try {
				checkInputFiles();
			} catch (InvalidArgumentException.InvalidInputFileException e) {
				MAIN_CONSOLE.println(e.getMessage());
				finalize(-1);
			}

			// 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 if (options.getSubstTypeCode() == 3)
				options.setNumModels(11);
			else
				options.setNumModels(203);

			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);
			options.setCandidateModels();
		}
		// build set of models

		// calculate likelihoods with phyml in the command line
		RunPhyml runPhyml;
		if (MPJ_RUN) {
			if (options.threadScheduling && options.getNumberOfThreads() > 0) {
				runPhyml = new RunPhymlHybrid(MPJ_ME, MPJ_SIZE,
						new ConsoleProgressObserver(options), options,
						getCandidateModels(), options.getNumberOfThreads());
			} else {
				runPhyml = new RunPhymlMPJ(
						new ConsoleProgressObserver(options), options,
						getCandidateModels());
			}
		} else {
			runPhyml = new RunPhymlThread(new ConsoleProgressObserver(options),
					options, getCandidateModels());
		}
		runPhyml.execute();

		if (MPJ_ME == 0) {

			List<Model> bestModels = new ArrayList<Model>();
			// do AIC if selected
			if (options.doAIC) {
				myAIC = new AIC(options.writePAUPblock, options.doImportances,
						options.doModelAveraging, options.confidenceInterval);
				myAIC.compute();
				minAIC = myAIC.getMinModel();
				AICwasCalculated = true;
				bestModels.add(minAIC);
			}

			// do AICc if selected
			if (options.doAICc) {
				myAICc = new AICc(options.writePAUPblock,
						options.doImportances, options.doModelAveraging,
						options.confidenceInterval);
				myAICc.compute();
				minAICc = myAICc.getMinModel();
				AICcwasCalculated = true;
				if (!bestModels.contains(minAICc)) {
					bestModels.add(minAICc);
				}
			}

			// do BIC if selected
			if (options.doBIC) {
				myBIC = new BIC(options.writePAUPblock, options.doImportances,
						options.doModelAveraging, options.confidenceInterval);
				myBIC.compute();
				minBIC = myBIC.getMinModel();
				BICwasCalculated = true;
				if (!bestModels.contains(minBIC)) {
					bestModels.add(minBIC);
				}
			}

			// do DT if selected
			if (options.doDT) {
				myDT = new DT(options.writePAUPblock, options.doImportances,
						options.doModelAveraging, options.confidenceInterval);
				myDT.compute();
				minDT = myDT.getMinModel();
				DTwasCalculated = true;
			}

			if (options.isAmbiguous() && options.isForceCheckULnL()) {
				runPhyml.executeIgnoreGaps(bestModels.toArray(new Model[0]));
			}

			if (options.doAIC) {
				myAIC.print(MAIN_CONSOLE);
				if (options.doAveragedPhylogeny) {
					consensusAIC = new RunConsense(myAIC,
							options.consensusType, options.confidenceInterval);
				}
			}
			if (options.doAICc) {
				myAICc.print(MAIN_CONSOLE);
				if (options.doAveragedPhylogeny) {
					consensusAICc = new RunConsense(myAICc,
							options.consensusType, options.confidenceInterval);
				}
			}
			if (options.doBIC) {
				myBIC.print(MAIN_CONSOLE);
				if (options.doAveragedPhylogeny) {
					consensusBIC = new RunConsense(myBIC,
							options.consensusType, options.confidenceInterval);
				}
			}
			if (options.doDT) {
				myDT.print(MAIN_CONSOLE);
				if (options.doAveragedPhylogeny) {
					consensusDT = new RunConsense(myDT, options.consensusType,
							options.confidenceInterval);
				}
			}
			// do hLRT if selected
			if (options.doHLRT) {
				myHLRT = new HLRT(options);
				myHLRT.compute(!options.backwardHLRTSelection,
						options.confidenceLevelHLRT, options.writePAUPblock);
			}

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

			Tree bestAIC = myAIC != null ? myAIC.getMinModel().getTree() : null;
			Tree bestAICc = myAICc != null ? myAICc.getMinModel().getTree()
					: null;
			Tree bestBIC = myBIC != null ? myBIC.getMinModel().getTree() : null;
			Tree bestDT = myDT != null ? myDT.getMinModel().getTree() : null;

			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE
					.println("---------------------------------------------------------------");
			MAIN_CONSOLE
					.println("*                                                             *");
			MAIN_CONSOLE
					.println("*                    SELECTION SUMMARY                        *");
			MAIN_CONSOLE
					.println("*                                                             *");
			MAIN_CONSOLE
					.println("---------------------------------------------------------------");
			MAIN_CONSOLE.println("");

			TreeSummary treeSummary = new TreeSummary(bestAIC, bestAICc,
					bestBIC, bestDT, candidateModels);

			treeSummary.print(MAIN_CONSOLE);

			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE.println(" ");
			MAIN_CONSOLE.println("::Best Models::");
			MAIN_CONSOLE.println(" ");
			if (myAIC == null && myAICc == null && myBIC == null
					&& myDT == null) {
				MAIN_CONSOLE.println("No information criterion was selected.");
			}
			MAIN_CONSOLE
					.println("\tModel \t\tf(a) \tf(c) \tf(g) \tf(t) \tkappa \ttitv "
							+ "\tRa\tRb\tRc\tRd\tRe\tRf\tpInv \tgamma");
			MAIN_CONSOLE
					.println("----------------------------------------------------------------------------------------------------------------------------------------");

			boolean icValid = true;
			if (myAIC != null) {
				Model minModel = myAIC.getMinModel();
				MAIN_CONSOLE.println("AIC \t" + getModelRow(minModel));
				icValid &= myAIC.isValid();
				// MAIN_CONSOLE.println("average \t" +
				// getModelAverageRow(myAIC));
			}
			if (myBIC != null) {
				Model minModel = myBIC.getMinModel();
				MAIN_CONSOLE.println("BIC \t" + getModelRow(minModel));
				icValid &= myBIC.isValid();
				// MAIN_CONSOLE.println("average \t" +
				// getModelAverageRow(myBIC));
			}
			if (myAICc != null) {
				Model minModel = myAICc.getMinModel();
				MAIN_CONSOLE.println("AICc \t" + getModelRow(minModel));
				icValid &= myAICc.isValid();
				// MAIN_CONSOLE.println("average \t" +
				// getModelAverageRow(myAICc));
			}
			if (myDT != null) {
				Model minModel = myDT.getMinModel();
				MAIN_CONSOLE.println("DT \t" + getModelRow(minModel));
				icValid &= myDT.isValid();
				// MAIN_CONSOLE.println("average \t" +
				// getModelAverageRow(myDT));

				MAIN_CONSOLE.println(" ");
				MAIN_CONSOLE.println("Program is done.");

			}
			if (!icValid) {
				MAIN_CONSOLE.println("* There might not be information enough to achieve reliable results. Please inspect the results carefully.");
			}
			
			if (ModelTestConfiguration.isHtmlLogEnabled()) {
				HtmlReporter.buildReport(options,
						ModelTest.getCandidateModels(), null, treeSummary);
			}
		} // end root

	} // end of runCommandLine

	private String getModelRow(Model model) {
		StringBuilder sb = new StringBuilder();
		sb.append(model.getName() + "\t");
		if (model.getName().length() < 8)
			sb.append("\t");
		sb.append(Utilities.format(model.getfA(), 4, 2, false) + "\t"
				+ Utilities.format(model.getfC(), 4, 2, false) + "\t"
				+ Utilities.format(model.getfG(), 4, 2, false) + "\t"
				+ Utilities.format(model.getfT(), 4, 2, false) + "\t"
				+ Utilities.format(model.getKappa(), 4, 2, false) + "\t"
				+ Utilities.format(model.getTitv(), 4, 2, false) + "\t"
				+ Utilities.format(model.getRa(), 7, 3, false) + " "
				+ Utilities.format(model.getRb(), 7, 3, false) + " "
				+ Utilities.format(model.getRc(), 7, 3, false) + " "
				+ Utilities.format(model.getRd(), 7, 3, false) + " "
				+ Utilities.format(model.getRe(), 7, 3, false) + " "
				+ Utilities.format(model.getRf(), 7, 3, false) + " ");
		if (model.ispI()) {
			sb.append(Utilities.format(model.getPinv(), 7, 2, false));
		} else {
			sb.append("N/A");
		}
		sb.append("\t");
		if (model.ispG()) {
			sb.append(Utilities.format(model.getShape(), 7, 2, false));
		} else {
			sb.append("N/A");
		}
		return sb.toString();
	}

	@SuppressWarnings("unused")
	private String getModelAverageRow(InformationCriterion ic) {
		StringBuilder sb = new StringBuilder();
		sb.append(Utilities.format(ic.getAfA(), 4, 2, false) + "\t");
		sb.append(Utilities.format(ic.getAfC(), 4, 2, false) + "\t");
		sb.append(Utilities.format(ic.getAfG(), 4, 2, false) + "\t");
		sb.append(Utilities.format(ic.getAfT(), 4, 2, false) + "\t");
		sb.append(Utilities.format(ic.getAkappa(), 4, 2, false) + "\t");
		sb.append(Utilities.format(ic.getAtitv(), 4, 2, false) + "\t");
		sb.append(Utilities.format(ic.getaRa(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getaRb(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getaRc(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getaRd(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getaRe(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getaRf(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getApinvI(), 7, 4, false) + "\t");
		sb.append(Utilities.format(ic.getApinvIG(), 7, 4, false) + "\t");
		sb.append(Utilities.format(ic.getAshapeG(), 7, 3, false) + "\t");
		sb.append(Utilities.format(ic.getAshapeIG(), 7, 3, false) + "\t");
		return sb.toString();
	}

	/****************************
	 * ParseArguments **************************** * Parses the command line for
	 * jModeltest * *
	 ************************************************************************/

	public void ParseArguments() {
		int i, j;
		String arg = "";
		String error = "\nCOMMAND LINE ERROR: ";
		File ckpFile = null;
		boolean isInputFile = false;
		boolean isIcForHcSet = false;
		boolean getPhylip = false;
		try {
			i = 0;
			while (i < arguments.length) {
				if (!arguments[i].startsWith("-")) {
					System.err
							.println(error
									+ "Arguments must start with \"-\". The ofending argument was: "
									+ arguments[i] + ".");
					System.err.print("      Arguments: ");
					for (String argument : arguments) {
						System.err.print(argument + " ");
					}
					System.err.println("");

					CommandLineError();
					System.exit(1);
				}

				arg = arguments[i++];

				if (arg.equals("-d")) {
					if (i < arguments.length) {
						options.setInputFile(new File(arguments[i++]));
						isInputFile = true;
					} else {
						System.err.println(error
								+ "-d option requires an input filename.");
						CommandLineError();
					}
				}

				else if (arg.equals("-o")) {
					String outFile = arguments[i++];
					try {
						MAIN_CONSOLE = new TextOutputStream(new PrintStream(
								outFile));
					} catch (FileNotFoundException e1) {
						System.err
								.println("An error has ocurred while trying to open the output file \""
										+ outFile + "\" for writing");
						e1.printStackTrace();
						System.exit(-1);
					}
				} else if (arg.equals("-s")) {
					if (i < arguments.length) {
						String type = arguments[i++];
						try {
							int number = Integer.parseInt(type);
							switch (number) {
							case 3:
								options.setSubstTypeCode(0);
								break;
							case 5:
								options.setSubstTypeCode(1);
								break;
							case 7:
								options.setSubstTypeCode(2);
								break;
							case 11:
								options.setSubstTypeCode(3);
								break;
							case 203:
								options.setSubstTypeCode(4);
								break;
							default:
								System.err
										.println(error
												+ "-s substitution types have to be 3,5,7,11 only.");
								CommandLineError();
							}
						} catch (NumberFormatException e) {
							System.err
									.println(error
											+ "-s option requires a number for the substitution types: 3,5,7,11.");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-s option requires a number for the substitution types: 3,5,7,11.");
						CommandLineError();
					}
				}

				else if (arg.equals("-ckp")) {
					if (i < arguments.length) {
						ckpFile = new File(arguments[i++]);
					} else {
						System.err.println(error
								+ "-ckp option requires a checkpoint filename.");
						CommandLineError();
					}
				}
				else if (arg.equals("-f")) {
					options.doF = true;
				}

				else if (arg.equals("-i")) {
					options.doI = true;
				}

				else if (arg.equals("-seed")) {
					if (i < arguments.length) {
						try {
							int rngSeed;
							String type = arguments[i++];
							rngSeed = Integer.parseInt(type);
              options.setRngSeed(rngSeed);
						} catch (NumberFormatException e) {
							System.err
									.println(error
											+ "-seed option requires a RNG seed number.");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-g option requires a number of gamma categories.");
						CommandLineError();
					}
				}
				else if (arg.equals("-g")) {
					if (i < arguments.length) {
						try {
							options.doG = true;
							String type = arguments[i++];
							options.numGammaCat = Integer.parseInt(type);
						} catch (NumberFormatException e) {
							System.err
									.println(error
											+ "-g option requires a number of gamma categories.");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-g option requires a number of gamma categories.");
						CommandLineError();
					}
				}

				else if (arg.equals("-G")) {
					if (i < arguments.length) {
						try {
							Double threshold = Double
									.parseDouble(arguments[i++]);
							options.setGuidedSearchThreshold(threshold);
						} catch (NumberFormatException e) {
							System.err.println(error
									+ "-G option requires a threshold.");
							CommandLineError();
						}
					} else {
						System.err.println(error
								+ "-G option requires a threshold.");
						CommandLineError();
					}
				}

				else if (arg.equals("-H")) {
					if (i < arguments.length) {
						String criterion = arguments[i++];
						if (criterion.equals("AIC")) {
							options.setHeuristicInformationCriterion(InformationCriterion.IC_AIC);
						} else if (criterion.equals("BIC")) {
							options.setHeuristicInformationCriterion(InformationCriterion.IC_BIC);
						} else if (criterion.equals("AICc")) {
							options.setHeuristicInformationCriterion(InformationCriterion.IC_AICc);
						} else {
							System.err
									.println(error
											+ "-H argument is invalid (AIC, BIC, AICc).");
							CommandLineError();
						}
						isIcForHcSet = true;
					} else {
						System.err
								.println(error
										+ "-H option requires an argument (AIC, BIC, AICc).");
						CommandLineError();
					}
				} else if (arg.equals("-t")) {
					if (i < arguments.length) {
						String type = arguments[i++];

						if (type.equalsIgnoreCase("fixed")) {
							options.fixedTopology = true;
							options.optimizeMLTopology = false;
							options.userTopologyExists = false;
						} else if (type.equalsIgnoreCase("BIONJ")) {
							options.fixedTopology = false;
							options.optimizeMLTopology = false;
							options.userTopologyExists = false;
						} else if (type.equalsIgnoreCase("ML")) {
							options.fixedTopology = false;
							options.optimizeMLTopology = true;
							options.userTopologyExists = false;
						} else {
							System.err
									.println(error
											+ "-t option requires a type of base tree for likelihod calculations: "
											+ "\"fixed\", \"BIONJ\" or \"ML\" only");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-t option requires a type of base tree for likelihod calculations: "
										+ "fixed, BIONJ or ML");
						CommandLineError();
					}
				} else if (arg.equals("-getPhylip")) {
					getPhylip = true;
				}

				else if (arg.equals("-u")) {
					if (i < arguments.length) {
						options.setInputTreeFile(new File(arguments[i++]));
						options.fixedTopology = false;
						options.optimizeMLTopology = false;
						options.userTopologyExists = true;
					} else {
						System.err
								.println(error
										+ "-u option requires an file name for the tree file");
						CommandLineError();
					}
				}

				else if (arg.equals("-n")) {
					if (i < arguments.length) {
						options.setExecutionName(arguments[i++]);
					} else {
						System.err
								.println(error
										+ "-n option requires a name for the execution");
						CommandLineError();
					}
				}
				
				else if (arg.equals("-S")) {
					if (i < arguments.length) {
						String type = arguments[i++];

						if (type.equalsIgnoreCase("NNI")) {
							options.treeSearchOperations = ApplicationOptions.TreeSearch.NNI;
						} else if (type.equalsIgnoreCase("SPR")) {
							options.treeSearchOperations = ApplicationOptions.TreeSearch.SPR;
						} else if (type.equalsIgnoreCase("BEST")) {
							options.treeSearchOperations = ApplicationOptions.TreeSearch.BEST;
						} else {
							System.err
									.println(error
											+ "-S option requires a type of tree topology search operation: "
											+ "\"NNI\", \"SPR\" or \"BEST\" only");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-S option requires a type of tree topology search operation: "
										+ "\"NNI\", \"SPR\", \"BEST\"");
						CommandLineError();
					}
				}

				else if (arg.equals("-AIC")) {
					options.doAIC = true;
				}

				else if (arg.equals("-AICc")) {
					options.doAICc = true;
				}

				else if (arg.equals("-BIC")) {
					options.doBIC = true;
				}

				else if (arg.equals("-DT")) {
					options.doDT = true;
				}

				else if (arg.equals("-uLnL")) {
					options.setForceCheckULnL(true);
				}

				else if (arg.equals("-p")) {
					options.doImportances = true;
				}

				else if (arg.equals("-v")) {
					options.doImportances = true;
					options.doModelAveraging = true;
				}

				else if (arg.equals("-w")) {
					options.writePAUPblock = true;
				}

				else if (arg.equals("-c")) {
					if (i < arguments.length) {
						try {
							String type = arguments[i++];
							options.confidenceInterval = Double
									.parseDouble(type);
						} catch (NumberFormatException e) {
							System.err
									.println(error
											+ "-c option requires a number (0-1) for the model selection confidence interval.");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-c option requires a number (0-1) for the model selection confidence interval.");
						CommandLineError();
					}
				}

				else if (arg.equals("-help")) {
					PrintUsage();
				}

				else if (arg.equals("-hLRT")) {
					options.doHLRT = true;
				}

				else if (arg.equals("-O")) {
					if (i < arguments.length) {
						String type = arguments[i++];
						char[] array = type.toCharArray();
						Arrays.sort(array);

						String validString = "";

						if (type.length() == 5) {
							validString = "fgptv";
						} else if (type.length() == 6) {
							validString = "fgptvw";
						} else if (type.length() == 7) {
							validString = "fgptvwx";
						} else {
							System.err
									.println(error
											+ "-O option requires a 5, 6 or 7 specific letter string with the order of tests (default is ftvwxgp)"
											+ "\n            f=freq, t=titvi, v=2ti4tv(subst=3)/2ti(subst>3), w=2tv, x=4tv, g=gamma, p=pinv"
											+ "\n            this argument is used only if -hLRT argument is set"
											+ "\n            'f','t','v','g','p' are mandatory in any order. 'w' is optional, and 'x' requires 'w' to be present"
											+ "\n            thus, length should be 5, 6 *including 'w') or 7 (including both 'w' and 'x')"
											+ "\n            e.g., -hLRT -O gpfvwxt"
											);
							CommandLineError();
						}

						char[] valid = validString.toCharArray();
						if (!Arrays.equals(array, valid)) {
							System.err
							.println(String.valueOf(array) + " " + String.valueOf(valid));
							System.err
							.println(error
									+ "-O option requires a 5, 6 or 7 specific letter string with the order of tests (default is ftvwxgp)"
									+ "\n            f=freq, t=titvi, v=2ti4tv(subst=3)/2ti(subst>3), w=2tv, x=4tv, g=gamma, p=pinv"
									+ "\n            this argument is used only if -hLRT argument is set"
									+ "\n            'f','t','v','g','p' are mandatory in any order. 'w' is optional, and 'x' requires 'w' to be present"
									+ "\n            thus, length should be 5, 6 *including 'w') or 7 (including both 'w' and 'x')"
									+ "\n            e.g., -hLRT -O gpfvwxt"
									);
							CommandLineError();
						} else {
							testingOrder = new Vector<String>();
							for (j = 0; j < type.length(); j++) {
								if (type.charAt(j) == 'f')
									testingOrder.addElement("freq");
								else if (type.charAt(j) == 't')
									testingOrder.addElement("titv");
								else if (type.charAt(j) == 'v') {
									if (options.getSubstTypeCode() == 0)
										testingOrder.addElement("2ti4tv");
									else if (options.getSubstTypeCode() >= 1)
										testingOrder.addElement("2ti");
								} else if (type.charAt(j) == 'w') {
									if (options.getSubstTypeCode() >= 1)
										testingOrder.addElement("2tv");
								} else if (type.charAt(j) == 'x') {
									if (options.getSubstTypeCode() > 1)
										testingOrder.addElement("4tv");
								} else if (type.charAt(j) == 'g')
									testingOrder.addElement("gamma");
								else if (type.charAt(j) == 'p')
									testingOrder.addElement("pinv");
							}
						}
					}
				}

				else if (arg.equals("-dLRT")) {
					options.doDLRT = true;
				}

				else if (arg.equals("-r")) {
					options.backwardHLRTSelection = true;
				}

				else if (arg.equals("-h")) {
					if (i < arguments.length) {
						try {
							String type = arguments[i++];
							options.confidenceLevelHLRT = Double
									.parseDouble(type);
						} catch (NumberFormatException e) {
							System.err
									.println(error
											+ "-h option requires a number (0-1) for the hLRT confidence interval.");
							CommandLineError();
						}
					} else {
						System.err
								.println(error
										+ "-h option requires a number (0-1) for the hLRT confidence interval.");
						CommandLineError();
					}
				}

				else if (arg.equals("-a")) {
					options.doAveragedPhylogeny = true;
				}

				else if (arg.equals("-z")) {
					options.consensusType = "strict";
				}

				else if (arg.equals("-sims")) {
					if (i < arguments.length) {
						options.simulationsName = arguments[i++];
						options.doingSimulations = true;
					} else {
						System.err
								.println(error
										+ "-sims option requires a name for the simulations files");
						CommandLineError();
					}

				} else if (arg.equals("-machinesfile")) {
					if (i < arguments.length) {

						File machinesFile = new File(arguments[i++]);
						if (!(machinesFile.exists() && machinesFile.canRead())) {
							if (MPJ_ME == 0) {
								System.err
										.println(error
												+ "Machines file does not exists or it is not readable");
							}
							CommandLineError();
						}

						boolean hostsError = false;
						try {
							TextInputStream machinesInputStream = new TextInputStream(
									machinesFile.getAbsolutePath());
							String line;

							HOSTS_TABLE = new Hashtable<String, Integer>();
							while ((line = machinesInputStream.readLine()) != null) {
								String hostProcs[] = line.split(":");
								if (hostProcs.length == 2) {
									try {

										int numberOfThreads = Integer
												.parseInt(hostProcs[1]);

										HOSTS_TABLE.put(hostProcs[0],
												new Integer(numberOfThreads));

									} catch (NumberFormatException e) {
										hostsError = true;
										break;
										// Warn and continue with 1 processor
									}
								} else {
									hostsError = true;
								}

								if (hostsError) {
									if (MPJ_ME == 0) {
										System.err.println("");
										System.err
												.println("WARNING: Machines File format is wrong.");
										System.err
												.println("         Each line should have the following format:");
										System.err
												.println("         HOSTNAME:NUMBER_OF_PROCESORS");
										System.err
												.println("Using a single thread");
										System.err.println("");
									}
									HOSTS_TABLE = null;
									options.setNumberOfThreads(1);
								}
							}
							options.setMachinesFile(machinesFile);

						} catch (FileNotFoundException e) {
							System.err
									.println(error
											+ "Machines file does not exists or it is not readable");
							CommandLineError();
						}
					} else {
						System.err.println(error
								+ "-machinesfile option requires a filename");
						CommandLineError();
					}
				} else if (arg.equals("-tr")) {
					if (HOSTS_TABLE != null) {
						String type = arguments[i++];
						System.err
								.println("WARNING: Machines File has been specified. -tr "
										+ type + " argument will be ignored.");
					} else {
						if (i < arguments.length) {
							try {
								String type = arguments[i++];
								options.setNumberOfThreads(Integer
										.parseInt(type));
								options.threadScheduling = true;
							} catch (NumberFormatException e) {
								System.err
										.println(error
												+ "-tr option requires the number of processors to compute.");
								CommandLineError();
							}
						} else {
							System.err
									.println(error
											+ "-tr option requires the number of processors to compute.");
							CommandLineError();
						}
					}

				} else if (arg.equals("--set-local-config")) {
					if (i < arguments.length) {
						String localConfigFile = arguments[i++];
						File testFile = new File(localConfigFile);
						if (!(testFile.exists() && testFile.canRead())) {
							System.err.println(error
									+ "Config file " + localConfigFile +" does not exist or cannot be read.");
							CommandLineError();
						}
						ModelTest.CONFIG_FILE = testFile.getAbsolutePath();
					} else {
						System.err.println(error
								+ "--set-local-config option requires a config file.");
						CommandLineError();
					}
				}  else if (arg.equals("--set-property")) {
					if (i < arguments.length) {
						String propertyDef = arguments[i++];
						if (propertyDef.indexOf("=") == -1) {
							System.err.println(error
									+ "--set-property option requires a property definition with the following format:");
							System.err.println("                    --set-property <property_name>=<property_value>");
							System.err.println("              e.g., --set-property log-dir=myHome/myLogDir");
							CommandLineError();
						}
						String propertyName = propertyDef.split("=")[0].trim();
						String propertyValue = propertyDef.split("=")[1].trim();
						USERDEF_PROPERTIES.setProperty(propertyName, propertyValue);
					} else {
						System.err.println(error
								+ "--set-property option requires a property definition.");
						System.err.println(error
								+ "e.g., --set-property log-dir=myHome/myLogDir");
						CommandLineError();
					}
				}
				else {
					System.err.println(error + "the argument \" " + arg
							+ "\" is unknown. Check its syntax.");
					CommandLineError();
				}
			} // while
			if (!isInputFile) {
				System.err.println(error
						+ "Input File is required (-d argument)");
				CommandLineError();
			}
			if (getPhylip) {
				MAIN_CONSOLE.print("\n\nReading data file \"" + options.getInputFile().getName()
						+ "\"...");

				if (options.getInputFile().exists()) {

					try {
						File outputFile = new File(options.getInputFile().getAbsolutePath() + ".phy");
						ModelTestService.readAlignment(options.getInputFile(),
								outputFile, false);
						MAIN_CONSOLE.println(" OK.");
						MAIN_CONSOLE.println("Result written into " + outputFile.getPath());
						MAIN_CONSOLE.println("");
					} catch (Exception e)// file cannot be read correctly
					{
						System.err.println("\nThe specified file \""
								+ options.getInputFile().getAbsolutePath()
								+ "\" cannot be read as an alignment");
						MAIN_CONSOLE.println(" failed.\n" + e.getMessage());
						throw new InvalidArgumentException.InvalidAlignmentFileException(
								options.getInputFile());
					}
				} else // file does not exist
				{
					System.err.println("\nThe specified file \""
							+ options.getInputFile().getAbsolutePath() + "\" cannot be found");
					throw new InvalidArgumentException.UnexistentAlignmentFileException(
							options.getInputFile());
				}
				finalize(0);
			}
			
			if (ckpFile != null) {
				if (!loadCheckpoint(ckpFile)) {
					System.err.println("\nThe specified checkpoint file \""
							+ options.getInputFile().getAbsolutePath()
							+ "\" cannot be read");
					finalize(0);
				}
			}
			
			if (testingOrder != null && !options.doHLRT) {
				System.err.println("\nWARNING: Hypothesis testing order has been set, but hierarchical Likelihood Ratio Test was not selected with '-hLRT' parameter");
				System.err.println("         Thus, this option has no effect\n");
			}
			
			if (isIcForHcSet && options.getSubstTypeCode() != 4) {
				System.err.println("\nWARNING: Information Criterion for hierarchical clustering has been specified, but hierarchical clustering applies only for 203 substitution schemes");
				System.err.println("         Thus, this option has no effect\n");
			}
		} catch (Exception e) {
			e.printStackTrace();
		}

		// add default IC
		options.doBIC |= !(options.doAIC || options.doAICc || options.doDT);
	}

	static public void CommandLineError() {
		System.err.println("\nPlease run with -help argument for info about the usage\n");
		System.exit(1);
	}
	
	/****************************
	 * PrintUsage **************************** * Prints command line usage * *
	 ************************************************************************/
	static public void PrintUsage() {
		if (MPJ_ME == 0) {
			System.err.println("\njModelTest command usage");
			System.err.println(
			    "java -jar jModelTest.jar -d sequenceFileName"
			+ "\n                        [-getPhylip]"
			+ "\n                        [-ckp checkpointFileName.ckp]"
			+ "\n                        [-n executionName]"
			+ "\n                        [-t fixed|BIONJ|ML] [-u userTreeFileName] [-o outputFileName]"
			+ "\n                        [-S NNI|SPR|BEST]"
			+ "\n                        [-AIC] [-AICc] [-BIC] [-DT] [-c confidenceInterval]"
			+ "\n                        [-s 3|5|7|11|203]"
			+ "\n                        [-f] [-i] [-g numberOfCategories]"
			+ "\n                        [-uLNL]"
			+ "\n                        [-dLRT] [-h confidenceInterval] [-hLRT] [-O {ftvwxgp}]"
			+ "\n                        [-a] [-z] [-p] [-v] [-w]"
			+ "\n                        [-tr numberOfThreads] [-machinesfile machinesFileName]"
			);
			String usage = 
					  "\n     -a"
					+ "\n         estimate model-averaged phylogeny for each active criterion (e.g., -a) (default is false)"
					+ "\n\n     -AIC"
					+ "\n         calculate the Akaike Information Criterion (e.g., -AIC) (default is false)"
					+ "\n\n     -AICc"
					+ "\n         calculate the corrected Akaike Information Criterion (e.g., -AICc) (default is false)"
					+ "\n\n     -BIC"
					+ "\n         calculate the Bayesian Information Criterion (e.g., -BIC) (default is false)"
					+ "\n\n     -DT"
					+ "\n         calculate the decision theory criterion (e.g., -DT) (default is false)"
					+ "\n\n     -c confidenceInterval"
					+ "\n         confidence interval (e.g., -c 90) (default is 100)"
					+ "\n\n     -ckp checkpointFileName"
					+ "\n         Loads a checkpointing file"
					+ "\n\n     -d sequenceFileName"
					+ "\n         input data file (e.g., -d data.phy)"
					+ "\n\n     -dLRT"
					+ "\n         do dynamical likelihood ratio tests (e.g., -dLRT)(default is false)"
					+ "\n\n     -f"
					+ "\n         include models with unequals base frecuencies (e.g., -f) (default is false)"
					+ "\n\n     -g numberOfCategories"
					+ "\n         include models with rate variation among sites and number of categories (e.g., -g 8) (default is false & 4 categories)"
					+ "\n\n     -G threshold"
					+ "\n         heuristic search. Requires a threshold > 0 (e.g., -G 0.1)"
					+ "\n\n     -getPhylip"
					+ "\n         converts the input file into phylip format"
					+ "\n\n     -h confidenceInterval"
					+ "\n         confidence level for the hLRTs (e.g., -a0.002) (default is 0.01)"
					+ "\n\n     -H informationCriterion"
					+ "\n         information criterion for clustering search (AIC, AICc, BIC). (default is BIC)"
					+ "\n         this argument applies only for 203 substitution schemes (e.g., -s 203 -H AIC)"
					+ "\n\n     -help"
					+ "\n         displays this help message"
					+ "\n\n     -hLRT"
					+ "\n         do hierarchical likelihood ratio tests (default is false)"
					+ "\n         hypothesis testing order can be specified with -O argument"
					+ "\n\n     -i"
					+ "\n         include models with a proportion invariable sites (e.g., -i) (default is false)"
					+ "\n\n     -machinesfile manchinesFileName"
					+ "\n         gets the processors per host from a machines file"
					+ "\n\n     -n executionName"
					+ "\n         execution name for appending to the log filenames (default: current time yyyyMMddhhmmss)"
					+ "\n\n     -o outputFileName"
					+ "\n         set output file (e.g., -o jmodeltest.out)"
					+ "\n\n     -O hypothesisOrder"
					+ "\n         hypothesis order for the hLRTs (e.g., -hLRT -O gpftv) (default is ftvwxgp)"
					+ "\n            f=freq, t=titvi, v=2ti4tv(subst=3)/2ti(subst>3), w=2tv, x=4tv, g=gamma, p=pinv"
					+ "\n            this argument is used only if -hLRT argument is set"
					+ "\n            'f','t','v','g','p' are mandatory in any order. 'w' is optional, and 'x' requires 'w' to be present"
					+ "\n            thus, length should be 5, 6 *including 'w') or 7 (including both 'w' and 'x')"
					+ "\n            e.g., -hLRT -O gpfvwxt"
					+ "\n\n     -p"
					+ "\n         calculate parameter importances (e.g., -p) (default is false)"
					+ "\n\n     -r"
					+ "\n         backward selection for the hLRT (e.g., -r) (default is forward)"
					+ "\n\n     -s numberOfSubstitutionSchemes"
					+ "\n         number of substitution schemes (e.g., -s 11) (it has to be 3,5,7,11,203; default is 3)"
					+ "\n\n     -seed RNG_seed"
					+ "\n         set the seed for the random number generator (default is 12345; use -1 for a random one)"
					+ "\n\n     --set-local-config localConfigurationFile"
					+ "\n         set a local configuration file in replacement of conf/jmodeltest.conf"
					+ "\n\n     --set-property propertyName=propertyValue"
					+ "\n         set a new value for a property contained in the configuration file (conf/jmodeltest.conf)"
					+ "\n\n     -S NNI|SPR|BEST"
					+ "\n         tree topology search operation option (NNI (fast), SPR (a bit slower), BEST (best of NNI and SPR)) (default is BEST)"
					+ "\n\n     -t fixed|BIONJ|ML"
					+ "\n             base tree for likelihood calculations (e.g., -t BIONJ)"
					+ "\n             fixed  (common BIONJ-JC topology)"
					+ "\n             BIONJ  (Neighbor-Joining topology)"
					+ "\n             ML     (Maximum Likelihood topology) (default)"
					+ "\n\n     -tr numberOfThreads"
					+ "\n         number of threads to execute (default is "
					+ Runtime.getRuntime().availableProcessors()
					+ ")"
					+ "\n\n     -u treeFileName"
					+ "\n         user tree for likelihood calculations  (e.g., -u data.tre)"
					+ "\n\n     -uLnL"
					+ "\n         calculate delta AIC,AICc,BIC against unconstrained likelihood (e.g., -uLnL)"
					+ "\n\n        (default is false if the input alignment has gaps or ambiguous characters)"
					+ "\n\n     -v"
					+ "\n         do model averaging and parameter importances (e.g., -v) (default is false)"
					+ "\n\n     -w"
					+ "\n         write PAUP block (e.g., -w) (default is false)"
					+ "\n\n     -z"
					+ "\n         strict consensus type for model-averaged phylogeny (e.g., -z) (default is majority rule)"
					+ "\n\n Command line: java -jar jModeltest.jar -d sequenceFileName [arguments]"
					+ "\n\n Example: java -jar jModeltest.jar -d sequenceFileName -i -f -g 4 -BIC -AIC -AICc -DT -v -a -w";
			System.err.println(usage);
			System.err.println(" ");
		}
		System.exit(0);
	}

	/****************************
	 * printHeader ****************************** * Prints header information at
	 * start * * *
	 ***********************************************************************/

	static public void printHeader(TextOutputStream stream) {
		// we can set styles using the editor pane
		// I am using doc to stream ....
		stream.print("-------------------------- ");
		stream.print("jModeltest " + CURRENT_VERSION);
		stream.println(" --------------------------");
		stream.println("(c) 2011-onwards D. Darriba, G.L. Taboada, R. Doallo and D. Posada,");
		stream.println("(1) Department of Biochemistry, Genetics and Immunology");
		stream.println("    University of Vigo, 36310 Vigo, Spain.");
		stream.println("(2) Department of Electronics and Systems");
		stream.println("    University of A Coruna, 15071 A Coruna, Spain.");
		stream.println("e-mail: ddarriba@udc.es, dposada@uvigo.es");
		stream.println("--------------------------------------------------------------------------------");
		stream.println(" ");
		java.util.Date current_time = new java.util.Date();
		stream.println(current_time.toString());
		stream.println(System.getProperty("os.name") + " "
				+ System.getProperty("os.version") + ", arch: "
				+ System.getProperty("os.arch") + ", bits: "
				+ System.getProperty("sun.arch.data.model") + ", numcores: "
				+ Runtime.getRuntime().availableProcessors());
		stream.println(" ");
	}

	/****************************
	 * printNotice ****************************** * Prints notice information at
	 * start up * * *
	 ***********************************************************************/

	static public void printNotice(TextOutputStream stream) {
		// stream.println("\n******************************* NOTICE ************************************");
		stream.println("jModelTest " + CURRENT_VERSION);
		stream.println("Copyright (C) 2011 D. Darriba, G.L. Taboada, R. Doallo and D. Posada");
		stream.println("This program comes with ABSOLUTELY NO WARRANTY");
		stream.println("This is free software, and you are welcome to redistribute it under certain");
		stream.println("conditions");
		stream.println(" ");
		stream.println("Notice: This program may contain errors. Please inspect results carefully.");
		stream.println(" ");
		// stream.println("***************************************************************************\n");

	}

	/****************************
	 * printCitation ****************************** * Prints citation
	 * information at start up * * *
	 ***********************************************************************/

	static public void printCitation(TextOutputStream stream) {
		// stream.println("\n******************************* CITATION *********************************");
		stream.println("--------------------------------------------------------------------------------");
		stream.println("Citation: Darriba D, Taboada GL, Doallo R and Posada D. 2012.");
		stream.println("          \"jModelTest 2: more models, new heuristics and parallel computing\".");
		stream.println("          Nature Methods 9(8), 772.");
		stream.println("--------------------------------------------------------------------------------");
		stream.println(" ");
		// stream.println("***************************************************************************\n");

	}

	/****************************
	 * writePaupBlockk *************************** * Prints a block of PAUP
	 * commands for the best model * *
	 ************************************************************************/
	static public void WritePaupBlock(TextOutputStream stream,
			String criterion, Model model) {
		try {
			stream.println("\n--\nPAUP* Commands Block:");
			stream.println(" If you want to load the selected model and associated estimates in PAUP*,");
			stream.println(" attach the next block of commands after the data in your PAUP file:");
			stream.print("\n[!\nLikelihood settings from best-fit model (");
			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());
			} catch (IOException e) {
				System.err.println("\nThe specified tree file \""
						+ treefile.getName() + "\" cannot be found");
				MAIN_CONSOLE.println(" failed.\n");
				throw new InvalidArgumentException.UnexistentTreeFileException(
						treefile.getAbsolutePath());
			} catch (TreeParseException e) {
				System.err.println("\nCannot parse tree file \""
						+ treefile.getName() + "\"");
				MAIN_CONSOLE.println(" failed.\n");
				throw new InvalidArgumentException.InvalidTreeFileException(
						treefile.getAbsolutePath());
			}
			if (tree != null) {
				options.setUserTree(TreeUtilities.toNewick(tree, true, false,
						false));
				TextOutputStream out = new TextOutputStream(options
						.getTreeFile().getAbsolutePath());
				out.print(options.getUserTree());
				out.close();
				MAIN_CONSOLE.println(" OK.");
			} else // tree is not valid
			{
				System.err.println("\nUnexpected error parsing \""
						+ treefile.getName() + "\"");
				MAIN_CONSOLE.println(" failed.\n");
				throw new InvalidArgumentException.InvalidTreeFileException(
						treefile.getAbsolutePath());
			}

		}
	}

	public static TextOutputStream setMainConsole(TextOutputStream mainConsole) {
		ModelTest.MAIN_CONSOLE = mainConsole;
		return mainConsole;
	}

	public static TextOutputStream getMainConsole() {
		return MAIN_CONSOLE;
	}

	public static TextOutputStream setPhymlConsole(TextOutputStream phymlConsole) {
		ModelTest.PHYML_CONSOLE = phymlConsole;
		return phymlConsole;
	}

	public static TextOutputStream getPhymlConsole() {
		return PHYML_CONSOLE;
	}
	
	public static TextOutputStream getCurrentOutStream() {
		return CURRENT_OUT_STREAM;
	}

	public static void setCurrentOutStream(TextOutputStream currentOutStream) {
		CURRENT_OUT_STREAM = currentOutStream;
	}

	public static AIC getMyAIC() {
		if (!AICwasCalculated)
			throw new WeakStateException.UninitializedCriterionException("AIC");
		return myAIC;
	}

	public static void setMyAIC(AIC myAIC) {
		ModelTest.myAIC = myAIC;
		ModelTest.minAIC = myAIC != null ? myAIC.getMinModel() : null;
		AICwasCalculated = (myAIC != null);
	}

	public static boolean testAIC() {
		return AICwasCalculated;
	}

	public static AICc getMyAICc() {
		if (!AICcwasCalculated)
			throw new WeakStateException.UninitializedCriterionException("AICc");
		return myAICc;
	}

	public static void setMyAICc(AICc myAICc) {
		ModelTest.myAICc = myAICc;
		ModelTest.minAICc = myAICc != null ? myAICc.getMinModel() : null;
		AICcwasCalculated = (myAICc != null);
	}

	public static boolean testAICc() {
		return AICcwasCalculated;
	}

	public static BIC getMyBIC() {
		if (!BICwasCalculated)
			throw new WeakStateException.UninitializedCriterionException("BIC");
		return myBIC;
	}

	public static void setMyBIC(BIC myBIC) {
		ModelTest.myBIC = myBIC;
		ModelTest.minBIC = myBIC != null ? myBIC.getMinModel() : null;
		BICwasCalculated = (myBIC != null);
	}

	public static boolean testBIC() {
		return BICwasCalculated;
	}

	public static DT getMyDT() {
		if (!DTwasCalculated)
			throw new WeakStateException.UninitializedCriterionException("DT");
		return myDT;
	}

	public static void setMyDT(DT myDT) {
		ModelTest.myDT = myDT;
		ModelTest.minDT = myDT != null ? myDT.getMinModel() : null;
		DTwasCalculated = (myDT != null);
	}

	public static boolean testDT() {
		return DTwasCalculated;
	}

	public static RunConsense getConsensusAIC() {
		return consensusAIC;
	}

	public static RunConsense getConsensusAICc() {
		return consensusAICc;
	}

	public static RunConsense getConsensusBIC() {
		return consensusBIC;
	}

	public static RunConsense getConsensusDT() {
		return consensusDT;
	}

	public static void setConsensusAIC(RunConsense pConsensusAIC) {
		consensusAIC = pConsensusAIC;
	}

	public static void setConsensusAICc(RunConsense pConsensusAICc) {
		consensusAICc = pConsensusAICc;
	}

	public static void setConsensusBIC(RunConsense pConsensusBIC) {
		consensusBIC = pConsensusBIC;
	}

	public static void setConsensusDT(RunConsense pConsensusDT) {
		consensusDT = pConsensusDT;
	}

	/**
	 * Finalizes the MPJ runtime environment. When an error occurs, it aborts
	 * the execution of every other processes.
	 * 
	 * @param status
	 *            the finalization status
	 */
	public static void finalize(int status) {

		if (status != 0) {
			if (MPJ_RUN) {
				MPI.COMM_WORLD.Abort(status);
			}
		}

		if (MPJ_RUN) {
			MPI.Finalize();
		}

		System.exit(status);

	}

	/**
	 * @param minDLRT
	 *            the minDLRT to set
	 */
	public static void setMinDLRT(Model minDLRT) {
		ModelTest.minDLRT = minDLRT;
	}

	/**
	 * @return the minDLRT
	 */
	public static Model getMinDLRT() {
		return minDLRT;
	}

	/**
	 * @param minHLRT
	 *            the minHLRT to set
	 */
	public static void setMinHLRT(Model minHLRT) {
		ModelTest.minHLRT = minHLRT;
	}

	/**
	 * @return the minHLRT
	 */
	public static Model getMinHLRT() {
		return minHLRT;
	}

	/**
	 * @return the minDT
	 */
	public static Model getMinDT() {
		return minDT;
	}

	/**
	 * @return the minBIC
	 */
	public static Model getMinBIC() {
		return minBIC;
	}

	/**
	 * @return the minAICc
	 */
	public static Model getMinAICc() {
		return minAICc;
	}

	/**
	 * @return the minAIC
	 */
	public static Model getMinAIC() {
		return minAIC;
	}

	/**
	 * @param candidateModels
	 *            the candidateModels to set
	 */
	public static void setCandidateModels(Model[] candidateModels) {
		ApplicationOptions.getInstance().setNumModels(candidateModels.length);
		ModelTest.candidateModels = candidateModels;
	}

	/**
	 * @return the loadedModels
	 */
	public static Model[] getLoadedModels() {
		return loadedModels;
	}
	
	/**
	 * @return the candidateModels
	 */
	public static Model[] getCandidateModels() {
		return candidateModels;
	}

	/**
	 * @return a single candidate model
	 */
	public static Model getCandidateModel(int index) {
		return candidateModels[index];
	}

	public static String getHostname() {
		return hostname;
	}

	public static void purgeModels() {
		List<Model> modelList = new ArrayList<Model>();
		for (Model model : candidateModels) {
			if (model.getLnL() > 0) {
				modelList.add(model);
			}
		}
		candidateModels = modelList.toArray(new Model[0]);
		ApplicationOptions.getInstance().setNumModels(candidateModels.length);
	}

	public class NullPrinter extends OutputStream {

		@Override
		public void write(int arg0) throws IOException {
			// DO NOTHING

		}

	}

} // class ModelTest public class ModelTestService {

	/**
     * Read alignment into a temporary file.
     * 
     * @param inputFile the file to be read
     * @param outputFile the file to write in
     * 
     * @throws AlignmentParseException the alignment parse exception.
     * @throws FileNotFoundException Signals that the input filename does not exist.
     * @throws IOException Signals that an I/O exception has occured.
     */
	public static String readAlignment(File inputFile, File outputFile)
            throws AlignmentParseException, FileNotFoundException, IOException {
		return readAlignment(inputFile, outputFile, true);
	}
	
	/**
     * Read alignment.
     * 
     * @param inputFile the file to be read
     * @param outputFile the file to write in
     * 
     * @throws AlignmentParseException the alignment parse exception.
     * @throws FileNotFoundException Signals that the input filename does not exist.
     * @throws IOException Signals that an I/O exception has occured.
     */
    public static String readAlignment(File inputFile, File outputFile, boolean deleteOnExit)
            throws AlignmentParseException, FileNotFoundException, IOException {

        StringBuilder text = new StringBuilder();

        BufferedReader br = new BufferedReader(
        		new FileReader(inputFile.getAbsolutePath()));
        String s;
        while ((s = br.readLine()) != null) {
            text.append(s).append("\r\n");
        }
        br.close();


        //Get options
        String in = text.toString();
        String inO = "";
        String inP = "";
        String inF = "";
        boolean autodetect = true;
        boolean collapse = false;
        boolean gaps = false;
        boolean missing = false;
        int limit = 0;
        String out = "";
        String outO = "Linux";
        switch (Utilities.findCurrentOS()) {
        case Utilities.OS_WINDOWS:
            outO = "Windows";
	        break;
        case Utilities.OS_OSX:
        case Utilities.OS_LINUX:
            outO = "Linux";
            break;
        }
        String outP = "PhyML";
        String outF = "PHYLIP";
        boolean lower = false;
        boolean numbers = false;
        boolean sequential = false;
        boolean match = false;

        //Get converter and convert MSA
        Factory factory = new DefaultFactory();
        Converter converter;

		PrintStream outPS = System.err;
		System.setErr(new PrintStream(new NullPrintStream()));
        try {
            converter = factory.getConverter(inO, inP, inF, autodetect,
                    collapse, gaps, missing, limit,
                    outO, outP, outF, lower, numbers, sequential, match, "ModelTestService");
            out = converter.convert(in);
        } catch (UnsupportedOperationException ex) {
            throw new AlignmentParseException("There's some error in your data: " + ex.getMessage());
        } catch (ParseException ex) {
            throw new AlignmentParseException("There's some error in your data: " + ex.getMessage());
        }
        System.setErr(outPS);
        
        if (outputFile.exists()) {
        	outputFile.delete();
        }
        outputFile.createNewFile();
        if (deleteOnExit) {
        	outputFile.deleteOnExit();
        }
        	
        FileWriter fw = new FileWriter(outputFile);
        fw.write(out);
        fw.close();
        
        return out;
    }
    
	public InformationCriterion doIC(int ic, boolean writePAUPblock, boolean doImportances,
			boolean doModelAveraging, double confidenceInterval) {
		InformationCriterion criterion;
		switch (ic) {
		case InformationCriterion.IC_AIC:
			criterion = new AIC(writePAUPblock, doImportances,
					doModelAveraging, confidenceInterval);
			break;
		case InformationCriterion.IC_AICc:
			criterion = new AICc(writePAUPblock, doImportances,
					doModelAveraging, confidenceInterval);
			break;
		case InformationCriterion.IC_BIC:
			criterion = new BIC(writePAUPblock, doImportances,
					doModelAveraging, confidenceInterval);
			break;
		case InformationCriterion.IC_DT:
			criterion = new DT(writePAUPblock, doImportances, doModelAveraging,
					confidenceInterval);
			break;
		default:
			throw new InvalidArgumentException.UnexistentCriterionException(ic);
		}
		criterion.compute();
		
		return criterion;
	}
	
} 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 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 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 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
import java.io.File;
import java.io.IOException;
import java.io.PrintWriter;
import java.io.PushbackReader;
import java.io.StringReader;
import java.util.Arrays;
import java.util.Locale;
import pal.tree.Tree;
import pal.tree.TreeParseException;
import es.uvigo.darwin.jmodeltest.ApplicationOptions;
import es.uvigo.darwin.jmodeltest.ModelTest;
import es.uvigo.darwin.jmodeltest.ModelTestService;
import es.uvigo.darwin.jmodeltest.exe.RunConsense;
import es.uvigo.darwin.jmodeltest.exe.RunPhyml;
import es.uvigo.darwin.jmodeltest.exe.RunPhymlThread;
import es.uvigo.darwin.jmodeltest.io.AlignmentReader;
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.TreeUtilities;
			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", 