Examples with SphericalPolarCoordinates dr.geo.math.SphericalPolarCoordinates used on opensource projects

Example 1 with SphericalPolarCoordinates

use of dr.geo.math.SphericalPolarCoordinates in project beast-mcmc by beast-dev.

the class DiscreteRatePriorGenerator method getKilometerGreatCircleDistance.

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));
}

Example 2 with SphericalPolarCoordinates

use of dr.geo.math.SphericalPolarCoordinates in project beast-mcmc by beast-dev.

the class DiffusionRateStatistic method getStatisticValue.

public double getStatisticValue(int dim) {
    String traitName = traitLikelihoods.get(0).getTraitName();
    double treelength = 0;
    double treeDistance = 0;
    double maxDistanceFromRoot = 0;
    double maxDistanceOverTimeFromRoot = 0;
    //double[] rates =  null;
    List<Double> rates = new ArrayList<Double>();
    //double[] diffusionCoefficients =  null;
    List<Double> diffusionCoefficients = new ArrayList<Double>();
    double waDiffusionCoefficient = 0;
    double lowerHeight = heightLowers[dim];
    double upperHeight = Double.MAX_VALUE;
    if (heightLowers.length == 1) {
        upperHeight = heightUpper;
    } else {
        if (dim > 0) {
            if (!cumulative) {
                upperHeight = heightLowers[dim - 1];
            }
        }
    }
    for (AbstractMultivariateTraitLikelihood traitLikelihood : traitLikelihoods) {
        MultivariateTraitTree tree = traitLikelihood.getTreeModel();
        BranchRateModel branchRates = traitLikelihood.getBranchRateModel();
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef node = tree.getNode(i);
            if (node != tree.getRoot()) {
                NodeRef parentNode = tree.getParent(node);
                if ((tree.getNodeHeight(parentNode) > lowerHeight) && (tree.getNodeHeight(node) < upperHeight)) {
                    double[] trait = traitLikelihood.getTraitForNode(tree, node, traitName);
                    double[] parentTrait = traitLikelihood.getTraitForNode(tree, parentNode, traitName);
                    double[] traitUp = parentTrait;
                    double[] traitLow = trait;
                    double timeUp = tree.getNodeHeight(parentNode);
                    double timeLow = tree.getNodeHeight(node);
                    double rate = (branchRates != null ? branchRates.getBranchRate(tree, node) : 1.0);
                    MultivariateDiffusionModel diffModel = traitLikelihood.diffusionModel;
                    double[] precision = diffModel.getPrecisionParameter().getParameterValues();
                    if (tree.getNodeHeight(parentNode) > upperHeight) {
                        timeUp = upperHeight;
                        //TODO: implement TrueNoise??
                        traitUp = imputeValue(trait, parentTrait, upperHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, false);
                    }
                    if (tree.getNodeHeight(node) < lowerHeight) {
                        timeLow = lowerHeight;
                        traitLow = imputeValue(trait, parentTrait, lowerHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, false);
                    }
                    double time = timeUp - timeLow;
                    treelength += time;
                    double[] rootTrait = traitLikelihood.getTraitForNode(tree, tree.getRoot(), traitName);
                    if (useGreatCircleDistances && (trait.length == 2)) {
                        // Great Circle distance
                        SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(traitLow[0], traitLow[1]);
                        SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(traitUp[0], traitUp[1]);
                        double distance = coord1.distance(coord2);
                        treeDistance += distance;
                        double dc = Math.pow(distance, 2) / (4 * time);
                        diffusionCoefficients.add(dc);
                        waDiffusionCoefficient += dc * time;
                        rates.add(distance / time);
                        SphericalPolarCoordinates rootCoord = new SphericalPolarCoordinates(rootTrait[0], rootTrait[1]);
                        double tempDistanceFromRoot = rootCoord.distance(coord2);
                        if (tempDistanceFromRoot > maxDistanceFromRoot) {
                            maxDistanceFromRoot = tempDistanceFromRoot;
                            maxDistanceOverTimeFromRoot = tempDistanceFromRoot / (tree.getNodeHeight(tree.getRoot()) - timeLow);
                            //distance between traitLow and traitUp for maxDistanceFromRoot
                            if (timeUp == upperHeight) {
                                maxDistanceFromRoot = distance;
                                maxDistanceOverTimeFromRoot = distance / time;
                            }
                        }
                    } else {
                        double distance = getNativeDistance(traitLow, traitUp);
                        treeDistance += distance;
                        double dc = Math.pow(distance, 2) / (4 * time);
                        diffusionCoefficients.add(dc);
                        waDiffusionCoefficient += dc * time;
                        rates.add(distance / time);
                        double tempDistanceFromRoot = getNativeDistance(traitLow, rootTrait);
                        if (tempDistanceFromRoot > maxDistanceFromRoot) {
                            maxDistanceFromRoot = tempDistanceFromRoot;
                            maxDistanceOverTimeFromRoot = tempDistanceFromRoot / (tree.getNodeHeight(tree.getRoot()) - timeLow);
                            //distance between traitLow and traitUp for maxDistanceFromRoot
                            if (timeUp == upperHeight) {
                                maxDistanceFromRoot = distance;
                                maxDistanceOverTimeFromRoot = distance / time;
                            }
                        }
                    }
                }
            }
        }
    }
    if (summaryStat == summaryStatistic.DIFFUSION_RATE) {
        if (summaryMode == Mode.AVERAGE) {
            return DiscreteStatistics.mean(toArray(rates));
        } else if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(rates));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(rates));
            return Math.sqrt(DiscreteStatistics.variance(toArray(rates), mean)) / mean;
        } else {
            return treeDistance / treelength;
        }
    } else if (summaryStat == summaryStatistic.DIFFUSION_COEFFICIENT) {
        if (summaryMode == Mode.AVERAGE) {
            return DiscreteStatistics.mean(toArray(diffusionCoefficients));
        } else if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(diffusionCoefficients));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(diffusionCoefficients));
            return Math.sqrt(DiscreteStatistics.variance(toArray(diffusionCoefficients), mean)) / mean;
        } else {
            return waDiffusionCoefficient / treelength;
        }
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE) {
        return maxDistanceFromRoot;
    } else {
        return maxDistanceOverTimeFromRoot;
    }
}

Example 3 with SphericalPolarCoordinates

use of dr.geo.math.SphericalPolarCoordinates in project beast-mcmc by beast-dev.

the class GreatCircleDiffusionModel method calculateLogDensity.

protected double calculateLogDensity(double[] start, double[] stop, double time) {
    SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(start[0], start[1]);
    SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(stop[0], stop[1]);
    double distance = coord1.distance(coord2);
    double inverseVariance = precision.getParameterValue(0) / time;
    // in the normalization constant
    if (coefficient == null)
        return -LOG2PI + Math.log(inverseVariance) - 0.5 * (distance * distance * inverseVariance);
    double coef = -coefficient.getParameterValue(0);
    return -LOG2PI + coef * Math.log(inverseVariance) - 0.5 * (distance * distance * Math.pow(inverseVariance, coef));
}

Example 4 with SphericalPolarCoordinates

use of dr.geo.math.SphericalPolarCoordinates in project beast-mcmc by beast-dev.

the class HyperSphereDistribution method latLongToCartesianInnerProduct.

private static double latLongToCartesianInnerProduct(double[] x, double[] y, double radius) {
    // x[] and y[] should be in the form (lat, long)
    if (x.length != 2 || y.length != 2) {
        throw new RuntimeException("Wrong dimensions");
    }
    final SphericalPolarCoordinates coordX = new SphericalPolarCoordinates(x[0], x[1], radius);
    final SphericalPolarCoordinates coordY = new SphericalPolarCoordinates(y[0], y[1], radius);
    return coordX.getCartesianCoordinates().dot(coordY.getCartesianCoordinates());
}

Example 5 with SphericalPolarCoordinates

use of dr.geo.math.SphericalPolarCoordinates in project beast-mcmc by beast-dev.

the class TimeSlicer method getKilometerGreatCircleDistance.

private static double getKilometerGreatCircleDistance(double[] location1, double[] location2) {
    SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(location1[0], location1[1]);
    SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(location2[0], location2[1]);
    return (coord1.distance(coord2));
}