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

/**
 * Write segment file based on maximum-likelihood estimates of the minor allele fraction at SNP sites,
 * assuming the specified allelic bias.  These estimates are converted to target coverages,
 * which are written to a temporary file and then passed to {@link RCBSSegmenter}.
 * @param snps                  TargetCollection of allelic counts at SNP sites
 * @param sampleName            sample name
 * @param outputFile            segment file to write to and return
 * @param allelicBias           allelic bias to use in estimate of minor allele fraction
 */
public static void writeSegmentFile(final TargetCollection<AllelicCount> snps, final String sampleName,
                                    final File outputFile, final double allelicBias) {
    Utils.validateArg(snps.totalSize() > 0, "Must have a positive number of SNPs to perform SNP segmentation.");
    try {
        final File targetsFromSNPCountsFile = File.createTempFile("targets-from-snps", ".tsv");

        final List<Target> targets = snps.targets().stream()
                .map(ac -> new Target(name(ac), ac.getInterval())).collect(Collectors.toList());

        final RealMatrix minorAlleleFractions = new Array2DRowRealMatrix(snps.targetCount(), 1);
        minorAlleleFractions.setColumn(0, snps.targets().stream()
                .mapToDouble(ac -> ac.estimateMinorAlleleFraction(allelicBias)).toArray());

        ReadCountCollectionUtils.write(targetsFromSNPCountsFile, new ReadCountCollection(targets, Collections.singletonList(sampleName), minorAlleleFractions));

        //segment SNPs based on observed log_2 minor allele fraction (log_2 is applied in CBS.R)
        RCBSSegmenter.writeSegmentFile(sampleName, targetsFromSNPCountsFile.getAbsolutePath(),
                outputFile.getAbsolutePath(), false);
    } catch (final IOException e) {
        throw new UserException.CouldNotCreateOutputFile("Could not create temporary output file during " +
                "SNP segmentation.", e);
    }
}
 

/**
 * Subsets the count columns in the read-count collection.
 *
 * <p>
 *     Creates a brand-new read-count collection. Changes in the new instance won't affect this one and vice-versa.
 * </p>
 *
 * @param columnsToKeep column names to keep in the result read-count collection.
 * @return never {@code null}.
 */
public ReadCountCollection subsetColumns(final Set<String> columnsToKeep) {
    Utils.nonNull(columnsToKeep, "the set of input columns to keep cannot be null.");
    Utils.nonEmpty(columnsToKeep, "the number of columns to keep must be greater than 0");
    if (!new HashSet<>(columnNames).containsAll(columnsToKeep)) {
        throw unknownColumnToKeepNames(columnsToKeep);
    }

    if (columnsToKeep.size() == columnNames.size())  {
        return new ReadCountCollection(targets, columnNames, counts.copy(), false);
    }

    final int[] columnsToKeepIndices = IntStream.range(0, columnNames.size())
            .filter(i -> columnsToKeep.contains(columnNames.get(i))).toArray();
    final List<String> resultColumnNames = Arrays.stream(columnsToKeepIndices).mapToObj(columnNames::get).collect(Collectors.toList());

    final RealMatrix resultCountsM = new Array2DRowRealMatrix(counts.getRowDimension(), columnsToKeepIndices.length);
    for (int i = 0; i < columnsToKeepIndices.length; i++) {
        resultCountsM.setColumn(i, counts.getColumn(columnsToKeepIndices[i]));
    }

    return new ReadCountCollection(targets, Collections.unmodifiableList(resultColumnNames), resultCountsM, false);
}
 

@Override
public void update(@Nonnull final Object arg, @Nonnull final double[] outScores)
        throws HiveException {
    double x = PrimitiveObjectInspectorUtils.getDouble(arg, oi);
    xRing.add(x).toArray(xSeries, true /* FIFO */);

    // need to wait until the buffer is filled
    if (!xRing.isFull()) {
        outScores[0] = 0.d;
    } else {
        // create past trajectory matrix and find its left singular vectors
        RealMatrix H = new Array2DRowRealMatrix(window, nPastWindow);
        for (int i = 0; i < nPastWindow; i++) {
            H.setColumn(i, Arrays.copyOfRange(xSeries, i, i + window));
        }

        // create current trajectory matrix and find its left singular vectors
        RealMatrix G = new Array2DRowRealMatrix(window, nCurrentWindow);
        int currentHead = pastSize + currentOffset;
        for (int i = 0; i < nCurrentWindow; i++) {
            G.setColumn(i,
                Arrays.copyOfRange(xSeries, currentHead + i, currentHead + i + window));
        }

        switch (scoreFunc) {
            case svd:
                outScores[0] = computeScoreSVD(H, G);
                break;
            case ika:
                outScores[0] = computeScoreIKA(H, G);
                break;
            default:
                throw new IllegalStateException("Unexpected score function: " + scoreFunc);
        }
    }
}
 

private static Genome makeGenome(final List<List<ReadCountRecord.SingleSampleRecord>> targetCoverages,
                                 final List<List<AllelicCount>> snpCounts) {
    final List<ReadCountRecord.SingleSampleRecord> records = new ArrayList<>();
    final List<AllelicCount> snps = new ArrayList<>();
    targetCoverages.stream().forEach(records::addAll);
    snpCounts.stream().forEach(snps::addAll);

    final List<Target> targets = records.stream().map(ReadCountRecord::getTarget).collect(Collectors.toList());
    final RealMatrix counts = new Array2DRowRealMatrix(targets.size(), 1);
    counts.setColumn(0, records.stream().mapToDouble(r -> r.getCount()).toArray());
    return new Genome(new ReadCountCollection(targets, Arrays.asList("sample"), counts), snps);
}
 

@Test
public void testBasicSoundness() {
    for (final RealMatrix transitionMatrix : TRANSITION_MATRICES) {
        final IntegerCopyNumberTransitionProbabilityCache cache = new IntegerCopyNumberTransitionProbabilityCache(
                new IntegerCopyNumberTransitionMatrix(transitionMatrix, 0));
        for (final int dist : DISTANCES) {
            final RealMatrix transitionMatrixExponentiated = cache.getTransitionProbabilityMatrix(dist);

            /* assert positivity */
            Assert.assertTrue(Arrays.stream(transitionMatrixExponentiated.getData())
                    .flatMapToDouble(Arrays::stream)
                    .allMatch(d -> d >= 0));

            /* assert conservation of probability */
            for (int c = 0; c < transitionMatrix.getColumnDimension(); c++) {
                Assert.assertEquals(Arrays.stream(transitionMatrixExponentiated.getColumn(c)).sum(), 1.0, EPSILON);
            }

            /* assert correctness, T(2*d) = T(d)*T(d) */
            assertEqualMatrices(cache.getTransitionProbabilityMatrix(2*dist),
                    transitionMatrixExponentiated.multiply(transitionMatrixExponentiated));
        }

        /* assert loss of initial state over long distances, i.e. all columns must be equal */
        final RealMatrix longRangeTransitionMatrix = cache.getTransitionProbabilityMatrix(Integer.MAX_VALUE);
        final double[] firstColumn = longRangeTransitionMatrix.getColumn(0);
        final RealMatrix syntheticLongRangeTransitionMatrix = new Array2DRowRealMatrix(firstColumn.length,
                firstColumn.length);
        for (int i = 0; i < firstColumn.length; i++) {
            syntheticLongRangeTransitionMatrix.setColumn(i, firstColumn);
        }
        assertEqualMatrices(longRangeTransitionMatrix, syntheticLongRangeTransitionMatrix);

        final double[] stationary = cache.getStationaryProbabilityVector().toArray();
        ArrayAsserts.assertArrayEquals(stationary, firstColumn, EPSILON);
    }
}
 

private void computeADFStatistics() {
	double[] y = diff(ts);
	RealMatrix designMatrix = null;
	int k = lag+1;
	int n = ts.length - 1;
	
	RealMatrix z = MatrixUtils.createRealMatrix(laggedMatrix(y, k)); //has rows length(ts) - 1 - k + 1
	RealVector zcol1 = z.getColumnVector(0); //has length length(ts) - 1 - k + 1
	double[] xt1 = subsetArray(ts, k-1, n-1);  //ts[k:(length(ts) - 1)], has length length(ts) - 1 - k + 1
	double[] trend = sequence(k,n); //trend k:n, has length length(ts) - 1 - k + 1
	if (k > 1) {
		RealMatrix yt1 = z.getSubMatrix(0, ts.length - 1 - k, 1, k-1); //same as z but skips first column
		//build design matrix as cbind(xt1, 1, trend, yt1)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3 + k - 1);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
		designMatrix.setSubMatrix(yt1.getData(), 0, 3);
		
	} else {
		//build design matrix as cbind(xt1, 1, tt)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
	}
	/*OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
	regression.setNoIntercept(true);
	regression.newSampleData(zcol1.toArray(), designMatrix.getData());
	double[] beta = regression.estimateRegressionParameters();
	double[] sd = regression.estimateRegressionParametersStandardErrors();
	*/
	RidgeRegression regression = new RidgeRegression(designMatrix.getData(), zcol1.toArray());
	regression.updateCoefficients(.0001);
	double[] beta = regression.getCoefficients();
	double[] sd = regression.getStandarderrors();
	
	double t = beta[0] / sd[0];
	if (t <= PVALUE_THRESHOLD) {
		this.needsDiff = true;
	} else {
		this.needsDiff = false;
	}	
}
 

/**
 * Applies rank transform to each of the columns of <code>matrix</code>
 * using the current <code>rankingAlgorithm</code>
 *
 * @param matrix matrix to transform
 */
private void rankTransform(RealMatrix matrix) {
    for (int i = 0; i < matrix.getColumnDimension(); i++) {
        matrix.setColumn(i, rankingAlgorithm.rank(matrix.getColumn(i)));
    }
}
 

/**
 * Applies rank transform to each of the columns of <code>matrix</code>
 * using the current <code>rankingAlgorithm</code>
 *
 * @param matrix matrix to transform
 */
private void rankTransform(RealMatrix matrix) {
    for (int i = 0; i < matrix.getColumnDimension(); i++) {
        matrix.setColumn(i, rankingAlgorithm.rank(matrix.getColumn(i)));
    }
}
 

/**
 * Applies rank transform to each of the columns of <code>matrix</code>
 * using the current <code>rankingAlgorithm</code>
 *
 * @param matrix matrix to transform
 */
private void rankTransform(RealMatrix matrix) {
    for (int i = 0; i < matrix.getColumnDimension(); i++) {
        matrix.setColumn(i, rankingAlgorithm.rank(matrix.getColumn(i)));
    }
}
 

/**
 * Applies rank transform to each of the columns of <code>matrix</code>
 * using the current <code>rankingAlgorithm</code>
 *
 * @param matrix matrix to transform
 */
private void rankTransform(RealMatrix matrix) {
    for (int i = 0; i < matrix.getColumnDimension(); i++) {
        matrix.setColumn(i, rankingAlgorithm.rank(matrix.getColumn(i)));
    }
}