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

Example 6 with WrappedVector

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

the class AbstractParticleOperator method updateAction.

void updateAction(WrappedVector action, ReadableVector velocity, int eventIndex) {
    if (TEST_CRITICAL_REGION) {
        if (nativeZigZag.inCriticalRegion()) {
            nativeZigZag.exitCriticalRegion();
        }
    }
    WrappedVector column = getPrecisionColumn(eventIndex);
    if (TIMING) {
        timer.startTimer("updateAction");
    }
    final double[] a = action.getBuffer();
    final double[] c = column.getBuffer();
    final double twoV = 2 * velocity.get(eventIndex);
    for (int i = 0, len = a.length; i < len; ++i) {
        a[i] += twoV * c[i];
    }
    if (TIMING) {
        timer.stopTimer("updateAction");
    }
    if (mask != null) {
        applyMask(a);
    }
}

Example 7 with WrappedVector

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

the class BouncyParticleOperator method integrateTrajectory.

@Override
double integrateTrajectory(WrappedVector position) {
    WrappedVector velocity = drawInitialVelocity();
    WrappedVector gradient = getInitialGradient();
    WrappedVector action = getPrecisionProduct(velocity);
    BounceState bounceState = new BounceState(drawTotalTravelTime());
    while (bounceState.remainingTime > 0) {
        if (bounceState.type == Type.BOUNDARY) {
            updateAction(action, velocity, bounceState.index);
        } else {
            action = getPrecisionProduct(velocity);
        }
        double v_Phi_x = -innerProduct(velocity, gradient);
        double v_Phi_v = innerProduct(velocity, action);
        double tMin = Math.max(0.0, -v_Phi_x / v_Phi_v);
        double U_min = tMin * tMin / 2 * v_Phi_v + tMin * v_Phi_x;
        double bounceTime = getBounceTime(v_Phi_v, v_Phi_x, U_min);
        MinimumTravelInformation travelInfo = getTimeToBoundary(position, velocity);
        double refreshTime = getRefreshTime();
        bounceState = doBounce(bounceState.remainingTime, bounceTime, travelInfo, refreshTime, position, velocity, gradient, action);
    }
    storedVelocity = velocity;
    return 0.0;
}

Example 8 with WrappedVector

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

the class IrreversibleZigZagOperator method integrateTrajectory.

double integrateTrajectory(WrappedVector position) {
    WrappedVector momentum = drawInitialMomentum();
    WrappedVector velocity = drawInitialVelocity(momentum);
    WrappedVector gradient = getInitialGradient();
    WrappedVector action = getPrecisionProduct(velocity);
    BounceState bounceState = new BounceState(drawTotalTravelTime());
    int count = 0;
    if (TIMING) {
        timer.startTimer("integrateTrajectory");
    }
    while (bounceState.isTimeRemaining()) {
        final MinimumTravelInformation firstBounce;
        if (TIMING) {
            timer.startTimer("getNext");
        }
        firstBounce = getNextBounce(position, velocity, action, gradient, momentum);
        if (TIMING) {
            timer.stopTimer("getNext");
        }
        if (CPP_NEXT_BOUNCE) {
            MinimumTravelInformation test = testNative(position, velocity, action, gradient);
            System.exit(-1);
        }
        bounceState = doBounce(bounceState, firstBounce, position, velocity, action, gradient, momentum);
        ++count;
    }
    if (TIMING) {
        timer.stopTimer("integrateTrajectory");
    }
    return 0.0;
}

Example 9 with WrappedVector

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

the class NoUTurnOperator method takeOneStep.

private double[] takeOneStep(long m, double[] initialPosition) {
    double[] endPosition = Arrays.copyOf(initialPosition, initialPosition.length);
    // final double[][] mass = massProvider.getMass();
    final WrappedVector initialMomentum = mask(preconditioning.drawInitialMomentum(), mask);
    final double initialJointDensity = getJointProbability(gradientProvider, initialMomentum);
    double logSliceU = Math.log(MathUtils.nextDouble()) + initialJointDensity;
    TreeState trajectoryTree = new TreeState(initialPosition, initialMomentum.getBuffer(), 1, true);
    // Trajectory of Hamiltonian dynamics endowed with a binary tree structure.
    int height = 0;
    while (trajectoryTree.flagContinue) {
        double[] tmp = updateTrajectoryTree(trajectoryTree, height, logSliceU, initialJointDensity);
        if (tmp != null) {
            endPosition = tmp;
        }
        height++;
        if (height > options.maxHeight) {
            trajectoryTree.flagContinue = false;
        }
    }
    stepSizeInformation.update(m, trajectoryTree.cumAcceptProb, trajectoryTree.numAcceptProbStates);
    return endPosition;
}

Example 10 with WrappedVector

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

the class MultivariateConditionalOnTipsRealizedDelegate method simulateTraitForInternalNode.

// private ReadableVector getMeanWithDrift(final ReadableVector mean,
// final ReadableVector drift) {
// return new ReadableVector.Sum(mean, drift);
// }
// private ReadableVector getMeanWithDrift(double[] mean, int offsetMean, double[] drift, int dim) {
// for (int i = 0;i < dim; ++i) {
// tmpDrift[i] = mean[offsetMean + i] + drift[i];
// }
// return new WrappedVector.Raw(tmpDrift, 0, dimTrait);
// }
private void simulateTraitForInternalNode(final int offsetSample, final int offsetParent, final int offsetPartial, final double branchPrecision) {
    if (!Double.isInfinite(branchPrecision)) {
        // Here we simulate X_j | X_pa(j), Y
        final WrappedVector M0 = new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait);
        final DenseMatrix64F P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
        // final ReadableVector parentSample = new WrappedVector.Raw(sample, offsetParent, dimTrait);
        final ReadableVector M1;
        final DenseMatrix64F P1;
        M1 = getMeanBranch(offsetParent);
        P1 = getPrecisionBranch(branchPrecision);
        // if (hasNoDrift) {
        // M1 = parentSample; // new WrappedVector.Raw(sample, offsetParent, dimTrait);
        // P1 = new DenseMatrix64F(dimTrait, dimTrait);
        // CommonOps.scale(branchPrecision, Pd, P1);
        // } else {
        // M1 = getMeanWithDrift(parentSample,
        // new WrappedVector.Raw(displacementBuffer, 0, dimTrait)); //getMeanWithDrift(sample, offsetParent, displacementBuffer, dimTrait);
        // P1 = DenseMatrix64F.wrap(dimTrait, dimTrait, precisionBuffer);
        // }
        final WrappedVector M2 = new WrappedVector.Raw(tmpMean, 0, dimTrait);
        final DenseMatrix64F P2 = new DenseMatrix64F(dimTrait, dimTrait);
        final DenseMatrix64F V2 = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.add(P0, P1, P2);
        safeInvert2(P2, V2, false);
        weightedAverage(M0, P0, M1, P1, M2, V2, dimTrait);
        final WrappedMatrix C2;
        if (NEW_CHOLESKY) {
            DenseMatrix64F tC2 = getCholeskyOfVariance(V2, dimTrait);
            C2 = new WrappedMatrix.WrappedDenseMatrix(tC2);
            MultivariateNormalDistribution.nextMultivariateNormalCholesky(// input mean
            M2, // input variance
            C2, // input variance
            1.0, new WrappedVector.Raw(sample, offsetSample, dimTrait), tmpEpsilon);
        } else {
            double[][] tC2 = getCholeskyOfVariance(V2.getData(), dimTrait);
            C2 = new WrappedMatrix.ArrayOfArray(tC2);
            MultivariateNormalDistribution.nextMultivariateNormalCholesky(// input mean
            M2, // input variance
            C2, // input variance
            1.0, new WrappedVector.Raw(sample, offsetSample, dimTrait), tmpEpsilon);
        }
        if (DEBUG) {
            System.err.println("sT F I N");
            System.err.println("M0: " + M0);
            System.err.println("P0: " + P0);
            System.err.println("M1: " + M1);
            System.err.println("P1: " + P1);
            System.err.println("M2: " + M2);
            System.err.println("V2: " + V2);
            System.err.println("C2: " + C2);
            System.err.println("SS: " + new WrappedVector.Raw(sample, offsetSample, dimTrait));
            System.err.println("");
            if (!check(M2)) {
                System.exit(-1);
            }
        }
    } else {
        System.arraycopy(sample, offsetParent, sample, offsetSample, dimTrait);
        if (DEBUG) {
            System.err.println("sT F I N infinite branch precision");
            System.err.println("SS: " + new WrappedVector.Raw(sample, offsetSample, dimTrait));
        }
    }
}