/*
 * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
 * Copyright (C) $$Year-Rel$$ The Jalview Authors
 * 
 * This file is part of Jalview.
 * 
 * Jalview 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.
 *  
 * Jalview 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.
 * 
 * You should have received a copy of the GNU General Public License
 * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
 * The Jalview Authors are detailed in the 'AUTHORS' file.
 */
package jalview.analysis;

import jalview.api.analysis.DistanceScoreModelI;
import jalview.api.analysis.ScoreModelI;
import jalview.api.analysis.SimilarityScoreModelI;
import jalview.datamodel.AlignmentView;
import jalview.math.MatrixI;

import java.io.PrintStream;

/**
 * Performs Principal Component Analysis on given sequences
 */
public class PCA implements Runnable
{
  boolean jvCalcMode = true;

  MatrixI symm;

  double[] eigenvalue;

  MatrixI eigenvector;

  StringBuilder details = new StringBuilder(1024);

  private AlignmentView seqs;

  private ScoreModelI scoreModel;

  public PCA(AlignmentView s, ScoreModelI sm)
  {
    this.seqs = s;

    scoreModel = sm;
    details.append("PCA calculation using " + sm.getName()
            + " sequence similarity matrix\n========\n\n");
  }

  /**
   * Returns Eigenvalue
   * 
   * @param i
   *          Index of diagonal within matrix
   * 
   * @return Returns value of diagonal from matrix
   */
  public double getEigenvalue(int i)
  {
    return eigenvector.getD()[i];
  }

  /**
   * DOCUMENT ME!
   * 
   * @param l
   *          DOCUMENT ME!
   * @param n
   *          DOCUMENT ME!
   * @param mm
   *          DOCUMENT ME!
   * @param factor
   *          DOCUMENT ME!
   * 
   * @return DOCUMENT ME!
   */
  public float[][] getComponents(int l, int n, int mm, float factor)
  {
    float[][] out = new float[getHeight()][3];

    for (int i = 0; i < getHeight(); i++)
    {
      out[i][0] = (float) component(i, l) * factor;
      out[i][1] = (float) component(i, n) * factor;
      out[i][2] = (float) component(i, mm) * factor;
    }

    return out;
  }

  /**
   * DOCUMENT ME!
   * 
   * @param n
   *          DOCUMENT ME!
   * 
   * @return DOCUMENT ME!
   */
  public double[] component(int n)
  {
    // n = index of eigenvector
    double[] out = new double[getHeight()];

    for (int i = 0; i < out.length; i++)
    {
      out[i] = component(i, n);
    }

    return out;
  }

  /**
   * DOCUMENT ME!
   * 
   * @param row
   *          DOCUMENT ME!
   * @param n
   *          DOCUMENT ME!
   * 
   * @return DOCUMENT ME!
   */
  double component(int row, int n)
  {
    double out = 0.0;

    for (int i = 0; i < symm.width(); i++)
    {
      out += (symm.getValue(row, i) * eigenvector.getValue(i, n));
    }

    return out / eigenvector.getD()[n];
  }

  public String getDetails()
  {
    return details.toString();
  }

  /**
   * DOCUMENT ME!
   */
  @Override
  public void run()
  {
    PrintStream ps = new PrintStream(System.out)
    {
      @Override
      public void print(String x)
      {
        details.append(x);
      }

      @Override
      public void println()
      {
        details.append("\n");
      }
    };

    // long now = System.currentTimeMillis();
    try
    {
      details.append("PCA Calculation Mode is "
              + (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
              + "\n");

      eigenvector = computeSimilarity(seqs);

      details.append(" --- OrigT * Orig ---- \n");
      eigenvector.print(ps, "%8.2f");

      symm = eigenvector.copy();

      eigenvector.tred();

      details.append(" ---Tridiag transform matrix ---\n");
      details.append(" --- D vector ---\n");
      eigenvector.printD(ps, "%15.4e");
      ps.println();
      details.append("--- E vector ---\n");
      eigenvector.printE(ps, "%15.4e");
      ps.println();

      // Now produce the diagonalization matrix
      eigenvector.tqli();
    } catch (Exception q)
    {
      q.printStackTrace();
      details.append("\n*** Unexpected exception when performing PCA ***\n"
              + q.getLocalizedMessage());
      details.append("*** Matrices below may not be fully diagonalised. ***\n");
    }

    details.append(" --- New diagonalization matrix ---\n");
    eigenvector.print(ps, "%8.2f");
    details.append(" --- Eigenvalues ---\n");
    eigenvector.printD(ps, "%15.4e");
    ps.println();
    /*
     * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
     * (int ev=0;ev<symm.rows; ev++) {
     * 
     * ps.print(","+component(seq, ev)); } ps.println(); }
     */
    // System.out.println(("PCA.run() took "
    // + (System.currentTimeMillis() - now) + "ms"));
  }

  /**
   * Computes a pairwise similarity matrix for the given sequence regions using
   * the configured score model. If the score model is a similarity model, then
   * it computes the result directly. If it is a distance model, then use it to
   * compute pairwise distances, and convert these to similarity scores by
   * substracting from the maximum value.
   * 
   * @param av
   * @return
   */
  MatrixI computeSimilarity(AlignmentView av)
  {
    MatrixI result = null;
    if (scoreModel instanceof SimilarityScoreModelI)
    {
      result = ((SimilarityScoreModelI) scoreModel).findSimilarities(av);
    }
    else if (scoreModel instanceof DistanceScoreModelI)
    {
      result = ((DistanceScoreModelI) scoreModel).findDistances(av);
      result.reverseRange(true);
    }
    else
    {
      System.err
              .println("Unexpected type of score model, cannot calculate similarity");
    }

    return result;
  }

  public void setJvCalcMode(boolean calcMode)
  {
    this.jvCalcMode = calcMode;
  }

  /**
   * Answers the N dimensions of the NxN PCA matrix. This is the number of
   * sequences involved in the pairwise score calculation.
   * 
   * @return
   */
  public int getHeight()
  {
    // TODO can any of seqs[] be null?
    return seqs.getSequences().length;
  }
}