use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk-protected by broadinstitute.
the class NormalizeSomaticReadCountsIntegrationTest method assertBetaHatsRobustToOutliers.
/**
* Asserts that the calculation of beta hats is not significantly affected by zero-coverage outlier counts
* We perform this check by randomly setting some coverages to zero in copy ratio space (-infinity in log space).
* betaHats imputes 0 in log space (1 in copy ratio space) whenever coverage is below a certain low threshold
* and should thus be robust to this type of noise.
*/
private void assertBetaHatsRobustToOutliers(final ReadCountCollection preTangentNormalized, final File ponFile) {
try (final HDF5File ponReader = new HDF5File(ponFile)) {
final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
final List<String> ponTargets = pon.getPanelTargetNames();
final RealMatrix input = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);
// randomly set some entries to zero in copy-ratio space (-infinity in log space)
final Random random = new Random(13);
final double noiseProportion = 0.01;
final RealMatrix noisyInput = input.copy();
noisyInput.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(final int row, final int column, final double value) {
return random.nextDouble() < noiseProportion ? Double.NEGATIVE_INFINITY : value;
}
});
final RealMatrix betaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), input, PCATangentNormalizationUtils.EPSILON);
final RealMatrix noisyBetaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), noisyInput, PCATangentNormalizationUtils.EPSILON);
final RealMatrix difference = betaHats.subtract(noisyBetaHats);
difference.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
@Override
public void visit(final int row, int column, double value) {
Assert.assertEquals(value, 0, 0.01);
}
});
}
}
use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk by broadinstitute.
the class CopyNumberTriStateTransitionProbabilityCacheUnitTest method markovianPropertiesTest.
//test various properties of a transition matrix
@Test(dataProvider = "meanEventSizeAndEventStartProbability")
public void markovianPropertiesTest(final double meanEventSize, final double eventStartProbability) {
final CopyNumberTriStateTransitionProbabilityCache cache = new CopyNumberTriStateTransitionProbabilityCache(meanEventSize, eventStartProbability);
for (final int d : DISTANCES) {
//check symmetries -- these are part of the model and need not be true in the future
final RealMatrix transitionMatrix = cache.getAsMatrixInProbabilitySpace(d);
assertSymmetries(transitionMatrix);
//check that columns sums equal 1
for (int column = 0; column < transitionMatrix.getColumnDimension(); column++) {
Assert.assertEquals(MathUtils.sum(transitionMatrix.getColumn(column)), 1, EPSILON);
}
//check that all elements are positive
transitionMatrix.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
@Override
public void visit(int row, int column, double value) {
Assert.assertTrue(value >= 0);
}
});
//check that T(2d) = T(d)*T(d)
assertEqualMatrices(cache.getAsMatrixInProbabilitySpace(2 * d), transitionMatrix.multiply(transitionMatrix));
//check that the largest eigenvalue of the transition matrix is 1 (this corresponds to the asymptotic stationary state)
Assert.assertEquals(MathUtils.arrayMax(new EigenDecomposition(transitionMatrix).getRealEigenvalues()), 1, EPSILON);
}
// check that at long distances memory of the initial state is lost and all initial distributions tend toward
// the same asymptotic stationary distribution. That is, all columns of the large-distance transition matrix are equal
final RealMatrix asymptoticMatrix = cache.getAsMatrixInProbabilitySpace(HUGE_DISTANCE);
for (int column = 1; column < asymptoticMatrix.getColumnDimension(); column++) {
Assert.assertEquals(asymptoticMatrix.getColumnVector(0).subtract(asymptoticMatrix.getColumnVector(column)).getL1Norm(), 0, EPSILON);
}
}
use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk by broadinstitute.
the class NormalizeSomaticReadCountsIntegrationTest method assertBetaHatsRobustToOutliers.
/**
* Asserts that the calculation of beta hats is not significantly affected by zero-coverage outlier counts
* We perform this check by randomly setting some coverages to zero in copy ratio space (-infinity in log space).
* betaHats imputes 0 in log space (1 in copy ratio space) whenever coverage is below a certain low threshold
* and should thus be robust to this type of noise.
*/
private void assertBetaHatsRobustToOutliers(final ReadCountCollection preTangentNormalized, final File ponFile) {
try (final HDF5File ponReader = new HDF5File(ponFile)) {
final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
final List<String> ponTargets = pon.getPanelTargetNames();
final RealMatrix input = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);
// randomly set some entries to zero in copy-ratio space (-infinity in log space)
final Random random = new Random(13);
final double noiseProportion = 0.01;
final RealMatrix noisyInput = input.copy();
noisyInput.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(final int row, final int column, final double value) {
return random.nextDouble() < noiseProportion ? Double.NEGATIVE_INFINITY : value;
}
});
final RealMatrix betaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), input, PCATangentNormalizationUtils.EPSILON);
final RealMatrix noisyBetaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), noisyInput, PCATangentNormalizationUtils.EPSILON);
final RealMatrix difference = betaHats.subtract(noisyBetaHats);
difference.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
@Override
public void visit(final int row, int column, double value) {
Assert.assertEquals(value, 0, 0.01);
}
});
}
}
use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk-protected by broadinstitute.
the class CopyNumberTriStateTransitionProbabilityCacheUnitTest method markovianPropertiesTest.
//test various properties of a transition matrix
@Test(dataProvider = "meanEventSizeAndEventStartProbability")
public void markovianPropertiesTest(final double meanEventSize, final double eventStartProbability) {
final CopyNumberTriStateTransitionProbabilityCache cache = new CopyNumberTriStateTransitionProbabilityCache(meanEventSize, eventStartProbability);
for (final int d : DISTANCES) {
//check symmetries -- these are part of the model and need not be true in the future
final RealMatrix transitionMatrix = cache.getAsMatrixInProbabilitySpace(d);
assertSymmetries(transitionMatrix);
//check that columns sums equal 1
for (int column = 0; column < transitionMatrix.getColumnDimension(); column++) {
Assert.assertEquals(MathUtils.sum(transitionMatrix.getColumn(column)), 1, EPSILON);
}
//check that all elements are positive
transitionMatrix.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
@Override
public void visit(int row, int column, double value) {
Assert.assertTrue(value >= 0);
}
});
//check that T(2d) = T(d)*T(d)
assertEqualMatrices(cache.getAsMatrixInProbabilitySpace(2 * d), transitionMatrix.multiply(transitionMatrix));
//check that the largest eigenvalue of the transition matrix is 1 (this corresponds to the asymptotic stationary state)
Assert.assertEquals(MathUtils.arrayMax(new EigenDecomposition(transitionMatrix).getRealEigenvalues()), 1, EPSILON);
}
// check that at long distances memory of the initial state is lost and all initial distributions tend toward
// the same asymptotic stationary distribution. That is, all columns of the large-distance transition matrix are equal
final RealMatrix asymptoticMatrix = cache.getAsMatrixInProbabilitySpace(HUGE_DISTANCE);
for (int column = 1; column < asymptoticMatrix.getColumnDimension(); column++) {
Assert.assertEquals(asymptoticMatrix.getColumnVector(0).subtract(asymptoticMatrix.getColumnVector(column)).getL1Norm(), 0, EPSILON);
}
}
Aggregations