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

Example 6 with SphericalPolarCoordinates

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

the class Branch2dRateToGrid 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));
}

Example 7 with SphericalPolarCoordinates

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

the class ContinuousDiffusionStatistic method getStatisticValue.

public double getStatisticValue(int dim) {
    double treeLength = 0;
    double treeDistance = 0;
    double totalMaxDistanceFromRoot = 0;
    // can only be used when cumulative and not associated with discrete state (not based on the distances on the branches from the root up that point)
    double maxDistanceFromRootCumulative = 0;
    double maxBranchDistanceFromRoot = 0;
    // can only be used when cumulative and not associated with discrete state (not based on the distances on the branches from the root up that point)
    double maxDistanceOverTimeFromRootWA = 0;
    double maxBranchDistanceOverTimeFromRootWA = 0;
    List<Double> rates = new ArrayList<Double>();
    List<Double> distances = new ArrayList<Double>();
    List<Double> times = new ArrayList<Double>();
    List<Double> traits = new ArrayList<Double>();
    List<double[]> traits2D = 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();
        String traitName = traitLikelihood.getTraitName();
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef node = tree.getNode(i);
            if (node != tree.getRoot()) {
                NodeRef parentNode = tree.getParent(node);
                boolean testNode = true;
                if (branchset.equals(BranchSet.CLADE)) {
                    try {
                        testNode = inClade(tree, node, taxonList);
                    } catch (TreeUtils.MissingTaxonException mte) {
                        throw new RuntimeException(mte.toString());
                    }
                } else if (branchset.equals(BranchSet.BACKBONE)) {
                    if (backboneTime > 0) {
                        testNode = onAncestralPathTime(tree, node, backboneTime);
                    } else {
                        try {
                            testNode = onAncestralPathTaxa(tree, node, taxonList);
                        } catch (TreeUtils.MissingTaxonException mte) {
                            throw new RuntimeException(mte.toString());
                        }
                    }
                }
                if (testNode) {
                    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);
                        //                        System.out.println(rate);
                        MultivariateDiffusionModel diffModel = traitLikelihood.diffusionModel;
                        double[] precision = diffModel.getPrecisionParameter().getParameterValues();
                        History history = null;
                        if (stateString != null) {
                            history = setUpHistory(markovJumpLikelihood.getHistoryForNode(tree, node, SITE), markovJumpLikelihood.getStatesForNode(tree, node)[SITE], markovJumpLikelihood.getStatesForNode(tree, parentNode)[SITE], timeLow, timeUp);
                        }
                        if (tree.getNodeHeight(parentNode) > upperHeight) {
                            timeUp = upperHeight;
                            traitUp = imputeValue(trait, parentTrait, upperHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, trueNoise);
                            if (stateString != null) {
                                history.truncateUpper(timeUp);
                            }
                        }
                        if (tree.getNodeHeight(node) < lowerHeight) {
                            timeLow = lowerHeight;
                            traitLow = imputeValue(trait, parentTrait, lowerHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, trueNoise);
                            if (stateString != null) {
                                history.truncateLower(timeLow);
                            }
                        }
                        if (dimension > traitLow.length) {
                            System.err.println("specified trait dimension for continuous trait summary, " + dimension + ", is > dimensionality of trait, " + traitLow.length + ". No trait summarized.");
                        } else {
                            traits.add(traitLow[(dimension - 1)]);
                        }
                        if (traitLow.length == 2) {
                            traits2D.add(traitLow);
                        }
                        double time;
                        if (stateString != null) {
                            time = history.getStateTime(stateString);
                        //                            System.out.println("tine before = "+(timeUp - timeLow)+", time after= "+time);
                        } else {
                            time = timeUp - timeLow;
                        }
                        treeLength += time;
                        times.add(time);
                        //setting up continuous trait values for heights in discrete trait history
                        if (stateString != null) {
                            history.setTraitsforHeights(traitUp, traitLow, precision, rate, trueNoise);
                        }
                        double[] rootTrait = traitLikelihood.getTraitForNode(tree, tree.getRoot(), traitName);
                        double timeFromRoot = (tree.getNodeHeight(tree.getRoot()) - timeLow);
                        if (useGreatCircleDistances && (trait.length == 2)) {
                            // Great Circle distance
                            double distance;
                            if (stateString != null) {
                                distance = history.getStateGreatCircleDistance(stateString);
                            } else {
                                distance = getGreatCircleDistance(traitLow, traitUp);
                            }
                            distances.add(distance);
                            if (time > 0) {
                                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 tempDistanceFromRootLow = rootCoord.distance(new SphericalPolarCoordinates(traitUp[0], traitUp[1]));
                            if (tempDistanceFromRootLow > totalMaxDistanceFromRoot) {
                                totalMaxDistanceFromRoot = tempDistanceFromRootLow;
                                if (stateString != null) {
                                    double[] stateTimeDistance = getStateTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, precision, branchRates, true);
                                    if (stateTimeDistance[0] > 0) {
                                        maxDistanceFromRootCumulative = tempDistanceFromRootLow * (stateTimeDistance[0] / timeFromRoot);
                                        maxDistanceOverTimeFromRootWA = maxDistanceFromRootCumulative / stateTimeDistance[0];
                                        maxBranchDistanceFromRoot = stateTimeDistance[1];
                                        maxBranchDistanceOverTimeFromRootWA = stateTimeDistance[1] / stateTimeDistance[0];
                                    }
                                } else {
                                    maxDistanceFromRootCumulative = tempDistanceFromRootLow;
                                    maxDistanceOverTimeFromRootWA = tempDistanceFromRootLow / timeFromRoot;
                                    double[] timeDistance = getTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, true);
                                    maxBranchDistanceFromRoot = timeDistance[1];
                                    maxBranchDistanceOverTimeFromRootWA = timeDistance[1] / timeDistance[0];
                                }
                                //distance between traitLow and traitUp for maxDistanceFromRootCumulative
                                if (timeUp == upperHeight) {
                                    if (time > 0) {
                                        maxDistanceFromRootCumulative = distance;
                                        maxDistanceOverTimeFromRootWA = distance / time;
                                        maxBranchDistanceFromRoot = distance;
                                        maxBranchDistanceOverTimeFromRootWA = distance / time;
                                    }
                                }
                            }
                        } else {
                            double distance;
                            if (stateString != null) {
                                distance = history.getStateNativeDistance(stateString);
                            } else {
                                distance = getNativeDistance(traitLow, traitUp);
                            }
                            distances.add(distance);
                            if (time > 0) {
                                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 > totalMaxDistanceFromRoot) {
                                totalMaxDistanceFromRoot = tempDistanceFromRoot;
                                if (stateString != null) {
                                    double[] stateTimeDistance = getStateTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, precision, branchRates, false);
                                    if (stateTimeDistance[0] > 0) {
                                        maxDistanceFromRootCumulative = tempDistanceFromRoot * (stateTimeDistance[0] / timeFromRoot);
                                        maxDistanceOverTimeFromRootWA = maxDistanceFromRootCumulative / stateTimeDistance[0];
                                        maxBranchDistanceFromRoot = stateTimeDistance[1];
                                        maxBranchDistanceOverTimeFromRootWA = stateTimeDistance[1] / stateTimeDistance[0];
                                    }
                                } else {
                                    maxDistanceFromRootCumulative = tempDistanceFromRoot;
                                    maxDistanceOverTimeFromRootWA = tempDistanceFromRoot / timeFromRoot;
                                    double[] timeDistance = getTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, false);
                                    maxBranchDistanceFromRoot = timeDistance[1];
                                    maxBranchDistanceOverTimeFromRootWA = timeDistance[1] / timeDistance[0];
                                }
                                //distance between traitLow and traitUp for maxDistanceFromRootCumulative
                                if (timeUp == upperHeight) {
                                    if (time > 0) {
                                        maxDistanceFromRootCumulative = distance;
                                        maxDistanceOverTimeFromRootWA = distance / time;
                                        maxBranchDistanceFromRoot = distance;
                                        maxBranchDistanceOverTimeFromRootWA = 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) {
            final double mean = DiscreteStatistics.mean(toArray(rates));
            return Math.sqrt(DiscreteStatistics.variance(toArray(rates), mean)) / mean;
        //weighted average
        } else {
            return treeDistance / treeLength;
        }
    } else if (summaryStat == summaryStatistic.TRAIT) {
        if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(traits));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(traits));
            return Math.sqrt(DiscreteStatistics.variance(toArray(traits), mean)) / mean;
        // default is average. A warning is thrown by the parser when trying to use WEIGHTED_AVERAGE
        } else {
            return DiscreteStatistics.mean(toArray(traits));
        }
    } else if (summaryStat == summaryStatistic.TRAIT2DAREA) {
        double area = getAreaFrom2Dtraits(traits2D, 0.99);
        return area;
    } 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;
        }
    //wavefront distance
    //TODO: restrict to non state-specific wavefrontDistance/rate
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE) {
        return maxDistanceFromRootCumulative;
    //            return maxBranchDistanceFromRoot;
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE_PHYLO) {
        return maxBranchDistanceFromRoot;
    //wavefront rate, only weighted average TODO: extend for average, median, COEFFICIENT_OF_VARIATION?
    } else if (summaryStat == summaryStatistic.WAVEFRONT_RATE) {
        return maxDistanceOverTimeFromRootWA;
    //            return maxBranchDistanceOverTimeFromRootWA;
    } else if (summaryStat == summaryStatistic.DIFFUSION_DISTANCE) {
        return treeDistance;
    //DIFFUSION_TIME
    } else if (summaryStat == summaryStatistic.DISTANCE_TIME_CORRELATION) {
        if (summaryMode == Mode.SPEARMAN) {
            return getSpearmanRho(convertDoubles(times), convertDoubles(distances));
        } else if (summaryMode == Mode.R_SQUARED) {
            Regression r = new Regression(convertDoubles(times), convertDoubles(distances));
            return r.getRSquared();
        } else {
            Regression r = new Regression(convertDoubles(times), convertDoubles(distances));
            return r.getCorrelationCoefficient();
        }
    } else {
        return treeLength;
    }
}

Example 8 with SphericalPolarCoordinates

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

the class ContinuousDiffusionStatistic method getGreatCircleDistance.

public double getGreatCircleDistance(double[] loc1, double[] loc2) {
    SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(loc1[0], loc1[1]);
    SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(loc2[0], loc2[1]);
    return coord1.distance(coord2);
}

Example 9 with SphericalPolarCoordinates

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

the class DiffusionRateCovarianceStatistic method getStatisticValue.

/**
     * @return whatever Philippe wants
     */
public double getStatisticValue(int dim) {
    String traitName = traitLikelihoods.get(0).getTraitName();
    for (AbstractMultivariateTraitLikelihood traitLikelihood : traitLikelihoods) {
        MultivariateTraitTree tree = traitLikelihood.getTreeModel();
        int counter = 0;
        int index = 0;
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef child = tree.getNode(i);
            NodeRef parent = tree.getParent(child);
            if (parent != null & !tree.isRoot(parent)) {
                double[] childTrait = traitLikelihood.getTraitForNode(tree, child, traitName);
                double[] parentTrait = traitLikelihood.getTraitForNode(tree, parent, traitName);
                double childTime = tree.getBranchLength(child);
                double parentTime = tree.getBranchLength(parent);
                NodeRef grandParent = tree.getParent(parent);
                double[] grandParentTrait = traitLikelihood.getTraitForNode(tree, grandParent, traitName);
                if (useGreatCircleDistances && (childTrait.length == 2)) {
                    // Great Circle distance
                    SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(childTrait[0], childTrait[1]);
                    SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(parentTrait[0], parentTrait[1]);
                    double childDistance = coord1.distance(coord2);
                    SphericalPolarCoordinates coord3 = new SphericalPolarCoordinates(grandParentTrait[0], grandParentTrait[1]);
                    double parentDistance = coord2.distance(coord3);
                    if (!diffusionCoefficient) {
                        childRate[index] = childDistance / childTime;
                        parentRate[index] = parentDistance / parentTime;
                    } else {
                        childRate[index] = Math.pow(childDistance, 2) / (4 * childTime);
                        parentRate[index] = Math.pow(parentDistance, 2) / (4 * parentTime);
                    }
                } else {
                    double childDistance = getNativeDistance(childTrait, parentTrait);
                    double parentDistance = getNativeDistance(parentTrait, grandParentTrait);
                    if (!diffusionCoefficient) {
                        childRate[index] = childDistance / childTime;
                        parentRate[index] = parentDistance / parentTime;
                    } else {
                        childRate[index] = Math.pow(childDistance, 2) / (4 * childTime);
                        parentRate[index] = Math.pow(parentDistance, 2) / (4 * parentTime);
                    }
                }
                index += 1;
            }
        }
    }
    return DiscreteStatistics.covariance(childRate, parentRate);
}