Java Code Examples for org.apache.commons.math3.stat.descriptive.rank.Median#evaluate()

/**
 * 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;
}
 

@Test(dataProvider = "readCountOnlyData")
public void testSubtractMedianOfMedians(final ReadCountCollection readCounts) {
    final RealMatrix counts = readCounts.counts();
    final Median median = new Median();
    final double[] columnMedians = IntStream.range(0, counts.getColumnDimension())
            .mapToDouble(i -> median.evaluate(counts.getColumn(i))).toArray();
    final double center = median.evaluate(columnMedians);
    final double[][] expected = new double[counts.getRowDimension()][];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = counts.getRow(i).clone();
        for (int j = 0; j < expected[i].length; j++) {
            expected[i][j] -= center;
        }
    }
    HDF5PCACoveragePoNCreationUtils.subtractMedianOfMedians(readCounts, NULL_LOGGER);
    final RealMatrix newCounts = readCounts.counts();
    Assert.assertEquals(newCounts.getColumnDimension(), expected[0].length);
    Assert.assertEquals(newCounts.getRowDimension(), expected.length);
    for (int i = 0; i < expected.length; i++) {
        for (int j = 0; j < expected[i].length; j++) {
            Assert.assertEquals(newCounts.getEntry(i, j), expected[i][j], 0.000001);
        }
    }
}
 

public static void main(String[] args) {
    final double[] values = new java.util.Random().doubles(5000).toArray();
    final Percentile stat1 = new Percentile(0.5);
    final Median stat2 = new Median();
    for (double value : values) stat1.add(value);
    final double result1 = stat1.getValue();
    final double result2 = stat2.evaluate(values);
    if (result1 != result2) {
        throw new RuntimeException("Error: " + result1 + " != " + result2);
    }
}
 

public static double getMedian(float[] values) {
    double[] array = new double[values.length];
    for (int k = 0; k < values.length; k++) {
        array[k] = values[k];
    }
    Median median = new Median();
    return median.evaluate(array);
}
 

public double mad(double [] autoCorrelationValues){
    double [] tempTable = new double[autoCorrelationValues.length];
    Median m = new Median();
    double medianValue = m.evaluate(autoCorrelationValues);
    for(int i=0 ; i<autoCorrelationValues.length ;i++){
            tempTable[i] = Math.abs(autoCorrelationValues[i] - medianValue);
    }
    return m.evaluate(tempTable); //return the median of tempTable, the equation (13) in the paper
}
 

/**
 * Updates the likelihood of the NonRef allele (if present) based on the likelihoods of a set of non-symbolic
 * <p>
 *     This method does
 * </p>
 *
 * @param allelesToConsider
 */
@SuppressWarnings("unchecked")  // for the cast (A) Allele.NON_REF_ALLELE below
public void updateNonRefAlleleLikelihoods(final AlleleList<A> allelesToConsider) {
    final int nonRefAlleleIndex = indexOfAllele((A) Allele.NON_REF_ALLELE);
    if ( nonRefAlleleIndex < 0) {
        return;
    }
    final int alleleCount = alleles.numberOfAlleles();
    final int nonSymbolicAlleleCount = alleleCount - 1;
    // likelihood buffer reused across evidence:
    final double[] qualifiedAlleleLikelihoods = new double[nonSymbolicAlleleCount];
    final Median medianCalculator = new Median();
    for (int s = 0; s < samples.numberOfSamples(); s++) {
        final double[][] sampleValues = valuesBySampleIndex[s];
        final int evidenceCount = evidenceBySampleIndex.get(s).size();
        for (int r = 0; r < evidenceCount; r++) {
            final BestAllele bestAllele = searchBestAllele(s, r, true);
            int numberOfQualifiedAlleleLikelihoods = 0;
            for (int i = 0; i < alleleCount; i++) {
                final double alleleLikelihood = sampleValues[i][r];
                if (i != nonRefAlleleIndex && alleleLikelihood < bestAllele.likelihood
                        && !Double.isNaN(alleleLikelihood) && allelesToConsider.indexOfAllele(alleles.getAllele(i)) != MISSING_INDEX) {
                    qualifiedAlleleLikelihoods[numberOfQualifiedAlleleLikelihoods++] = alleleLikelihood;
                }
            }
            final double nonRefLikelihood = medianCalculator.evaluate(qualifiedAlleleLikelihoods, 0, numberOfQualifiedAlleleLikelihoods);
            // when the median is NaN that means that all applicable likekihoods are the same as the best
            // so the evidence is not informative at all given the existing alleles. Unless there is only one (or zero) concrete
            // alleles with give the same (the best) likelihood to the NON-REF. When there is only one (or zero) concrete
            // alleles we set the NON-REF likelihood to NaN.
            sampleValues[nonRefAlleleIndex][r] = !Double.isNaN(nonRefLikelihood) ? nonRefLikelihood
                    : nonSymbolicAlleleCount <= 1 ? Double.NaN : bestAllele.likelihood;
        }
    }
}
 

public double getMedian() {
    Median median = new Median();
    return median.evaluate(doubleArray.getElements());
}
 

public IntendedSequence build (final String intendedSequence, final BarcodeNeighborGroup neighbors) {
	IntendedSequence result = null;
	String root = neighbors.getRootSequence();

	Integer deletedBasePos=null;
	Character deletedBase=null;
	Integer umiCountsIntended = null;
	Double umiBiasIntended = null;


	// need a non-null intendedSequence to check for deletions.
	if (intendedSequence!=null) {
		umiCountsIntended=umiCounts.getCountForKey(intendedSequence);
		umiBiasIntended = umiBias.get(intendedSequence);

		LevenshteinDistanceResult r= LevenshteinDistance.computeLevenshteinDistanceResult(intendedSequence, root, 1, 1, 2);
		String [] ops  = r.getOperations();
		// any position before the last is D, and last is I.

		// gather up the deleted base and position.
		for (int i=0; i<ops.length-2; i++)
			// if a deletion at some base, or a substitution at the last base with an intended sequence, then you can get deletion base/position/rate.
			if (ops[i].equals("D") && ops[ops.length-1].equals("I")) {
				deletedBasePos=i+1; // the position is one based.
				deletedBase = intendedSequence.charAt(deletedBasePos-1); // accessing the array 0 based.
				break;
			}

		// Special case: substitution at the last base only.
		if (ops[ops.length-1].equals("S")) {
			deletedBasePos=ops.length; // the position is one based.
			deletedBase = intendedSequence.charAt(deletedBasePos-1); // accessing the array 0 based.
		}

	}
	int totalRelatedUMIs =neighbors.getNeighborCellBarcodes().stream().mapToInt(x -> this.umiCounts.getCountForKey(x)).sum();
	double [] neighborUMIs = neighbors.getNeighborCellBarcodes().stream().mapToDouble(x -> this.umiCounts.getCountForKey(x)).toArray();
	double [] neighborUMIBias = neighbors.getNeighborCellBarcodes().stream().mapToDouble(x -> this.umiBias.get(x)).toArray();

	Median m = new Median();
	double medianNeighborUMIBias = m.evaluate(neighborUMIBias);
	double medianRelatedSequenceUMIs = m.evaluate(neighborUMIs);

	List<String> neighborBC = neighbors.getNeighborCellBarcodes();

	result = new IntendedSequence(intendedSequence, neighborBC, deletedBase, deletedBasePos,
			umiCountsIntended, medianRelatedSequenceUMIs, totalRelatedUMIs, umiBiasIntended, medianNeighborUMIBias);

	return result;

}
 

public double theta1(double [] autoCorrelationValues, double gama){
    Median m = new Median();
    return  (m.evaluate(autoCorrelationValues) + (gama * mad(autoCorrelationValues)));				
}
 

public double theta2(double [][] crossCorrelationValues, double gama){
    double[] vecCrossCorrelation = transformMatrix(crossCorrelationValues);
    Median m = new Median();		
    return (m.evaluate(vecCrossCorrelation) + (gama * mad(vecCrossCorrelation)));		
}
 

@Test(dataProvider = "dataSets")
public void testAddNonRefAllele(final String[] samples, final Allele[] alleles, final Map<String,List<GATKRead>> reads) {
    final AlleleLikelihoods<GATKRead, Allele> original = new AlleleLikelihoods<>(new IndexedSampleList(samples), new IndexedAlleleList<>(alleles), reads);
    final AlleleLikelihoods<GATKRead, Allele> result = new AlleleLikelihoods<>(new IndexedSampleList(samples), new IndexedAlleleList<>(alleles), reads);
    final double[][][] originalLikelihoods = fillWithRandomLikelihoods(samples,alleles,original, result);

    result.addNonReferenceAllele(Allele.NON_REF_ALLELE);
    Assert.assertEquals(result.numberOfAlleles(), original.numberOfAlleles() + 1);
    Assert.assertEquals(result.indexOfAllele(Allele.NON_REF_ALLELE), result.numberOfAlleles() - 1);
    final double[][][] newLikelihoods = new double[originalLikelihoods.length][][];
    for (int s = 0; s < samples.length; s++) {
        newLikelihoods[s] = Arrays.copyOf(originalLikelihoods[s], originalLikelihoods[s].length + 1);
        final int sampleReadCount = original.sampleEvidenceCount(s);
        final int ordinarynumberOfAlleles = originalLikelihoods[s].length;
        newLikelihoods[s][ordinarynumberOfAlleles] = new double[sampleReadCount];
        for (int r = 0; r < sampleReadCount; r++) {

            //TODO secondBestLk is totaly irrelevant, and this code is really just a MathUtils.max to get bestLk
            double bestLk = newLikelihoods[s][0][r];
            double secondBestLk = Double.NEGATIVE_INFINITY;
            for (int a = 1; a < ordinarynumberOfAlleles; a++) {
                final double lk = originalLikelihoods[s][a][r];
                if (lk > bestLk) {
                    secondBestLk = bestLk;
                    bestLk = lk;
                } else if (lk > secondBestLk) {
                    secondBestLk = lk;
                }
            }
            final Median median = new Median();
            final List<Double> qualifylingLikelihoods = new ArrayList<>();
            for (int a = 0; a < ordinarynumberOfAlleles; a++) {
                if (originalLikelihoods[s][a][r] >= bestLk) continue;
                qualifylingLikelihoods.add(originalLikelihoods[s][a][r]);
            }
            final double medianLikelihood = median.evaluate(qualifylingLikelihoods.stream().mapToDouble(d -> d).toArray());
            // NaN is returned in cases whether there is no elements in qualifyingLikelihoods.
            // In such case we set the NON-REF likelihood to -Inf.
            final double expectedNonRefLk = !Double.isNaN(medianLikelihood) ? medianLikelihood
                    : ordinarynumberOfAlleles <= 1 ? Double.NaN : bestLk;
            newLikelihoods[s][ordinarynumberOfAlleles][r] = expectedNonRefLk;
        }
    }
    testLikelihoodMatrixQueries(samples,result,newLikelihoods);
}