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;
/**
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++)
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
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
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))
*
* @return
*/
- public MatrixI run ()
+ public MatrixI run () throws Exception
{
+ //initialse eigenMatrix and repMatrix
MatrixI eigenMatrix = scoresOld.copy();
MatrixI repMatrix = scoresOld.copy();
try
* 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
//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.
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);
{
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;
}
//-----------------------------
//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<Double, double[]> eigenPairs = new TreeMap<>(Comparator.reverseOrder());
- for (int i = 0; i < eigenVals.length; i++)
- {
- eigenPairs.put(eigenVals[i], eigenMatrix.getColumn(i));
- }
- */
TreeMap<Double, Integer> eigenPairs = new TreeMap<>(Comparator.reverseOrder());
for (int i = 0; i < eigenVals.length; i++)
{
}
// 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<Double, double[]> pair : eigenPairs.entrySet())
int l = 0;
for (Entry<Double, Integer> pair : eigenPairs.entrySet())
{
{
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;
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));
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
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++)
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);
repMatrix = repMatrixOLD.copy();
break;
}
- break;
}
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;
}
*
* @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);
//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++)
jTransposed[i][j] = df[j] / dx;
}
}
- MatrixI J = new Matrix(jTransposed).transpose(); // inverted
+ MatrixI J = new Matrix(jTransposed).transpose();
return J;
}
* 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
}
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);
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;
}
}
/**
- * 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
*/
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);
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;
}
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;
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;
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 {
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
*/
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;
}
git://source.jalview.org / jalview.git / commitdiff