use of org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor in project gatk by broadinstitute.
the class XHMMSegmentCallerBase method standardizeBySample.
/**
* Standardize read counts (per-sample).
* Note: modification is done in-place.
*
* @param counts original read counts
*/
private void standardizeBySample(final RealMatrix counts) {
final double[] columnMeans = GATKProtectedMathUtils.columnMeans(counts);
final double[] columnStdDev = GATKProtectedMathUtils.columnStdDevs(counts);
counts.walkInColumnOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(int row, int column, double value) {
return (value - columnMeans[column]) / columnStdDev[column];
}
});
}
use of org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor in project gatk by broadinstitute.
the class HDF5PCACoveragePoNCreationUtils method normalizeAndLogReadCounts.
/**
* Final pre-panel normalization that consists of dividing all counts by the median of
* its column and log it with base 2.
* <p>
* The normalization occurs in-place.
* </p>
*
* @param readCounts the input counts to normalize.
*/
@VisibleForTesting
static void normalizeAndLogReadCounts(final ReadCountCollection readCounts, final Logger logger) {
final RealMatrix counts = readCounts.counts();
final Median medianCalculator = new Median();
final double[] medians = IntStream.range(0, counts.getColumnDimension()).mapToDouble(col -> medianCalculator.evaluate(counts.getColumn(col))).toArray();
counts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(final int row, final int column, final double value) {
return Math.log(Math.max(EPSILON, value / medians[column])) * INV_LN_2;
}
});
logger.info("Counts normalized by the column median and log2'd.");
}
use of org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor in project gatk by broadinstitute.
the class HDF5PCACoveragePoNCreationUtils method subtractMedianOfMedians.
/**
* Calculates the median of column medians and subtract it from all counts.
* @param readCounts the input counts to center.
* @return the median of medians that has been subtracted from all counts.
*/
@VisibleForTesting
static double subtractMedianOfMedians(final ReadCountCollection readCounts, final Logger logger) {
final RealMatrix counts = readCounts.counts();
final Median medianCalculator = new Median();
final double[] columnMedians = MatrixSummaryUtils.getColumnMedians(counts);
final double medianOfMedians = medianCalculator.evaluate(columnMedians);
counts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(final int row, final int column, final double value) {
return value - medianOfMedians;
}
});
logger.info(String.format("Counts centered around the median of medians %.2f", medianOfMedians));
return medianOfMedians;
}
use of org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor in project gatk by broadinstitute.
the class HDF5PCACoveragePoNCreationUtils method calculateReducedPanelAndPInverses.
/**
* SVD and Pseudo inverse calculation.
*
* @param logNormalized the input counts for the SVD and reduction steps, fully normalized and already logged.
* @param requestedNumberOfEigensamples user requested number of eigensamples for the reduced panel.
* @return never {@code null}.
*/
@VisibleForTesting
static ReductionResult calculateReducedPanelAndPInverses(final ReadCountCollection logNormalized, final OptionalInt requestedNumberOfEigensamples, final Logger logger, final JavaSparkContext ctx) {
if (ctx == null) {
logger.warn("No Spark context provided, not going to use Spark...");
}
final RealMatrix logNormalizedCounts = logNormalized.counts();
final int numberOfCountColumns = logNormalizedCounts.getColumnDimension();
logger.info("Starting the SVD decomposition of the log-normalized counts ...");
final long svdStartTime = System.currentTimeMillis();
final SVD logNormalizedSVD = SVDFactory.createSVD(logNormalized.counts(), ctx);
final long svdEndTime = System.currentTimeMillis();
logger.info(String.format("Finished the SVD decomposition of the log-normal counts. Elapse of %d seconds", (svdEndTime - svdStartTime) / 1000));
final int numberOfEigensamples = determineNumberOfEigensamples(requestedNumberOfEigensamples, numberOfCountColumns, logNormalizedSVD, logger);
logger.info(String.format("Including %d eigensamples in the reduced PoN", numberOfEigensamples));
final double[] singularValues = logNormalizedSVD.getSingularValues();
final RealMatrix reducedCounts = logNormalizedSVD.getU().getSubMatrix(0, logNormalizedCounts.getRowDimension() - 1, 0, numberOfEigensamples - 1).copy();
reducedCounts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(final int row, final int column, final double value) {
return singularValues[column] * value;
}
});
// The Pseudo inverse comes nearly for free from having run the SVD decomposition.
final RealMatrix logNormalizedPseudoInverse = logNormalizedSVD.getPinv();
logger.info("Calculating the reduced PoN inverse matrix...");
final long riStartTime = System.currentTimeMillis();
final RealMatrix reducedCountsPseudoInverse = SVDFactory.createSVD(reducedCounts, ctx).getPinv();
final long riEndTime = System.currentTimeMillis();
logger.info(String.format("Finished calculating the reduced PoN inverse matrix. Elapse of %d seconds", (riEndTime - riStartTime) / 1000));
return new ReductionResult(logNormalizedPseudoInverse, reducedCounts, reducedCountsPseudoInverse, logNormalizedSVD.getSingularValues());
}
use of org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor in project gatk by broadinstitute.
the class PCATangentNormalizationUtils method factorNormalize.
/**
* Target-factor normalizes a {@link RealMatrix} in-place given target factors..
*/
static void factorNormalize(final RealMatrix input, final double[] targetFactors) {
Utils.nonNull(input, "Input matrix cannot be null.");
Utils.nonNull(targetFactors, "Target factors cannot be null.");
Utils.validateArg(targetFactors.length == input.getRowDimension(), "Number of target factors does not correspond to the number of rows.");
// Divide all counts by the target factor for the row.
input.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(final int row, final int column, final double value) {
return value / targetFactors[row];
}
});
}
Aggregations