use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project GDSC-SMLM by aherbert.
the class ApacheLVMFitter method computeFit.
public FitStatus computeFit(double[] y, final double[] y_fit, double[] a, double[] a_dev) {
int n = y.length;
try {
// Different convergence thresholds seem to have no effect on the resulting fit, only the number of
// iterations for convergence
final double initialStepBoundFactor = 100;
final double costRelativeTolerance = 1e-10;
final double parRelativeTolerance = 1e-10;
final double orthoTolerance = 1e-10;
final double threshold = Precision.SAFE_MIN;
// Extract the parameters to be fitted
final double[] initialSolution = getInitialSolution(a);
// TODO - Pass in more advanced stopping criteria.
// Create the target and weight arrays
final double[] yd = new double[n];
final double[] w = new double[n];
for (int i = 0; i < n; i++) {
yd[i] = y[i];
w[i] = 1;
}
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
//@formatter:off
LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd).weight(new DiagonalMatrix(w));
if (f instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) f).canComputeValuesAndJacobian()) {
// Compute together, or each individually
builder.model(new ValueAndJacobianFunction() {
final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) f;
public Pair<RealVector, RealMatrix> value(RealVector point) {
final double[] p = point.toArray();
final Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
return new Pair<RealVector, RealMatrix>(new ArrayRealVector(result.getFirst(), false), new Array2DRowRealMatrix(result.getSecond(), false));
}
public RealVector computeValue(double[] params) {
return new ArrayRealVector(fun.computeValues(params), false);
}
public RealMatrix computeJacobian(double[] params) {
return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
}
});
} else {
// Compute separately
builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) f, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) f, a, n));
}
LeastSquaresProblem problem = builder.build();
Optimum optimum = optimizer.optimize(problem);
final double[] parameters = optimum.getPoint().toArray();
setSolution(a, parameters);
iterations = optimum.getIterations();
evaluations = optimum.getEvaluations();
if (a_dev != null) {
try {
double[][] covar = optimum.getCovariances(threshold).getData();
setDeviationsFromMatrix(a_dev, covar);
} catch (SingularMatrixException e) {
// Matrix inversion failed. In order to return a solution
// return the reciprocal of the diagonal of the Fisher information
// for a loose bound on the limit
final int[] gradientIndices = f.gradientIndices();
final int nparams = gradientIndices.length;
GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams);
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
calculator.findLinearised(nparams, y, a, alpha, beta, (NonLinearFunction) f);
FisherInformationMatrix m = new FisherInformationMatrix(alpha);
setDeviations(a_dev, m.crlb(true));
}
}
// Compute function value
if (y_fit != null) {
Gaussian2DFunction f = (Gaussian2DFunction) this.f;
f.initialise0(a);
f.forEach(new ValueProcedure() {
int i = 0;
public void execute(double value) {
y_fit[i] = value;
}
});
}
// As this is unweighted then we can do this to get the sum of squared residuals
// This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
// just computes the dot product anyway.
value = optimum.getResiduals().dotProduct(optimum.getResiduals());
} catch (TooManyEvaluationsException e) {
return FitStatus.TOO_MANY_EVALUATIONS;
} catch (TooManyIterationsException e) {
return FitStatus.TOO_MANY_ITERATIONS;
} catch (ConvergenceException e) {
// Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
return FitStatus.SINGULAR_NON_LINEAR_MODEL;
} catch (Exception e) {
// TODO - Find out the other exceptions from the fitter and add return values to match.
return FitStatus.UNKNOWN;
}
return FitStatus.OK;
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class ReCapSegCallerUnitTest method testMakeCalls.
@Test
public void testMakeCalls() {
final List<Target> targets = new ArrayList<>();
final List<String> columnNames = Arrays.asList("Sample");
final List<Double> coverage = new ArrayList<>();
//add amplification targets
for (int i = 0; i < 10; i++) {
final SimpleInterval interval = new SimpleInterval("chr", 100 + 2 * i, 101 + 2 * i);
targets.add(new Target(interval));
coverage.add(ParamUtils.log2(2.0));
}
//add deletion targets
for (int i = 0; i < 10; i++) {
final SimpleInterval interval = new SimpleInterval("chr", 300 + 2 * i, 301 + 2 * i);
targets.add(new Target(interval));
coverage.add(ParamUtils.log2(0.5));
}
//add targets that don't belong to a segment
for (int i = 1; i < 10; i++) {
final SimpleInterval interval = new SimpleInterval("chr", 400 + 2 * i, 401 + 2 * i);
targets.add(new Target(interval));
coverage.add(ParamUtils.log2(1.0));
}
//add obviously neutral targets with some small spread
for (int i = -5; i < 6; i++) {
final SimpleInterval interval = new SimpleInterval("chr", 500 + 2 * i, 501 + 2 * i);
targets.add(new Target(interval));
coverage.add(ParamUtils.log2(0.01 * i + 1));
}
//add spread-out targets to a neutral segment (mean near zero)
for (int i = -5; i < 6; i++) {
final SimpleInterval interval = new SimpleInterval("chr", 700 + 2 * i, 701 + 2 * i);
targets.add(new Target(interval));
coverage.add(ParamUtils.log2(0.1 * i + 1));
}
final RealMatrix coverageMatrix = new Array2DRowRealMatrix(targets.size(), 1);
coverageMatrix.setColumn(0, coverage.stream().mapToDouble(x -> x).toArray());
final int n = targets.size();
final int m = coverageMatrix.getRowDimension();
final ReadCountCollection counts = new ReadCountCollection(targets, columnNames, coverageMatrix);
List<ModeledSegment> segments = new ArrayList<>();
//amplification
segments.add(new ModeledSegment(new SimpleInterval("chr", 100, 200), 100, ParamUtils.log2(2.0)));
//deletion
segments.add(new ModeledSegment(new SimpleInterval("chr", 300, 400), 100, ParamUtils.log2(0.5)));
//neutral
segments.add(new ModeledSegment(new SimpleInterval("chr", 450, 550), 100, ParamUtils.log2(1)));
//neutral
segments.add(new ModeledSegment(new SimpleInterval("chr", 650, 750), 100, ParamUtils.log2(1)));
List<ModeledSegment> calls = ReCapSegCaller.makeCalls(counts, segments);
Assert.assertEquals(calls.get(0).getCall(), ReCapSegCaller.AMPLIFICATION_CALL);
Assert.assertEquals(calls.get(1).getCall(), ReCapSegCaller.DELETION_CALL);
Assert.assertEquals(calls.get(2).getCall(), ReCapSegCaller.NEUTRAL_CALL);
Assert.assertEquals(calls.get(3).getCall(), ReCapSegCaller.NEUTRAL_CALL);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class ReadCountCollectionUtilsUnitTest method tooManyZerosData.
@DataProvider(name = "tooManyZerosData")
public Object[][] tooManyZerosData() {
final double[] zeroProbabilities = new double[] { .001, .01, .02, 0.1 };
final List<Object[]> result = new ArrayList<>();
final Random rdn = new Random(13);
final int columnCount = 100;
final int targetCount = 100;
final List<String> columnNames = IntStream.range(0, columnCount).mapToObj(i -> "sample_" + (i + 1)).collect(Collectors.toList());
final List<Target> targets = IntStream.range(0, targetCount).mapToObj(i -> new Target("target_" + (i + 1))).collect(Collectors.toList());
for (final double zeroProbability : zeroProbabilities) {
final double[][] counts = new double[columnCount][targetCount];
for (int i = 0; i < counts.length; i++) {
for (int j = 0; j < counts[0].length; j++) {
counts[i][j] = rdn.nextDouble() <= zeroProbability ? 0.0 : rdn.nextDouble();
}
}
final ReadCountCollection readCounts = new ReadCountCollection(targets, columnNames, new Array2DRowRealMatrix(counts, false));
result.add(new Object[] { readCounts });
}
return result.toArray(new Object[result.size()][]);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk 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-protected by broadinstitute.
the class SomaticGenotypingEngine method getAsRealMatrix.
//convert a likelihood matrix of alleles x reads into a RealMatrix
public static RealMatrix getAsRealMatrix(final LikelihoodMatrix<Allele> matrix) {
final RealMatrix result = new Array2DRowRealMatrix(matrix.numberOfAlleles(), matrix.numberOfReads());
result.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
@Override
public double visit(int row, int column, double value) {
return matrix.get(row, column);
}
});
return result;
}
Aggregations