Java Code Examples for org.apache.commons.math3.linear.RealMatrix#copy()

@Test
public void testGCCorrection() {
    final int numSamples = 5;
    final int numIntervals = 10000;

    final Pair<RealMatrix, double[]> data = simulateData(numSamples, numIntervals);
    final RealMatrix readCounts = data.getLeft();
    final double[] intervalGCContent = data.getRight();

    final RealMatrix correctedCoverage = readCounts.copy();
    GCBiasCorrector.correctGCBias(correctedCoverage, intervalGCContent);
    Utils.nonNull(correctedCoverage);
    final StandardDeviation stdevEvaluator = new StandardDeviation();
    final double[] stdDevsBySample = IntStream.range(0, correctedCoverage.getRowDimension())
            .mapToDouble(r -> stdevEvaluator.evaluate(correctedCoverage.getRow(r))).toArray();
    Arrays.stream(stdDevsBySample).forEach(x -> Assert.assertTrue(x < NON_GC_BIAS_NOISE_LEVEL * MEAN_READ_DEPTH));

    //check that GC-bias correction is approximately idempotent -- if you correct again, very little should happen
    final RealMatrix recorrectedCoverage = correctedCoverage.copy();
    GCBiasCorrector.correctGCBias(recorrectedCoverage, intervalGCContent);
    final double correctedChange = correctedCoverage.subtract(readCounts).getFrobeniusNorm();
    final double recorrectedChange = recorrectedCoverage.subtract(correctedCoverage).getFrobeniusNorm();
    Assert.assertTrue(recorrectedChange < correctedChange / 10.);
}
 

/**
 * Calculate the beta-hats that best fit case read counts given the panel of normals.
 *
 * @param normalsPseudoinverse the log-normalized or reduced-panel pseudoinverse from a panel of normals
 * @param input a {@code TxS} matrix where {@code T} is the number of targets and {@code S} the number of count groups (e.g. case samples).
 * @return never {@code null} an {@code NxS} matrix, where N is the number of samples in
 *  the panel and S the original name of count groups.
 */
@VisibleForTesting
public static RealMatrix calculateBetaHats(final RealMatrix normalsPseudoinverse,
                                           final RealMatrix input,
                                           final double epsilon) {
    Utils.nonNull(normalsPseudoinverse, "Normals inverse matrix cannot be null.");
    Utils.nonNull(input, "Input counts cannot be null.");
    Utils.validateArg(epsilon > 0, String.format("Invalid epsilon value, must be > 0: %f", epsilon));
    final double targetThreshold = (Math.log(epsilon) / Math.log(2)) + 1;

    // copy case samples in order to mask targets in-place and mask (set to zero) targets with coverage below threshold
    final RealMatrix maskedInput = input.copy();
    maskedInput.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int row, final int column, final double value) {
            return value > targetThreshold ? value : 0;
        }
    });

    return normalsPseudoinverse.multiply(maskedInput);
}
 

/**
 * @param weight Weight matrix.
 * @throws NonSquareMatrixException if the argument is not
 * a square matrix.
 */
public Weight(RealMatrix weight) {
    if (weight.getColumnDimension() != weight.getRowDimension()) {
        throw new NonSquareMatrixException(weight.getColumnDimension(),
                                           weight.getRowDimension());
    }

    weightMatrix = weight.copy();
}
 

/**
 * @param weight Weight matrix.
 * @throws NonSquareMatrixException if the argument is not
 * a square matrix.
 */
public Weight(RealMatrix weight) {
    if (weight.getColumnDimension() != weight.getRowDimension()) {
        throw new NonSquareMatrixException(weight.getColumnDimension(),
                                           weight.getRowDimension());
    }

    weightMatrix = weight.copy();
}
 

/**
 * @param weight Weight matrix.
 * @throws NonSquareMatrixException if the argument is not
 * a square matrix.
 */
public Weight(RealMatrix weight) {
    if (weight.getColumnDimension() != weight.getRowDimension()) {
        throw new NonSquareMatrixException(weight.getColumnDimension(),
                                           weight.getRowDimension());
    }

    weightMatrix = weight.copy();
}
 

/**
 * @param weight Weight matrix.
 * @throws NonSquareMatrixException if the argument is not
 * a square matrix.
 */
public Weight(RealMatrix weight) {
    if (weight.getColumnDimension() != weight.getRowDimension()) {
        throw new NonSquareMatrixException(weight.getColumnDimension(),
                                           weight.getRowDimension());
    }

    weightMatrix = weight.copy();
}
 

/**
 * @param weight Weight matrix.
 * @throws NonSquareMatrixException if the argument is not
 * a square matrix.
 */
public Weight(RealMatrix weight) {
    if (weight.getColumnDimension() != weight.getRowDimension()) {
        throw new NonSquareMatrixException(weight.getColumnDimension(),
                                           weight.getRowDimension());
    }

    weightMatrix = weight.copy();
}
 

/**
 * Create a SpearmansCorrelation with the given input data matrix
 * and ranking algorithm.
 *
 * @param dataMatrix matrix of data with columns representing
 * variables to correlate
 * @param rankingAlgorithm ranking algorithm
 */
public SpearmansCorrelation(final RealMatrix dataMatrix, final RankingAlgorithm rankingAlgorithm) {
    this.data = dataMatrix.copy();
    this.rankingAlgorithm = rankingAlgorithm;
    rankTransform(data);
    rankCorrelation = new PearsonsCorrelation(data);
}
 

/**
 * Builds a simple converter for correlated residuals with the
 * specified weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 *
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations,
                             final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/** Build a simple converter for correlated residuals with the specific weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations, final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/** Build a simple converter for correlated residuals with the specific weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations, final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/** Build a simple converter for correlated residuals with the specific weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations, final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/** Build a simple converter for correlated residuals with the specific weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations, final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/**
 * Create a SpearmansCorrelation with the given input data matrix
 * and ranking algorithm.
 *
 * @param dataMatrix matrix of data with columns representing
 * variables to correlate
 * @param rankingAlgorithm ranking algorithm
 */
public SpearmansCorrelation(final RealMatrix dataMatrix, final RankingAlgorithm rankingAlgorithm) {
    this.data = dataMatrix.copy();
    this.rankingAlgorithm = rankingAlgorithm;
    rankTransform(data);
    rankCorrelation = new PearsonsCorrelation(data);
}
 

/**
 * Builds a simple converter for correlated residuals with the
 * specified weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 *
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations,
                             final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/** Build a simple converter for correlated residuals with the specific weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations, final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/**
 * Builds a simple converter for correlated residuals with the
 * specified weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
 * will be triggered while computing the scalar objective.
 * </p>
 *
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @throws DimensionMismatchException if the observations vector and the scale
 * matrix dimensions do not match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorFunction function,
                             final double[] observations,
                             final RealMatrix scale) {
    if (observations.length != scale.getColumnDimension()) {
        throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}
 

/**
 * Computes the Spearman's rank correlation matrix for the columns of the
 * input matrix.
 *
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    RealMatrix matrixCopy = matrix.copy();
    rankTransform(matrixCopy);
    return new PearsonsCorrelation().computeCorrelationMatrix(matrixCopy);
}
 

/**
 * Computes the Spearman's rank correlation matrix for the columns of the
 * input matrix.
 *
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    RealMatrix matrixCopy = matrix.copy();
    rankTransform(matrixCopy);
    return new PearsonsCorrelation().computeCorrelationMatrix(matrixCopy);
}
 

/**
 * Computes the Spearman's rank correlation matrix for the columns of the
 * input matrix.
 *
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    RealMatrix matrixCopy = matrix.copy();
    rankTransform(matrixCopy);
    return new PearsonsCorrelation().computeCorrelationMatrix(matrixCopy);
}