From: MorellThomas Date: Sat, 16 Sep 2023 10:56:12 +0000 (+0200) Subject: Almost working X-Git-Tag: Release_2_11_4_0~31^2~36 X-Git-Url: http://source.jalview.org/gitweb/?a=commitdiff_plain;h=f55d511de05a6b47542053778b19bc6c5dd40b4e;p=jalview.git Almost working --- diff --git a/src/jalview/analysis/AlignSeq.java b/src/jalview/analysis/AlignSeq.java index 4ec0457..7165096 100755 --- a/src/jalview/analysis/AlignSeq.java +++ b/src/jalview/analysis/AlignSeq.java @@ -58,11 +58,11 @@ public class AlignSeq private static final int MAX_NAME_LENGTH = 30; //&! - private static final int GAP_OPEN_COST = 120; - //private static final int GAP_OPEN_COST = 100; + //private static final int GAP_OPEN_COST = 120; + private static final int GAP_OPEN_COST = 100; - private static final int GAP_EXTEND_COST = 20; - //private static final int GAP_EXTEND_COST = 5; + //private static final int GAP_EXTEND_COST = 20; + private static final int GAP_EXTEND_COST = 5; private static final int GAP_INDEX = -1; @@ -508,7 +508,7 @@ public class AlignSeq int trace; maxscore = score[i][j] / 10f; - //&! get trailing gaps + //prepare trailing gaps while ((i < seq1.length - 1) || (j < seq2.length - 1)) { i++; @@ -526,7 +526,7 @@ public class AlignSeq count = (seq1.length + seq2.length) - 1; - //&! get trailing gaps + //get trailing gaps while ((i >= seq1.length) || (j >= seq2.length)) { if (i >= seq1.length) @@ -591,18 +591,18 @@ public class AlignSeq sb2.append(s2str.charAt(j)); } - //&! get initial gaps + //get initial gaps while (j > 0 || i > 0) { if (j > 0) { + j--; sb1.append("-"); sb2.append(s2str.charAt(j)); - j--; } else if (i > 0) { + i--; sb1.append(s1str.charAt(i)); sb2.append("-"); - i--; } } @@ -1436,7 +1436,8 @@ public class AlignSeq float prescore = score; // only debug score *= coverage; - System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %f, mean: %f, max: %d", preprescore, prescore, _max[1], score, this.meanScore, this.hypotheticMaxScore)); - this.alignmentScore = (preprescore < 1) ? Float.NaN : score; + System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %1.16f, mean: %f, max: %d", preprescore, prescore, _max[1], score, this.meanScore, this.hypotheticMaxScore)); + float minScore = 1f; + this.alignmentScore = (preprescore < minScore) ? Float.NaN : score; } } diff --git a/src/jalview/analysis/PaSiMap.java b/src/jalview/analysis/PaSiMap.java index 77c28be..4235779 100755 --- a/src/jalview/analysis/PaSiMap.java +++ b/src/jalview/analysis/PaSiMap.java @@ -59,8 +59,6 @@ public class PaSiMap implements Runnable */ private MatrixI pairwiseScores; - private MatrixI tridiagonal; - private MatrixI eigenMatrix; /** @@ -71,8 +69,6 @@ public class PaSiMap implements Runnable * @param sm * @param options */ - //public PaSiMap(AlignmentView sequences, ScoreModelI sm, - //&! viewport or panel? public PaSiMap(AlignmentViewport sequences, ScoreModelI sm, SimilarityParamsI options) { @@ -95,7 +91,7 @@ public class PaSiMap implements Runnable } /** - * DOCUMENT ME! + * Returns coordinates for each datapoint * * @param l * DOCUMENT ME! @@ -103,8 +99,7 @@ public class PaSiMap implements Runnable * DOCUMENT ME! * @param mm * DOCUMENT ME! - * @param factor - * DOCUMENT ME! + * @param factor ~ is 1 * * @return DOCUMENT ME! */ @@ -162,19 +157,6 @@ public class PaSiMap implements Runnable */ double component(int row, int n) { - /* - double out = 0.0; - - for (int i = 0; i < pairwiseScores.width(); i++) - { - double pairwiseScore = pairwiseScores.getValue(row, i); - double eigenScore = eigenMatrix.getValue(i, n); - out += (pairwiseScores.getValue(row, i) * eigenMatrix.getValue(i, n)); - } - - return out / eigenMatrix.getD()[n]; - */ - System.out.println(String.format("row %d, col %d", row, n)); return eigenMatrix.getValue(row, n); } @@ -198,19 +180,6 @@ public class PaSiMap implements Runnable pairwiseScores.print(ps, "%8.2f"); /* - * tridiagonal matrix, with D and E vectors - */ - /* - sb.append(" ---Tridiag transform matrix ---\n"); - sb.append(" --- D vector ---\n"); - tridiagonal.printD(ps, "%15.4e"); - ps.println(); - sb.append("--- E vector ---\n"); - tridiagonal.printE(ps, "%15.4e"); - ps.println(); - */ - - /* * eigenvalues matrix, with D vector */ sb.append(" --- New diagonalization matrix ---\n"); @@ -225,10 +194,10 @@ public class PaSiMap implements Runnable /** * Performs the PaSiMap calculation * - * creates a new gui/PairwiseAlignPanel with the input sequences (<++>/AlignmentViewport) + * creates a new gui/PairwiseAlignPanel with the input sequences (AlignmentViewport) * uses analysis/AlignSeq to creatue the pairwise alignments and calculate the AlignmentScores (float for each pair) * gets all float[][] scores from the gui/PairwiseAlignPanel - * checks the connections for each sequence with <++>/AlignmentViewport seqs.calculateConnectivity(float[][] scores, int dim) (from analysis/Connectivity) -- throws an Exception if insufficient + * checks the connections for each sequence with AlignmentViewport seqs.calculateConnectivity(float[][] scores, int dim) (from analysis/Connectivity) -- throws an Exception if insufficient * creates a math/MatrixI pairwiseScores of the float[][] scores * copys the scores and fills the diagonal to create a symmetric matrix using math/Matrix.fillDiagonal() * performs the analysis/ccAnalysis with the symmetric matrix @@ -239,8 +208,6 @@ public class PaSiMap implements Runnable { try { - // run needleman regardless if aligned or not - // gui.PairwiseAlignPanel <++> PairwiseAlignPanel alignment = new PairwiseAlignPanel(seqs, true); float[][] scores = alignment.getAlignmentScores(); //bigger index first -- eg scores[14][13] @@ -251,13 +218,10 @@ public class PaSiMap implements Runnable ccAnalysis cc = new ccAnalysis(pairwiseScores, dim); pairwiseScores = cc.run(); - tridiagonal = pairwiseScores.copy(); - //tridiagonal.tred(); /** perform the eigendecomposition for the plot */ - eigenMatrix = tridiagonal.copy(); + eigenMatrix = pairwiseScores.copy(); eigenMatrix.setD(pairwiseScores.getD()); - //eigenMatrix.tqli(); System.out.println(getDetails()); @@ -343,14 +307,4 @@ public class PaSiMap implements Runnable { eigenMatrix = m; } - - public MatrixI getTridiagonal() - { - return tridiagonal; - } - - public void setTridiagonal(MatrixI tridiagonal) - { - this.tridiagonal = tridiagonal; - } } diff --git a/src/jalview/analysis/ccAnalysis.java b/src/jalview/analysis/ccAnalysis.java index 5e149b8..f7d5c58 100755 --- a/src/jalview/analysis/ccAnalysis.java +++ b/src/jalview/analysis/ccAnalysis.java @@ -46,20 +46,10 @@ import java.lang.System; import java.util.Arrays; import java.util.ArrayList; import java.util.Comparator; -import java.util.HashSet; import java.util.Map.Entry; import java.util.TreeMap; -import org.apache.commons.math3.analysis.MultivariateVectorFunction; -import org.apache.commons.math3.analysis.MultivariateMatrixFunction; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder; -import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer; -import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.Array2DRowRealMatrix; -import org.apache.commons.math3.linear.RealVector; -import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.SingularValueDecomposition; /** @@ -73,7 +63,7 @@ public class ccAnalysis public ccAnalysis(MatrixI scores, byte dim) { - //&! round matrix to .4f to be same as in pasimap + // round matrix to .4f to be same as in pasimap for (int i = 0; i < scores.height(); i++) { for (int j = 0; j < scores.width(); j++) @@ -88,8 +78,9 @@ public class ccAnalysis this.dim = dim; } - /** TODO - * DOCUMENT ME + /** + * Initialise a distrust-score for each hypothesis (h) of hSigns + * distrust = conHypNum - proHypNum * * @param hSigns ~ hypothesis signs (+/-) for each sequence * @param scores ~ input score matrix @@ -127,8 +118,8 @@ public class ccAnalysis return distrustScores; } - /** TODO - * DOCUMENT ME + /** + * Optemise hypothesis concerning the sign of the hypothetical value for each hSigns by interpreting the pairwise correlation coefficients as scalar products * * @param hSigns ~ hypothesis signs (+/-) * @param distrustScores @@ -158,7 +149,7 @@ public class ccAnalysis if (NOTmaxDistrustIndex != maxDistrustIndex) { byte hBSign = hSigns[NOTmaxDistrustIndex]; - double cell = scores.getRow(maxDistrustIndex)[NOTmaxDistrustIndex]; + double cell = scores.getValue(maxDistrustIndex, NOTmaxDistrustIndex); // distrust only changed if not NaN if (!Double.isNaN(cell)) @@ -193,8 +184,9 @@ public class ccAnalysis * * @return */ - public MatrixI run () + public MatrixI run () throws Exception { + //initialse eigenMatrix and repMatrix MatrixI eigenMatrix = scoresOld.copy(); MatrixI repMatrix = scoresOld.copy(); try @@ -207,9 +199,9 @@ public class ccAnalysis * products. */ - //&! debug - System.out.println("input:"); + System.out.println("Input correlation matrix:"); eigenMatrix.print(System.out, "%1.4f "); + int matrixWidth = eigenMatrix.width(); // square matrix, so width == height int matrixElementsTotal = (int) Math.pow(matrixWidth, 2); //total number of elemts @@ -238,14 +230,11 @@ public class ccAnalysis //Estimate absolute values for each h_i by determining sqrt of mean of //non-NaN absolute values for every row. - //<++> hAbs = //np.sqrt(np.nanmean(np.absolute(eigenMatrix), axis=1)) - MiscMath.print(eigenMatrix.absolute().meanRow(), "%1.8f"); - double[] hAbs = MiscMath.sqrt(eigenMatrix.absolute().meanRow()); //np.sqrt(np.nanmean(np.absolute(eigenMatrix), axis=1)) + double[] hAbs = MiscMath.sqrt(eigenMatrix.absolute().meanRow()); //Combine estimated signs with absolute values in obtain total value for //each h_i. double[] hValues = MiscMath.elementwiseMultiply(hSigns, hAbs); - //<++>hValues.reshape((1,matrixWidth)); // doesnt it already look like this /*Complement symmetric matrix by using the scalar products of estimated *values of h_i to replace NaN-cells. @@ -265,7 +254,7 @@ public class ccAnalysis if (Double.isNaN(cell)) { //calculate scalar product as new cell value - cell = hValues[rowIndex] * hValues[columnIndex]; // something is wrong with hAbs andhValues and everything!!!!!!!!!!!! + cell = hValues[rowIndex] * hValues[columnIndex]; //fill in new value in cell and symmetric partner eigenMatrix.setValue(rowIndex, columnIndex, cell); eigenMatrix.setValue(columnIndex, rowIndex, cell); @@ -304,7 +293,6 @@ public class ccAnalysis { for (int hIndex = 0; hIndex < matrixWidth; hIndex++) { - //@python newH = np.sum(hValues * eigenMatrix[hIndex]) / np.sum(hValues ** 2) double newH = Arrays.stream(MiscMath.elementwiseMultiply(hValues, eigenMatrix.getRow(hIndex))).sum() / Arrays.stream(MiscMath.elementwiseMultiply(hValues, hValues)).sum(); hValues[hIndex] = newH; } @@ -339,31 +327,19 @@ public class ccAnalysis //----------------------------- //Use complemented symmetric matrix to calculate final representative //vectors. - //&! debug - System.out.println("after estimating:"); - eigenMatrix.print(System.out, "%1.8f "); //Eigendecomposition. eigenMatrix.tred(); - System.out.println("tred"); - eigenMatrix.print(System.out, "%8.2f"); - eigenMatrix.tqli(); - System.out.println("eigenvals"); - eigenMatrix.printD(System.out, "%2.4f "); - System.out.println(); - System.out.println("tqli"); + + System.out.println("eigenmatrix"); eigenMatrix.print(System.out, "%8.2f"); + System.out.println(); + System.out.println("uncorrected eigenvalues"); + eigenMatrix.printD(System.out, "%2.4f "); double[] eigenVals = eigenMatrix.getD(); - /* - TreeMap eigenPairs = new TreeMap<>(Comparator.reverseOrder()); - for (int i = 0; i < eigenVals.length; i++) - { - eigenPairs.put(eigenVals[i], eigenMatrix.getColumn(i)); - } - */ TreeMap eigenPairs = new TreeMap<>(Comparator.reverseOrder()); for (int i = 0; i < eigenVals.length; i++) { @@ -371,11 +347,10 @@ public class ccAnalysis } // matrix of representative eigenvectors (each row is a vector) - double[][] _repMatrix = new double[eigenVals.length][dim]; //last ones were dim + double[][] _repMatrix = new double[eigenVals.length][dim]; double[][] _oldMatrix = new double[eigenVals.length][dim]; double[] correctedEigenValues = new double[dim]; - //for (Entry pair : eigenPairs.entrySet()) int l = 0; for (Entry pair : eigenPairs.entrySet()) { @@ -387,12 +362,9 @@ public class ccAnalysis { eigenValue /= correctionFactor; } - //l++; - //correctedEigenValues[dim - l] = eigenValue; correctedEigenValues[l] = eigenValue; for (int j = 0; j < eigenVector.length; j++) { - //_repMatrix[j][dim - l] = (eigenValue < 0) ? 0.0 : Math.sqrt(eigenValue) * eigenVector[j]; _repMatrix[j][l] = (eigenValue < 0) ? 0.0 : - Math.sqrt(eigenValue) * eigenVector[j]; double tmpOldScore = scoresOld.getColumn(column)[j]; _oldMatrix[j][dim - l - 1] = (Double.isNaN(tmpOldScore)) ? 0.0 : tmpOldScore; @@ -409,23 +381,11 @@ public class ccAnalysis repMatrix = new Matrix(_repMatrix); repMatrix.setD(correctedEigenValues); - MatrixI oldMatrix = new Matrix(_oldMatrix); //TODO do i even need it anymore? - - System.out.println("old matrix"); - oldMatrix.print(System.out, "%8.2f"); - - System.out.println("scoresOld"); - scoresOld.print(System.out, "%1.4f "); - - System.out.println("rep matrix"); - repMatrix.print(System.out, "%1.8f "); + MatrixI oldMatrix = new Matrix(_oldMatrix); MatrixI dotMatrix = repMatrix.postMultiply(repMatrix.transpose()); - System.out.println("dot matrix"); - dotMatrix.print(System.out, "%1.8f "); - double rmsd = scoresOld.rmsd(dotMatrix); //TODO do i need this here? - System.out.println(rmsd); + double rmsd = scoresOld.rmsd(dotMatrix); System.out.println("iteration, rmsd, maxDiff, rmsdDiff"); System.out.println(String.format("0, %8.5f, -, -", rmsd)); @@ -439,20 +399,15 @@ public class ccAnalysis for (int hAIndex = 0; hAIndex < oldMatrix.height(); hAIndex++) { double[] row = oldMatrix.getRow(hAIndex); - double[] hA = repMatrix.getRow(hAIndex); // inverted + double[] hA = repMatrix.getRow(hAIndex); hAIndex = hAIndex; //find least-squares-solution fo rdifferences between original scores and representative vectors - //--> originalToEquasionSystem(hA, hAIndex, repMatrix, row) --> double[] - System.out.println(String.format("||||||||||||||||||||||||||||||||||||||||||||\nIteration: %d", iteration)); - //repMatrix = new Matrix( new double[][]{{ 0.92894902, -0.25013783, -0.0051076 }, { 0.91955135, -0.25024707, -0.00516568}, { 0.90957348, -0.21717002, -0.14259899}, { 0.90298063, -0.21678816, -0.14697814}, { 0.9157065, -0.04646437, 0.10105454}, { 0.9050301, -0.04689785, 0.18964432}, { 0.92498545, 0.03881933, 0.14771523}, { 0.88008842, 0.0142395, 0.11356242}, { 0.94276528, 0.25591474, -0.07190911}, { 0.93939976, 0.23375136, -0.07530558}, { 0.93550782, 0.24635804, -0.07317563}, { 0.92497731, 0.21729928, -0.02361162}}); - //scoresOld = new Matrix( new double[][]{{ Double.NaN, 0.9914, 0.8879, 0.8803, 0.8528, 0.8481, 0.8434, 0.8174, 0.8174, 0.8232, 0.8148, 0.8114}, {0.9914, Double.NaN, 0.8792, 0.8717, 0.8441, 0.8395, 0.8347, 0.8087, 0.8085, 0.8143, 0.8059, 0.8024}, {0.8879, 0.8792, Double.NaN, 0.9578, 0.8289, 0.8075, 0.8268, 0.7815, 0.8165, 0.8168, 0.8084, 0.7974}, {0.8803, 0.8717, 0.9578, Double.NaN, 0.838, 0.7974, 0.8058, 0.7798, 0.8094, 0.8098, 0.8067, 0.7905}, {0.8528, 0.8441, 0.8289, 0.838, Double.NaN, 0.8954, 0.8412, 0.7879, 0.8418, 0.8384, 0.8389, 0.8556}, {0.8481, 0.8395, 0.8075, 0.7974, 0.8954, Double.NaN, 0.8699, 0.8106, 0.8267, 0.8234, 0.8222, 0.8241}, {0.8434, 0.8347, 0.8268, 0.8058, 0.8412, 0.8699, Double.NaN, 0.869, 0.8745, 0.8583, 0.8661, 0.8593}, {0.8174, 0.8087, 0.7815, 0.7798, 0.7879, 0.8106, 0.869, Double.NaN, 0.8273, 0.8331, 0.8137, 0.8029}, {0.8174, 0.8085, 0.8165, 0.8094, 0.8418, 0.8267, 0.8745, 0.8273, Double.NaN, 0.967, 0.978, 0.9373}, {0.8232, 0.8143, 0.8168, 0.8098, 0.8384, 0.8234, 0.8583, 0.8331, 0.967, Double.NaN, 0.9561, 0.9337}, {0.8148, 0.8059, 0.8084, 0.8067, 0.8389, 0.8222, 0.8661, 0.8137, 0.978, 0.9561, Double.NaN, 0.9263}, {0.8114, 0.8024, 0.7974, 0.7905, 0.8556, 0.8241, 0.8593, 0.8029, 0.9373, 0.9337, 0.9263, Double.NaN}}); double[] hAlsm = leastSquaresOptimisation(repMatrix, scoresOld, hAIndex); // update repMatrix with new hAlsm for (int j = 0; j < repMatrix.width(); j++) { repMatrix.setValue(hAIndex, j, hAlsm[j]); } - break; } // dot product of representative vecotrs yields a matrix with values approximating the correlation matrix @@ -461,7 +416,6 @@ public class ccAnalysis rmsd = scoresOld.rmsd(dotMatrix); // calculate maximum change of representative vectors of current iteration - //repMatrix.subtract(repMatrixOLD).print(System.out, "%8.2f"); MatrixI diff = repMatrix.subtract(repMatrixOLD).absolute(); double maxDiff = 0.0; for (int i = 0; i < diff.height(); i++) @@ -471,7 +425,6 @@ public class ccAnalysis maxDiff = (diff.getValue(i, j) > maxDiff) ? diff.getValue(i, j) : maxDiff; } } - System.out.println(String.format("maxDiff: %f", maxDiff)); // calculate rmsd between current and previous estimation double rmsdDiff = dotMatrix.rmsd(dotMatrixOLD); @@ -483,7 +436,6 @@ public class ccAnalysis repMatrix = repMatrixOLD.copy(); break; } - break; } @@ -492,8 +444,7 @@ public class ccAnalysis Console.error("Error computing cc_analysis: " + q.getMessage()); q.printStackTrace(); } - //repMatrix = repMatrix.transpose(); - System.out.println("final repMatrix"); + System.out.println("final coordinates:"); repMatrix.print(System.out, "%8.2f"); return repMatrix; } @@ -521,7 +472,6 @@ public class ccAnalysis * * @return */ - //private double[] originalToEquasionSystem(double[] x, int hAIndex, MatrixI h, MatrixI originalScores) private double[] originalToEquasionSystem(double[] hA, MatrixI repMatrix, MatrixI scoresOld, int hAIndex) { double[] originalRow = scoresOld.getRow(hAIndex); @@ -556,11 +506,6 @@ public class ccAnalysis //hA = x0 double[] hA = repMatrix.getRow(hAIndex); double[] f0 = originalToEquasionSystem(hA, repMatrix, scoresOld, hAIndex); - System.out.println("Approximate derivative with "); - System.out.print("hA (x): "); - MiscMath.print(hA, "%1.8f"); - System.out.print("f0: "); - MiscMath.print(f0, "%1.8f"); double[] signX0 = new double[hA.length]; double[] xAbs = new double[hA.length]; for (int i = 0; i < hA.length; i++) @@ -592,7 +537,7 @@ public class ccAnalysis jTransposed[i][j] = df[j] / dx; } } - MatrixI J = new Matrix(jTransposed).transpose(); // inverted + MatrixI J = new Matrix(jTransposed).transpose(); return J; } @@ -651,7 +596,7 @@ public class ccAnalysis * solve a trust-region problem arising in least-squares optimisation * @param n ~ number of variables * @param m ~ number of residuals - * @param uf ~ <++> + * @param uf * @param s ~ singular values of J * @param V ~ transpose of VT * @param Delta ~ radius of a trust region @@ -672,7 +617,7 @@ public class ccAnalysis } if (fullRank) { - double[] p = MiscMath.elementwiseMultiply(V.sumProduct(MiscMath.elementwiseDivide(uf, s)), -1); // inverted and roughly fine + double[] p = MiscMath.elementwiseMultiply(V.sumProduct(MiscMath.elementwiseDivide(uf, s)), -1); if (MiscMath.norm(p) <= Delta) { TrustRegion result = new TrustRegion(p, 0.0, 0); @@ -735,25 +680,10 @@ public class ccAnalysis private double evaluateQuadratic(MatrixI J, double[] g, double[] s) { - //TODO s (-> stepH) is slightly different - double[] Js = J.sumProduct(s); //TODO completely wromg + double[] Js = J.sumProduct(s); double q = MiscMath.dot(Js, Js); double l = MiscMath.dot(s, g); - /* - System.out.println("doing evaluateQuadratic"); - System.out.println("inputs"); - System.out.print("J"); - J.print(System.out, "%f "); - System.out.print("g"); - MiscMath.print(g, "%f"); - System.out.print("s"); - MiscMath.print(s, "%f"); - System.out.print("\nJs"); - MiscMath.print(Js, "%f"); - System.out.println(String.format("0.5 * %f + %f", q, l)); - */ - return 0.5 * q + l; } @@ -791,37 +721,7 @@ public class ccAnalysis } /** - * check the termination condition for nonlinear least squares - * TODO can be removed and added just as one line in trf (: terminationStatus = (ftolSatisfied condition || xtolSatisfied condition) ? 1 : 0;) doesnt matter as long as distinguished between 0 and rest - * - * @param actualReduction - * @param cost - * @param stepNorm - * @param xNorm - * @param ratio - * - * @return - */ - private byte checkTermination(double actualReduction, double cost, double stepNorm, double xNorm, double ratio) - { - // default ftol and xtol = 1e-8 - boolean ftolSatisfied = actualReduction < (1e-8 * cost) && ratio > 0.25; - boolean xtolSatisfied = stepNorm < (1e-8 * (1e-8 + xNorm)); - - if (ftolSatisfied && xtolSatisfied) - { - return (byte) 4; - } else if (ftolSatisfied) { - return (byte) 2; - } else if (xtolSatisfied) { - return (byte) 3; - } else { - return (byte) 0; - } - } - - /** - * TODO DOCUMENT ME! + * trust region reflective algorithm * @param repMatrix ~ Matrix containing representative vectors * @param scoresOld ~ Matrix containing initial observations * @param index ~ current row index @@ -831,75 +731,41 @@ public class ccAnalysis */ private double[] trf(MatrixI repMatrix, MatrixI scoresOld, int index, MatrixI J) { - System.out.println("-----------------\nStart of trf"); //hA = x0 - double[] hA = repMatrix.getRow(index); //inverted + double[] hA = repMatrix.getRow(index); double[] f0 = originalToEquasionSystem(hA, repMatrix, scoresOld, index); - int nfev = 1; // ?? - int njev = 1; // ?? + int nfev = 1; int m = J.height(); int n = J.width(); double cost = 0.5 * MiscMath.dot(f0, f0); - double[] g = J.transpose().sumProduct(f0); // inverted - double[] scale = new double[hA.length]; - Arrays.fill(scale, 1); // ?? + double[] g = J.transpose().sumProduct(f0); double Delta = MiscMath.norm(hA); - int maxNfev = hA.length * 100; // ?? - double alpha = 0.0; // ?? "Levenberg-Marquardt" parameter - - System.out.println("Checking initial values:"); - System.out.print("hA (x): "); - MiscMath.print(hA, "%1.8f"); - System.out.print("f0: "); - MiscMath.print(f0, "%1.8f"); - //System.out.println(String.format("nfev: %d, njev: %d, maxNfev: %d", nfev, njev, maxNfev)); - //System.out.println(String.format("m: %d, n: %d", m, n)); - System.out.println(String.format("cost: %1.8f, Delta: %1.8f, alpha: %1.8f", cost, Delta, alpha)); - System.out.print("g: "); - MiscMath.print(g, "%1.8f"); + int maxNfev = hA.length * 100; + double alpha = 0.0; // "Levenberg-Marquardt" parameter double gNorm = 0; byte terminationStatus = 0; int iteration = 0; - System.out.println("outer while loop starts"); while (true) { - System.out.println(String.format("iteration: %d", iteration)); - gNorm = MiscMath.norm(g); if (terminationStatus != 0 || nfev == maxNfev) { - System.out.println(String.format("outer loop broken with terminationStatus: %d and nfev: %d", terminationStatus, nfev)); break; } - // d = scale SingularValueDecomposition svd = new SingularValueDecomposition(new Array2DRowRealMatrix(J.asArray())); - // svd not 100% correct -> origin of problems - MatrixI U = new Matrix(svd.getU().getData()); // inverted - double[] s = svd.getSingularValues(); // ?? - MatrixI V = new Matrix(svd.getV().getData()).transpose(); //TODO inverted origin of probelm + MatrixI U = new Matrix(svd.getU().getData()); + double[] s = svd.getSingularValues(); + MatrixI V = new Matrix(svd.getV().getData()).transpose(); double[] uf = U.transpose().sumProduct(f0); - System.out.println("After SVD"); - System.out.println(String.format("gNorm: %1.8f", gNorm)); - System.out.print("U: "); - U.print(System.out, "%1.8f "); - System.out.print("s: "); - MiscMath.print(s, "%1.8f"); - System.out.print("V: "); - V.print(System.out, "%1.8f "); - System.out.print("uf: "); - MiscMath.print(uf, "%1.8f"); - double actualReduction = -1; double[] xNew = new double[hA.length]; double[] fNew = new double[f0.length]; double costNew = 0; double stepHnorm = 0; - System.out.println("Inner while loop starts"); - while (actualReduction <= 0 && nfev < maxNfev) { TrustRegion trustRegion = solveLsqTrustRegion(n, m, uf, s, V, Delta, alpha); @@ -914,19 +780,9 @@ public class ccAnalysis stepHnorm = MiscMath.norm(stepH); - System.out.println("After TrustRegion"); - System.out.print("stepH: "); - MiscMath.print(stepH, "%1.8f "); - System.out.println(String.format("alpha: %1.8f, nIterations: %d, predictedReduction: %1.8f, nfev: %d, stepHnorm: %1.8f", alpha, nIterations, predictedReduction, nfev, stepHnorm)); - System.out.print("xNew: "); - MiscMath.print(xNew, "%1.8f "); - System.out.print("fNew: "); - MiscMath.print(fNew, "%1.8f "); - if (MiscMath.countNaN(fNew) > 0) { Delta = 0.25 * stepHnorm; - System.out.println(String.format("Loop continued with %d NaNs and Delta: %1.8f", MiscMath.countNaN(fNew), Delta)); continue; } @@ -938,7 +794,10 @@ public class ccAnalysis double DeltaNew = updatedTrustRegion[0]; double ratio = updatedTrustRegion[1]; - terminationStatus = checkTermination(actualReduction, cost, stepHnorm, MiscMath.norm(hA), ratio); + // default ftol and xtol = 1e-8 + boolean ftolSatisfied = actualReduction < (1e-8 * cost) && ratio > 0.25; + boolean xtolSatisfied = stepHnorm < (1e-8 * (1e-8 + xNorm)); + terminationStatus = (ftolSatisfied || xtolSatisfied) ? (byte) 1 : (byte) 0; if (terminationStatus != 0) { break; @@ -947,10 +806,7 @@ public class ccAnalysis alpha *= Delta / DeltaNew; Delta = DeltaNew; - System.out.println(String.format("actualReduction: %1.8f, alpha: %1.8f, Delta: %1.8f, termination_status: %d, cost_new: %1.8f", actualReduction, alpha, Delta, terminationStatus, costNew)); - //break; } - System.out.println(String.format("actualReduction before check: %1.8f", actualReduction)); if (actualReduction > 0) { hA = xNew; @@ -958,9 +814,6 @@ public class ccAnalysis cost = costNew; J = approximateDerivative(repMatrix, scoresOld, index); - System.out.println("J in the end"); - J.print(System.out, "%1.8f "); - njev++; g = J.transpose().sumProduct(f0); } else { @@ -970,26 +823,13 @@ public class ccAnalysis iteration++; } - System.out.println("into OptimizeResult"); - System.out.println("x (hA)"); - MiscMath.print(hA, "%1.8f"); - System.out.println(String.format("cost: %1.8f", cost)); - System.out.println("f0"); - MiscMath.print(f0, "%1.8f"); - System.out.println("J"); - J.print(System.out, "%1.8f "); - System.out.println("g"); - MiscMath.print(g, "%1.8f"); - System.out.println(String.format("gNorm: %1.8f", gNorm)); - System.out.println(String.format("nfev: %d", nfev)); - System.out.println(String.format("njev: %d", njev)); - System.out.println(String.format("terminationStatus: %d", terminationStatus)); - // OptimizeResult(x, cost, f0, J, g, gNorm, 0 in shape x, nfev, njev, terminationStatus) return hA; } /** - * TODO DOCUMENT ME! + * performs the least squares optimisation + * adapted from https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html#scipy.optimize.least_squares + * * @param repMatrix ~ Matrix containing representative vectors * @param scoresOld ~ Matrix containing initial observations * @param index ~ current row index @@ -998,10 +838,7 @@ public class ccAnalysis */ private double[] leastSquaresOptimisation(MatrixI repMatrix, MatrixI scoresOld, int index) { - System.out.println("lsq starts!!!"); MatrixI J = approximateDerivative(repMatrix, scoresOld, index); - System.out.println("J"); - J.print(System.out, "%1.8f "); double[] result = trf(repMatrix, scoresOld, index, J); return result; } diff --git a/src/jalview/gui/CalculationChooser.java b/src/jalview/gui/CalculationChooser.java index 6a3aa0b..21d985c 100644 --- a/src/jalview/gui/CalculationChooser.java +++ b/src/jalview/gui/CalculationChooser.java @@ -89,13 +89,13 @@ public class CalculationChooser extends JPanel private static final int MIN_PCA_SELECTION = 4; - private static final int MIN_PASIMAP_SELECTION = 8; //&! + private static final int MIN_PASIMAP_SELECTION = 8; AlignFrame af; JRadioButton pca; - JRadioButton pasimap; //&! initialize JRadioButton object for pasimap + JRadioButton pasimap; JRadioButton neighbourJoining; @@ -124,7 +124,6 @@ public class CalculationChooser extends JPanel */ private PCAPanel pcaPanel; - //&! private PaSiMapPanel pasimapPanel; /** @@ -207,7 +206,7 @@ public class CalculationChooser extends JPanel pcaBorderless.add(pca, FlowLayout.LEFT); calcChoicePanel.add(pcaBorderless, FlowLayout.LEFT); - //&! create pasimap panel + // create pasimap panel JPanel pasimapBorderless = new JPanel(new FlowLayout(FlowLayout.LEFT)); // create new JPanel (button) for pasimap pasimapBorderless.setBorder( BorderFactory.createEmptyBorder(2, b.left, 2, b.right)); // set border (margin) for button (same as treePanel and pca) @@ -222,7 +221,7 @@ public class CalculationChooser extends JPanel ButtonGroup calcTypes = new ButtonGroup(); calcTypes.add(pca); - calcTypes.add(pasimap); //&! add pasimap to the calculation types + calcTypes.add(pasimap); calcTypes.add(neighbourJoining); calcTypes.add(averageDistance); diff --git a/src/jalview/gui/PaSiMapPanel.java b/src/jalview/gui/PaSiMapPanel.java index f3a1f52..cd1dd5a 100644 --- a/src/jalview/gui/PaSiMapPanel.java +++ b/src/jalview/gui/PaSiMapPanel.java @@ -91,7 +91,7 @@ public class PaSiMapPanel extends GPCAPanel public PaSiMapPanel(AlignmentPanel alignPanel, String modelName, SimilarityParamsI params) { - super(); + super(3); // dim = 3 this.av = alignPanel.av; this.ap = alignPanel; boolean nucleotide = av.getAlignment().isNucleotide(); @@ -109,8 +109,6 @@ public class PaSiMapPanel extends GPCAPanel boolean selected = av.getSelectionGroup() != null && av.getSelectionGroup().getSize() > 0; - //&! do i need seqstrings? -> no - //AlignmentView seqstrings = av.getAlignmentView(selected); SequenceI[] seqs; if (!selected) { @@ -124,8 +122,6 @@ public class PaSiMapPanel extends GPCAPanel ScoreModelI scoreModel = ScoreModels.getInstance() .getScoreModel(modelName, ap); setPasimapModel( - //&! - //new PaSiMapModel(seqstrings, seqs, nucleotide, scoreModel, params)); new PaSiMapModel(av, seqs, nucleotide, scoreModel, params)); PaintRefresher.Register(this, av.getSequenceSetId()); @@ -715,8 +711,6 @@ public class PaSiMapPanel extends GPCAPanel * * @param data */ - //public void setInputData(AlignmentView data) - //&! viewport or panel? public void setInputData(AlignmentViewport data) { getPasimapModel().setInputData(data); diff --git a/src/jalview/gui/PairwiseAlignPanel.java b/src/jalview/gui/PairwiseAlignPanel.java index e51024e..fafee94 100755 --- a/src/jalview/gui/PairwiseAlignPanel.java +++ b/src/jalview/gui/PairwiseAlignPanel.java @@ -56,7 +56,12 @@ public class PairwiseAlignPanel extends GPairwiseAlignPanel * * @param viewport * DOCUMENT ME! + * @param endGaps ~ toggle gaps and the beginning and end of sequences */ + public PairwiseAlignPanel(AlignmentViewport viewport) + { + this(viewport, false); + } public PairwiseAlignPanel(AlignmentViewport viewport, boolean endGaps) { super(); @@ -109,12 +114,10 @@ public class PairwiseAlignPanel extends GPairwiseAlignPanel as.calcScoreMatrix(); if (endGaps) { - //&! as.traceAlignmentWithEndGaps(); } else { as.traceAlignment(); } - //&! as.scoreAlignment(); if (!first) @@ -125,7 +128,6 @@ public class PairwiseAlignPanel extends GPairwiseAlignPanel first = false; as.printAlignment(System.out); scores[i][j] = as.getMaxScore() / as.getASeq1().length; - //&! alignmentScores[i][j] = as.getAlignmentScore(); totscore = totscore + scores[i][j]; diff --git a/src/jalview/jbgui/GPCAPanel.java b/src/jalview/jbgui/GPCAPanel.java index cf7bc4e..ec5b209 100755 --- a/src/jalview/jbgui/GPCAPanel.java +++ b/src/jalview/jbgui/GPCAPanel.java @@ -85,6 +85,23 @@ public class GPCAPanel extends JInternalFrame zCombobox.addItem("dim " + i); } } + public GPCAPanel(int dim) + { + try + { + jbInit(); + } catch (Exception e) + { + e.printStackTrace(); + } + + for (int i = 1; i <= dim; i++) + { + xCombobox.addItem("dim " + i); + yCombobox.addItem("dim " + i); + zCombobox.addItem("dim " + i); + } + } private void jbInit() throws Exception { diff --git a/src/jalview/viewmodel/PaSiMapModel.java b/src/jalview/viewmodel/PaSiMapModel.java index 1f5e2dc..a97277f 100644 --- a/src/jalview/viewmodel/PaSiMapModel.java +++ b/src/jalview/viewmodel/PaSiMapModel.java @@ -38,8 +38,6 @@ public class PaSiMapModel /* * inputs */ - //private AlignmentView inputData; - //&! private AlignmentViewport inputData; private final SequenceI[] seqs; @@ -72,8 +70,6 @@ public class PaSiMapModel * @param modelName * @param params */ - //public PaSiMapModel(AlignmentView seqData, SequenceI[] sqs, boolean nuc, - //&! viewport or panel? public PaSiMapModel(AlignmentViewport seqData, SequenceI[] sqs, boolean nuc, ScoreModelI modelName, SimilarityParamsI params) { @@ -103,14 +99,9 @@ public class PaSiMapModel int width = pasimap.getWidth(); int height = pasimap.getHeight(); - // top = pasimap.getM().height() - 1; - //&! - //top = height - 1; top = height; points = new Vector<>(); - //&! - //Point[] scores = pasimap.getComponents(top - 1, top - 2, top - 3, 100); Point[] scores = pasimap.getComponents(width - 1, width - 2, width - 3, 1); for (int i = 0; i < top; i++) @@ -163,8 +154,6 @@ public class PaSiMapModel public void updateRcView(int dim1, int dim2, int dim3) { // note: actual indices for components are dim1-1, etc (patch for JAL-1123) - //&! - //Point[] scores = pasimap.getComponents(dim1 - 1, dim2 - 1, dim3 - 1, 100); Point[] scores = pasimap.getComponents(dim1 - 1, dim2 - 1, dim3 - 1, 1); for (int i = 0; i < pasimap.getWidth(); i++) @@ -178,15 +167,11 @@ public class PaSiMapModel return pasimap.getDetails(); } - //&! - //public AlignmentView getInputData() public AlignmentViewport getInputData() { return inputData; } - //&! - //public void setInputData(AlignmentView data) public void setInputData(AlignmentViewport data) { inputData = data;