1 |
|
|
2 |
|
|
3 |
|
|
4 |
|
|
5 |
|
|
6 |
|
|
7 |
|
|
8 |
|
|
9 |
|
|
10 |
|
|
11 |
|
|
12 |
|
|
13 |
|
|
14 |
|
|
15 |
|
|
16 |
|
|
17 |
|
|
18 |
|
|
19 |
|
|
20 |
|
|
21 |
|
|
22 |
|
|
23 |
|
|
24 |
|
|
25 |
|
|
26 |
|
|
27 |
|
|
28 |
|
|
29 |
|
|
30 |
|
|
31 |
|
|
32 |
|
|
33 |
|
|
34 |
|
|
35 |
|
|
36 |
|
|
37 |
|
package jalview.analysis; |
38 |
|
|
39 |
|
import jalview.bin.Console; |
40 |
|
import jalview.math.MatrixI; |
41 |
|
import jalview.math.Matrix; |
42 |
|
import jalview.math.MiscMath; |
43 |
|
|
44 |
|
import java.lang.Math; |
45 |
|
import java.lang.System; |
46 |
|
import java.util.Arrays; |
47 |
|
import java.util.ArrayList; |
48 |
|
import java.util.Comparator; |
49 |
|
import java.util.Map.Entry; |
50 |
|
import java.util.TreeMap; |
51 |
|
|
52 |
|
import org.apache.commons.math3.linear.Array2DRowRealMatrix; |
53 |
|
import org.apache.commons.math3.linear.SingularValueDecomposition; |
54 |
|
|
55 |
|
|
56 |
|
|
57 |
|
|
|
|
| 0% |
Uncovered Elements: 433 (433) |
Complexity: 76 |
Complexity Density: 0.25 |
|
58 |
|
public class ccAnalysis |
59 |
|
{ |
60 |
|
private byte dim = 0; |
61 |
|
|
62 |
|
private MatrixI scoresOld; |
63 |
|
|
|
|
| 0% |
Uncovered Elements: 12 (12) |
Complexity: 4 |
Complexity Density: 0.67 |
|
64 |
0 |
public ccAnalysis(MatrixI scores, byte dim)... |
65 |
|
{ |
66 |
|
|
67 |
0 |
for (int i = 0; i < scores.height(); i++) |
68 |
|
{ |
69 |
0 |
for (int j = 0; j < scores.width(); j++) |
70 |
|
{ |
71 |
0 |
if (!Double.isNaN(scores.getValue(i, j))) |
72 |
|
{ |
73 |
0 |
scores.setValue(i, j, |
74 |
|
(double) Math.round(scores.getValue(i, j) * (int) 10000) |
75 |
|
/ 10000); |
76 |
|
} |
77 |
|
} |
78 |
|
} |
79 |
0 |
this.scoresOld = scores; |
80 |
0 |
this.dim = dim; |
81 |
|
} |
82 |
|
|
83 |
|
|
84 |
|
|
85 |
|
|
86 |
|
|
87 |
|
@param |
88 |
|
|
89 |
|
@param |
90 |
|
|
91 |
|
|
92 |
|
@return |
93 |
|
|
|
|
| 0% |
Uncovered Elements: 23 (23) |
Complexity: 5 |
Complexity Density: 0.33 |
|
94 |
0 |
private int[] initialiseDistrusts(byte[] hSigns, MatrixI scores)... |
95 |
|
{ |
96 |
0 |
int[] distrustScores = new int[scores.width()]; |
97 |
|
|
98 |
|
|
99 |
0 |
for (int i = 0; i < scores.width(); i++) |
100 |
|
{ |
101 |
0 |
byte hASign = hSigns[i]; |
102 |
0 |
int conHypNum = 0; |
103 |
0 |
int proHypNum = 0; |
104 |
|
|
105 |
0 |
for (int j = 0; j < scores.width(); j++) |
106 |
|
{ |
107 |
0 |
double cell = scores.getRow(i)[j]; |
108 |
0 |
byte hBSign = hSigns[j]; |
109 |
0 |
if (!Double.isNaN(cell)) |
110 |
|
{ |
111 |
0 |
byte cellSign = (byte) Math.signum(cell); |
112 |
|
|
113 |
0 |
if (cellSign == hASign * hBSign) |
114 |
|
{ |
115 |
0 |
proHypNum++; |
116 |
|
} |
117 |
|
else |
118 |
|
{ |
119 |
0 |
conHypNum++; |
120 |
|
} |
121 |
|
} |
122 |
|
} |
123 |
0 |
distrustScores[i] = conHypNum - proHypNum; |
124 |
|
|
125 |
|
} |
126 |
0 |
return distrustScores; |
127 |
|
} |
128 |
|
|
129 |
|
|
130 |
|
|
131 |
|
|
132 |
|
|
133 |
|
|
134 |
|
@param |
135 |
|
|
136 |
|
@param |
137 |
|
@param |
138 |
|
|
139 |
|
|
140 |
|
@return |
141 |
|
|
|
|
| 0% |
Uncovered Elements: 28 (28) |
Complexity: 6 |
Complexity Density: 0.33 |
|
142 |
0 |
private byte[] optimiseHypothesis(byte[] hSigns, int[] distrustScores,... |
143 |
|
MatrixI scores) |
144 |
|
{ |
145 |
|
|
146 |
0 |
int[] maxes = MiscMath.findMax(distrustScores); |
147 |
0 |
int maxDistrustIndex = maxes[0]; |
148 |
0 |
int maxDistrust = maxes[1]; |
149 |
|
|
150 |
|
|
151 |
0 |
if (maxDistrust > 0) |
152 |
|
{ |
153 |
|
|
154 |
0 |
hSigns[maxDistrustIndex] *= -1; |
155 |
|
|
156 |
0 |
distrustScores[maxDistrustIndex] *= -1; |
157 |
|
|
158 |
|
|
159 |
0 |
byte hASign = hSigns[maxDistrustIndex]; |
160 |
0 |
for (int NOTmaxDistrustIndex = 0; NOTmaxDistrustIndex < distrustScores.length; NOTmaxDistrustIndex++) |
161 |
|
{ |
162 |
0 |
if (NOTmaxDistrustIndex != maxDistrustIndex) |
163 |
|
{ |
164 |
0 |
byte hBSign = hSigns[NOTmaxDistrustIndex]; |
165 |
0 |
double cell = scores.getValue(maxDistrustIndex, |
166 |
|
NOTmaxDistrustIndex); |
167 |
|
|
168 |
|
|
169 |
0 |
if (!Double.isNaN(cell)) |
170 |
|
{ |
171 |
0 |
byte cellSign = (byte) Math.signum(cell); |
172 |
|
|
173 |
|
|
174 |
|
|
175 |
0 |
if (cellSign == hASign * hBSign) |
176 |
|
{ |
177 |
0 |
distrustScores[NOTmaxDistrustIndex] -= 2; |
178 |
|
} |
179 |
|
else |
180 |
|
{ |
181 |
0 |
distrustScores[NOTmaxDistrustIndex] += 2; |
182 |
|
} |
183 |
|
} |
184 |
|
} |
185 |
|
} |
186 |
|
|
187 |
0 |
return optimiseHypothesis(hSigns, distrustScores, scores); |
188 |
|
|
189 |
|
} |
190 |
|
else |
191 |
|
{ |
192 |
0 |
return hSigns; |
193 |
|
} |
194 |
|
} |
195 |
|
|
196 |
|
|
197 |
|
|
198 |
|
|
199 |
|
|
200 |
|
|
201 |
|
|
202 |
|
@param |
203 |
|
|
204 |
|
|
205 |
|
@return |
206 |
|
|
|
|
| 0% |
Uncovered Elements: 152 (152) |
Complexity: 21 |
Complexity Density: 0.18 |
|
207 |
0 |
public MatrixI run() throws Exception... |
208 |
|
{ |
209 |
|
|
210 |
0 |
MatrixI eigenMatrix = scoresOld.copy(); |
211 |
0 |
MatrixI repMatrix = scoresOld.copy(); |
212 |
0 |
try |
213 |
|
{ |
214 |
|
|
215 |
|
|
216 |
|
|
217 |
|
|
218 |
|
|
219 |
|
|
220 |
|
|
221 |
|
|
222 |
0 |
System.out.println("Input correlation matrix:"); |
223 |
0 |
eigenMatrix.print(System.out, "%1.4f "); |
224 |
|
|
225 |
0 |
int matrixWidth = eigenMatrix.width(); |
226 |
|
|
227 |
0 |
int matrixElementsTotal = (int) Math.pow(matrixWidth, 2); |
228 |
|
|
229 |
|
|
230 |
0 |
float correctionFactor = (float) (matrixElementsTotal |
231 |
|
- eigenMatrix.countNaN()) / (float) matrixElementsTotal; |
232 |
|
|
233 |
|
|
234 |
|
|
235 |
|
|
236 |
|
|
237 |
|
|
238 |
|
|
239 |
|
|
240 |
|
|
241 |
|
|
242 |
|
|
243 |
|
|
244 |
|
|
245 |
|
|
246 |
|
|
247 |
0 |
byte[] hSigns = new byte[matrixWidth]; |
248 |
0 |
Arrays.fill(hSigns, (byte) 1); |
249 |
|
|
250 |
|
|
251 |
0 |
hSigns = optimiseHypothesis(hSigns, |
252 |
|
initialiseDistrusts(hSigns, eigenMatrix), eigenMatrix); |
253 |
|
|
254 |
|
|
255 |
|
|
256 |
0 |
double[] hAbs = MiscMath.sqrt(eigenMatrix.absolute().meanRow()); |
257 |
|
|
258 |
|
|
259 |
|
|
260 |
0 |
double[] hValues = MiscMath.elementwiseMultiply(hSigns, hAbs); |
261 |
|
|
262 |
|
|
263 |
|
|
264 |
|
|
265 |
|
|
266 |
|
|
267 |
|
|
268 |
|
|
269 |
0 |
ArrayList<int[]> estimatedPositions = new ArrayList<int[]>(); |
270 |
|
|
271 |
|
|
272 |
0 |
for (int rowIndex = 0; rowIndex < matrixWidth - 1; rowIndex++) |
273 |
|
{ |
274 |
0 |
for (int columnIndex = rowIndex |
275 |
0 |
+ 1; columnIndex < matrixWidth; columnIndex++) |
276 |
|
{ |
277 |
0 |
double cell = eigenMatrix.getValue(rowIndex, columnIndex); |
278 |
0 |
if (Double.isNaN(cell)) |
279 |
|
{ |
280 |
|
|
281 |
0 |
cell = hValues[rowIndex] * hValues[columnIndex]; |
282 |
|
|
283 |
0 |
eigenMatrix.setValue(rowIndex, columnIndex, cell); |
284 |
0 |
eigenMatrix.setValue(columnIndex, rowIndex, cell); |
285 |
|
|
286 |
0 |
estimatedPositions.add(new int[] { rowIndex, columnIndex }); |
287 |
|
} |
288 |
|
} |
289 |
|
} |
290 |
|
|
291 |
|
|
292 |
0 |
for (int diagonalIndex = 0; diagonalIndex < matrixWidth; diagonalIndex++) |
293 |
|
{ |
294 |
0 |
double cell = Math.pow(hValues[diagonalIndex], 2); |
295 |
0 |
eigenMatrix.setValue(diagonalIndex, diagonalIndex, cell); |
296 |
0 |
estimatedPositions.add(new int[] { diagonalIndex, diagonalIndex }); |
297 |
|
} |
298 |
|
|
299 |
|
|
300 |
|
|
301 |
|
|
302 |
|
|
303 |
0 |
System.out.print("initial values: [ "); |
304 |
0 |
for (double h : hValues) |
305 |
|
{ |
306 |
0 |
System.out.print(String.format("%1.4f, ", h)); |
307 |
|
} |
308 |
0 |
System.out.println(" ]"); |
309 |
|
|
310 |
0 |
double[] hValuesOld = new double[matrixWidth]; |
311 |
|
|
312 |
0 |
int iterationCount = 0; |
313 |
|
|
314 |
|
|
315 |
|
|
316 |
|
|
317 |
0 |
while (true) |
318 |
|
{ |
319 |
0 |
for (int hIndex = 0; hIndex < matrixWidth; hIndex++) |
320 |
|
{ |
321 |
0 |
double newH = Arrays |
322 |
|
.stream(MiscMath.elementwiseMultiply(hValues, |
323 |
|
eigenMatrix.getRow(hIndex))) |
324 |
|
.sum() |
325 |
|
/ Arrays.stream( |
326 |
|
MiscMath.elementwiseMultiply(hValues, hValues)) |
327 |
|
.sum(); |
328 |
0 |
hValues[hIndex] = newH; |
329 |
|
} |
330 |
|
|
331 |
0 |
System.out.print(String.format("iteration %d: [ ", iterationCount)); |
332 |
0 |
for (double h : hValues) |
333 |
|
{ |
334 |
0 |
System.out.print(String.format("%1.4f, ", h)); |
335 |
|
} |
336 |
0 |
System.out.println(" ]"); |
337 |
|
|
338 |
|
|
339 |
0 |
for (int[] pair : estimatedPositions) |
340 |
|
{ |
341 |
0 |
double newVal = hValues[pair[0]] * hValues[pair[1]]; |
342 |
0 |
eigenMatrix.setValue(pair[0], pair[1], newVal); |
343 |
0 |
eigenMatrix.setValue(pair[1], pair[0], newVal); |
344 |
|
} |
345 |
|
|
346 |
0 |
iterationCount++; |
347 |
|
|
348 |
|
|
349 |
0 |
if (MiscMath.allClose(hValues, hValuesOld, 0d, 1e-05d, false)) |
350 |
|
{ |
351 |
0 |
break; |
352 |
|
} |
353 |
|
|
354 |
|
|
355 |
0 |
System.arraycopy(hValues, 0, hValuesOld, 0, hValues.length); |
356 |
|
} |
357 |
|
|
358 |
|
|
359 |
|
|
360 |
|
|
361 |
|
|
362 |
|
|
363 |
0 |
eigenMatrix.tred(); |
364 |
0 |
eigenMatrix.tqli(); |
365 |
|
|
366 |
0 |
System.out.println("eigenmatrix"); |
367 |
0 |
eigenMatrix.print(System.out, "%8.2f"); |
368 |
0 |
System.out.println(); |
369 |
0 |
System.out.println("uncorrected eigenvalues"); |
370 |
0 |
eigenMatrix.printD(System.out, "%2.4f "); |
371 |
0 |
System.out.println(); |
372 |
|
|
373 |
0 |
double[] eigenVals = eigenMatrix.getD(); |
374 |
|
|
375 |
0 |
TreeMap<Double, Integer> eigenPairs = new TreeMap<>( |
376 |
|
Comparator.reverseOrder()); |
377 |
0 |
for (int i = 0; i < eigenVals.length; i++) |
378 |
|
{ |
379 |
0 |
eigenPairs.put(eigenVals[i], i); |
380 |
|
} |
381 |
|
|
382 |
|
|
383 |
0 |
double[][] _repMatrix = new double[eigenVals.length][dim]; |
384 |
0 |
double[][] _oldMatrix = new double[eigenVals.length][dim]; |
385 |
0 |
double[] correctedEigenValues = new double[dim]; |
386 |
|
|
387 |
0 |
int l = 0; |
388 |
0 |
for (Entry<Double, Integer> pair : eigenPairs.entrySet()) |
389 |
|
{ |
390 |
0 |
double eigenValue = pair.getKey(); |
391 |
0 |
int column = pair.getValue(); |
392 |
0 |
double[] eigenVector = eigenMatrix.getColumn(column); |
393 |
|
|
394 |
0 |
if (l >= 1) |
395 |
|
{ |
396 |
0 |
eigenValue /= correctionFactor; |
397 |
|
} |
398 |
0 |
correctedEigenValues[l] = eigenValue; |
399 |
0 |
for (int j = 0; j < eigenVector.length; j++) |
400 |
|
{ |
401 |
0 |
_repMatrix[j][l] = (eigenValue < 0) ? 0.0 |
402 |
|
: -Math.sqrt(eigenValue) * eigenVector[j]; |
403 |
0 |
double tmpOldScore = scoresOld.getColumn(column)[j]; |
404 |
0 |
_oldMatrix[j][dim - l - 1] = (Double.isNaN(tmpOldScore)) ? 0.0 |
405 |
|
: tmpOldScore; |
406 |
|
} |
407 |
0 |
l++; |
408 |
0 |
if (l >= dim) |
409 |
|
{ |
410 |
0 |
break; |
411 |
|
} |
412 |
|
} |
413 |
|
|
414 |
0 |
System.out.println("correctedEigenValues"); |
415 |
0 |
MiscMath.print(correctedEigenValues, "%2.4f "); |
416 |
|
|
417 |
0 |
repMatrix = new Matrix(_repMatrix); |
418 |
0 |
repMatrix.setD(correctedEigenValues); |
419 |
0 |
MatrixI oldMatrix = new Matrix(_oldMatrix); |
420 |
|
|
421 |
0 |
MatrixI dotMatrix = repMatrix.postMultiply(repMatrix.transpose()); |
422 |
|
|
423 |
0 |
double rmsd = scoresOld.rmsd(dotMatrix); |
424 |
|
|
425 |
0 |
System.out.println("iteration, rmsd, maxDiff, rmsdDiff"); |
426 |
0 |
System.out.println(String.format("0, %8.5f, -, -", rmsd)); |
427 |
|
|
428 |
|
|
429 |
0 |
for (int iteration = 1; iteration < 21; iteration++) |
430 |
|
|
431 |
|
{ |
432 |
0 |
MatrixI repMatrixOLD = repMatrix.copy(); |
433 |
0 |
MatrixI dotMatrixOLD = dotMatrix.copy(); |
434 |
|
|
435 |
|
|
436 |
0 |
for (int hAIndex = 0; hAIndex < oldMatrix.height(); hAIndex++) |
437 |
|
{ |
438 |
0 |
double[] row = oldMatrix.getRow(hAIndex); |
439 |
0 |
double[] hA = repMatrix.getRow(hAIndex); |
440 |
0 |
hAIndex = hAIndex; |
441 |
|
|
442 |
|
|
443 |
0 |
double[] hAlsm = leastSquaresOptimisation(repMatrix, scoresOld, |
444 |
|
hAIndex); |
445 |
|
|
446 |
0 |
for (int j = 0; j < repMatrix.width(); j++) |
447 |
|
{ |
448 |
0 |
repMatrix.setValue(hAIndex, j, hAlsm[j]); |
449 |
|
} |
450 |
|
} |
451 |
|
|
452 |
|
|
453 |
|
|
454 |
0 |
dotMatrix = repMatrix.postMultiply(repMatrix.transpose()); |
455 |
|
|
456 |
0 |
rmsd = scoresOld.rmsd(dotMatrix); |
457 |
|
|
458 |
|
|
459 |
|
|
460 |
0 |
MatrixI diff = repMatrix.subtract(repMatrixOLD).absolute(); |
461 |
0 |
double maxDiff = 0.0; |
462 |
0 |
for (int i = 0; i < diff.height(); i++) |
463 |
|
{ |
464 |
0 |
for (int j = 0; j < diff.width(); j++) |
465 |
|
{ |
466 |
0 |
maxDiff = (diff.getValue(i, j) > maxDiff) ? diff.getValue(i, j) |
467 |
|
: maxDiff; |
468 |
|
} |
469 |
|
} |
470 |
|
|
471 |
|
|
472 |
0 |
double rmsdDiff = dotMatrix.rmsd(dotMatrixOLD); |
473 |
|
|
474 |
0 |
System.out.println(String.format("%d, %8.5f, %8.5f, %8.5f", |
475 |
|
iteration, rmsd, maxDiff, rmsdDiff)); |
476 |
|
|
477 |
0 |
if (!(Math.abs(maxDiff) > 1e-06)) |
478 |
|
{ |
479 |
0 |
repMatrix = repMatrixOLD.copy(); |
480 |
0 |
break; |
481 |
|
} |
482 |
|
} |
483 |
|
|
484 |
|
} catch (Exception q) |
485 |
|
{ |
486 |
0 |
Console.error("Error computing cc_analysis: " + q.getMessage()); |
487 |
0 |
q.printStackTrace(); |
488 |
|
} |
489 |
0 |
System.out.println("final coordinates:"); |
490 |
0 |
repMatrix.print(System.out, "%1.8f "); |
491 |
0 |
return repMatrix; |
492 |
|
} |
493 |
|
|
494 |
|
|
495 |
|
|
496 |
|
|
497 |
|
|
498 |
|
|
499 |
|
|
500 |
|
|
501 |
|
|
502 |
|
|
503 |
|
|
504 |
|
|
505 |
|
|
506 |
|
|
507 |
|
|
508 |
|
|
509 |
|
|
510 |
|
|
511 |
|
|
512 |
|
|
513 |
|
|
514 |
|
|
515 |
|
|
516 |
|
|
517 |
|
@param |
518 |
|
|
519 |
|
@param |
520 |
|
|
521 |
|
@param |
522 |
|
|
523 |
|
@param |
524 |
|
|
525 |
|
|
526 |
|
|
527 |
|
@return |
528 |
|
|
|
|
| 0% |
Uncovered Elements: 14 (14) |
Complexity: 3 |
Complexity Density: 0.3 |
|
529 |
0 |
private double[] originalToEquasionSystem(double[] hA, MatrixI repMatrix,... |
530 |
|
MatrixI scoresOld, int hAIndex) |
531 |
|
{ |
532 |
0 |
double[] originalRow = scoresOld.getRow(hAIndex); |
533 |
0 |
int nans = MiscMath.countNaN(originalRow); |
534 |
0 |
double[] result = new double[originalRow.length - nans]; |
535 |
|
|
536 |
|
|
537 |
0 |
int resultIndex = 0; |
538 |
0 |
for (int hBIndex = 0; hBIndex < originalRow.length; hBIndex++) |
539 |
|
{ |
540 |
0 |
double pairwiseCC = originalRow[hBIndex]; |
541 |
|
|
542 |
0 |
if (!Double.isNaN(pairwiseCC)) |
543 |
|
{ |
544 |
0 |
double[] hB = repMatrix.getRow(hBIndex); |
545 |
0 |
result[resultIndex++] = MiscMath |
546 |
|
.sum(MiscMath.elementwiseMultiply(hA, hB)) - pairwiseCC; |
547 |
|
} |
548 |
|
else |
549 |
|
{ |
550 |
|
} |
551 |
|
} |
552 |
0 |
return result; |
553 |
|
} |
554 |
|
|
555 |
|
|
556 |
|
|
557 |
|
|
558 |
|
@param |
559 |
|
|
560 |
|
@param |
561 |
|
|
562 |
|
|
563 |
|
@return |
564 |
|
|
|
|
| 0% |
Uncovered Elements: 37 (37) |
Complexity: 7 |
Complexity Density: 0.28 |
|
565 |
0 |
private MatrixI approximateDerivative(MatrixI repMatrix,... |
566 |
|
MatrixI scoresOld, int hAIndex) |
567 |
|
{ |
568 |
|
|
569 |
0 |
double[] hA = repMatrix.getRow(hAIndex); |
570 |
0 |
double[] f0 = originalToEquasionSystem(hA, repMatrix, scoresOld, |
571 |
|
hAIndex); |
572 |
0 |
double[] signX0 = new double[hA.length]; |
573 |
0 |
double[] xAbs = new double[hA.length]; |
574 |
0 |
for (int i = 0; i < hA.length; i++) |
575 |
|
{ |
576 |
0 |
signX0[i] = (hA[i] >= 0) ? 1 : -1; |
577 |
0 |
xAbs[i] = (Math.abs(hA[i]) >= 1.0) ? Math.abs(hA[i]) : 1.0; |
578 |
|
} |
579 |
0 |
double rstep = Math.pow(Math.ulp(1.0), 0.5); |
580 |
|
|
581 |
0 |
double[] h = new double[hA.length]; |
582 |
0 |
for (int i = 0; i < hA.length; i++) |
583 |
|
{ |
584 |
0 |
h[i] = rstep * signX0[i] * xAbs[i]; |
585 |
|
} |
586 |
|
|
587 |
0 |
int m = f0.length; |
588 |
0 |
int n = hA.length; |
589 |
0 |
double[][] jTransposed = new double[n][m]; |
590 |
0 |
for (int i = 0; i < h.length; i++) |
591 |
|
{ |
592 |
0 |
double[] x = new double[h.length]; |
593 |
0 |
System.arraycopy(hA, 0, x, 0, h.length); |
594 |
0 |
x[i] += h[i]; |
595 |
0 |
double dx = x[i] - hA[i]; |
596 |
0 |
double[] df = originalToEquasionSystem(x, repMatrix, scoresOld, |
597 |
|
hAIndex); |
598 |
0 |
for (int j = 0; j < df.length; j++) |
599 |
|
{ |
600 |
0 |
df[j] -= f0[j]; |
601 |
0 |
jTransposed[i][j] = df[j] / dx; |
602 |
|
} |
603 |
|
} |
604 |
0 |
MatrixI J = new Matrix(jTransposed).transpose(); |
605 |
0 |
return J; |
606 |
|
} |
607 |
|
|
608 |
|
|
609 |
|
|
610 |
|
|
611 |
|
@param |
612 |
|
@param |
613 |
|
@param |
614 |
|
@param |
615 |
|
|
616 |
|
@return |
617 |
|
|
|
|
| 0% |
Uncovered Elements: 5 (5) |
Complexity: 1 |
Complexity Density: 0.2 |
|
618 |
0 |
private double[] phiAndDerivative(double alpha, double[] suf, double[] s,... |
619 |
|
double Delta) |
620 |
|
{ |
621 |
0 |
double[] denom = MiscMath |
622 |
|
.elementwiseAdd(MiscMath.elementwiseMultiply(s, s), alpha); |
623 |
0 |
double pNorm = MiscMath.norm(MiscMath.elementwiseDivide(suf, denom)); |
624 |
0 |
double phi = pNorm - Delta; |
625 |
|
|
626 |
0 |
double phiPrime = -MiscMath.sum(MiscMath.elementwiseDivide( |
627 |
|
MiscMath.elementwiseMultiply(suf, suf), |
628 |
|
MiscMath.elementwiseMultiply( |
629 |
|
MiscMath.elementwiseMultiply(denom, denom), denom))) |
630 |
|
/ pNorm; |
631 |
0 |
return new double[] { phi, phiPrime }; |
632 |
|
} |
633 |
|
|
634 |
|
|
635 |
|
|
636 |
|
|
|
|
| 0% |
Uncovered Elements: 10 (10) |
Complexity: 4 |
Complexity Density: 0.67 |
|
637 |
|
private class TrustRegion |
638 |
|
{ |
639 |
|
private double[] step; |
640 |
|
|
641 |
|
private double alpha; |
642 |
|
|
643 |
|
private int iteration; |
644 |
|
|
|
|
| 0% |
Uncovered Elements: 3 (3) |
Complexity: 1 |
Complexity Density: 0.33 |
|
645 |
0 |
public TrustRegion(double[] step, double alpha, int iteration)... |
646 |
|
{ |
647 |
0 |
this.step = step; |
648 |
0 |
this.alpha = alpha; |
649 |
0 |
this.iteration = iteration; |
650 |
|
} |
651 |
|
|
|
|
| 0% |
Uncovered Elements: 1 (1) |
Complexity: 1 |
Complexity Density: 1 |
|
652 |
0 |
public double[] getStep()... |
653 |
|
{ |
654 |
0 |
return this.step; |
655 |
|
} |
656 |
|
|
|
|
| 0% |
Uncovered Elements: 1 (1) |
Complexity: 1 |
Complexity Density: 1 |
|
657 |
0 |
public double getAlpha()... |
658 |
|
{ |
659 |
0 |
return this.alpha; |
660 |
|
} |
661 |
|
|
|
|
| 0% |
Uncovered Elements: 1 (1) |
Complexity: 1 |
Complexity Density: 1 |
|
662 |
0 |
public int getIteration()... |
663 |
|
{ |
664 |
0 |
return this.iteration; |
665 |
|
} |
666 |
|
} |
667 |
|
|
668 |
|
|
669 |
|
|
670 |
|
|
671 |
|
@param |
672 |
|
|
673 |
|
@param |
674 |
|
|
675 |
|
@param |
676 |
|
@param |
677 |
|
|
678 |
|
@param |
679 |
|
|
680 |
|
@param |
681 |
|
|
682 |
|
@param |
683 |
|
|
684 |
|
|
685 |
|
@return |
686 |
|
|
|
|
| 0% |
Uncovered Elements: 52 (52) |
Complexity: 12 |
Complexity Density: 0.35 |
|
687 |
0 |
private TrustRegion solveLsqTrustRegion(int n, int m, double[] uf,... |
688 |
|
double[] s, MatrixI V, double Delta, double alpha) |
689 |
|
{ |
690 |
0 |
double[] suf = MiscMath.elementwiseMultiply(s, uf); |
691 |
|
|
692 |
|
|
693 |
0 |
boolean fullRank = false; |
694 |
0 |
if (m >= n) |
695 |
|
{ |
696 |
0 |
double threshold = s[0] * Math.ulp(1.0) * m; |
697 |
0 |
fullRank = s[s.length - 1] > threshold; |
698 |
|
} |
699 |
0 |
if (fullRank) |
700 |
|
{ |
701 |
0 |
double[] p = MiscMath.elementwiseMultiply( |
702 |
|
V.sumProduct(MiscMath.elementwiseDivide(uf, s)), -1); |
703 |
0 |
if (MiscMath.norm(p) <= Delta) |
704 |
|
{ |
705 |
0 |
TrustRegion result = new TrustRegion(p, 0.0, 0); |
706 |
0 |
return result; |
707 |
|
} |
708 |
|
} |
709 |
|
|
710 |
0 |
double alphaUpper = MiscMath.norm(suf) / Delta; |
711 |
0 |
double alphaLower = 0.0; |
712 |
0 |
if (fullRank) |
713 |
|
{ |
714 |
0 |
double[] phiAndPrime = phiAndDerivative(0.0, suf, s, Delta); |
715 |
0 |
alphaLower = -phiAndPrime[0] / phiAndPrime[1]; |
716 |
|
} |
717 |
|
|
718 |
0 |
alpha = (!fullRank && alpha == 0.0) |
719 |
|
? alpha = Math.max(0.001 * alphaUpper, |
720 |
|
Math.pow(alphaLower * alphaUpper, 0.5)) |
721 |
|
: alpha; |
722 |
|
|
723 |
0 |
int iteration = 0; |
724 |
0 |
while (iteration < 10) |
725 |
|
{ |
726 |
0 |
alpha = (alpha < alphaLower || alpha > alphaUpper) |
727 |
|
? alpha = Math.max(0.001 * alphaUpper, |
728 |
|
Math.pow(alphaLower * alphaUpper, 0.5)) |
729 |
|
: alpha; |
730 |
0 |
double[] phiAndPrime = phiAndDerivative(alpha, suf, s, Delta); |
731 |
0 |
double phi = phiAndPrime[0]; |
732 |
0 |
double phiPrime = phiAndPrime[1]; |
733 |
|
|
734 |
0 |
alphaUpper = (phi < 0) ? alpha : alphaUpper; |
735 |
0 |
double ratio = phi / phiPrime; |
736 |
0 |
alphaLower = Math.max(alphaLower, alpha - ratio); |
737 |
0 |
alpha -= (phi + Delta) * ratio / Delta; |
738 |
|
|
739 |
0 |
if (Math.abs(phi) < 0.01 * Delta) |
740 |
|
{ |
741 |
0 |
break; |
742 |
|
} |
743 |
0 |
iteration++; |
744 |
|
} |
745 |
|
|
746 |
|
|
747 |
0 |
double[] tmp = MiscMath.elementwiseDivide(suf, MiscMath |
748 |
|
.elementwiseAdd(MiscMath.elementwiseMultiply(s, s), alpha)); |
749 |
0 |
double[] p = MiscMath.elementwiseMultiply(V.sumProduct(tmp), -1); |
750 |
|
|
751 |
|
|
752 |
|
|
753 |
|
|
754 |
0 |
p = MiscMath.elementwiseMultiply(p, Delta / MiscMath.norm(p)); |
755 |
|
|
756 |
0 |
TrustRegion result = new TrustRegion(p, alpha, iteration + 1); |
757 |
0 |
return result; |
758 |
|
} |
759 |
|
|
760 |
|
|
761 |
|
|
762 |
|
|
763 |
|
|
764 |
|
@param |
765 |
|
|
766 |
|
@param |
767 |
|
|
768 |
|
@param |
769 |
|
|
770 |
|
|
771 |
|
@return |
772 |
|
|
|
|
| 0% |
Uncovered Elements: 4 (4) |
Complexity: 1 |
Complexity Density: 0.25 |
|
773 |
0 |
private double evaluateQuadratic(MatrixI J, double[] g, double[] s)... |
774 |
|
{ |
775 |
|
|
776 |
0 |
double[] Js = J.sumProduct(s); |
777 |
0 |
double q = MiscMath.dot(Js, Js); |
778 |
0 |
double l = MiscMath.dot(s, g); |
779 |
|
|
780 |
0 |
return 0.5 * q + l; |
781 |
|
} |
782 |
|
|
783 |
|
|
784 |
|
|
785 |
|
|
786 |
|
@param |
787 |
|
@param |
788 |
|
@param |
789 |
|
@param |
790 |
|
@param |
791 |
|
|
792 |
|
@return |
793 |
|
|
|
|
| 0% |
Uncovered Elements: 19 (19) |
Complexity: 7 |
Complexity Density: 0.64 |
|
794 |
0 |
private double[] updateTrustRegionRadius(double Delta,... |
795 |
|
double actualReduction, double predictedReduction, |
796 |
|
double stepNorm, boolean boundHit) |
797 |
|
{ |
798 |
0 |
double ratio = 0; |
799 |
0 |
if (predictedReduction > 0) |
800 |
|
{ |
801 |
0 |
ratio = actualReduction / predictedReduction; |
802 |
|
} |
803 |
0 |
else if (predictedReduction == 0 && actualReduction == 0) |
804 |
|
{ |
805 |
0 |
ratio = 1; |
806 |
|
} |
807 |
|
else |
808 |
|
{ |
809 |
0 |
ratio = 0; |
810 |
|
} |
811 |
|
|
812 |
0 |
if (ratio < 0.25) |
813 |
|
{ |
814 |
0 |
Delta = 0.25 * stepNorm; |
815 |
|
} |
816 |
0 |
else if (ratio > 0.75 && boundHit) |
817 |
|
{ |
818 |
0 |
Delta *= 2.0; |
819 |
|
} |
820 |
|
|
821 |
0 |
return new double[] { Delta, ratio }; |
822 |
|
} |
823 |
|
|
824 |
|
|
825 |
|
|
826 |
|
|
827 |
|
@param |
828 |
|
|
829 |
|
@param |
830 |
|
|
831 |
|
@param |
832 |
|
|
833 |
|
@param |
834 |
|
|
835 |
|
|
836 |
|
@return |
837 |
|
|
|
|
| 0% |
Uncovered Elements: 72 (72) |
Complexity: 8 |
Complexity Density: 0.13 |
|
838 |
0 |
private double[] trf(MatrixI repMatrix, MatrixI scoresOld, int index,... |
839 |
|
MatrixI J) |
840 |
|
{ |
841 |
|
|
842 |
0 |
double[] hA = repMatrix.getRow(index); |
843 |
0 |
double[] f0 = originalToEquasionSystem(hA, repMatrix, scoresOld, index); |
844 |
0 |
int nfev = 1; |
845 |
0 |
int m = J.height(); |
846 |
0 |
int n = J.width(); |
847 |
0 |
double cost = 0.5 * MiscMath.dot(f0, f0); |
848 |
0 |
double[] g = J.transpose().sumProduct(f0); |
849 |
0 |
double Delta = MiscMath.norm(hA); |
850 |
0 |
int maxNfev = hA.length * 100; |
851 |
0 |
double alpha = 0.0; |
852 |
|
|
853 |
0 |
double gNorm = 0; |
854 |
0 |
boolean terminationStatus = false; |
855 |
0 |
int iteration = 0; |
856 |
|
|
857 |
0 |
while (true) |
858 |
|
{ |
859 |
0 |
gNorm = MiscMath.norm(g); |
860 |
0 |
if (terminationStatus || nfev == maxNfev) |
861 |
|
{ |
862 |
0 |
break; |
863 |
|
} |
864 |
0 |
SingularValueDecomposition svd = new SingularValueDecomposition( |
865 |
|
new Array2DRowRealMatrix(J.asArray())); |
866 |
0 |
MatrixI U = new Matrix(svd.getU().getData()); |
867 |
0 |
double[] s = svd.getSingularValues(); |
868 |
0 |
MatrixI V = new Matrix(svd.getV().getData()).transpose(); |
869 |
0 |
double[] uf = U.transpose().sumProduct(f0); |
870 |
|
|
871 |
0 |
double actualReduction = -1; |
872 |
0 |
double[] xNew = new double[hA.length]; |
873 |
0 |
double[] fNew = new double[f0.length]; |
874 |
0 |
double costNew = 0; |
875 |
0 |
double stepHnorm = 0; |
876 |
|
|
877 |
0 |
while (actualReduction <= 0 && nfev < maxNfev) |
878 |
|
{ |
879 |
0 |
TrustRegion trustRegion = solveLsqTrustRegion(n, m, uf, s, V, Delta, |
880 |
|
alpha); |
881 |
0 |
double[] stepH = trustRegion.getStep(); |
882 |
0 |
alpha = trustRegion.getAlpha(); |
883 |
0 |
int nIterations = trustRegion.getIteration(); |
884 |
0 |
double predictedReduction = -(evaluateQuadratic(J, g, stepH)); |
885 |
|
|
886 |
0 |
xNew = MiscMath.elementwiseAdd(hA, stepH); |
887 |
0 |
fNew = originalToEquasionSystem(xNew, repMatrix, scoresOld, index); |
888 |
0 |
nfev++; |
889 |
|
|
890 |
0 |
stepHnorm = MiscMath.norm(stepH); |
891 |
|
|
892 |
0 |
if (MiscMath.countNaN(fNew) > 0) |
893 |
|
{ |
894 |
0 |
Delta = 0.25 * stepHnorm; |
895 |
0 |
continue; |
896 |
|
} |
897 |
|
|
898 |
|
|
899 |
0 |
costNew = 0.5 * MiscMath.dot(fNew, fNew); |
900 |
0 |
actualReduction = cost - costNew; |
901 |
|
|
902 |
0 |
double[] updatedTrustRegion = updateTrustRegionRadius(Delta, |
903 |
|
actualReduction, predictedReduction, stepHnorm, |
904 |
|
stepHnorm > (0.95 * Delta)); |
905 |
0 |
double DeltaNew = updatedTrustRegion[0]; |
906 |
0 |
double ratio = updatedTrustRegion[1]; |
907 |
|
|
908 |
|
|
909 |
0 |
boolean ftolSatisfied = actualReduction < (1e-8 * cost) |
910 |
|
&& ratio > 0.25; |
911 |
0 |
boolean xtolSatisfied = stepHnorm < (1e-8 |
912 |
|
* (1e-8 + MiscMath.norm(hA))); |
913 |
0 |
terminationStatus = ftolSatisfied || xtolSatisfied; |
914 |
0 |
if (terminationStatus) |
915 |
|
{ |
916 |
0 |
break; |
917 |
|
} |
918 |
|
|
919 |
0 |
alpha *= Delta / DeltaNew; |
920 |
0 |
Delta = DeltaNew; |
921 |
|
|
922 |
|
} |
923 |
0 |
if (actualReduction > 0) |
924 |
|
{ |
925 |
0 |
hA = xNew; |
926 |
0 |
f0 = fNew; |
927 |
0 |
cost = costNew; |
928 |
|
|
929 |
0 |
J = approximateDerivative(repMatrix, scoresOld, index); |
930 |
|
|
931 |
0 |
g = J.transpose().sumProduct(f0); |
932 |
|
} |
933 |
|
else |
934 |
|
{ |
935 |
0 |
stepHnorm = 0; |
936 |
0 |
actualReduction = 0; |
937 |
|
} |
938 |
0 |
iteration++; |
939 |
|
} |
940 |
|
|
941 |
0 |
return hA; |
942 |
|
} |
943 |
|
|
944 |
|
|
945 |
|
|
946 |
|
|
947 |
|
|
948 |
|
@param |
949 |
|
|
950 |
|
@param |
951 |
|
|
952 |
|
@param |
953 |
|
|
954 |
|
|
955 |
|
@return |
956 |
|
|
|
|
| 0% |
Uncovered Elements: 3 (3) |
Complexity: 1 |
Complexity Density: 0.33 |
|
957 |
0 |
private double[] leastSquaresOptimisation(MatrixI repMatrix,... |
958 |
|
MatrixI scoresOld, int index) |
959 |
|
{ |
960 |
0 |
MatrixI J = approximateDerivative(repMatrix, scoresOld, index); |
961 |
0 |
double[] result = trf(repMatrix, scoresOld, index, J); |
962 |
0 |
return result; |
963 |
|
} |
964 |
|
|
965 |
|
} |