Test version to switch branches
[jalview.git] / src / jalview / analysis / PaSiMap.java
index 77c28be..81d44f5 100755 (executable)
@@ -26,6 +26,7 @@ import jalview.bin.Console;
 import jalview.datamodel.AlignmentView;
 import jalview.datamodel.Point;
 import jalview.datamodel.SequenceI;
+import jalview.datamodel.SequenceGroup;
 import jalview.gui.PairwiseAlignPanel;
 import jalview.gui.PaSiMapPanel;
 import jalview.math.Matrix;
@@ -35,6 +36,8 @@ import jalview.viewmodel.AlignmentViewport;
 
 import java.io.PrintStream;
 import java.util.Hashtable;
+import java.util.HashMap;
+import java.util.HashSet;
 import java.util.Enumeration;
 
 /**
@@ -57,9 +60,9 @@ public class PaSiMap implements Runnable
   /*
    * outputs
    */
-  private MatrixI pairwiseScores;
+  private PairwiseAlignPanel alignment;
 
-  private MatrixI tridiagonal;
+  private MatrixI pairwiseScores;
 
   private MatrixI eigenMatrix;
 
@@ -71,8 +74,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 +96,7 @@ public class PaSiMap implements Runnable
   }
 
   /**
-   * DOCUMENT ME!
+   * Returns coordinates for each datapoint
    * 
    * @param l
    *          DOCUMENT ME!
@@ -103,8 +104,7 @@ public class PaSiMap implements Runnable
    *          DOCUMENT ME!
    * @param mm
    *          DOCUMENT ME!
-   * @param factor
-   *          DOCUMENT ME!
+   * @param factor ~ is 1
    * 
    * @return DOCUMENT ME!
    */
@@ -119,12 +119,6 @@ public class PaSiMap implements Runnable
       float z = (float) component(i, mm) * factor;
       out[i] = new Point(x, y, z);
     }
-    //&!
-    System.out.println("Points:");
-    for (Point point : out)
-    {
-      System.out.println(point.toString());
-    }
 
     return out;
   }
@@ -140,7 +134,7 @@ public class PaSiMap implements Runnable
   public double[] component(int n)
   {
     // n = index of eigenvector
-    double[] out = new double[getHeight()];
+    double[] out = new double[getWidth()];
 
     for (int i = 0; i < out.length; i++)
     {
@@ -162,19 +156,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);
   }
 
@@ -192,43 +173,30 @@ public class PaSiMap implements Runnable
     PrintStream ps = wrapOutputBuffer(sb);
 
     /*
-     * pairwise similarity scores
+     * coordinates matrix, with D vector
      */
-    sb.append(" --- OrigT * Orig ---- \n");
-    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");
+    sb.append(" --- Pairwise correlation coefficients ---\n");
+    pairwiseScores.print(ps, "%8.6f ");
     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");
-    eigenMatrix.print(ps, "%8.2f");
     sb.append(" --- Eigenvalues ---\n");
     eigenMatrix.printD(ps, "%15.4e");
     ps.println();
 
+    sb.append(" --- Coordinates ---\n");
+    eigenMatrix.print(ps, "%8.6f ");
+    ps.println();
+
     return sb.toString();
   }
 
   /**
    * 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,27 +207,52 @@ public class PaSiMap implements Runnable
   {
     try
     {
-      // run needleman regardless if aligned or not
-      // gui.PairwiseAlignPanel <++>
-      PairwiseAlignPanel alignment = new PairwiseAlignPanel(seqs, true);
+for (SequenceI see : seqs.getAlignment().getSequencesArray())
+{
+System.out.println(see.getName());
+}
+
+      int nSeqs = seqs.getAlignment().getHeight();
+      float[][] scores = new float[nSeqs][nSeqs];      // rows, cols
+
+      int nSplits = 1;
+      //while ((nSeqs / nSplits) > 300)                // heap full at 341
+      while (((float) nSeqs / nSplits) > 5f)           // heap full at 341
+       nSplits++;
+      int splitSeqs = (int) Math.ceil((float) nSeqs / nSplits);
+System.out.println(String.format("%d -> %d splits into %d seqs", nSeqs, nSplits, splitSeqs));
+
+      int[] splitIndices = new int[nSplits];
+      for (int i = 0; i < nSplits; i++)
+      {
+       splitIndices[i] = splitSeqs * (i + 1);  //exclusive!!
+      }
+
+      HashMap<int[], Float> valuesForScores = splitCombineAndAlign(seqs.getAlignment().getSequencesArray(), splitIndices, splitSeqs);
+
+      for (int[] coords : valuesForScores.keySet())
+      {
+       scores[coords[0]][coords[1]] = valuesForScores.get(coords);
+      }
+pairwiseScores = new Matrix(scores);
+pairwiseScores.print(System.out, "%1.4f ");
+
+/*
+      alignment = new PairwiseAlignPanel(seqs, true, 100, 5);
       float[][] scores = alignment.getAlignmentScores();       //bigger index first -- eg scores[14][13]
 
-      Hashtable<SequenceI, Integer> connectivity = seqs.calculateConnectivity(scores, dim);
+      //Hashtable<SequenceI, Integer> connectivity = seqs.calculateConnectivity(scores, dim);
 
       pairwiseScores = new Matrix(scores);
+pairwiseScores.print(System.out, "%1.4f ");
+/*
       pairwiseScores.fillDiagonal();
 
-      ccAnalysis cc = new ccAnalysis(pairwiseScores, dim);
-      pairwiseScores = cc.run();
-      tridiagonal = pairwiseScores.copy();
-      //tridiagonal.tred();
+      eigenMatrix = pairwiseScores.copy();
 
-      /** perform the eigendecomposition for the plot */
-      eigenMatrix = tridiagonal.copy();
-      eigenMatrix.setD(pairwiseScores.getD());
-      //eigenMatrix.tqli();
-      System.out.println(getDetails());
-      
+      ccAnalysis cc = new ccAnalysis(pairwiseScores, dim);
+      eigenMatrix = cc.run();
+*/
 
     } catch (Exception q)
     {
@@ -269,6 +262,163 @@ public class PaSiMap implements Runnable
   }
 
   /**
+  * aligns sequences in splits
+  * Splits each split into halves and aligns them agains halves of other splits
+  *
+  * @param seqs
+  * @param i ~ indices of split
+  * @param s ~ sequences per split
+  *
+  * @return  a map of where to put in scores, value ~ scores[n][m] = v
+  **/
+  protected HashMap<int[], Float> splitCombineAndAlign(SequenceI[] seqArray, int[] i, int s)
+  {
+    HashMap<int[], Float> result = new HashMap<int[], Float>();
+
+    int[][] allGroups = new int[i.length][s];
+    for (int g = 0; g < i.length; g++) // group g
+    {
+      int e = 0;       // index going through allGroups[][e]
+      for (int j = g * s; j < i[g]; j++)  // goes through all numbers in one group
+      {
+        allGroups[g][e++] = j >= seqArray.length ? -1 : j;
+      }
+    }
+    
+    int g = 0; // group count
+    for (int[] group : allGroups)
+    {
+      HashSet<SequenceI> sg = new HashSet<SequenceI>();
+      //SequenceGroup sg = new SequenceGroup();
+      for (int index : group)
+      {
+       if (index == -1)
+         continue;
+       //sg.addSequence(seqArray[index], false);
+       sg.add(seqArray[index]);
+      }
+      SequenceI[] sgArray = new SequenceI[sg.size()];
+      int k = 0;
+      for (SequenceI seq : sg)
+      {
+       sgArray[k++] = seq;
+      }
+      //seqs.setSelectionGroup(sg);
+      //PairwiseAlignPanel pap = new PairwiseAlignPanel(seqs, true, 100, 5);
+      //float[][] scores = pap.getAlignmentScores();   //bigger index first -- eg scores[14][13]
+      float[][] scores = simulateAlignment(sgArray);
+      for (int s1 = 0; s1 < scores.length; s1++)       // row
+      {
+       result.put(new int[]{s1 + g * s, s1 + g * s}, Float.NaN);       // self score = Float.NaN
+       for (int s2 = 0; s2 < s1; s2++)                 // col
+       {
+         result.put(new int[]{s1 + g * s, s2 + g * s}, scores[s1][s2]);
+       }
+      }   
+      g++;
+    }
+
+    int smallS = (int) Math.ceil((float) s/2);
+    int[][] newGroups = new int[i.length * 2][smallS];
+    
+    g = 0;
+    for (int[] group : allGroups)
+    {
+      int[] split1 = new int[smallS];
+      int[] split2 = new int[smallS];
+      for (int k = 0; k < group.length; k++)   
+      {
+       if (k < smallS)
+         split1[k] = group[k];
+       else
+         split2[k - smallS] = group[k];
+      }
+      newGroups[g++] = split1;
+      newGroups[g++] = split2;
+    }
+
+    // align each subsplit with subsplits from other split groups
+    for (int subsplitN = 0; subsplitN < newGroups.length; subsplitN++)
+    {
+      int c = 1;               // current split block
+      while (newGroups[subsplitN][0] > smallS * c)
+      {
+       c++;
+      }
+      for (int nextSplit = subsplitN + 1; nextSplit < newGroups.length; nextSplit++)
+      {
+       if (newGroups[nextSplit][0] >= s * c)   // if next subsplit of next split group -> align seqs
+       {
+          HashSet<SequenceI> sg = new HashSet<SequenceI>();
+          //SequenceGroup sg = new SequenceGroup();
+          for (int index : newGroups[subsplitN])
+          {
+           if (index == -1)
+             continue;
+           //sg.addSequence(seqArray[index], false);
+           sg.add(seqArray[index]);
+          }
+          for (int index : newGroups[nextSplit])
+          {
+           if (index == -1)
+             continue;
+           //sg.addSequence(seqArray[index], false);
+           sg.add(seqArray[index]);
+          }
+          SequenceI[] sgArray = new SequenceI[sg.size()];
+          int k = 0;
+          for (SequenceI seq : sg)
+          {
+           sgArray[k++] = seq;
+          }
+          //seqs.setSelectionGroup(sg);
+          //PairwiseAlignPanel pap = new PairwiseAlignPanel(seqs, true, 100, 5);
+          //float[][] scores = pap.getAlignmentScores();       //bigger index first -- eg scores[14][13]
+          float[][] scores = simulateAlignment(sgArray);
+          for (int s1 = 0; s1 < scores.length; s1++)   // row
+          {
+           for (int s2 = 0; s2 < s1; s2++)                     // col
+           {
+             if (s1 >= smallS && s2 < smallS)
+               result.put(new int[]{s1 + (nextSplit-1) * smallS, s2 + subsplitN * smallS}, scores[s1][s2]);
+           }
+          }   
+       }
+      }
+    }
+
+    return result;
+  }
+
+  /**
+  * simulate the alignment of a PairwiseAlignPanel
+  *
+  * @param seqs
+  * @return alignment scores
+  */
+  protected float[][] simulateAlignment(SequenceI[] seqs)
+  {
+    float[][] result = new float[seqs.length][seqs.length];
+    for (int i = 1; i < seqs.length; i++)
+    {
+      for (int j = 0; j < i; j++)
+      {
+       String[] seqStrings = new String[2];
+       seqStrings[0] = seqs[i].getSequenceAsString();
+       seqStrings[1] = seqs[j].getSequenceAsString();
+
+       AlignSeq as = new AlignSeq(seqs[i], seqStrings[0], seqs[j], seqStrings[1], AlignSeq.PEP);
+       as.calcScoreMatrix();
+       as.traceAlignmentWithEndGaps();
+       as.scoreAlignment();
+       as.printAlignment(System.out);
+       result[i][j] = as.getAlignmentScore();
+      }
+    }
+    return result;
+  }
+
+  /**
    * Returns a PrintStream that wraps (appends its output to) the given
    * StringBuilder
    * 
@@ -303,7 +453,7 @@ public class PaSiMap implements Runnable
   public int getHeight()
   {
     // TODO can any of seqs[] be null?
-    return pairwiseScores.height();// seqs.getSequences().length;
+    return eigenMatrix.height();// seqs.getSequences().length;
   }
 
   /**
@@ -315,7 +465,7 @@ public class PaSiMap implements Runnable
   public int getWidth()
   {
     // TODO can any of seqs[] be null?
-    return pairwiseScores.width();// seqs.getSequences().length;
+    return eigenMatrix.width();// seqs.getSequences().length;
   }
 
   /**
@@ -344,13 +494,18 @@ public class PaSiMap implements Runnable
     eigenMatrix = m;
   }
 
-  public MatrixI getTridiagonal()
+  public PairwiseAlignPanel getAlignments()
+  {
+    return alignment;
+  }
+
+  public String getAlignmentOutput()
   {
-    return tridiagonal;
+    return alignment.getAlignmentOutput();
   }
 
-  public void setTridiagonal(MatrixI tridiagonal)
+  public byte getDim()
   {
-    this.tridiagonal = tridiagonal;
+    return dim;
   }
 }