use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class IntegerCopyNumberTransitionMatrixUnitTest method testPadding.
@Test
public void testPadding() {
final IntegerCopyNumberTransitionMatrix data = new IntegerCopyNumberTransitionMatrix(new Array2DRowRealMatrix(new double[][] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } }), 2);
final RealMatrix expected = new Array2DRowRealMatrix(new double[][] { { 1.0 / 12, 2.0 / 15, 3.0 / 18, 0, 0 }, { 4.0 / 12, 5.0 / 15, 6.0 / 18, 0, 0 }, { 7.0 / 12, 8.0 / 15, 9.0 / 18, 0, 0 }, { 0, 0, 0, 1, 0 }, { 0, 0, 0, 0, 1 } });
Assert.assertEquals(data.getTransitionMatrix().subtract(expected).getNorm(), 0, 1e-12);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk-protected by broadinstitute.
the class SegmentUtilsUnitTest method testUnionSegments.
/**
* Test for {@link SegmentUtils#unionSegments}. Expected behavior:
* <p>
* On chr1 {@link SegmentUtils#collectBreakpointsByContig} gives:
* </p>
* <p>
* 1, 5, 10, 20, 40, 40, 42, 90, 91, 115, 125, 140.
* </p>
* <p>
* Then {@link SegmentUtils#constructUntrimmedSegments} finds the segments:
* </p>
* <p>
* [1, 4], [5, 10], [11, 19], [20, 40], [41, 41], [42, 89], [90, 91], [92, 114], [115, 125], [126, 140].
* </p>
* <p>
* and returns the non-empty segments:
* </p>
* <p>
* [1, 4], [5, 10], [20, 40], [42, 89], [90, 91], [92, 114], [115, 125], [126, 140].
* </p>
* <p>
* Then {@link SegmentUtils#mergeSpuriousStartsAndEnds} merges the last segment left to form [115, 140],
* and {@link SegmentMergeUtils#mergeSpuriousMiddles} randomly merges segment [92, 114] left or right.
* </p>
* <p>
* Finally, {@link SegmentUtils#trimInterval} gives:
* </p>
* <p>
* [1, 10], [20, 40], [42, 42], [90, 114], [115, 140] (if [92, 114] merged left) or
* </p>
* <p>
* [1, 10], [20, 40], [42, 42], [90, 91], [92, 140] (if [92, 114] merged right)
* </p>
* <p>
* The remaining empty segment on chr2 is retained.
*/
@Test
public void testUnionSegments() {
final String sampleName = "placeholder_sample_name";
final List<Target> targets = new ArrayList<Target>();
targets.add(new Target("t1", new SimpleInterval("chr1", 1, 10)));
targets.add(new Target("t2", new SimpleInterval("chr1", 20, 30)));
targets.add(new Target("t3", new SimpleInterval("chr1", 31, 40)));
targets.add(new Target("t4", new SimpleInterval("chr1", 90, 100)));
targets.add(new Target("t5", new SimpleInterval("chr1", 110, 120)));
targets.add(new Target("t6", new SimpleInterval("chr1", 130, 140)));
final RealMatrix zeroCoverageMatrix = new Array2DRowRealMatrix(targets.size(), 1);
final ReadCountCollection counts = new ReadCountCollection(targets, Collections.singletonList(sampleName), zeroCoverageMatrix);
final AllelicCount snp1 = new AllelicCount(new SimpleInterval("chr1", 5, 5), 0, 1);
final AllelicCount snp2 = new AllelicCount(new SimpleInterval("chr1", 40, 40), 0, 1);
final AllelicCount snp3 = new AllelicCount(new SimpleInterval("chr1", 42, 42), 0, 1);
final AllelicCount snp4 = new AllelicCount(new SimpleInterval("chr1", 91, 91), 0, 1);
final AllelicCount snp5 = new AllelicCount(new SimpleInterval("chr1", 115, 115), 0, 1);
final AllelicCount snp6 = new AllelicCount(new SimpleInterval("chr1", 125, 125), 0, 1);
final AllelicCount snp7 = new AllelicCount(new SimpleInterval("chr2", 10, 10), 0, 1);
final List<AllelicCount> snps = Arrays.asList(snp1, snp2, snp3, snp4, snp5, snp6, snp7);
final List<SimpleInterval> targetSegments = Arrays.asList(new SimpleInterval("chr1", 1, 10), new SimpleInterval("chr1", 20, 40), new SimpleInterval("chr1", 90, 140));
final List<SimpleInterval> snpSegments = Arrays.asList(new SimpleInterval("chr1", 5, 40), new SimpleInterval("chr1", 42, 91), new SimpleInterval("chr1", 115, 125), new SimpleInterval("chr2", 10, 10));
final List<SimpleInterval> unionedSegments = SegmentUtils.unionSegments(targetSegments, snpSegments, new Genome(counts, snps));
final List<SimpleInterval> expectedLeft = Arrays.asList(new SimpleInterval("chr1", 1, 10), new SimpleInterval("chr1", 20, 40), new SimpleInterval("chr1", 42, 42), new SimpleInterval("chr1", 90, 114), new SimpleInterval("chr1", 115, 140), new SimpleInterval("chr2", 10, 10));
final List<SimpleInterval> expectedRight = Arrays.asList(new SimpleInterval("chr1", 1, 10), new SimpleInterval("chr1", 20, 40), new SimpleInterval("chr1", 42, 42), new SimpleInterval("chr1", 90, 91), new SimpleInterval("chr1", 92, 140), new SimpleInterval("chr2", 10, 10));
Assert.assertTrue(unionedSegments.equals(expectedLeft) || unionedSegments.equals(expectedRight));
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class SparkConverter method convertSparkMatrixToRealMatrix.
/**
* Convert a local (not distributed) Spark Matrix to an Apache Commons matrix.
*
* @param r Never {@code null}
* @return Not {@code null}
*/
public static RealMatrix convertSparkMatrixToRealMatrix(final Matrix r) {
final RealMatrix result = new Array2DRowRealMatrix(r.numRows(), r.numCols());
final double[] columnMajorMat = r.toArray();
for (int i = 0; i < r.numRows(); i++) {
result.setRow(i, Arrays.copyOfRange(columnMajorMat, i * r.numCols(), i * r.numCols() + r.numCols()));
}
return result;
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class OjAlgoSingularValueDecomposer method createSVD.
/**
* Create a SVD instance using ojAlgo.
*
* @param m matrix that is not {@code null}
* @return SVD instance that is never {@code null}
*/
@Override
public SVD createSVD(final RealMatrix m) {
Utils.nonNull(m, "Cannot create SVD on a null matrix.");
logger.info("Converting Apache RealMatrix to ojAlgo Matrix...");
final OjAlgoAdapter pds = OjAlgoAdapter.of(m);
logger.info("Calculating SVD...");
final SingularValue<Double> svd = SingularValue.make(pds);
svd.compute(pds);
// If this code is too slow, conversion could be implemented by backing a RealMatrix with an ojAlgo Array2D
logger.info("Converting ojAlgo Matrix to Apache RealMatrix...");
final RealMatrix pinv = new Array2DRowRealMatrix(svd.getInverse().toRawCopy2D());
final RealMatrix u = new Array2DRowRealMatrix(svd.getQ1().toRawCopy2D());
final RealMatrix v = new Array2DRowRealMatrix(svd.getQ2().toRawCopy2D());
final double[] s = svd.getSingularValues().toRawCopy1D();
return new SimpleSVD(u, s, v, pinv);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class CopyRatioSegmenterUnitTest method testSegmentation.
@Test
public void testSegmentation() {
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
final List<Double> trueLog2CopyRatios = Arrays.asList(-2.0, 0.0, 1.4);
final double trueMemoryLength = 1e5;
final double trueStandardDeviation = 0.2;
final CopyRatioHMM trueModel = new CopyRatioHMM(trueLog2CopyRatios, trueWeights, trueMemoryLength, trueStandardDeviation);
final int chainLength = 10000;
final List<SimpleInterval> positions = randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
final List<Integer> trueStates = trueModel.generateHiddenStateChain(positions);
final List<Double> trueLog2CopyRatioSequence = trueStates.stream().map(n -> trueLog2CopyRatios.get(n)).collect(Collectors.toList());
final List<Double> data = trueLog2CopyRatioSequence.stream().map(cr -> generateData(trueStandardDeviation, cr, rng)).collect(Collectors.toList());
final List<Target> targets = positions.stream().map(Target::new).collect(Collectors.toList());
final ReadCountCollection rcc = new ReadCountCollection(targets, Arrays.asList("SAMPLE"), new Array2DRowRealMatrix(data.stream().mapToDouble(x -> x).toArray()));
final CopyRatioSegmenter segmenter = new CopyRatioSegmenter(10, rcc);
final List<ModeledSegment> segments = segmenter.getModeledSegments();
final double[] segmentCopyRatios = segments.stream().flatMap(s -> Collections.nCopies((int) s.getTargetCount(), s.getSegmentMeanInLog2CRSpace()).stream()).mapToDouble(x -> x).toArray();
final double averageCopyRatioError = IntStream.range(0, trueLog2CopyRatioSequence.size()).mapToDouble(n -> Math.abs(segmentCopyRatios[n] - trueLog2CopyRatioSequence.get(n))).average().getAsDouble();
Assert.assertEquals(averageCopyRatioError, 0, 0.025);
}
Aggregations