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

Example 11 with WrappedVector

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

the class AbstractValuesViaFullConditionalDelegate method getTraitForNode.

@Override
protected double[] getTraitForNode(NodeRef node) {
    assert simulationProcess != null;
    assert node != null;
    final int nodeBuffer = likelihoodDelegate.getActiveNodeIndex(node.getNumber());
    if (node.getNumber() >= tree.getExternalNodeCount()) {
        // Not external node
        return new double[0];
    // return new MeanAndVariance(new double[0]);
    }
    // new double[dimPartial * numTraits];
    double[] conditionalNodeBuffer = null;
    likelihoodDelegate.getPostOrderPartial(node.getNumber(), partialNodeBuffer);
    final double[] sample = new double[dimTrait * numTraits];
    int partialOffset = 0;
    int sampleOffset = 0;
    for (int trait = 0; trait < numTraits; ++trait) {
        if (missingInformation.isPartiallyMissing(node.getNumber(), trait)) {
            if (conditionalNodeBuffer == null) {
                conditionalNodeBuffer = new double[dimPartial * numTraits];
                simulationProcess.cacheSimulatedTraits(node);
                cdi.getPreOrderPartial(nodeBuffer, conditionalNodeBuffer);
            }
            System.err.println("Missing tip = " + node.getNumber() + " (" + nodeBuffer + "), trait = " + trait);
            final WrappedVector preMean = new WrappedVector.Raw(conditionalNodeBuffer, partialOffset, dimTrait);
            final DenseMatrix64F preVar = wrap(conditionalNodeBuffer, partialOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
            final WrappedVector postObs = new WrappedVector.Raw(partialNodeBuffer, partialOffset, dimTrait);
            System.err.println("post: " + postObs);
            System.err.println("pre : " + preMean);
            System.err.println("V: " + preVar);
            if (!missingInformation.isCompletelyMissing(node.getNumber(), trait)) {
                final PartiallyMissingInformation.HashedIntArray intArray = missingInformation.getMissingIndices(node.getNumber(), trait);
                final int[] missing = intArray.getArray();
                final int[] observed = intArray.getComplement();
                ConditionalVarianceAndTransform2 transform = new ConditionalVarianceAndTransform2(preVar, missing, observed);
                final WrappedVector cM = transform.getConditionalMean(// Tip value
                partialNodeBuffer, // Tip value
                partialOffset, conditionalNodeBuffer, // Mean value
                partialOffset);
                computeValueWithMissing(// input mean
                cM, // input variance,
                transform.getConditionalCholesky(), // output sample
                new WrappedVector.Indexed(sample, sampleOffset, missing, missing.length), transform.getTemporaryStorage());
                System.err.println("cM: " + cM);
                System.err.println("CV: " + transform.getConditionalVariance());
                System.err.println("value: " + new WrappedVector.Raw(sample, sampleOffset, dimTrait));
            }
        } else {
            computeValueWithNoMissing(partialNodeBuffer, partialOffset, sample, sampleOffset, dimTrait);
        }
        partialOffset += dimPartial;
        sampleOffset += dimTrait;
    }
    return sample;
// return new MeanAndVariance(sample);
}

Example 12 with WrappedVector

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

the class SafeMultivariateWithDriftIntegrator method partialMean.

@Override
void partialMean(int ibo, int jbo, int kbo, int ido, int jdo) {
    if (TIMING) {
        startTime("peel4");
    }
    final double[] displacementi = vectorDispi;
    final double[] displacementj = vectorDispj;
    for (int g = 0; g < dimTrait; ++g) {
        displacementi[g] = partials[ibo + g] - displacements[ido + g];
        displacementj[g] = partials[jbo + g] - displacements[jdo + g];
    }
    final double[] tmp = vectorPMk;
    computeWeightedSum(displacementi, displacementj, dimTrait, tmp);
    final WrappedVector kPartials = new WrappedVector.Raw(partials, kbo, dimTrait);
    final WrappedVector wrapTmp = new WrappedVector.Raw(tmp, 0, dimTrait);
    safeSolve(matrixPk, wrapTmp, kPartials, false);
    if (TIMING) {
        endTime("peel4");
        startTime("peel5");
    }
    if (DEBUG) {
        System.err.print("\t\tdisp i:");
        for (int e = 0; e < dimTrait; ++e) {
            System.err.print(" " + displacements[ido + e]);
        }
        System.err.println("");
        System.err.print("\t\tdisp j:");
        for (int e = 0; e < dimTrait; ++e) {
            System.err.print(" " + displacements[jdo + e]);
        }
    }
// return ck;
}

Example 13 with WrappedVector

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

the class AbstractParticleOperator method doOperation.

@Override
public double doOperation() {
    if (shouldUpdatePreconditioning()) {
        preconditioning = setupPreconditioning();
    }
    WrappedVector position = getInitialPosition();
    double hastingsRatio = integrateTrajectory(position);
    setParameter(position, parameter);
    if (CHECK_MATRIX_ILL_CONDITIONED & getCount() % 100 == 0) {
        productProvider.getTimeScaleEigen();
    }
    return hastingsRatio;
}

Example 14 with WrappedVector

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

the class AbstractZigZagOperator method integrateTrajectory.

@Override
double integrateTrajectory(WrappedVector position) {
    String signString;
    if (DEBUG_SIGN) {
        signString = printSign(position);
        System.err.println(signString);
    }
    if (TIMING) {
        timer.startTimer("warmUp");
    }
    WrappedVector momentum = drawInitialMomentum();
    WrappedVector velocity = drawInitialVelocity(momentum);
    WrappedVector gradient = getInitialGradient();
    WrappedVector action = getPrecisionProduct(velocity);
    BounceState bounceState = new BounceState(drawTotalTravelTime());
    if (TIMING) {
        timer.stopTimer("warmUp");
    }
    int count = 0;
    double[] p, v, a, g, m;
    if (TEST_NATIVE_OPERATOR) {
        p = position.getBuffer().clone();
        v = velocity.getBuffer().clone();
        a = action.getBuffer().clone();
        g = gradient.getBuffer().clone();
        m = momentum.getBuffer().clone();
        nativeZigZag.operate(columnProvider, p, v, a, g, m, bounceState.remainingTime);
    }
    if (TEST_CRITICAL_REGION) {
        nativeZigZag.enterCriticalRegion(position.getBuffer(), velocity.getBuffer(), action.getBuffer(), gradient.getBuffer(), momentum.getBuffer());
    }
    if (TIMING) {
        timer.startTimer("integrateTrajectory");
    }
    while (bounceState.isTimeRemaining()) {
        if (DEBUG) {
            debugBefore(position, count);
        }
        final MinimumTravelInformation firstBounce;
        if (taskPool != null) {
            if (FUSE) {
                if (TIMING) {
                    timer.startTimer("getNext");
                }
                firstBounce = getNextBounceParallel(position, velocity, action, gradient, momentum);
                if (TIMING) {
                    timer.stopTimer("getNext");
                }
                if (TEST_NATIVE_BOUNCE) {
                    testNative(firstBounce, position, velocity, action, gradient, momentum);
                }
            } else {
                MinimumTravelInformation boundaryBounce = getNextBoundaryBounce(position, velocity);
                MinimumTravelInformation gradientBounce = getNextGradientBounceParallel(action, gradient, momentum);
                firstBounce = (boundaryBounce.time < gradientBounce.time) ? new MinimumTravelInformation(boundaryBounce.time, boundaryBounce.index, Type.BOUNDARY) : new MinimumTravelInformation(gradientBounce.time, gradientBounce.index, Type.GRADIENT);
            }
        } else {
            if (FUSE) {
                if (TIMING) {
                    timer.startTimer("getNext");
                }
                firstBounce = getNextBounce(position, velocity, action, gradient, momentum);
                if (TIMING) {
                    timer.stopTimer("getNext");
                }
                if (TEST_NATIVE_BOUNCE) {
                    testNative(firstBounce, position, velocity, action, gradient, momentum);
                }
            } else {
                if (TIMING) {
                    timer.startTimer("getNextBoundary");
                }
                MinimumTravelInformation boundaryBounce = getNextBoundaryBounce(position, velocity);
                if (TIMING) {
                    timer.stopTimer("getNextBoundary");
                    timer.startTimer("getNextGradient");
                }
                MinimumTravelInformation gradientBounce = getNextGradientBounce(action, gradient, momentum);
                if (TIMING) {
                    timer.stopTimer("getNextGradient");
                }
                firstBounce = (boundaryBounce.time < gradientBounce.time) ? new MinimumTravelInformation(boundaryBounce.time, boundaryBounce.index, Type.BOUNDARY) : new MinimumTravelInformation(gradientBounce.time, gradientBounce.index, Type.GRADIENT);
            }
        }
        bounceState = doBounce(bounceState, firstBounce, position, velocity, action, gradient, momentum);
        if (DEBUG) {
            debugAfter(bounceState, position);
            String newSignString = printSign(position);
            System.err.println(newSignString);
            if (bounceState.type != Type.BOUNDARY && signString.compareTo(newSignString) != 0) {
                System.err.println("Sign error");
            }
        }
        ++count;
    }
    if (TIMING) {
        timer.stopTimer("integrateTrajectory");
    }
    if (TEST_CRITICAL_REGION) {
        nativeZigZag.exitCriticalRegion();
    }
    if (TEST_NATIVE_OPERATOR) {
        if (!close(p, position.getBuffer())) {
            System.err.println("c: " + new WrappedVector.Raw(p, 0, 10));
            System.err.println("c: " + new WrappedVector.Raw(position.getBuffer(), 0, 10));
        } else {
            System.err.println("close");
        }
    }
    if (DEBUG_SIGN) {
        printSign(position);
    }
    return 0.0;
}

Example 15 with WrappedVector

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

the class HamiltonianMonteCarloOperator method leapFrog.

private double leapFrog() throws NumericInstabilityException {
    if (DEBUG) {
        System.err.println("HMC step size: " + stepSize);
    }
    final double[] position = leapFrogEngine.getInitialPosition();
    final WrappedVector momentum = mask(preconditioning.drawInitialMomentum(), mask);
    final double prop = getKineticEnergy(momentum) + leapFrogEngine.getParameterLogJacobian();
    leapFrogEngine.updateMomentum(position, momentum.getBuffer(), mask(gradientProvider.getGradientLogDensity(), mask), stepSize / 2);
    int nStepsThisLeap = getNumberOfSteps();
    for (int i = 0; i < nStepsThisLeap; i++) {
        try {
            leapFrogEngine.updatePosition(position, momentum, stepSize);
        } catch (ArithmeticException e) {
            throw new NumericInstabilityException();
        }
        if (i < (nStepsThisLeap - 1)) {
            try {
                leapFrogEngine.updateMomentum(position, momentum.getBuffer(), mask(gradientProvider.getGradientLogDensity(), mask), stepSize);
            } catch (ArithmeticException e) {
                throw new NumericInstabilityException();
            }
        }
    }
    leapFrogEngine.updateMomentum(position, momentum.getBuffer(), mask(gradientProvider.getGradientLogDensity(), mask), stepSize / 2);
    final double res = getKineticEnergy(momentum) + leapFrogEngine.getParameterLogJacobian();
    // hasting ratio
    return prop - res;
}