/*
* DiscreteRatePriorGenerator.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package dr.app.phylogeography.tools;
import dr.app.tools.TimeSlicer;
import dr.app.util.Arguments;
import dr.geo.math.SphericalPolarCoordinates;
import dr.stats.DiscreteStatistics;
import org.jdom.Element;
import org.jdom.output.Format;
import org.jdom.output.XMLOutputter;
import java.io.*;
import java.util.ArrayList;
import java.util.List;
import java.util.StringTokenizer;
/**
* @author Philippe Lemey
* @author Andrew Rambaut
* @author Marc A. Suchard
*/
public class DiscreteRatePriorGenerator {
public static final String HELP = "help";
public static final String COORDINATES = "coordinates";
public static final String DENSITIES = "densities";
public static final String GENERICS = "generics";
public static final String FORMAT = "format";
public static final String MODEL = "model";
public DiscreteRatePriorGenerator(String[] locations, Double[] latitudes, Double[] longitudes, Double[] densities) {
if (locations == null) {
//System.out.println(locations[0]);
System.err.println("no locations specified!");
} else {
this.locations = locations;
}
this.latitudes = latitudes;
this.longitudes = longitudes;
this.densities = densities;
if ((latitudes == null)||(longitudes == null)) {
progressStream.println("no latitudes or longitudes specified!");
} else {
distances = getUpperTriangleDistanceMatrix(latitudes,longitudes);
}
if (densities != null) {
densityDonorMatrix = getDensityMatrix(densities, true);
densityRecipientMatrix = getDensityMatrix(densities, false);
}
}
// for the time being, locations, latitudes and longitudes are not required
private String[] locations;
private Double[] latitudes;
private Double[] longitudes;
private final Double[] densities;
private double[] distances;
private double[] densityDonorMatrix;
private double[] densityRecipientMatrix;
private double[] getUpperTriangleDistanceMatrix(Double[] latitudes, Double[] longitudes) {
double[] distances = new double[(latitudes.length * (latitudes.length - 1))/2];
int distanceCounter = 0;
int pairwiseCounter1 = 0;
for (int c = 0; c < latitudes.length; c++) {
pairwiseCounter1 ++;
for (int d = pairwiseCounter1; d < latitudes.length; d++) {
distances[distanceCounter] = getKilometerGreatCircleDistance(latitudes[c],longitudes[c],latitudes[d],longitudes[d]);
distanceCounter++;
}
}
return distances;
}
private double[] getFullDistanceMatrix(Double[] latitudes, Double[] longitudes) {
double[] distances = new double[latitudes.length * latitudes.length];
int distanceCounter = 0;
for (int a = 0; a < latitudes.length; a++) {
//resultsStream.print(locations[a]+"\t");
for (int b = 0; b < latitudes.length; b++) {
distances[distanceCounter] = getKilometerGreatCircleDistance(latitudes[a],longitudes[a],latitudes[b],longitudes[b]);
distanceCounter++;
}
}
return distances;
}
private double[] getFullCoordDiffMatrix(Double[] coordinates, boolean positive, boolean negative) {
double[] differences = new double[coordinates.length * coordinates.length];
int differenceCounter = 0;
for (int a = 0; a < coordinates.length; a++) {
for (int b = 0; b < coordinates.length; b++) {
double difference = coordinates[b]-coordinates[a];
if (difference > 0 && positive) {
differences[differenceCounter] = difference;
} else if (difference < 0 && negative) {
differences[differenceCounter] = difference;
} else {
differences[differenceCounter] = 0;
}
differenceCounter++;
}
}
return differences;
}
private static double getKilometerGreatCircleDistance(double lat1, double long1, double lat2, double long2) {
SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(lat1, long1);
SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(lat2, long2);
return (coord1.distance(coord2));
}
private double[] getDensityMatrix(Double[] densities, boolean donor) {
double[] returnMatrix = new double[densities.length * (densities.length - 1)];
int distanceCounter = 0;
int matrixEntry;
for (int c = 0; c < densities.length; c++) {
for (int d = 0; d < densities.length; d++) {
if (c == d) {
continue;
}
if (donor) {
matrixEntry = c;
} else {
matrixEntry = d;
}
returnMatrix[distanceCounter] = densities[matrixEntry];
distanceCounter++;
}
}
return returnMatrix;
}
private void printLocations(String[] locs, boolean upper) {
int pairwiseCounter1 = 0;
for (int c = 0; c < locs.length; c++) {
pairwiseCounter1 ++;
for (int d = pairwiseCounter1; d < locs.length; d++) {
if (upper) {
System.out.println(locs[c]+"\t"+locs[d]);
} else {
System.out.println(locs[d]+"\t"+locs[c]);
}
}
}
}
private void printLocationsDensitiesDistances (String[] locs, boolean upper) {
double[] firstDensity = getSingleElementFromFullMatrix(densities,false);
double[] secondDensity = getSingleElementFromFullMatrix(densities,true);
System.out.println("location1\tlocation2\tdistance\tdensity1\tdensity2");
int pairwiseCounter1 = 0;
int arrayCounter = 0;
for (int c = 0; c < locs.length; c++) {
pairwiseCounter1 ++;
for (int d = pairwiseCounter1; d < locs.length; d++) {
if (upper) {
System.out.println(locs[c]+"\t"+locs[d]+"\t"+distances[arrayCounter]+"\t"+firstDensity[arrayCounter]+"\t"+secondDensity[arrayCounter]);
arrayCounter ++;
} else {
System.out.println(locs[d]+"\t"+locs[c]+"\t"+distances[arrayCounter]+"\t"+firstDensity[(firstDensity.length/2)+arrayCounter]+"\t"+secondDensity[(secondDensity.length/2)+arrayCounter]);
arrayCounter ++;
}
}
}
}
private void printFullDistanceMatrix(String name, boolean locationNames) {
try {
PrintWriter outFile = new PrintWriter(new FileWriter(name), true);
outFile.print("location");
for (int a = 0; a < locations.length; a++) {
outFile.print(locations[a]+"\t");
for (int b = 0; b < locations.length; b++) {
double lat1 = latitudes[a];
double lat2 = latitudes[b];
double long1 = longitudes[a];
double long2 = longitudes[b];
double distance = getKilometerGreatCircleDistance(lat1,long1,lat2,long2);
outFile.print(distance+"\t");
}
outFile.print("\r");
}
outFile.close();
} catch(IOException io) {
System.err.print("Error writing to file: " + name);
}
}
enum OutputFormat {
TAB,
SPACE,
XML
}
enum Model {
PRIOR,
GLM,
JUMP
}
public void output(String outputFileName, OutputFormat outputFormat, Model model) {
//printLocations(locations,false);
//printLocationsDensitiesDistances(locations,true);
//printLocationsDensitiesDistances(locations,false);
//printFullDistanceMatrix("test.txt", false);
resultsStream = System.out;
if (outputFileName != null) {
try {
resultsStream = new PrintStream(new File(outputFileName));
} catch (IOException e) {
System.err.println("Error opening file: "+outputFileName);
System.exit(1);
}
}
if (model == model.PRIOR) {
int predictor = 0;
if (distances != null) {
outputStringLine("reversible priors: distances",outputFormat);
resultsStream.print("\r");
outputArray(distances, outputFormat, model, predictor, false);
resultsStream.print("\r");
outputStringLine("reversible priors: normalized inverse distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(distances,true, true, false, false), outputFormat, model, predictor, false);
resultsStream.print("\r");
outputStringLine("reversible priors: normalized inverse log distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(distances,true, true, false, true), outputFormat, model, predictor, false);
resultsStream.print("\r");
}
if (densities != null) {
outputStringLine("reversible priors: densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getUpperTrianglePairwiseProductMatrix(densities),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
outputStringLine("reversible priors: normalized densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getUpperTrianglePairwiseProductMatrix(densities),false,true,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
outputStringLine("reversible priors: normalized log densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getUpperTrianglePairwiseProductMatrix(densities),false,true,false,true), outputFormat, model, predictor, false);
resultsStream.print("\r");
if (distances != null) {
outputStringLine("reversible priors: product of normalized densities divided by normalized distances",outputFormat);
resultsStream.print("\r");
outputArray(productOfArrays(transform(getUpperTrianglePairwiseProductMatrix(densities),false,true,false,false),transform(distances,false, true,false,false),true), outputFormat, model, predictor, false);
resultsStream.print("\r");
outputStringLine("reversible priors: normalized product of densities divided by distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(productOfArrays(getUpperTrianglePairwiseProductMatrix(densities),distances, true),false,true,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
outputStringLine("reversible priors: normalized log product of densities divided by distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(productOfArrays(getUpperTrianglePairwiseProductMatrix(densities),distances, true),false,true,false,true), outputFormat, model, predictor, false);
resultsStream.print("\r");
}
}
if (distances != null) {
outputStringLine("nonreversible priors: distances",outputFormat);
resultsStream.print("\r");
outputArray(distances, outputFormat, model, predictor, true);
resultsStream.print("\r");
outputStringLine("nonreversible priors: normalized inverse distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(distances,true, true, false, false), outputFormat, model, predictor, true);
resultsStream.print("\r");
outputStringLine("nonreversible priors: normalized inverse log distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(distances,true, true, false, true), outputFormat, model,predictor, true);
resultsStream.print("\r");
}
if (densities != null) {
outputStringLine("nonreversible priors: densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getUpperTrianglePairwiseProductMatrix(densities),false,false,false,false), outputFormat, model, predictor, true);
resultsStream.print("\r");
outputStringLine("nonreversible priors: normalized densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getUpperTrianglePairwiseProductMatrix(densities),false,true,false,false), outputFormat, model, predictor, true);
resultsStream.print("\r");
outputStringLine("nonreversible priors: normalized log densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getUpperTrianglePairwiseProductMatrix(densities),false,true,false,true), outputFormat, model, predictor, true);
resultsStream.print("\r");
if (distances != null) {
outputStringLine("nonreversible priors: product of normalized densities divided by normalized distances",outputFormat);
resultsStream.print("\r");
outputArray(productOfArrays(transform(getUpperTrianglePairwiseProductMatrix(densities),false,true,false,false),transform(distances,false, true, false, false),true), outputFormat, model, predictor, true);
resultsStream.print("\r");
outputStringLine("nonreversible priors: normalized product of densities divided by distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(productOfArrays(getUpperTrianglePairwiseProductMatrix(densities),distances, true),false,true,false,false), outputFormat, model, predictor, true);
resultsStream.print("\r");
outputStringLine("nonreversible priors: normalized log product of densities divided by distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(productOfArrays(getUpperTrianglePairwiseProductMatrix(densities),distances, true),false,true,false,true), outputFormat, model, predictor, true);
resultsStream.print("\r");
}
}
} else if (model == model.GLM) {
int predictor = 1;
//TODO: fully implement glm output
if (distances != null) {
//printLocations(locations, true);
//printLocations(locations, false);
outputStringLine("predictor: distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(distances,false, false, true, false), outputFormat, model, predictor, true);
predictor++;
outputStringLine("predictor: standardized log distances",outputFormat);
resultsStream.print("\r");
outputArray(transform(distances,false, false, true, true), outputFormat, model, predictor, true);
}
if (densities != null) {
predictor++;
outputStringLine("predictor: standardized donor densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getSingleElementFromFullMatrix(densities, false),false,false,true,true), outputFormat, model, predictor, false);
predictor++;
outputStringLine("predictor: standardized recipient densities",outputFormat);
resultsStream.print("\r");
outputArray(transform(getSingleElementFromFullMatrix(densities, true),false,false,true,true), outputFormat, model, predictor, false);
// products not necessary, and shouldn't be normalized but standardized, also log
// outputStringLine("predictor: normalized product of densities divided by distances",outputFormat);
// resultsStream.print("\r");
// outputArray(transform(productOfArrays(getUpperTrianglePairwiseProductMatrix(densities),distances, true),false,true,false,false), outputFormat, predictor, true);
// resultsStream.print("\r");
// predictor++;
// outputStringLine("predictor: normalized log product of densities divided by distances",outputFormat);
// resultsStream.print("\r");
// outputArray(transform(productOfArrays(getUpperTrianglePairwiseProductMatrix(densities),distances, true),false,true,false,true), outputFormat, predictor, true);
// resultsStream.print("\r");
}
} else if (model == model.JUMP) {
int predictor = 0;
outputStringLine("great circle distance jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullDistanceMatrix(latitudes, longitudes),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
outputStringLine("latitude jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullCoordDiffMatrix(latitudes, true, true),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
outputStringLine("longitude jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullCoordDiffMatrix(longitudes, true, true),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
outputStringLine("westward jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullCoordDiffMatrix(longitudes, false, true),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
outputStringLine("eastward jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullCoordDiffMatrix(longitudes, true, false),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
outputStringLine("northward jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullCoordDiffMatrix(latitudes, true, false),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
outputStringLine("southward jump matrix",outputFormat);
resultsStream.print("\r");
outputArray(transform(getFullCoordDiffMatrix(latitudes, false, true),false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
// outputStringLine("latitude reward parameter",outputFormat);
// resultsStream.print("\r");
// outputArray(transform(objectToPrimitiveArray(latitudes),false,false,false,false), outputFormat, model, predictor, false);
// resultsStream.print("\r");
// predictor ++;
// outputStringLine("longitude reward parameter",outputFormat);
// resultsStream.print("\r");
// outputArray(transform(objectToPrimitiveArray(longitudes),false,false,false,false), outputFormat, model, predictor, false);
// resultsStream.print("\r");
// predictor ++;
for (int i = 0; i < locations.length; i++) {
double[] locationReward = getLocationReward(i,locations.length);
outputStringLine(locations[i]+" reward parameter",outputFormat);
resultsStream.print("\r");
outputArray(transform(locationReward,false,false,false,false), outputFormat, model, predictor, false);
resultsStream.print("\r");
predictor ++;
}
}
}
private double[] getUpperTrianglePairwiseProductMatrix(Double[] matrixValues) {
double[] pairwiseProducts = new double[(matrixValues.length * (matrixValues.length - 1))/2];
int counter = 0;
int pairwiseCounter1 = 0;
for(Double matrixValue : matrixValues) {
pairwiseCounter1++;
for(int d = pairwiseCounter1; d < matrixValues.length; d++) {
pairwiseProducts[counter] = matrixValue * matrixValues[d];
counter++;
}
}
return pairwiseProducts;
}
private double[] getSingleElementFromFullMatrix(Double[] matrixValues, boolean secondElement) {
double[] singleElements = new double[(matrixValues.length * (matrixValues.length - 1))];
//System.out.println("matrixsize "+(matrixValues.length * (matrixValues.length - 1)));
//get upper matrix
int counter = 0;
int pairwiseCounter1 = 0;
for(Double matrixValue : matrixValues) {
pairwiseCounter1++;
for(int d = pairwiseCounter1; d < matrixValues.length; d++) {
if (secondElement) {
singleElements[counter] = matrixValues[d];
} else {
singleElements[counter] = matrixValue;
}
counter++;
}
}
//get lower matrix
int pairwiseCounter2 = 0;
for(Double matrixValue : matrixValues) {
pairwiseCounter2++;
for(int d = pairwiseCounter2; d < matrixValues.length; d++) {
if (secondElement) {
singleElements[counter] = matrixValue;
} else {
singleElements[counter] = matrixValues[d];
}
counter++;
}
}
//System.out.println("counter "+counter);
return singleElements;
}
private double[] productOfArrays(double[] inputArray1, double[] inputArray2, boolean devision) {
if (inputArray1.length != inputArray2.length) {
System.err.println("trying to get a product of arrays of unequals size!");
System.exit(1);
}
double[] productOfArray = new double[inputArray1.length];
for (int i = 0; i < inputArray1.length; i++) {
if (devision) {
productOfArray[i] = inputArray1[i]/inputArray2[i];
} else {
}
productOfArray[i] = inputArray1[i]*inputArray2[i];
}
return productOfArray;
}
// normalize is really rescaling the vector to have a mean = 1 here, standardize is rescaling it to have mean = 0 and variance =1
private double[] transform(double[] inputArray, boolean inverse, boolean normalize, boolean standardize, boolean log) {
double[] transformedDistances = new double[inputArray.length];
for(int u=0; u<inputArray.length; u++) {
double distance = inputArray[u];
if (log) {
distance = Math.log(distance);
}
if (inverse) {
distance = 1/distance;
}
transformedDistances[u] = distance;
}
double meanDistance = 1;
double stdev = 0;
if (normalize || standardize) {
meanDistance = DiscreteStatistics.mean(transformedDistances);
}
if (standardize) {
stdev = Math.sqrt(DiscreteStatistics.variance(transformedDistances));
}
for(int v=0; v<inputArray.length; v++) {
if (normalize) {
transformedDistances[v] = (transformedDistances[v]/meanDistance);
} else if (standardize) {
transformedDistances[v] = ((transformedDistances[v] - meanDistance)/stdev);
}
}
return transformedDistances;
}
private void outputArray(double[] array, OutputFormat outputFormat, Model model, int predictor, boolean nonreversible) {
StringBuilder sb1 = new StringBuilder();
StringBuilder sb2 = new StringBuilder();
String sep;
if (outputFormat == OutputFormat.TAB ) {
sep = "\t";
} else {
sep = " ";
}
//String newLine = "\n";
//here we break up the jump parameter output in a lenght*by*length matrix
// int entryCounter = 1;
// double length = Math.sqrt(array.length);
// System.out.println(length);
// for(double anArray : array) {
// if (model == model.JUMP) {
// if ( (entryCounter % length) > 0 ) {
// sb1.append(anArray + sep);
// sb2.append(1 + sep);
// entryCounter ++;
// } else {
// System.out.println("return");
// sb1.append(newLine+anArray + sep);
// sb2.append(newLine+1 + sep);
// entryCounter ++;
// }
//
// } else {
// sb1.append(anArray + sep);
// sb2.append(1 + sep);
// }
// }
for(double anArray : array) {
sb1.append(anArray + sep);
sb2.append(1 + sep);
}
if (outputFormat == OutputFormat.XML ) {
if (model == model.GLM) {
Element parameter = new Element("parameter");
parameter.setAttribute("id","predictor"+predictor);
if (nonreversible) {
parameter.setAttribute("value",(sb1.toString()+sb1.toString()));
} else {
parameter.setAttribute("value",sb1.toString());
}
XMLOutputter xmlOutputter = new XMLOutputter(Format.getPrettyFormat().setTextMode(Format.TextMode.PRESERVE));
try {
xmlOutputter.output(parameter,resultsStream);
} catch (IOException e) {
System.err.println("IO Exception encountered: "+e.getMessage());
System.exit(-1);
}
resultsStream.print("\r");
} else if (model == model.PRIOR) {
Element priorElement = new Element("multivariateGammaPrior");
Element data = new Element("data");
Element parameter1 = new Element("parameter");
parameter1.setAttribute("idref","rates");
data.addContent(parameter1);
Element meanParameter = new Element("meanParameter");
Element parameter2 = new Element("parameter");
if (nonreversible) {
parameter2.setAttribute("value",(sb1.toString()+sb1.toString()));
} else {
parameter2.setAttribute("value",sb1.toString());
}
meanParameter.addContent(parameter2);
Element coefficientOfVariation = new Element("coefficientOfVariation");
Element parameter3 = new Element("parameter");
if (nonreversible) {
parameter3.setAttribute("value",sb2.toString()+sb2.toString());
} else {
parameter3.setAttribute("value",sb2.toString());
}
coefficientOfVariation.addContent(parameter3);
priorElement.addContent(data);
priorElement.addContent(meanParameter);
priorElement.addContent(coefficientOfVariation);
XMLOutputter xmlOutputter = new XMLOutputter(Format.getPrettyFormat().setTextMode(Format.TextMode.PRESERVE));
try {
xmlOutputter.output(priorElement,resultsStream);
} catch (IOException e) {
System.err.println("IO Exception encountered: "+e.getMessage());
System.exit(-1);
}
resultsStream.print("\r");
} else if (model == model.JUMP) {
Element parameter = new Element("parameter");
parameter.setAttribute("id","jump"+predictor);
parameter.setAttribute("value",sb1.toString());
XMLOutputter xmlOutputter = new XMLOutputter(Format.getPrettyFormat().setTextMode(Format.TextMode.PRESERVE));
try {
xmlOutputter.output(parameter,resultsStream);
} catch (IOException e) {
System.err.println("IO Exception encountered: "+e.getMessage());
System.exit(-1);
}
resultsStream.print("\r");
}
} else {
if (nonreversible) {
resultsStream.print(sb1.toString()+sb1.toString()+"\r");
} else {
resultsStream.print(sb1+"\r");
}
}
}
private void outputStringLine(String outputString, OutputFormat outputFormat) {
if (outputFormat == OutputFormat.XML ) {
resultsStream.print("<!-- ");
}
resultsStream.print(outputString);
if (outputFormat == OutputFormat.XML ) {
resultsStream.print(" -->");
} else {
resultsStream.print("\r");
}
}
private double[] objectToPrimitiveArray(Double[] array) {
double[] returnArray = new double[array.length];
int counter = 0;
for(Double anArray : array) {
returnArray[counter] = anArray.doubleValue();
counter ++;
}
return returnArray;
}
private double[] getLocationReward(int indicator, int length) {
double[] returnArray = new double[length];
for (int i = 0; i < length; i++) {
if (i == indicator) {
returnArray[i] = 1.0;
} else {
returnArray[i] = 0.0;
}
}
return returnArray;
}
// Messages to stderr, output to stdout
private static final PrintStream progressStream = System.err;
private PrintStream resultsStream;
private static final String commandName = "discreteRatePriorGenerator";
public static void printUsage(Arguments arguments) {
arguments.printUsage(commandName, "[<output-file-name>]");
progressStream.println();
progressStream.println(" Example: " + commandName + " coordinates.txt ratePriors.txt");
progressStream.println();
}
private static ArrayList parseCoordinatesFile(String inputFile, String[] locations, Double[] latitudes, Double[] longitudes) {
ArrayList<Object[]> returnList = new ArrayList<Object[]>();
List<String> countList = new ArrayList<String>();
try{
BufferedReader reader = new BufferedReader(new FileReader(inputFile));
String line = reader.readLine();
while (line != null && !line.equals("")) {
countList.add(line);
line = reader.readLine();
}
} catch (IOException e) {
System.err.println("Error reading " + inputFile);
System.exit(1);
}
if (countList.size()>0) {
locations = new String[countList.size()];
latitudes = new Double[countList.size()];
longitudes = new Double[countList.size()];
for(int i=0; i<countList.size(); i++) {
StringTokenizer tokens = new StringTokenizer(countList.get(i));
locations[i] = tokens.nextToken("\t");
latitudes[i] = Double.parseDouble(tokens.nextToken("\t"));
longitudes[i] = Double.parseDouble(tokens.nextToken("\t"));
//System.out.println(locations[i]+"\t"+latitudes[i]+"\t"+longitudes[i]);
}
}
returnList.add(locations);
returnList.add(latitudes);
returnList.add(longitudes);
return returnList;
}
private static ArrayList parseSingleMeasureFile(String inputFile, String[] locations, Double[] densities) {
ArrayList<Object[]> returnList = new ArrayList<Object[]>();
boolean locationsSpecified = true;
if (locations == null) {
locationsSpecified = false;
}
List<String> countList = new ArrayList<String>();
try{
BufferedReader reader = new BufferedReader(new FileReader(inputFile));
String line = reader.readLine();
while (line != null && !line.equals("")) {
countList.add(line);
line = reader.readLine();
}
} catch (IOException e) {
System.err.println("Error reading " + inputFile);
System.exit(1);
}
if (countList.size()>0) {
if (!locationsSpecified) {
locations = new String[countList.size()];
}
densities = new Double[countList.size()];
for(int i=0; i<countList.size(); i++) {
StringTokenizer tokens = new StringTokenizer(countList.get(i));
String location = tokens.nextToken("\t");
if (!locationsSpecified) {
locations[i] = location;
} else {
if (!(locations[i].equals(location))) {
System.err.println("Error in location specification in different files: " + locations[i] + "is not = " + location);
}
}
densities[i] = Double.parseDouble(tokens.nextToken("\t"));
}
}
returnList.add(locations);
returnList.add(densities);
return returnList;
}
public static void main(String[] args) throws IOException {
String outputFileName = null;
String[] locations = null;
Double[] latitudes = null;
Double[] longitudes = null;
Double[] densities = null;
OutputFormat outputFormat = OutputFormat.XML;
Model model = Model.PRIOR;
Arguments arguments = new Arguments(
new Arguments.Option[]{
new Arguments.StringOption(COORDINATES, "coordinate file", "specifies a tab-delimited file with coordinates for the locations"),
new Arguments.StringOption(DENSITIES, "density file", "specifies a tab-delimited file with densities for the locations"),
new Arguments.StringOption(GENERICS, "generics file", "specifies a tab-delimited file-list to use as measures for the locations"),
new Arguments.StringOption(FORMAT, TimeSlicer.enumNamesToStringArray(OutputFormat.values()),false,
"prior output format [default = XML]"),
new Arguments.StringOption(MODEL, TimeSlicer.enumNamesToStringArray(Model.values()),false,
"model output [default = rate priors]"),
//example: new Arguments.RealOption(MRSD,"specifies the most recent sampling data in fractional years to rescale time [default=0]"),
new Arguments.Option(HELP, "option to print this message"),
});
try {
arguments.parseArguments(args);
} catch (Arguments.ArgumentException ae) {
progressStream.println(ae);
printUsage(arguments);
System.exit(1);
}
if (arguments.hasOption(HELP)) {
printUsage(arguments);
System.exit(0);
}
String coordinatesFileString = arguments.getStringOption(COORDINATES);
if (coordinatesFileString != null) {
ArrayList LocsLatsLongs = parseCoordinatesFile(coordinatesFileString,locations,latitudes,longitudes);
locations = (String[]) LocsLatsLongs.get(0);
latitudes = (Double[]) LocsLatsLongs.get(1);
longitudes = (Double[]) LocsLatsLongs.get(2);
}
String densitiesFileString = arguments.getStringOption(DENSITIES);
if (densitiesFileString != null) {
ArrayList LocsDens = parseSingleMeasureFile(densitiesFileString,locations,densities);
locations = (String[]) LocsDens.get(0);
densities = (Double[]) LocsDens.get(1);
}
//TODO: support reading any measure (GENERICS)
String summaryFormat = arguments.getStringOption(FORMAT);
if (summaryFormat != null) {
outputFormat = OutputFormat.valueOf(summaryFormat.toUpperCase());
}
String modelComponent = arguments.getStringOption(MODEL);
if (modelComponent != null) {
model = Model.valueOf(modelComponent.toUpperCase());
}
final String[] args2 = arguments.getLeftoverArguments();
switch (args2.length) {
case 0:
printUsage(arguments);
System.exit(1);
case 1:
outputFileName = args2[0];
break;
default: {
System.err.println("Unknown option: " + args2[1]);
System.err.println();
printUsage(arguments);
System.exit(1);
}
}
DiscreteRatePriorGenerator rates = new DiscreteRatePriorGenerator(locations, latitudes, longitudes, densities);
rates.output(outputFileName, outputFormat, model);
System.exit(0);
}
private static void printArray(Double[] array, String name) {
try {
PrintWriter outFile = new PrintWriter(new FileWriter(name), true);
for(Double anArray : array) {
outFile.print(anArray + "\t");
}
outFile.close();
} catch(IOException io) {
System.err.print("Error writing to file: " + name);
}
}
private static void printArray(String[] array, String name) {
try {
PrintWriter outFile = new PrintWriter(new FileWriter(name), true);
for(String anArray : array) {
outFile.print(anArray + "\t");
}
outFile.close();
} catch(IOException io) {
System.err.print("Error writing to file: " + name);
}
}
}