Examples with Vector dr.math.matrixAlgebra.Vector used on opensource projects

Example 6 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class TopographicalMap method mapToXML.

//	static
public static Element mapToXML(double[][] map, int numLevels) {
    Element root = new Element(TERRAIN);
    Element tileTypes = new Element(TYPE_TYPE);
    Element size = new Element("size");
    Element content = new Element("content");
    root.addContent(tileTypes);
    root.addContent(size);
    root.addContent(content);
    // Make sizes
    size.addContent(new Element("rows").addContent(Integer.toString(map.length)));
    size.addContent(new Element("columns").addContent(Integer.toString(map[0].length)));
    size.addContent(new Element("height").addContent(Integer.toString(map.length)));
    size.addContent(new Element("width").addContent(Integer.toString(map[0].length)));
    // Make cutpoints
    double min = Double.MAX_VALUE;
    double max = Double.MIN_VALUE;
    for (int i = 0; i < map.length; i++) {
        for (int j = 0; j < map[0].length; j++) {
            if (!Double.isNaN(map[i][j])) {
                //				}
                if (map[i][j] < min)
                    min = map[i][j];
                if (map[i][j] > max)
                    max = map[i][j];
            }
        }
    }
    //		System.err.println("min = ");
    double[] cuts = new double[numLevels];
    //		cuts[numLevels-1] = max;
    double range = max - min;
    double delta = range / numLevels;
    for (int i = 0; i < numLevels; i++) cuts[i] = min + delta * (i + 1);
    int deltaColor = 255 / numLevels;
    // Make tile types
    double average = min + delta / 2;
    int grey = 255 - deltaColor;
    for (int i = 0; i < numLevels; i++) {
        Element type = new Element("type");
        type.addContent(new Element("name").addContent("L" + Integer.toString(i)));
        type.addContent(new Element("cost").addContent(Integer.toString((int) average)));
        Element color = new Element("color");
        color.setAttribute("r", Integer.toString(grey));
        color.setAttribute("g", Integer.toString(grey));
        color.setAttribute("b", Integer.toString(grey));
        type.addContent(color);
        tileTypes.addContent(type);
        average += delta;
        grey -= deltaColor;
    }
    Element blocked = new Element("type");
    blocked.addContent(new Element("name").addContent("B"));
    blocked.addContent(new Element("blocked"));
    Element color = new Element("color");
    color.setAttribute("r", Integer.toString(255));
    color.setAttribute("g", Integer.toString(255));
    color.setAttribute("b", Integer.toString(255));
    blocked.addContent(color);
    tileTypes.addContent(blocked);
    int count = 0;
    // Make content
    content.addContent(new Element("default").addContent("B"));
    StringBuffer ascii = new StringBuffer();
    boolean newLine = true;
    for (int r = 0; r < map.length; r++) {
        for (int c = 0; c < map[0].length; c++) {
            if (!Double.isNaN(map[r][c])) {
                count++;
                double value = map[r][c];
                //				if( value != Double.NaN ) {
                int cut = 0;
                try {
                    while (value > cuts[cut]) cut++;
                } catch (Exception e) {
                    System.err.println("Error: " + e);
                    System.err.println("value = " + value);
                    System.err.println("cuts = " + new Vector(cuts));
                    System.err.println("min = " + min);
                    System.err.println("max = " + max);
                    System.exit(-1);
                }
                Element cell = new Element("column");
                cell.setAttribute("r", Integer.toString(r));
                cell.setAttribute("c", Integer.toString(c));
                cell.setAttribute("length", "1");
                cell.addContent("L" + Integer.toString(cut));
                //					if( count < 10)
                content.addContent(cell);
                ascii.append((cut + 1) + " ");
            //				}  else {
            //					System.err.println("about time");
            //					System.exit(-1);
            } else {
                ascii.append("NA ");
            }
        }
        ascii.append("\n");
    }
    System.out.println(ascii);
    return root;
}

Example 7 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class NonPhylogeneticMultivariateTraitLikelihood method calculateLogLikelihood.

public double calculateLogLikelihood() {
    double[][] traitPrecision = diffusionModel.getPrecisionmatrix();
    double logDetTraitPrecision = Math.log(diffusionModel.getDeterminantPrecisionMatrix());
    double[] marginalRoot = tmp2;
    if (computeWishartStatistics) {
        wishartStatistics = new WishartSufficientStatistics(dimTrait);
    }
    // Compute the contribution of each datum at the root
    SufficientStatistics stats = computeInnerProductsForTips(traitPrecision, tmp2);
    double conditionalSumWeight = stats.sumWeight;
    double conditionalProductWeight = stats.productWeight;
    double innerProducts = stats.innerProduct;
    int nonMissingTips = stats.nonMissingTips;
    // Add in prior and integrate
    double sumWeight = conditionalSumWeight + rootPriorSampleSize;
    double productWeight = conditionalProductWeight * rootPriorSampleSize;
    double rootPrecision = productWeight / sumWeight;
    final int rootIndex = treeModel.getRoot().getNumber();
    int rootOffset = dim * rootIndex;
    for (int datum = 0; datum < numData; ++datum) {
        // Determine marginal root (scaled) mean
        for (int d = 0; d < dimTrait; ++d) {
            marginalRoot[d] = conditionalSumWeight * meanCache[rootOffset + d] + rootPriorSampleSize * rootPriorMean[d];
        }
        // Compute outer product contribution from prior
        double yAy1 = computeWeightedAverageAndSumOfSquares(rootPriorMean, Ay, traitPrecision, dimTrait, rootPriorSampleSize);
        // TODO Only need to compute once
        innerProducts += yAy1;
        if (DEBUG_NO_TREE) {
            System.err.println("OP for root");
            System.err.println("Value  = " + new Vector(rootPriorMean));
            System.err.print("Prec   = \n" + new Matrix(traitPrecision));
            System.err.println("Weight = " + rootPriorSampleSize + "\n");
        }
        // Compute outer product differences to complete square
        double yAy2 = computeWeightedAverageAndSumOfSquares(marginalRoot, Ay, traitPrecision, dimTrait, 1.0 / sumWeight);
        innerProducts -= yAy2;
        // Add prior on root contribution
        if (computeWishartStatistics) {
            final double[] outerProducts = wishartStatistics.getScaleMatrix();
            final double weight = conditionalSumWeight * rootPriorSampleSize / sumWeight;
            for (int i = 0; i < dimTrait; i++) {
                final double diffi = meanCache[rootOffset + i] - rootPriorMean[i];
                for (int j = 0; j < dimTrait; j++) {
                    outerProducts[i * dimTrait + j] += diffi * weight * (meanCache[rootOffset + j] - rootPriorMean[j]);
                }
            }
            wishartStatistics.incrementDf(1);
        }
        rootOffset += dimTrait;
    }
    if (DEBUG_NO_TREE) {
        System.err.println("SumWeight    : " + sumWeight);
        System.err.println("ProductWeight: " + productWeight);
        System.err.println("Total OP     : " + innerProducts);
    }
    // Compute log likelihood
    double logLikelihood = -LOG_SQRT_2_PI * dimTrait * nonMissingTips * numData + 0.5 * logDetTraitPrecision * nonMissingTips * numData + 0.5 * Math.log(rootPrecision) * dimTrait * numData - 0.5 * innerProducts;
    if (DEBUG_NO_TREE) {
        System.err.println("logLikelihood (final) = " + logLikelihood);
        System.err.println("numData = " + numData);
    }
    // Should redraw internal node states when needed
    areStatesRedrawn = false;
    return logLikelihood;
}

Example 8 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class NonPhylogeneticMultivariateTraitLikelihood method computeInnerProductsForTips.

// Useful identity for computing outerproducts for Wishart statistics
// \sum (y_i - \bar{y}) (y_i - \bar{y})^{t} = \sum y_i y_i^{t} - n \bar{y} \bar{y}^t
private SufficientStatistics computeInnerProductsForTips(double[][] traitPrecision, double[] tmpVector) {
    // Compute the contribution of each datum at the root
    final int rootIndex = treeModel.getRoot().getNumber();
    final int meanOffset = dim * rootIndex;
    // Zero-out root mean
    for (int d = 0; d < dim; ++d) {
        meanCache[meanOffset + d] = 0;
    }
    double innerProducts = 0.0;
    // Compute the contribution of each datum at the root
    double productWeight = 1.0;
    double sumWeight = 0.0;
    int nonMissingTips = 0;
    for (int i = 0; i < treeModel.getExternalNodeCount(); ++i) {
        NodeRef tipNode = treeModel.getExternalNode(i);
        final int tipNumber = tipNode.getNumber();
        double tipWeight = 0.0;
        if (!missingTraits.isCompletelyMissing(tipNumber)) {
            tipWeight = 1.0 / getLengthToRoot(tipNode);
            int tipOffset = dim * tipNumber;
            int rootOffset = dim * rootIndex;
            for (int datum = 0; datum < numData; ++datum) {
                // Add weighted tip value
                for (int d = 0; d < dimTrait; ++d) {
                    meanCache[rootOffset + d] += tipWeight * meanCache[tipOffset + d];
                    tmpVector[d] = meanCache[tipOffset + d];
                }
                // Compute outer product
                double yAy = computeWeightedAverageAndSumOfSquares(tmpVector, Ay, traitPrecision, dimTrait, tipWeight);
                innerProducts += yAy;
                if (DEBUG_NO_TREE) {
                    System.err.println("OP for " + tipNumber + " = " + yAy);
                    System.err.println("Value  = " + new Vector(tmpVector));
                    System.err.print("Prec   =\n" + new Matrix(traitPrecision));
                    System.err.println("weight = " + tipWeight + "\n");
                }
                tipOffset += dimTrait;
                rootOffset += dimTrait;
            }
            if (computeWishartStatistics) {
                incrementOuterProducts(tipNumber, tipWeight);
            }
        }
        if (tipWeight > 0.0) {
            sumWeight += tipWeight;
            productWeight *= tipWeight;
            ++nonMissingTips;
        }
    }
    lowerPrecisionCache[rootIndex] = sumWeight;
    normalize(meanCache, meanOffset, dim, sumWeight);
    if (computeWishartStatistics) {
        incrementOuterProducts(rootIndex, -sumWeight);
        wishartStatistics.incrementDf(-1);
    }
    return new SufficientStatistics(sumWeight, productWeight, innerProducts, nonMissingTips);
}

Example 9 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class IntegratedMultivariateTraitLikelihood method postOrderTraverse.

void postOrderTraverse(MultivariateTraitTree treeModel, NodeRef node, double[][] precisionMatrix, double logDetPrecisionMatrix, boolean cacheOuterProducts) {
    final int thisNumber = node.getNumber();
    if (treeModel.isExternal(node)) {
        if (missingTraits.isCompletelyMissing(thisNumber)) {
            upperPrecisionCache[thisNumber] = 0;
            // Needed in the pre-order traversal
            lowerPrecisionCache[thisNumber] = 0;
        } else {
            // not missing tip trait
            // changeou
            //    upperPrecisionCache[thisNumber] = (1.0 / getRescaledBranchLengthForPrecision(node)) * Math.pow(cacheHelper.getOUFactor(node), 2);
            upperPrecisionCache[thisNumber] = cacheHelper.getUpperPrecFactor(node) * Math.pow(cacheHelper.getOUFactor(node), 2);
            lowerPrecisionCache[thisNumber] = Double.POSITIVE_INFINITY;
        }
        return;
    }
    final NodeRef childNode0 = treeModel.getChild(node, 0);
    final NodeRef childNode1 = treeModel.getChild(node, 1);
    postOrderTraverse(treeModel, childNode0, precisionMatrix, logDetPrecisionMatrix, cacheOuterProducts);
    postOrderTraverse(treeModel, childNode1, precisionMatrix, logDetPrecisionMatrix, cacheOuterProducts);
    final int childNumber0 = childNode0.getNumber();
    final int childNumber1 = childNode1.getNumber();
    final int meanOffset0 = dim * childNumber0;
    final int meanOffset1 = dim * childNumber1;
    final int meanThisOffset = dim * thisNumber;
    final double precision0 = upperPrecisionCache[childNumber0];
    final double precision1 = upperPrecisionCache[childNumber1];
    final double totalPrecision = precision0 + precision1;
    final double factorOU0 = cacheHelper.getOUFactor(childNode0);
    final double factorOU1 = cacheHelper.getOUFactor(childNode1);
    doPeel(thisNumber, meanThisOffset, meanOffset0, meanOffset1, totalPrecision, precision0, precision1, missingTraits, thisNumber, precisionMatrix, logDetPrecisionMatrix, factorOU0, factorOU1, cacheOuterProducts, // TODO arguments to remove
    node, // TODO arguments to remove
    childNode0, // TODO arguments to remove
    childNode1, true, false);
    final boolean DO_CLAMP = true;
    final boolean debug = false;
    if (DO_CLAMP && nodeToClampMap != null && nodeToClampMap.containsKey(node)) {
        // TODO precompute boolean contains for all nodes
        RestrictedPartials clamp = nodeToClampMap.get(node);
        final int clampIndex = clamp.getIndex();
        final int clampOffset = dim * clampIndex;
        final int spareOffset = dim * spareIndex;
        // Copy partial into meanCache // TODO Only when value changes
        for (int i = 0; i < dim; ++i) {
            meanCache[clampOffset + i] = clamp.getPartial(i);
        }
        if (debug) {
            System.err.println("BEFORE");
            System.err.println(new Vector(logRemainderDensityCache));
            System.err.println(new Vector(meanCache));
            System.err.println(new Vector(lowerPrecisionCache));
            System.err.println(new Vector(upperPrecisionCache));
            System.err.println("");
        }
        final double precisionThis = lowerPrecisionCache[thisNumber];
        final double precisionClamp = clamp.getPriorSampleSize() / rescaleLength(1.0);
        final double precisionNew = precisionThis + precisionClamp;
        doPeel(spareIndex, spareOffset, meanThisOffset, clampOffset, precisionNew, precisionThis, precisionClamp, missingTraits, clampIndex, precisionMatrix, logDetPrecisionMatrix, // TODO Do yet figured out for OU models
        1.0, // TODO Do yet figured out for OU models
        1.0, cacheOuterProducts, node, null, null, true, false);
        // Move values from clampIndex -> thisIndex
        lowerPrecisionCache[thisNumber] = lowerPrecisionCache[spareIndex];
        upperPrecisionCache[thisNumber] = upperPrecisionCache[spareIndex];
        for (int i = 0; i < dim; ++i) {
            meanCache[meanThisOffset + i] = meanCache[spareOffset + i];
        }
    }
    if (debug) {
        System.err.println(thisNumber);
        System.err.println(new Vector(logRemainderDensityCache));
        System.err.println(new Vector(meanCache));
        System.err.println(new Vector(lowerPrecisionCache));
        System.err.println(new Vector(upperPrecisionCache));
        System.err.println("");
    }
}

Example 10 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class KDEDistributionTest method testNormalKDE.

public void testNormalKDE() {
    int length = 100;
    Double[] sample = new Double[length];
    for (int i = 0; i < length; i++) {
        sample[i] = (double) (i + 1);
    }
    //        final int gridSize = 256;
    NormalKDEDistribution kde = new NormalKDEDistribution(sample, null, null, null);
    System.out.println("bw = " + kde.getBandWidth());
    assertEquals(rBandWidth[0], kde.getBandWidth(), tolerance);
    final int gridSize = NormalKDEDistribution.MINIMUM_GRID_SIZE;
    double[] testPoints = new double[gridSize];
    double x = kde.getFromPoint();
    double delta = (kde.getToPoint() - kde.getFromPoint()) / (gridSize - 1);
    for (int i = 0; i < gridSize; ++i) {
        testPoints[i] = x;
        x += delta;
    }
    System.err.println("Eval @ " + new Vector(testPoints));
    double[] testDensity = new double[gridSize];
    for (int i = 0; i < gridSize; ++i) {
        testDensity[i] = kde.pdf(testPoints[i]);
    }
    System.err.println("Den    " + new Vector(testDensity));
    System.err.println("den[0] = " + testDensity[0]);
    System.err.println("den[N] = " + testDensity[NormalKDEDistribution.MINIMUM_GRID_SIZE - 1]);
//  System.exit(-1);
}