/* Jalview - a java multiple alignment editor
 * Copyright (C) 1998  Michele Clamp
 *
 * 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 2
 * of the License, or (at your option) any later version.
 *
 * 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.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */
package jalview.analysis;


import java.util.*;
import jalview.gui.*;
import jalview.datamodel.*;


public class Conservation {
  Vector sequences;
  int    start;
  int    end;
  Vector seqNums; // vector of int vectors where first is sequence checksum
  int maxLength=0; //  used by quality calcs
  boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
  private void calcSeqNums() {
    for (int i=0; i<sequences.size(); i++) {
      calcSeqNum(i);
    }
  }
  private void calcSeqNum(int i) {
    String sq=null; // for dumb jbuilder not-inited exception warning
    int[] sqnum = null;
    if (i>-1 && i<sequences.size()) {
      sq = ((SequenceI)sequences.elementAt(i)).getSequence();
      if (seqNums.size()<=i) {
        seqNums.addElement(new int[sq.length()+1]);
      }
      if (sq.hashCode()!=((int[])seqNums.elementAt(i))[0])
       {
         int j, len;
         seqNumsChanged = true;
         sq = ((SequenceI) sequences.elementAt(i)).getSequence();
         len = sq.length();
         if (maxLength < len)
           maxLength = len;
         sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
         sqnum[0] = sq.hashCode();
         for (j = 1; j <= len; j++)
         {
           sqnum[j] = ( (Integer) jalview.schemes.ResidueProperties.aaHash.get(new
               String(sq.substring(j - 1, j)))).intValue(); // yuk
         }
         seqNums.setElementAt(sqnum, i);
       }
     } else {
       // JBPNote INFO level debug
       System.err.println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
     }
   }
   Vector total = new Vector();
   public Vector quality;
   public Double[] qualityRange = new Double[2];
   String consString = "";

   Sequence consSequence;
   Hashtable        propHash;
   int              threshold;

  String name = "";

  public Conservation(String name,Hashtable propHash, int threshold, Vector sequences, int start, int end) {
    this.name      = name;
    this.propHash  = propHash;
    this.threshold = threshold;
    this.sequences = sequences;
    this.start     = start;
    this.end       = end;
    seqNums = new Vector(sequences.size());
    calcSeqNums();
  }


  public void  calculate() {

    for (int i = start;i <= end; i++) {
      Hashtable resultHash  = null;
      Hashtable residueHash = null;

      resultHash  = new Hashtable();
      residueHash = new Hashtable();

      for (int j=0; j < sequences.size(); j++) {
        // JBPNote - have to make sure elements of the sequences vector
        //  are tested like this everywhere...
        if (sequences.elementAt(j) instanceof Sequence) {
          Sequence s = (Sequence)sequences.elementAt(j);

          if (s.getLength() > i) {
            String res = s.getSequence().substring(i,i+1);

            if (residueHash.containsKey(res)) {
              int count = ((Integer)residueHash.get(res)).intValue() ;
              count++;
              residueHash.put(res,new Integer(count));
            } else {
              residueHash.put(res,new Integer(1));
            }
          } else {
            if (residueHash.containsKey("-")) {
              int count = ((Integer)residueHash.get("-")).intValue() ;
              count++;
              residueHash.put("-",new Integer(count));
            } else {
              residueHash.put("-",new Integer(1));
            }
          }
        }
      }

      //What is the count threshold to count the residues in residueHash()
      int thresh = threshold*(sequences.size())/100;

      //loop over all the found residues
      Enumeration e = residueHash.keys();

      while (e.hasMoreElements()) {

        String res = (String)e.nextElement();
        if (((Integer)residueHash.get(res)).intValue() > thresh) {

          //Now loop over the properties
          Enumeration e2 = propHash.keys();

          while (e2.hasMoreElements()) {
            String    type = (String)e2.nextElement();
            Hashtable ht   = (Hashtable)propHash.get(type);

            //Have we ticked this before?
            if (! resultHash.containsKey(type)) {
              if (ht.containsKey(res)) {
                resultHash.put(type,ht.get(res));
              } else {
                resultHash.put(type,ht.get("-"));
              }
            } else if ( ((Integer)resultHash.get(type)).equals((Integer)ht.get(res)) == false) {
              resultHash.put(type,new Integer(-1));
            }
          }
        }
      }

      total.addElement(resultHash);
    }
  }

  public int countGaps(int j)
  {
    int count = 0;

    for (int i = 0; i < sequences.size();i++)
    {
      if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length())
      {  count++; continue;}

      char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);
      if (jalview.util.Comparison.isGap((c)))
        count++;

    }
    return count;
  }
  /***
   * countConsNGaps
   * returns gap count in int[0], and conserved residue count in int[1]
   */
  public int[] countConsNGaps(int j)
  {
    int count = 0;
    int cons=0;
    int nres = 0;
    int[] r = new int[2];
    char f='$';
    for (int i = 0; i < sequences.size();i++)
    {
      if( j >= ((Sequence)sequences.elementAt(i)).getSequence().length())
      {  count++;
      continue;}

      char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);
      if (jalview.util.Comparison.isGap((c)))
        count++;
      else {
        nres++;
        if (nres==1) {
          f = c;
          cons++;
        } else
          if (f == c) {
            cons++;
          }
      }
    }
    r[0] = (nres==cons) ? 1 : 0;
    r[1] = count;


    return r;
  }

  public  void  verdict(boolean consflag, float percentageGaps) {
    String consString = "";

    for (int i=start; i <= end; i++) {
      int[] gapcons = countConsNGaps(i);
      boolean cons = (gapcons[0]==1) ? true : false;
      int totGaps = gapcons[1];
      float pgaps = (float)totGaps*100/(float)sequences.size();
      //      System.out.println("percentage gaps = "+pgaps+"\n");
      if (percentageGaps > pgaps)
      {
        Hashtable resultHash = (Hashtable)total.elementAt(i-start);

        //Now find the verdict
        int         count = 0;
        Enumeration e3    = resultHash.keys();

        while (e3.hasMoreElements())
        {
          String type    = (String)e3.nextElement();
          Integer result = (Integer)resultHash.get(type);

          //Do we want to count +ve conservation or +ve and -ve cons.?

          if (consflag)
          {
            if (result.intValue() == 1)
              count++;
          }
          else
          {
            if (result.intValue() != -1)
              count++;
          }
        }

        if (count < 10)
          consString = consString + String.valueOf(count);// Conserved props!=Identity
        else
          consString = consString + ((gapcons[0]==1) ? "*" : "+");

      }
      else
      {
        consString = consString + "-";
      }
    }

    consSequence = new Sequence(name,consString,start,end);
  }

  public Sequence getConsSequence() {
    return consSequence;
  }

  // From Alignment.java in jalview118

  public void findQuality() {
    findQuality(0,maxLength-1);
  }

  int[][] cons2;

  private void percentIdentity2() {
    calcSeqNums(); // updates maxLength, too.
    if (cons2==null || seqNumsChanged) {
      cons2 = new int[maxLength][24];
      // Initialize the array
      for (int j=0;j<24;j++) {
        for (int i=0; i < maxLength;i++) {
          cons2[i][j] = 0;
        }
      }

      int sqnum[];
      int j = 0;
      while(j < sequences.size()) {
        sqnum=(int[])seqNums.elementAt(j);
        for (int i = 1; i < sqnum.length; i++) {
          cons2[i-1][sqnum[i]]++;
        }
        for (int i=sqnum.length-1; i<maxLength; i++) {
          cons2[i][23]++; // gap count
        }
        j++;
      }

      // unnecessary ?
      /* for (int i=start; i <= end; i++) {
           int max = -1000;
    int maxi = -1;
    int maxj = -1;

    for (int j=0;j<24;j++) {
        if (cons2[i][j] > max) {
        max = cons2[i][j];
        maxi = i;
        maxj = j;
      }

    }
  } */
    }

}


public void findQuality(int start, int end) {
    quality = new Vector();
    double max = -10000;
    String s = "";
    int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
    //Loop over columns // JBPNote Profiling info
    //    long ts = System.currentTimeMillis();
    //long te = System.currentTimeMillis();
    percentIdentity2();

    int size = seqNums.size();
    int[] lengths = new int[size];

    for (int l = 0; l < size; l++)
      lengths[l] = ((int[]) seqNums.elementAt(l)).length-1;

    for (int j=start; j <= end; j++) {
      double bigtot = 0;

      // First Xr = depends on column only
      double x[] = new double[24];

      for (int ii=0; ii < 24; ii++) {
	x[ii] = 0;
	try {
	  for (int i2=0; i2 < 24; i2++) {
	    x[ii]  += (double)cons2[j][i2] * BLOSUM62[ii][i2]+4;
	  }
	} catch (Exception e) {
	  System.err.println("Exception during quality calculation.");
          e.printStackTrace();
	}
	//System.out.println("X " + ii + " " + x[ii]);
	x[ii] /= (size);
	//System.out.println("X " + ii + " " + x[ii]);
      }
      // Now calculate D for each position and sum
      for (int k=0; k < size; k++) {
	double tot = 0;
	double[] xx = new double[24];
        int seqNum =
            (j<lengths[k])
            ? ((int[]) seqNums.elementAt(k))[j+1]
            : 23; // Sequence, or gap at the end

        // This is a loop over r
        for (int i=0; i < 23; i++) {
          double sr = 0;
          try {
            sr = (double)BLOSUM62[i][seqNum]+4;
          } catch (Exception e) {
            System.out.println("Exception in sr: " + e);
            e.printStackTrace();
          }
          //Calculate X with another loop over residues

          //  System.out.println("Xi " + i + " " + x[i] + " " + sr);
          xx[i] = x[i] - sr;

          tot += xx[i]*xx[i];
        }
        bigtot += Math.sqrt(tot);
      }
      // This is the quality for one column
      if (max < bigtot) {max = bigtot;}
      //      bigtot  = bigtot * (size-cons2[j][23])/size;

      quality.addElement(new Double(bigtot));

      s += "-";

      // Need to normalize by gaps
    }
    double newmax=-10000;
    for (int j=start; j <= end; j++) {
     double tmp =  ((Double)quality.elementAt(j)).doubleValue();
     tmp = (max - tmp)*(size-cons2[j][23])/size;
     //     System.out.println(tmp+ " " + j);
     quality.setElementAt(new Double(tmp),j);
     if (tmp>newmax)
       newmax = tmp;
    }
    //    System.out.println("Quality " + s);
    qualityRange[0] = new Double(0);
    qualityRange[1] = new Double(newmax);
  }


}