src/jalview/analysis/ccAnalysis.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21
  22 /*
  23 * Copyright 2018-2022 Kathy Su, Kay Diederichs
  24 *
  25 * This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
  26 *
  27 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
  28 *
  29 * You should have received a copy of the GNU General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
  30 */
  31
  32 /**
  33 * Ported from https://doi.org/10.1107/S2059798317000699 by
  34 * @AUTHOR MorellThomas
  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  * A class to model rectangular matrices of double values and operations on them
  57  */
  58 public class ccAnalysis
  59 {
  60   private byte dim = 0;         //dimensions
  61
  62   private MatrixI scoresOld;    //input scores
  63
  64   public ccAnalysis(MatrixI scores, byte dim)
  65   {
  66     // round matrix to .4f to be same as in pasimap
  67     for (int i = 0; i < scores.height(); i++)
  68     {
  69       for (int j = 0; j < scores.width(); j++)
  70       {
  71         if (!Double.isNaN(scores.getValue(i,j)))
  72         {
  73           scores.setValue(i, j, (double) Math.round(scores.getValue(i,j) * (int) 10000) / 10000);
  74         }
  75       }
  76     }
  77     this.scoresOld = scores;
  78     this.dim = dim;
  79   }
  80
  81   /**
  82   * Initialise a distrust-score for each hypothesis (h) of hSigns
  83   * distrust = conHypNum - proHypNum
  84   *
  85   * @param hSigns ~ hypothesis signs (+/-) for each sequence
  86   * @param scores ~ input score matrix
  87   *
  88   * @return distrustScores
  89   */
  90   private int[] initialiseDistrusts(byte[] hSigns, MatrixI scores)
  91   {
  92     int[] distrustScores = new int[scores.width()];
  93
  94     // loop over symmetric matrix
  95     for (int i = 0; i < scores.width(); i++)
  96     {
  97       byte hASign = hSigns[i];
  98       int conHypNum = 0;
  99       int proHypNum = 0;
 100
 101       for (int j = 0; j < scores.width(); j++)
 102       {
 103         double cell = scores.getRow(i)[j];      // value at [i][j] in scores
 104         byte hBSign = hSigns[j];
 105         if (!Double.isNaN(cell))
 106         {
 107           byte cellSign = (byte) Math.signum(cell);     //check if sign of matrix value fits hyptohesis
 108           if (cellSign == hASign * hBSign)
 109           {
 110             proHypNum++;
 111           } else {
 112             conHypNum++;
 113           }
 114         }
 115       }
 116       distrustScores[i] = conHypNum - proHypNum;        //create distrust score for each sequence
 117     }
 118     return distrustScores;
 119   }
 120
 121   /**
 122   * Optemise hypothesis concerning the sign of the hypothetical value for each hSigns by interpreting the pairwise correlation coefficients as scalar products
 123   *
 124   * @param hSigns ~ hypothesis signs (+/-)
 125   * @param distrustScores
 126   * @param scores ~ input score matrix
 127   *
 128   * @return hSigns
 129   */
 130   private byte[] optimiseHypothesis(byte[] hSigns, int[] distrustScores, MatrixI scores)
 131   {
 132     // get maximum distrust score
 133     int[] maxes = MiscMath.findMax(distrustScores);
 134     int maxDistrustIndex = maxes[0];
 135     int maxDistrust = maxes[1];
 136
 137     // if hypothesis is not optimal yet
 138     if (maxDistrust > 0)
 139     {
 140       //toggle sign for hI with maximum distrust
 141       hSigns[maxDistrustIndex] *= -1;
 142       // update distrust at same position
 143       distrustScores[maxDistrustIndex] *= -1;
 144
 145       // also update distrust scores for all hI that were not changed
 146       byte hASign = hSigns[maxDistrustIndex];
 147       for (int NOTmaxDistrustIndex = 0; NOTmaxDistrustIndex < distrustScores.length; NOTmaxDistrustIndex++)
 148       {
 149         if (NOTmaxDistrustIndex != maxDistrustIndex)
 150         {
 151           byte hBSign = hSigns[NOTmaxDistrustIndex];
 152           double cell = scores.getValue(maxDistrustIndex, NOTmaxDistrustIndex);
 153
 154           // distrust only changed if not NaN
 155           if (!Double.isNaN(cell))
 156           {
 157             byte cellSign = (byte) Math.signum(cell);
 158             // if sign of cell matches hypothesis decrease distrust by 2 because 1 more value supporting and 1 less contradicting
 159             // else increase by 2
 160             if (cellSign == hASign * hBSign)
 161             {
 162               distrustScores[NOTmaxDistrustIndex] -= 2;
 163             } else {
 164               distrustScores[NOTmaxDistrustIndex] += 2;
 165             }
 166           }
 167         }
 168       }
 169       //further optimisation necessary
 170       return optimiseHypothesis(hSigns, distrustScores, scores);
 171
 172     } else {
 173       return hSigns;
 174     }
 175   }
 176
 177   /**
 178   * takes the a symmetric MatrixI as input scores which may contain Double.NaN
 179   * approximate the missing values using hypothesis optimisation
 180   *
 181   * runs analysis
 182   *
 183   * @param scores ~ score matrix
 184   *
 185   * @return
 186   */
 187   public MatrixI run () throws Exception
 188   {
 189     //initialse eigenMatrix and repMatrix
 190     MatrixI eigenMatrix = scoresOld.copy();
 191     MatrixI repMatrix = scoresOld.copy();
 192     try
 193     {
 194     /*
 195     * Calculate correction factor for 2nd and higher eigenvalue(s).
 196     * This correction is NOT needed for the 1st eigenvalue, because the
 197     * unknown (=NaN) values of the matrix are approximated by presuming
 198     * 1-dimensional vectors as the basis of the matrix interpretation as dot
 199     * products.
 200     */
 201
 202     System.out.println("Input correlation matrix:");
 203     eigenMatrix.print(System.out, "%1.4f ");
 204
 205     int matrixWidth = eigenMatrix.width(); // square matrix, so width == height
 206     int matrixElementsTotal = (int) Math.pow(matrixWidth, 2);   //total number of elemts
 207
 208     float correctionFactor = (float) (matrixElementsTotal - eigenMatrix.countNaN()) / (float) matrixElementsTotal;
 209
 210     /*
 211     * Calculate hypothetical value (1-dimensional vector) h_i for each
 212     * dataset by interpreting the given correlation coefficients as scalar
 213     * products.
 214     */
 215
 216     /*
 217     * Memory for current hypothesis concerning sign of each h_i.
 218     * List of signs for all h_i in the encoding:
 219       * *  1: positive
 220       * *  0: zero
 221       * * -1: negative
 222     * Initial hypothesis: all signs are positive.
 223     */
 224     byte[] hSigns = new byte[matrixWidth];
 225     Arrays.fill(hSigns, (byte) 1);
 226
 227     //Estimate signs for each h_i by refining hypothesis on signs.
 228     hSigns = optimiseHypothesis(hSigns, initialiseDistrusts(hSigns, eigenMatrix), eigenMatrix);
 229
 230
 231     //Estimate absolute values for each h_i by determining sqrt of mean of
 232     //non-NaN absolute values for every row.
 233     double[] hAbs = MiscMath.sqrt(eigenMatrix.absolute().meanRow());
 234
 235     //Combine estimated signs with absolute values in obtain total value for
 236     //each h_i.
 237     double[] hValues = MiscMath.elementwiseMultiply(hSigns, hAbs);
 238
 239     /*Complement symmetric matrix by using the scalar products of estimated
 240     *values of h_i to replace NaN-cells.
 241     *Matrix positions that have estimated values
 242     *(only for diagonal and upper off-diagonal values, due to the symmetry
 243     *the positions of the lower-diagonal values can be inferred).
 244     List of tuples (row_idx, column_idx).*/
 245
 246     ArrayList<int[]> estimatedPositions = new ArrayList<int[]>();
 247
 248     // for off-diagonal cells
 249     for (int rowIndex = 0; rowIndex < matrixWidth - 1; rowIndex++)
 250     {
 251       for (int columnIndex = rowIndex + 1; columnIndex < matrixWidth; columnIndex++)
 252       {
 253         double cell = eigenMatrix.getValue(rowIndex, columnIndex);
 254         if (Double.isNaN(cell))
 255         {
 256           //calculate scalar product as new cell value
 257           cell = hValues[rowIndex] * hValues[columnIndex];
 258           //fill in new value in cell and symmetric partner
 259           eigenMatrix.setValue(rowIndex, columnIndex, cell);
 260           eigenMatrix.setValue(columnIndex, rowIndex, cell);
 261           //save positions of estimated values
 262           estimatedPositions.add(new int[]{rowIndex, columnIndex});
 263         }
 264       }
 265     }
 266
 267     // for diagonal cells
 268     for (int diagonalIndex = 0; diagonalIndex < matrixWidth; diagonalIndex++)
 269       {
 270         double cell = Math.pow(hValues[diagonalIndex], 2);
 271         eigenMatrix.setValue(diagonalIndex, diagonalIndex, cell);
 272         estimatedPositions.add(new int[]{diagonalIndex, diagonalIndex});
 273       }
 274
 275     /*Refine total values of each h_i:
 276     *Initialise h_values of the hypothetical non-existant previous iteration
 277     *with the correct format but with impossible values.
 278      Needed for exit condition of otherwise endless loop.*/
 279     System.out.print("initial values: [ ");
 280     for (double h : hValues)
 281     {
 282       System.out.print(String.format("%1.4f, ", h));
 283     }
 284     System.out.println(" ]");
 285
 286
 287     double[] hValuesOld = new double[matrixWidth];
 288
 289     int iterationCount = 0;
 290
 291     // repeat unitl values of h do not significantly change anymore
 292     while (true)
 293     {
 294       for (int hIndex = 0; hIndex < matrixWidth; hIndex++)
 295       {
 296         double newH = Arrays.stream(MiscMath.elementwiseMultiply(hValues, eigenMatrix.getRow(hIndex))).sum() / Arrays.stream(MiscMath.elementwiseMultiply(hValues, hValues)).sum();
 297         hValues[hIndex] = newH;
 298       }
 299
 300       System.out.print(String.format("iteration %d: [ ", iterationCount));
 301       for (double h : hValues)
 302       {
 303         System.out.print(String.format("%1.4f, ", h));
 304       }
 305       System.out.println(" ]");
 306
 307       //update values of estimated positions
 308       for (int[] pair : estimatedPositions)     // pair ~ row, col
 309       {
 310         double newVal = hValues[pair[0]] * hValues[pair[1]];
 311         eigenMatrix.setValue(pair[0], pair[1], newVal);
 312         eigenMatrix.setValue(pair[1], pair[0], newVal);
 313       }
 314
 315       iterationCount++;
 316
 317       //exit loop as soon as new values are similar to the last iteration
 318       if (MiscMath.allClose(hValues, hValuesOld, 0d, 1e-05d, false))
 319       {
 320         break;
 321       }
 322
 323       //save hValues for comparison in the next iteration
 324       System.arraycopy(hValues, 0, hValuesOld, 0, hValues.length);
 325     }
 326
 327     //-----------------------------
 328     //Use complemented symmetric matrix to calculate final representative
 329     //vectors.
 330
 331     //Eigendecomposition.
 332     eigenMatrix.tred();
 333     eigenMatrix.tqli();
 334
 335     System.out.println("eigenmatrix");
 336     eigenMatrix.print(System.out, "%8.2f");
 337     System.out.println();
 338     System.out.println("uncorrected eigenvalues");
 339     eigenMatrix.printD(System.out, "%2.4f ");
 340     System.out.println();
 341
 342     double[] eigenVals = eigenMatrix.getD();
 343
 344     TreeMap<Double, Integer> eigenPairs = new TreeMap<>(Comparator.reverseOrder());
 345     for (int i = 0; i < eigenVals.length; i++)
 346     {
 347       eigenPairs.put(eigenVals[i], i);
 348     }
 349
 350     // matrix of representative eigenvectors (each row is a vector)
 351     double[][] _repMatrix = new double[eigenVals.length][dim];
 352     double[][] _oldMatrix = new double[eigenVals.length][dim];
 353     double[] correctedEigenValues = new double[dim];
 354
 355     int l = 0;
 356     for (Entry<Double, Integer> pair : eigenPairs.entrySet())
 357     {
 358       double eigenValue = pair.getKey();
 359       int column = pair.getValue();
 360       double[] eigenVector = eigenMatrix.getColumn(column);
 361       //for 2nd and higher eigenvalues
 362       if (l >= 1)
 363       {
 364         eigenValue /= correctionFactor;
 365       }
 366       correctedEigenValues[l] = eigenValue;
 367       for (int j = 0; j < eigenVector.length; j++)
 368       {
 369         _repMatrix[j][l] = (eigenValue < 0) ? 0.0 : - Math.sqrt(eigenValue) * eigenVector[j];
 370         double tmpOldScore = scoresOld.getColumn(column)[j];
 371         _oldMatrix[j][dim - l - 1] = (Double.isNaN(tmpOldScore)) ? 0.0 : tmpOldScore;
 372       }
 373       l++;
 374       if (l >= dim)
 375       {
 376         break;
 377       }
 378     }
 379
 380     System.out.println("correctedEigenValues");
 381     MiscMath.print(correctedEigenValues, "%2.4f ");
 382
 383     repMatrix = new Matrix(_repMatrix);
 384     repMatrix.setD(correctedEigenValues);
 385     MatrixI oldMatrix = new Matrix(_oldMatrix);
 386
 387     MatrixI dotMatrix = repMatrix.postMultiply(repMatrix.transpose());
 388
 389     double rmsd = scoresOld.rmsd(dotMatrix);
 390
 391     System.out.println("iteration, rmsd, maxDiff, rmsdDiff");
 392     System.out.println(String.format("0, %8.5f, -, -", rmsd));
 393     // Refine representative vectors by minimising sum-of-squared deviates between dotMatrix and original  score matrix
 394     for (int iteration = 1; iteration < 21; iteration++)        // arbitrarily set to 20
 395     {
 396       MatrixI repMatrixOLD = repMatrix.copy();
 397       MatrixI dotMatrixOLD = dotMatrix.copy();
 398
 399       // for all rows/hA in the original matrix
 400       for (int hAIndex = 0; hAIndex < oldMatrix.height(); hAIndex++)
 401       {
 402         double[] row = oldMatrix.getRow(hAIndex);
 403         double[] hA = repMatrix.getRow(hAIndex);
 404         hAIndex = hAIndex;
 405         //find least-squares-solution fo rdifferences between original scores and representative vectors
 406         double[] hAlsm = leastSquaresOptimisation(repMatrix, scoresOld, hAIndex);
 407         // update repMatrix with new hAlsm
 408         for (int j = 0; j < repMatrix.width(); j++)
 409         {
 410           repMatrix.setValue(hAIndex, j, hAlsm[j]);
 411         }
 412       }
 413
 414       // dot product of representative vecotrs yields a matrix with values approximating the correlation matrix
 415       dotMatrix = repMatrix.postMultiply(repMatrix.transpose());
 416       // calculate rmsd between approximation and correlation matrix
 417       rmsd = scoresOld.rmsd(dotMatrix);
 418
 419       // calculate maximum change of representative vectors of current iteration
 420       MatrixI diff = repMatrix.subtract(repMatrixOLD).absolute();
 421       double maxDiff = 0.0;
 422       for (int i = 0; i < diff.height(); i++)
 423       {
 424         for (int j = 0; j < diff.width(); j++)
 425         {
 426           maxDiff = (diff.getValue(i, j) > maxDiff) ? diff.getValue(i, j) : maxDiff;
 427         }
 428       }
 429
 430       // calculate rmsd between current and previous estimation
 431       double rmsdDiff = dotMatrix.rmsd(dotMatrixOLD);
 432
 433       System.out.println(String.format("%d, %8.5f, %8.5f, %8.5f", iteration, rmsd, maxDiff, rmsdDiff));
 434
 435       if (!(Math.abs(maxDiff) > 1e-06))
 436       {
 437         repMatrix = repMatrixOLD.copy();
 438         break;
 439       }
 440     }
 441
 442
 443     } catch (Exception q)
 444     {
 445       Console.error("Error computing cc_analysis:  " + q.getMessage());
 446       q.printStackTrace();
 447     }
 448     System.out.println("final coordinates:");
 449     repMatrix.print(System.out, "%1.8f ");
 450     return repMatrix;
 451   }
 452
 453   /**
 454   * Create equations system using information on originally known
 455   * pairwise correlation coefficients (parsed from infile) and the
 456   * representative result vectors
 457   *
 458   * Each equation has the format:
 459   * hA * hA - pairwiseCC = 0
 460   * with:
 461   * hA: unknown variable
 462   * hB: known representative vector
 463   * pairwiseCC: known pairwise correlation coefficien
 464   *
 465   * The resulting equations system is overdetermined, if there are more
 466   * equations than unknown elements
 467   *
 468   * @param x ~ unknown n-dimensional column-vector
 469   * (needed for generating equations system, NOT to be specified by user).
 470   * @param hAIndex ~ index of currently optimised representative result vector.
 471   * @param h ~ matrix with row-wise listing of representative result vectors.
 472   * @param originalRow ~ matrix-row of originally parsed pairwise correlation coefficients.
 473   *
 474   * @return
 475   */
 476   private double[] originalToEquasionSystem(double[] hA, MatrixI repMatrix, MatrixI scoresOld, int hAIndex)
 477   {
 478     double[] originalRow = scoresOld.getRow(hAIndex);
 479     int nans = MiscMath.countNaN(originalRow);
 480     double[] result = new double[originalRow.length - nans];
 481
 482     //for all pairwiseCC in originalRow
 483     int resultIndex = 0;
 484     for (int hBIndex = 0; hBIndex < originalRow.length; hBIndex++)
 485     {
 486       double pairwiseCC = originalRow[hBIndex];
 487       // if not NaN -> create new equation and add it to the system
 488       if (!Double.isNaN(pairwiseCC))
 489       {
 490         double[] hB = repMatrix.getRow(hBIndex);
 491         result[resultIndex++] = MiscMath.sum(MiscMath.elementwiseMultiply(hA, hB)) - pairwiseCC;
 492       } else {
 493       }
 494     }
 495     return result;
 496   }
 497
 498   /**
 499   * returns the jacobian matrix
 500   * @param repMatrix ~ matrix of representative vectors
 501   * @param hAIndex ~ current row index
 502   *
 503   * @return
 504   */
 505   private MatrixI approximateDerivative(MatrixI repMatrix, MatrixI scoresOld, int hAIndex)
 506   {
 507     //hA = x0
 508     double[] hA = repMatrix.getRow(hAIndex);
 509     double[] f0 = originalToEquasionSystem(hA, repMatrix, scoresOld, hAIndex);
 510     double[] signX0 = new double[hA.length];
 511     double[] xAbs = new double[hA.length];
 512     for (int i = 0; i < hA.length; i++)
 513     {
 514       signX0[i] = (hA[i] >= 0) ? 1 : -1;
 515       xAbs[i] = (Math.abs(hA[i]) >= 1.0) ? Math.abs(hA[i]) : 1.0;
 516       }
 517     double rstep = Math.pow(Math.ulp(1.0), 0.5);
 518
 519     double[] h = new double [hA.length];
 520     for (int i = 0; i < hA.length; i++)
 521     {
 522       h[i] = rstep * signX0[i] * xAbs[i];
 523     }
 524
 525     int m = f0.length;
 526     int n = hA.length;
 527     double[][] jTransposed = new double[n][m];
 528     for (int i = 0; i < h.length; i++)
 529     {
 530       double[] x = new double[h.length];
 531       System.arraycopy(hA, 0, x, 0, h.length);
 532       x[i] += h[i];
 533       double dx = x[i] - hA[i];
 534       double[] df = originalToEquasionSystem(x, repMatrix, scoresOld, hAIndex);
 535       for (int j = 0; j < df.length; j++)
 536       {
 537         df[j] -= f0[j];
 538         jTransposed[i][j] = df[j] / dx;
 539       }
 540     }
 541     MatrixI J = new Matrix(jTransposed).transpose();
 542     return J;
 543   }
 544
 545   /**
 546   * norm of regularized (by alpha) least-squares solution minus Delta
 547   * @param alpha
 548   * @param suf
 549   * @param s
 550   * @param Delta
 551   *
 552   * @return
 553   */
 554   private double[] phiAndDerivative(double alpha, double[] suf, double[] s, double Delta)
 555   {
 556     double[] denom = MiscMath.elementwiseAdd(MiscMath.elementwiseMultiply(s, s), alpha);
 557     double pNorm = MiscMath.norm(MiscMath.elementwiseDivide(suf, denom));
 558     double phi = pNorm - Delta;
 559     // - sum ( suf**2 / denom**3) / pNorm
 560     double phiPrime = - MiscMath.sum(MiscMath.elementwiseDivide(MiscMath.elementwiseMultiply(suf, suf), MiscMath.elementwiseMultiply(MiscMath.elementwiseMultiply(denom, denom), denom))) / pNorm;
 561     return new double[]{phi, phiPrime};
 562   }
 563
 564   /**
 565   * class holding the result of solveLsqTrustRegion
 566   */
 567   private class TrustRegion
 568   {
 569     private double[] step;
 570     private double alpha;
 571     private int iteration;
 572
 573     public TrustRegion(double[] step, double alpha, int iteration)
 574     {
 575       this.step = step;
 576       this.alpha = alpha;
 577       this.iteration = iteration;
 578     }
 579
 580     public double[] getStep()
 581     {
 582       return this.step;
 583     }
 584
 585     public double getAlpha()
 586     {
 587       return this.alpha;
 588     }
 589
 590     public int getIteration()
 591     {
 592       return this.iteration;
 593     }
 594   }
 595
 596   /**
 597   * solve a trust-region problem arising in least-squares optimisation
 598   * @param n ~ number of variables
 599   * @param m ~ number of residuals
 600   * @param uf
 601   * @param s ~ singular values of J
 602   * @param V ~ transpose of VT
 603   * @param Delta ~ radius of a trust region
 604   * @param alpha ~ initial guess for alpha
 605   *
 606   * @return
 607   */
 608   private TrustRegion solveLsqTrustRegion(int n, int m, double[] uf, double[] s, MatrixI V, double Delta, double alpha)
 609   {
 610     double[] suf = MiscMath.elementwiseMultiply(s, uf);
 611
 612     //check if J has full rank and tr Gauss-Newton step
 613     boolean fullRank = false;
 614     if (m >= n)
 615     {
 616       double threshold = s[0] * Math.ulp(1.0) * m;
 617       fullRank = s[s.length - 1] > threshold;
 618     }
 619     if (fullRank)
 620     {
 621       double[] p = MiscMath.elementwiseMultiply(V.sumProduct(MiscMath.elementwiseDivide(uf, s)), -1);
 622       if (MiscMath.norm(p) <= Delta)
 623       {
 624         TrustRegion result = new TrustRegion(p, 0.0, 0);
 625         return result;
 626       }
 627     }
 628
 629     double alphaUpper = MiscMath.norm(suf) / Delta;
 630     double alphaLower = 0.0;
 631     if (fullRank)
 632     {
 633       double[] phiAndPrime = phiAndDerivative(0.0, suf, s, Delta);
 634       alphaLower = - phiAndPrime[0] / phiAndPrime[1];
 635     }
 636
 637     alpha = (!fullRank && alpha == 0.0) ? alpha = Math.max(0.001 * alphaUpper, Math.pow(alphaLower * alphaUpper, 0.5)) : alpha;
 638
 639     int iteration = 0;
 640     while (iteration < 10)      // 10 is default max_iter
 641     {
 642       alpha = (alpha < alphaLower || alpha > alphaUpper) ? alpha = Math.max(0.001 * alphaUpper, Math.pow(alphaLower * alphaUpper, 0.5)) : alpha;
 643       double[] phiAndPrime = phiAndDerivative(alpha, suf, s, Delta);
 644       double phi = phiAndPrime[0];
 645       double phiPrime = phiAndPrime[1];
 646
 647       alphaUpper = (phi < 0) ? alpha : alphaUpper;
 648       double ratio = phi / phiPrime;
 649       alphaLower = Math.max(alphaLower, alpha - ratio);
 650       alpha -= (phi + Delta) * ratio / Delta;
 651
 652       if (Math.abs(phi) < 0.01 * Delta) // default rtol set to 0.01
 653       {
 654         break;
 655       }
 656       iteration++;
 657     }
 658
 659     // p = - V.dot( suf / (s**2 + alpha))
 660     double[] tmp = MiscMath.elementwiseDivide(suf, MiscMath.elementwiseAdd(MiscMath.elementwiseMultiply(s, s), alpha));
 661     double[] p = MiscMath.elementwiseMultiply(V.sumProduct(tmp), -1);
 662
 663     // Make the norm of p equal to Delta, p is changed only slightly during this.
 664     // It is done to prevent p lie outside of the trust region
 665     p = MiscMath.elementwiseMultiply(p, Delta / MiscMath.norm(p));
 666
 667     TrustRegion result = new TrustRegion(p, alpha, iteration + 1);
 668     return result;
 669   }
 670
 671   /**
 672   * compute values of a quadratic function arising in least squares
 673   * function: 0.5 * s.T * (J.T * J + diag) * s + g.T * s
 674   *
 675   * @param J ~ jacobian matrix
 676   * @param g ~ gradient
 677   * @param s ~ steps and rows
 678   *
 679   * @return
 680   */
 681   private double evaluateQuadratic(MatrixI J, double[] g, double[] s)
 682   {
 683
 684     double[] Js = J.sumProduct(s);
 685     double q = MiscMath.dot(Js, Js);
 686     double l = MiscMath.dot(s, g);
 687
 688     return 0.5 * q + l;
 689   }
 690
 691   /**
 692   * update the radius of a trust region based on the cost reduction
 693   *
 694   * @param Delta
 695   * @param actualReduction
 696   * @param predictedReduction
 697   * @param stepNorm
 698   * @param boundHit
 699   *
 700   * @return
 701   */
 702   private double[] updateTrustRegionRadius(double Delta, double actualReduction, double predictedReduction, double stepNorm, boolean boundHit)
 703   {
 704     double ratio = 0;
 705     if (predictedReduction > 0)
 706     {
 707       ratio = actualReduction / predictedReduction;
 708     } else if (predictedReduction == 0 && actualReduction == 0) {
 709       ratio = 1;
 710     } else {
 711       ratio = 0;
 712     }
 713
 714     if (ratio < 0.25)
 715     {
 716       Delta = 0.25 * stepNorm;
 717     } else if (ratio > 0.75 && boundHit) {
 718       Delta *= 2.0;
 719     }
 720
 721     return new double[]{Delta, ratio};
 722   }
 723
 724   /**
 725   * trust region reflective algorithm
 726   * @param repMatrix ~ Matrix containing representative vectors
 727   * @param scoresOld ~ Matrix containing initial observations
 728   * @param index ~ current row index
 729   * @param J ~ jacobian matrix
 730   *
 731   * @return
 732   */
 733   private double[] trf(MatrixI repMatrix, MatrixI scoresOld, int index, MatrixI J)
 734   {
 735     //hA = x0
 736     double[] hA = repMatrix.getRow(index);
 737     double[] f0 = originalToEquasionSystem(hA, repMatrix, scoresOld, index);
 738     int nfev = 1;
 739     int m = J.height();
 740     int n = J.width();
 741     double cost = 0.5 * MiscMath.dot(f0, f0);
 742     double[] g = J.transpose().sumProduct(f0);
 743     double Delta = MiscMath.norm(hA);
 744     int maxNfev = hA.length * 100;
 745     double alpha = 0.0;         // "Levenberg-Marquardt" parameter
 746
 747     double gNorm = 0;
 748     boolean terminationStatus = false;
 749     int iteration = 0;
 750
 751     while (true)
 752     {
 753       gNorm = MiscMath.norm(g);
 754       if (terminationStatus || nfev == maxNfev)
 755       {
 756         break;
 757       }
 758       SingularValueDecomposition svd = new SingularValueDecomposition(new Array2DRowRealMatrix(J.asArray()));
 759       MatrixI U = new Matrix(svd.getU().getData());
 760       double[] s = svd.getSingularValues();
 761       MatrixI V = new Matrix(svd.getV().getData()).transpose();
 762       double[] uf = U.transpose().sumProduct(f0);
 763
 764       double actualReduction = -1;
 765       double[] xNew = new double[hA.length];
 766       double[] fNew = new double[f0.length];
 767       double costNew = 0;
 768       double stepHnorm = 0;
 769
 770       while (actualReduction <= 0 && nfev < maxNfev)
 771       {
 772         TrustRegion trustRegion = solveLsqTrustRegion(n, m, uf, s, V, Delta, alpha);
 773         double[] stepH = trustRegion.getStep();
 774         alpha = trustRegion.getAlpha();
 775         int nIterations = trustRegion.getIteration();
 776         double predictedReduction = - (evaluateQuadratic(J, g, stepH));
 777
 778         xNew = MiscMath.elementwiseAdd(hA, stepH);
 779         fNew = originalToEquasionSystem(xNew, repMatrix, scoresOld, index);
 780         nfev++;
 781
 782         stepHnorm = MiscMath.norm(stepH);
 783
 784         if (MiscMath.countNaN(fNew) > 0)
 785         {
 786           Delta = 0.25 * stepHnorm;
 787           continue;
 788         }
 789
 790         // usual trust-region step quality estimation
 791         costNew = 0.5 * MiscMath.dot(fNew, fNew);
 792         actualReduction = cost - costNew;
 793
 794         double[] updatedTrustRegion = updateTrustRegionRadius(Delta, actualReduction, predictedReduction, stepHnorm, stepHnorm > (0.95 * Delta));
 795         double DeltaNew = updatedTrustRegion[0];
 796         double ratio = updatedTrustRegion[1];
 797
 798         // default ftol and xtol = 1e-8
 799         boolean ftolSatisfied = actualReduction < (1e-8 * cost) && ratio > 0.25;
 800         boolean xtolSatisfied = stepHnorm < (1e-8 * (1e-8 + MiscMath.norm(hA)));
 801         terminationStatus = ftolSatisfied || xtolSatisfied;
 802         if (terminationStatus)
 803         {
 804           break;
 805         }
 806
 807         alpha *= Delta / DeltaNew;
 808         Delta = DeltaNew;
 809
 810       }
 811       if (actualReduction > 0)
 812       {
 813         hA = xNew;
 814         f0 = fNew;
 815         cost = costNew;
 816
 817         J = approximateDerivative(repMatrix, scoresOld, index);
 818
 819         g = J.transpose().sumProduct(f0);
 820       } else {
 821         stepHnorm = 0;
 822         actualReduction = 0;
 823       }
 824       iteration++;
 825     }
 826
 827     return hA;
 828   }
 829
 830   /**
 831   * performs the least squares optimisation
 832   * adapted from https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html#scipy.optimize.least_squares
 833   *
 834   * @param repMatrix ~ Matrix containing representative vectors
 835   * @param scoresOld ~ Matrix containing initial observations
 836   * @param index ~ current row index
 837   *
 838   * @return
 839   */
 840   private double[] leastSquaresOptimisation(MatrixI repMatrix, MatrixI scoresOld, int index)
 841   {
 842     MatrixI J = approximateDerivative(repMatrix, scoresOld, index);
 843     double[] result = trf(repMatrix, scoresOld, index, J);
 844     return result;
 845   }
 846
 847 }