Java Code Examples for org.apache.commons.math3.linear.RealMatrix#getRow()
The following examples show how to use
org.apache.commons.math3.linear.RealMatrix#getRow() .
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
| Source File: SpectralMethods.java From egads with GNU General Public License v3.0 | 6 votes |
|
public static RealMatrix createHankelMatrix(RealMatrix data, int windowSize) {
int n = data.getRowDimension();
int m = data.getColumnDimension();
int k = n - windowSize + 1;
RealMatrix res = MatrixUtils.createRealMatrix(k, m * windowSize);
double[] buffer = {};
for (int i = 0; i < n; ++i) {
double[] row = data.getRow(i);
buffer = ArrayUtils.addAll(buffer, row);
if (i >= windowSize - 1) {
RealMatrix mat = MatrixUtils.createRowRealMatrix(buffer);
res.setRowMatrix(i - windowSize + 1, mat);
buffer = ArrayUtils.subarray(buffer, m, buffer.length);
}
}
return res;
}
Example 2
| Source File: NormalizeSomaticReadCountsIntegrationTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 6 votes |
|
private void assertFactorNormalizedValues(final ReadCountCollection input, final ReadCountCollection factorNormalized) {
try (final HDF5File ponReader = new HDF5File(TEST_PON)) {
final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
final double[] targetFactors = pon.getTargetFactors();
final List<String> ponTargets = pon.getTargetNames();
final Map<String,Integer> ponTargetIndexes = new HashMap<>(ponTargets.size());
for (int i = 0; i < ponTargets.size(); i++) {
ponTargetIndexes.put(ponTargets.get(i),i);
}
final RealMatrix inputCounts = input.counts();
final RealMatrix factorNormalizedCounts = factorNormalized.counts();
for (int i = 0; i < factorNormalizedCounts.getRowDimension(); i++) {
final double factor = targetFactors[ponTargetIndexes.get(factorNormalized.targets().get(i).getName())];
final double[] inputValues = inputCounts.getRow(i);
final double[] outputValues = factorNormalizedCounts.getRow(i);
for (int j = 0; j < inputValues.length; j++) {
final double expected = inputValues[j] / factor;
Assert.assertEquals(outputValues[j],expected,0.0000001,"" + i + " , " + j);
}
}
}
}
Example 3
| Source File: NormalizeSomaticReadCountsIntegrationTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 6 votes |
|
private RealMatrix reorderTargetsToPoNOrder(final ReadCountCollection preTangentNormalized, final List<String> ponTargets) {
final RealMatrix preTangentNormalizedCounts = preTangentNormalized.counts();
final Map<String,Integer> ponTargetIndex = IntStream.range(0, ponTargets.size())
.boxed().collect(Collectors.toMap(ponTargets::get, Function.identity()));
// first we need to sort the input counts so that they match the
// target order in the PoN.
final double[][] ponPreparedInput = new double[ponTargets.size()][];
for (int i = 0; i < preTangentNormalizedCounts.getRowDimension(); i++) {
final Target target = preTangentNormalized.targets().get(i);
if (!ponTargetIndex.containsKey(target.getName()))
continue;
final int idx = ponTargetIndex.get(target.getName());
ponPreparedInput[idx] = preTangentNormalizedCounts.getRow(i);
}
// The actual input to create the beta-hats, sorted by the PoN targets:
return new Array2DRowRealMatrix(ponPreparedInput,false);
}
Example 4
| Source File: PoNTestUtils.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 6 votes |
|
/**
* Test whether two matrices are equal (within 1e-4)
* @param left never {@code null}
* @param right never {@code null}
* @param isAllowNegatedValues whether values that are just negated are still considered equal. True is useful for
* the outputs of some matrix operations, such as SVD.
*/
public static void assertEqualsMatrix(final RealMatrix left, final RealMatrix right, final boolean isAllowNegatedValues) {
Assert.assertEquals(left.getRowDimension(), right.getRowDimension());
Assert.assertEquals(left.getColumnDimension(), right.getColumnDimension());
for (int i = 0; i < left.getRowDimension(); i++) {
final double[] leftRow = left.getRow(i);
final double[] rightRow = right.getRow(i);
for (int j = 0; j < leftRow.length; j++) {
if (isAllowNegatedValues) {
Assert.assertEquals(Math.abs(leftRow[j]), Math.abs(rightRow[j]), DOUBLE_MATRIX_TOLERANCE);
} else {
Assert.assertEquals(leftRow[j], rightRow[j], DOUBLE_MATRIX_TOLERANCE);
}
}
}
}
Example 5
| Source File: ReadCountCollectionUnitTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
@Test(dataProvider="targetArrangeData")
public void testArrangeTargets(final ReadCountCollectionInfo info, final List<String> newOrder) {
final List<Target> targetsNewOrder = newOrder.stream()
.map(Target::new).collect(Collectors.toList());
final ReadCountCollection subject = info.newInstance();
final ReadCountCollection result = subject.arrangeTargets(targetsNewOrder);
final List<Target> afterTargets = result.targets();
Assert.assertEquals(afterTargets.size(), targetsNewOrder.size());
Assert.assertFalse(afterTargets.stream().anyMatch(t -> !targetsNewOrder.contains(t)));
final RealMatrix beforeCounts = subject.counts();
final RealMatrix afterCounts = result.counts();
Assert.assertEquals(beforeCounts.getColumnDimension(), afterCounts.getColumnDimension());
Assert.assertEquals(afterCounts.getRowDimension(), targetsNewOrder.size());
final int[] beforeIndexes = new int[targetsNewOrder.size()];
final int[] afterIndexes = new int[targetsNewOrder.size()];
int nextIdx = 0;
for (final Target target : targetsNewOrder) {
final int beforeIndex = subject.targets().indexOf(target);
final int afterIndex = result.targets().indexOf(target);
beforeIndexes[nextIdx] = beforeIndex;
afterIndexes[nextIdx++] = afterIndex;
}
// check that the counts are exactly the same.
for (int i = 0; i < beforeIndexes.length; i++) {
final double[] before = beforeCounts.getRow(beforeIndexes[i]);
final double[] after = afterCounts.getRow(afterIndexes[i]);
Assert.assertEquals(before, after);
}
}
Example 6
| Source File: ReadCountCollectionUnitTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
@Test(dataProvider="targetSubsetData")
public void testSubsetTargets(final ReadCountCollectionInfo info, final Set<String> targetNamesToKeep) {
final Set<Target> targetsToKeep = targetNamesToKeep.stream()
.map(Target::new).collect(Collectors.toSet());
final ReadCountCollection subject = info.newInstance();
final ReadCountCollection result = subject.subsetTargets(targetsToKeep);
final List<Target> afterTargets = result.targets();
Assert.assertEquals(afterTargets.size(), targetsToKeep.size());
Assert.assertFalse(afterTargets.stream().anyMatch(t -> !targetsToKeep.contains(t)));
final RealMatrix beforeCounts = subject.counts();
final RealMatrix afterCounts = result.counts();
Assert.assertEquals(beforeCounts.getColumnDimension(), afterCounts.getColumnDimension());
Assert.assertEquals(afterCounts.getRowDimension(), targetNamesToKeep.size());
final int[] beforeIndexes = new int[targetsToKeep.size()];
final int[] afterIndexes = new int[targetsToKeep.size()];
int nextIdx = 0;
for (final Target target : targetsToKeep) {
final int beforeIndex = subject.targets().indexOf(target);
final int afterIndex = result.targets().indexOf(target);
beforeIndexes[nextIdx] = beforeIndex;
afterIndexes[nextIdx++] = afterIndex;
}
// check that the order of targets in the output is kept to the original order.
for (int i = 0; i < beforeIndexes.length; i++) {
for (int j = 0; j < beforeIndexes.length; j++) {
Assert.assertEquals(Integer.compare(beforeIndexes[i], beforeIndexes[j]),
Integer.compare(afterIndexes[i], afterIndexes[j]));
}
}
// check that the counts are exactly the same.
for (int i = 0; i < beforeIndexes.length; i++) {
final double[] before = beforeCounts.getRow(beforeIndexes[i]);
final double[] after = afterCounts.getRow(afterIndexes[i]);
Assert.assertEquals(before, after);
}
}
Example 7
| Source File: MultivariateNormal.java From macrobase with Apache License 2.0 | 5 votes |
|
public MultivariateNormal(RealVector mean, RealMatrix sigma) {
double[][] arrayOfMatrix = new double[sigma.getColumnDimension()][sigma.getRowDimension()];
for (int i = 0; i < sigma.getColumnDimension(); i++) {
arrayOfMatrix[i] = sigma.getRow(i);
}
distribution = new MultivariateNormalDistribution(mean.toArray(), arrayOfMatrix);
}
Example 8
| Source File: HDF5LibraryUnitTest.java From gatk with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
/**
* Test whether two matrices are equal (within 1e-4)
* @param left never {@code null}
* @param right never {@code null}
*/
private static void assertEqualsMatrix(final RealMatrix left, final RealMatrix right) {
Assert.assertEquals(left.getRowDimension(), right.getRowDimension());
Assert.assertEquals(left.getColumnDimension(), right.getColumnDimension());
for (int i = 0; i < left.getRowDimension(); i++) {
final double[] leftRow = left.getRow(i);
final double[] rightRow = right.getRow(i);
for (int j = 0; j < leftRow.length; j++) {
Assert.assertEquals(leftRow[j], rightRow[j], DOUBLE_MATRIX_TOLERANCE);
}
}
}
Example 9
| Source File: HDF5UtilsUnitTest.java From gatk with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
/**
* Test whether two matrices are equal
* @param left never {@code null}
* @param right never {@code null}
*/
private static void assertEqualsMatrix(final RealMatrix left, final RealMatrix right, final double tolerance) {
Assert.assertEquals(left.getRowDimension(), right.getRowDimension());
Assert.assertEquals(left.getColumnDimension(), right.getColumnDimension());
for (int i = 0; i < left.getRowDimension(); i++) {
final double[] leftRow = left.getRow(i);
final double[] rightRow = right.getRow(i);
for (int j = 0; j < leftRow.length; j++) {
Assert.assertEquals(leftRow[j], rightRow[j], tolerance);
}
}
}
Example 10
| Source File: ReadCountCollectionUtils.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 4 votes |
|
/**
* Truncates the extreme count values in the input read-count collection.
* Values are forced to be bound by the percentile indicated with the input {@code percentile} which must be
* in the range [0 .. 50.0]. Values under that percentile and the complementary (1 - percentile) are set to the
* corresponding threshold value.
*
* <p>The imputation is done in-place, thus the input matrix is modified as a result of this call.</p>
*
* @param readCounts the input and output read-count matrix.
*/
public static void truncateExtremeCounts(final ReadCountCollection readCounts, final double percentile, final Logger logger) {
final RealMatrix counts = readCounts.counts();
final int targetCount = counts.getRowDimension();
final int columnCount = counts.getColumnDimension();
// Create a row major array of the counts.
final double[] values = Doubles.concat(counts.getData());
final Percentile bottomPercentileEvaluator = new Percentile(percentile);
final Percentile topPercentileEvaluator = new Percentile(100.0 - percentile);
final double bottomPercentileThreshold = bottomPercentileEvaluator.evaluate(values);
final double topPercentileThreshold = topPercentileEvaluator.evaluate(values);
long totalCounts = 0;
long bottomTruncatedCounts = 0;
long topTruncatedCounts = 0;
for (int i = 0; i < targetCount; i++) {
final double[] rowCounts = counts.getRow(i);
for (int j = 0; j < columnCount; j++) {
final double count = rowCounts[j];
totalCounts++;
if (count < bottomPercentileThreshold) {
counts.setEntry(i, j, bottomPercentileThreshold);
bottomTruncatedCounts++;
} else if (count > topPercentileThreshold) {
counts.setEntry(i, j, topPercentileThreshold);
topTruncatedCounts++;
}
}
}
if (topTruncatedCounts == 0 && bottomTruncatedCounts == 0) {
logger.info(String.format("None of the %d counts were truncated as they all fall in the non-extreme range " +
"[%.2f, %.2f]", totalCounts, bottomPercentileThreshold, topPercentileThreshold));
} else {
final double truncatedPercentage = ((double)(topTruncatedCounts + bottomTruncatedCounts) / totalCounts) * 100;
logger.info(String.format("Some counts (%d out of %d, %.2f%%) were truncated as they fall out of the " +
"non-extreme range [%.2f, %.2f]", topTruncatedCounts + bottomTruncatedCounts, totalCounts,
truncatedPercentage, bottomPercentileThreshold, topPercentileThreshold));
}
}
Example 11
| Source File: GaussNewtonOptimizer.java From astor with GNU General Public License v2.0 | 4 votes |
|
/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
final ConvergenceChecker<PointVectorValuePair> checker
= getConvergenceChecker();
// Computation will be useless without a checker (see "for-loop").
if (checker == null) {
throw new NullArgumentException();
}
final double[] targetValues = getTarget();
final int nR = targetValues.length; // Number of observed data.
final RealMatrix weightMatrix = getWeight();
// Diagonal of the weight matrix.
final double[] residualsWeights = new double[nR];
for (int i = 0; i < nR; i++) {
residualsWeights[i] = weightMatrix.getEntry(i, i);
}
final double[] currentPoint = getStartPoint();
final int nC = currentPoint.length;
// iterate until convergence is reached
PointVectorValuePair current = null;
int iter = 0;
for (boolean converged = false; !converged;) {
++iter;
// evaluate the objective function and its jacobian
PointVectorValuePair previous = current;
// Value of the objective function at "currentPoint".
final double[] currentObjective = computeObjectiveValue(currentPoint);
final double[] currentResiduals = computeResiduals(currentObjective);
final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
current = new PointVectorValuePair(currentPoint, currentObjective);
// build the linear problem
final double[] b = new double[nC];
final double[][] a = new double[nC][nC];
for (int i = 0; i < nR; ++i) {
final double[] grad = weightedJacobian.getRow(i);
final double weight = residualsWeights[i];
final double residual = currentResiduals[i];
// compute the normal equation
final double wr = weight * residual;
for (int j = 0; j < nC; ++j) {
b[j] += wr * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < nC; ++k) {
double[] ak = a[k];
double wgk = weight * grad[k];
for (int l = 0; l < nC; ++l) {
ak[l] += wgk * grad[l];
}
}
}
try {
// solve the linearized least squares problem
RealMatrix mA = new BlockRealMatrix(a);
DecompositionSolver solver = useLU ?
new LUDecomposition(mA).getSolver() :
new QRDecomposition(mA).getSolver();
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
// update the estimated parameters
for (int i = 0; i < nC; ++i) {
currentPoint[i] += dX[i];
}
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
// Check convergence.
if (previous != null) {
converged = checker.converged(iter, previous, current);
if (converged) {
cost = computeCost(currentResiduals);
// Update (deprecated) "point" field.
point = current.getPoint();
return current;
}
}
}
// Must never happen.
throw new MathInternalError();
}
Example 12
| Source File: GaussNewtonOptimizer.java From astor with GNU General Public License v2.0 | 4 votes |
|
/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
final ConvergenceChecker<PointVectorValuePair> checker
= getConvergenceChecker();
// Computation will be useless without a checker (see "for-loop").
if (checker == null) {
throw new NullArgumentException();
}
final double[] targetValues = getTarget();
final int nR = targetValues.length; // Number of observed data.
final RealMatrix weightMatrix = getWeight();
// Diagonal of the weight matrix.
final double[] residualsWeights = new double[nR];
for (int i = 0; i < nR; i++) {
residualsWeights[i] = weightMatrix.getEntry(i, i);
}
final double[] currentPoint = getStartPoint();
final int nC = currentPoint.length;
// iterate until convergence is reached
PointVectorValuePair current = null;
int iter = 0;
for (boolean converged = false; !converged;) {
++iter;
// evaluate the objective function and its jacobian
PointVectorValuePair previous = current;
// Value of the objective function at "currentPoint".
final double[] currentObjective = computeObjectiveValue(currentPoint);
final double[] currentResiduals = computeResiduals(currentObjective);
final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
current = new PointVectorValuePair(currentPoint, currentObjective);
// build the linear problem
final double[] b = new double[nC];
final double[][] a = new double[nC][nC];
for (int i = 0; i < nR; ++i) {
final double[] grad = weightedJacobian.getRow(i);
final double weight = residualsWeights[i];
final double residual = currentResiduals[i];
// compute the normal equation
final double wr = weight * residual;
for (int j = 0; j < nC; ++j) {
b[j] += wr * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < nC; ++k) {
double[] ak = a[k];
double wgk = weight * grad[k];
for (int l = 0; l < nC; ++l) {
ak[l] += wgk * grad[l];
}
}
}
try {
// solve the linearized least squares problem
RealMatrix mA = new BlockRealMatrix(a);
DecompositionSolver solver = useLU ?
new LUDecomposition(mA).getSolver() :
new QRDecomposition(mA).getSolver();
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
// update the estimated parameters
for (int i = 0; i < nC; ++i) {
currentPoint[i] += dX[i];
}
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
// Check convergence.
if (previous != null) {
converged = checker.converged(iter, previous, current);
if (converged) {
cost = computeCost(currentResiduals);
// Update (deprecated) "point" field.
point = current.getPoint();
return current;
}
}
}
// Must never happen.
throw new MathInternalError();
}
Example 13
| Source File: GaussNewtonOptimizer.java From astor with GNU General Public License v2.0 | 4 votes |
|
/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
final ConvergenceChecker<PointVectorValuePair> checker
= getConvergenceChecker();
// Computation will be useless without a checker (see "for-loop").
if (checker == null) {
throw new NullArgumentException();
}
final double[] targetValues = getTarget();
final int nR = targetValues.length; // Number of observed data.
final RealMatrix weightMatrix = getWeight();
// Diagonal of the weight matrix.
final double[] residualsWeights = new double[nR];
for (int i = 0; i < nR; i++) {
residualsWeights[i] = weightMatrix.getEntry(i, i);
}
final double[] currentPoint = getStartPoint();
final int nC = currentPoint.length;
// iterate until convergence is reached
PointVectorValuePair current = null;
int iter = 0;
for (boolean converged = false; !converged;) {
++iter;
// evaluate the objective function and its jacobian
PointVectorValuePair previous = current;
// Value of the objective function at "currentPoint".
final double[] currentObjective = computeObjectiveValue(currentPoint);
final double[] currentResiduals = computeResiduals(currentObjective);
final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
current = new PointVectorValuePair(currentPoint, currentObjective);
// build the linear problem
final double[] b = new double[nC];
final double[][] a = new double[nC][nC];
for (int i = 0; i < nR; ++i) {
final double[] grad = weightedJacobian.getRow(i);
final double weight = residualsWeights[i];
final double residual = currentResiduals[i];
// compute the normal equation
final double wr = weight * residual;
for (int j = 0; j < nC; ++j) {
b[j] += wr * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < nC; ++k) {
double[] ak = a[k];
double wgk = weight * grad[k];
for (int l = 0; l < nC; ++l) {
ak[l] += wgk * grad[l];
}
}
}
try {
// solve the linearized least squares problem
RealMatrix mA = new BlockRealMatrix(a);
DecompositionSolver solver = useLU ?
new LUDecomposition(mA).getSolver() :
new QRDecomposition(mA).getSolver();
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
// update the estimated parameters
for (int i = 0; i < nC; ++i) {
currentPoint[i] += dX[i];
}
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
// Check convergence.
if (previous != null) {
converged = checker.converged(iter, previous, current);
if (converged) {
cost = computeCost(currentResiduals);
// Update (deprecated) "point" field.
point = current.getPoint();
return current;
}
}
}
// Must never happen.
throw new MathInternalError();
}
Example 14
| Source File: GaussNewtonOptimizer.java From astor with GNU General Public License v2.0 | 4 votes |
|
/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
final ConvergenceChecker<PointVectorValuePair> checker
= getConvergenceChecker();
// Computation will be useless without a checker (see "for-loop").
if (checker == null) {
throw new NullArgumentException();
}
final double[] targetValues = getTarget();
final int nR = targetValues.length; // Number of observed data.
final RealMatrix weightMatrix = getWeight();
// Diagonal of the weight matrix.
final double[] residualsWeights = new double[nR];
for (int i = 0; i < nR; i++) {
residualsWeights[i] = weightMatrix.getEntry(i, i);
}
final double[] currentPoint = getStartPoint();
final int nC = currentPoint.length;
// iterate until convergence is reached
PointVectorValuePair current = null;
int iter = 0;
for (boolean converged = false; !converged;) {
++iter;
// evaluate the objective function and its jacobian
PointVectorValuePair previous = current;
// Value of the objective function at "currentPoint".
final double[] currentObjective = computeObjectiveValue(currentPoint);
final double[] currentResiduals = computeResiduals(currentObjective);
final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
current = new PointVectorValuePair(currentPoint, currentObjective);
// build the linear problem
final double[] b = new double[nC];
final double[][] a = new double[nC][nC];
for (int i = 0; i < nR; ++i) {
final double[] grad = weightedJacobian.getRow(i);
final double weight = residualsWeights[i];
final double residual = currentResiduals[i];
// compute the normal equation
final double wr = weight * residual;
for (int j = 0; j < nC; ++j) {
b[j] += wr * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < nC; ++k) {
double[] ak = a[k];
double wgk = weight * grad[k];
for (int l = 0; l < nC; ++l) {
ak[l] += wgk * grad[l];
}
}
}
try {
// solve the linearized least squares problem
RealMatrix mA = new BlockRealMatrix(a);
DecompositionSolver solver = useLU ?
new LUDecomposition(mA).getSolver() :
new QRDecomposition(mA).getSolver();
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
// update the estimated parameters
for (int i = 0; i < nC; ++i) {
currentPoint[i] += dX[i];
}
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
// Check convergence.
if (previous != null) {
converged = checker.converged(iter, previous, current);
if (converged) {
cost = computeCost(currentResiduals);
// Update (deprecated) "point" field.
point = current.getPoint();
return current;
}
}
}
// Must never happen.
throw new MathInternalError();
}
Example 15
| Source File: GaussNewtonOptimizer.java From astor with GNU General Public License v2.0 | 4 votes |
|
/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
final ConvergenceChecker<PointVectorValuePair> checker
= getConvergenceChecker();
// Computation will be useless without a checker (see "for-loop").
if (checker == null) {
throw new NullArgumentException();
}
final double[] targetValues = getTarget();
final int nR = targetValues.length; // Number of observed data.
final RealMatrix weightMatrix = getWeight();
// Diagonal of the weight matrix.
final double[] residualsWeights = new double[nR];
for (int i = 0; i < nR; i++) {
residualsWeights[i] = weightMatrix.getEntry(i, i);
}
final double[] currentPoint = getStartPoint();
final int nC = currentPoint.length;
// iterate until convergence is reached
PointVectorValuePair current = null;
int iter = 0;
for (boolean converged = false; !converged;) {
++iter;
// evaluate the objective function and its jacobian
PointVectorValuePair previous = current;
// Value of the objective function at "currentPoint".
final double[] currentObjective = computeObjectiveValue(currentPoint);
final double[] currentResiduals = computeResiduals(currentObjective);
final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
current = new PointVectorValuePair(currentPoint, currentObjective);
// build the linear problem
final double[] b = new double[nC];
final double[][] a = new double[nC][nC];
for (int i = 0; i < nR; ++i) {
final double[] grad = weightedJacobian.getRow(i);
final double weight = residualsWeights[i];
final double residual = currentResiduals[i];
// compute the normal equation
final double wr = weight * residual;
for (int j = 0; j < nC; ++j) {
b[j] += wr * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < nC; ++k) {
double[] ak = a[k];
double wgk = weight * grad[k];
for (int l = 0; l < nC; ++l) {
ak[l] += wgk * grad[l];
}
}
}
try {
// solve the linearized least squares problem
RealMatrix mA = new BlockRealMatrix(a);
DecompositionSolver solver = useLU ?
new LUDecomposition(mA).getSolver() :
new QRDecomposition(mA).getSolver();
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
// update the estimated parameters
for (int i = 0; i < nC; ++i) {
currentPoint[i] += dX[i];
}
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
// Check convergence.
if (previous != null) {
converged = checker.converged(iter, previous, current);
if (converged) {
cost = computeCost(currentResiduals);
// Update (deprecated) "point" field.
point = current.getPoint();
return current;
}
}
}
// Must never happen.
throw new MathInternalError();
}
