package bullet;

/**
 * Title:        Exploring Genetic Bullet
 * Description:
 * Copyright:    Copyright (c) 2002
 * Company:      IHES
 * @author Andrey Zinovyev
 * @version 1.0
 */

import vdaoengine.*;
import vdaoengine.data.*;
import vdaoengine.utils.*;
import vdaoengine.analysis.cluster.*;
import bioseqtools.*;
import org.biojava.bio.seq.*;
import org.biojava.bio.*;
import org.biojava.bio.symbol.*;
import org.biojava.bio.seq.io.*;
import bioseqtools.Utils;
import java.io.*;
import java.util.*;
import java.text.*;

public class GenBullet2 {

static  String datFile, genFile;
static  int taskType = 0, parameter = 3, windowSize = 120; // 0 - annotation
static  int sbeg = -1,send = -1;
static  int sAnnotbeg = -1,sAnnotend = -1;
static  VDataSet VD;
static Sequence seq;
static SequenceIterator seqI;
static VNormalization N;
static int optimalWindow;
static float multiplier = (float)1.2;
static int everyPosition = 21;
static int all=0,good1=0,good2=0,good3=0,alljunk1=0,alljunk2=0,alljunk3=0,goodjunk1=0,goodjunk2=0,goodjunk3=0,
goodexon1=0,goodexon2=0,goodexon3=0;
static int numberOfClasses = 7;
static int taskModifier = 0;
static float factor = 1.0f;
static float noiseLevel = 0f;


  public GenBullet2() {
  }
  public static void main(String[] args) {


  for(int i=0;i<args.length;i++)
      if (i<args.length)
        {
         System.out.println(args[i]);
         if (args[i].equals("-d")) datFile = new String(args[i+1]);
         if (args[i].equals("-g")) genFile = new String(args[i+1]);
         if (args[i].equals("-t")) taskType = Integer.parseInt(args[i+1]);
         if (args[i].equals("-p")) parameter = Integer.parseInt(args[i+1]);
         if (args[i].equals("-b")) sbeg = Integer.parseInt(args[i+1]);
         if (args[i].equals("-f")) send = Integer.parseInt(args[i+1]);
         if (args[i].equals("-ba")) sAnnotbeg = Integer.parseInt(args[i+1]);
         if (args[i].equals("-fa")) sAnnotend = Integer.parseInt(args[i+1]);
         if (args[i].equals("-e")) everyPosition = Integer.parseInt(args[i+1]);
         if (args[i].equals("-w")) windowSize = Integer.parseInt(args[i+1]);
         if (args[i].equals("-m")) multiplier = Float.parseFloat(args[i+1]);
         if (args[i].equals("-c")) numberOfClasses = Integer.parseInt(args[i+1]);
         if (args[i].equals("-k")) taskModifier = Integer.parseInt(args[i+1]);
         if (args[i].equals("-vf")) factor = Float.parseFloat(args[i+1]);
         if (args[i].equals("-nl")) noiseLevel = Float.parseFloat(args[i+1]);
        }

  for(int i=0;i<args.length;i++)
    System.out.print(args[i]+" ");
    System.out.println();

  try{
  switch (taskType) {
  case 0: LoadSeq(); ProduceAnnotation(true); break;
  case 1: CalculateOptimalWindow(); break;
  case 2: LoadSeq(); CalculateOptimalWindowUsingAnnot(); break;
  case 3: LoadSeq(); CalculateCodonUsage(); break;
  case 4: LoadSeq(); ProduceAnnotationWithKnownCodonUsage(); break;
  case 5: LoadSeq(); ProduceAnnotationMultiply(true); break;
  case 6: LoadSeq(); ProduceAnnotationBayesian(true); break;
  case 7: LoadSeq(); ExtractCodonUsage(); break;
  }
  }catch(Exception e){
    e.printStackTrace();
  }

  }


public static void LoadSeq(){
  System.out.print("Loading sequence... "+genFile);
  try{
  File GenBankFile = new File(genFile);
  BufferedReader eReader = new BufferedReader(
  new InputStreamReader(new FileInputStream(GenBankFile)));
  seqI = SeqIOTools.readGenbank(eReader);
  seq = seqI.nextSequence();
  System.out.println(" Seq("+seq.getURN()+"): "+seq.length()+" bp, Feat: "+seq.countFeatures());
  } catch (Exception e) {System.out.println(e);}
}

public static void ProduceAnnotation(boolean doOut){
try{

  if(sbeg==-1) sbeg=windowSize/2+2;
  if(send==-1) send=seq.length()-windowSize/2-2;
  if(sAnnotbeg==-1) sAnnotbeg=windowSize/2+2;
  if(sAnnotend==-1) sAnnotend=seq.length()-windowSize/2-2;


/*  VDataTable VT = new VDataTable();
  VT.LoadFromVDatFile(datFile);
  VD = new VDataSet();
  VTableSelector VS = new VTableSelector();
  VS.selectFirstNumericalColumns(VT,64);
  System.out.print("Loading dataset...");
  VD.loadFromDataTable(VT,VS);
  System.out.println(" Rec: "+VD.getPointCount());
  N = new VNormalization(VD.calcStatistics());
  VD.addToPreProcess(N);
  VD.preProcessData();*/

  if(taskModifier==1)
    VVectorCalc.distanceType=1;

  VD = new VDataSet();
  int n = (send-sbeg)/everyPosition;
  VD.massif = new float[n][64];
  VD.setPointCount(n);
  VD.setCoordCount(64);

  if (doOut) System.out.println("Calculating... "+n+" points");

  int veci[];

  String Sseq = seq.seqString();
  if(noiseLevel>0){
    Sseq = addNoise(Sseq,noiseLevel/3f,noiseLevel/3f,noiseLevel/3f);
  }

  int k=0;
  for(int i=sbeg;i<send;i+=everyPosition)
  {
  //System.out.print(i+" ");
  veci = bioseqtools.Utils.calcTripletFreq(Sseq,i-windowSize/2,i+windowSize/2);
  if(k<n)
    for(int j=0;j<64;j++) VD.massif[k][j]=veci[j];
  k++;
  }

  VStatistics stat = VD.calcStatistics();
//  System.out.println(stat);
  N = new VNormalization(stat);
  VD.clearPreProcess();
  VD.addToPreProcess(N);
  VD.preProcessData();

VKMeansClusterization KM = new VKMeansClusterization(numberOfClasses);
if (doOut) System.out.println("Clusterization... "+numberOfClasses+" classes");
KM.setClusterizationQuality(0.00001);
KM.doClusterization(VD);
if (doOut) System.out.println("IN: "+KM.iterationsNumber+" Disp: "+KM.calculateAverageIntraClassStdDev(VD));
KM.rearrangeClusterNumbers();

//int veci[];
float vecf[] = new float[64];
int cn;

all=0;good1=0;good2=0;good3=0;alljunk1=0;alljunk2=0;alljunk3=0;goodjunk1=0;goodjunk2=0;goodjunk3=0;
goodexon1=0;goodexon2=0;goodexon3=0;

if (doOut) System.out.println("Annotation from "+sAnnotbeg+" to "+sAnnotbeg+" every "+parameter+"....");
for(int i=sAnnotbeg;i<sAnnotend;i+=parameter)
if(i<seq.length()-windowSize/2-2)
  {
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  for(int j=0;j<vecf.length;j++) vecf[j]=veci[j];
  N.process(vecf);
  cn = KM.calculateClosestCentroid(vecf);
//  cn = (int)((numberOfClasses+1)*Math.random());
  String fn = new String(bioseqtools.Utils.GetFeatureName(seq,i));
  int cd = bioseqtools.Utils.CalcIsCoding(seq,i);
  int cd1 = 0;//bioseqtools.Utils.CalcIsCodingGene(seq,i);
  int cd2 = 0;//bioseqtools.Utils.CalcIsCoding2(seq,i);
  int compl = bioseqtools.Utils.CalcIsGeneComplement(seq,i);
  if (doOut) System.out.println(i+"\t"+cn+"\t"+fn+"\t"+cd+"\t"+cd1+"\t"+cd2+"\t"+(float)compl*0.5);
  all++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd==1))) good1++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd1==1))) good2++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd2==1))) good3++;
  if(cd==0) alljunk1++;
  if(cd1==0) alljunk2++;
  if(cd2==0) alljunk3++;
  if((cd==0)&&(cn==0)) goodjunk1++;
  if((cd1==0)&&(cn==0)) goodjunk2++;
  if((cd2==0)&&(cn==0)) goodjunk3++;
  if((cd==1)&&(cn>0)) goodexon1++;
  if((cd1==1)&&(cn>0)) goodexon2++;
  if((cd2==1)&&(cn>0)) goodexon3++;
  }

if (doOut) System.out.println("% of coding: "+(float)(all-alljunk1)/all);
if (doOut) System.out.println("Quality:");
if (doOut) System.out.println(" CDS only: "+(float)good1/all+" in junk: "+(float)goodjunk1/alljunk1+"(from "+alljunk1+") in exons(Sn): "+(float)(goodexon1)/(all-alljunk1)+" (from "+(all-alljunk1)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon1/(all-alljunk1)+" Sp: "+(float)goodexon1/(goodexon1+alljunk1-goodjunk1));
if (doOut) System.out.println(" Whole genes: "+(float)good2/all+" in junk: "+(float)goodjunk2/alljunk2+"(from "+alljunk2+") in genes: "+(float)(goodexon2)/(all-alljunk2)+" (from "+(all-alljunk2)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon2/(all-alljunk2)+" Sp: "+(float)goodexon2/(goodexon2+alljunk2-goodjunk2));
if (doOut) System.out.println(" CDS+tRNA+rRNA: "+(float)good3/all+" in junk: "+(float)goodjunk3/alljunk3+"(from "+alljunk3+") in exons: "+(float)(goodexon3)/(all-alljunk3)+" (from "+(all-alljunk3)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon3/(all-alljunk3)+" Sp: "+(float)goodexon3/(goodexon3+alljunk3-goodjunk3));

} catch (Exception e) {System.out.println(e);}
}

public static void CalculateOptimalWindow(){

  try
  {
  System.out.print("Loading sequence...");
  try{
  File GenBankFile = new File(genFile);
  BufferedReader eReader = new BufferedReader(
  new InputStreamReader(new FileInputStream(GenBankFile)));
  SequenceIterator seqI = SeqIOTools.readGenbank(eReader);
  seq = seqI.nextSequence();
  System.out.println(" Seq("+seq.getURN()+"): "+seq.length()+" bp, Feat: "+seq.countFeatures());
  } catch (Exception e) {System.out.println(e);}
  if(sbeg==-1) sbeg=windowSize/2+2;
  if(send==-1) send=seq.length()-windowSize/2-2;

  VD = new VDataSet();
  int n = (send-sbeg)/parameter;
  VD.massif = new float[n][64];
  VD.setPointCount(n);
  VD.setCoordCount(64);

  System.out.println("Calculating... "+n+" points");

  windowSize = 21;
  int veci[];

  while(windowSize<seq.length()/3)
  {
  System.out.print(windowSize+"\t");

  int k=0;
  sbeg=windowSize/2+2;
  send=seq.length()-windowSize/2-2;
  for(int i=sbeg;i<send;i+=parameter)
  {
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  if(k<n)
    for(int j=0;j<64;j++) VD.massif[k][j]=veci[j];
  k++;
  }

  VStatistics stat = VD.calcStatistics();
//  System.out.println(stat);
  N = new VNormalization(stat);
  VD.clearPreProcess();
  VD.addToPreProcess(N);
  VD.preProcessData();

  VKMeansClusterization KM = new VKMeansClusterization(numberOfClasses);
  KM.setClusterizationQuality(0.00001);
  KM.doClusterization(VD);

  System.out.println(KM.calculateAverageIntraClassStdDev(VD));
  windowSize*=multiplier;
  }

  } catch(Exception e){System.out.println(e);};
}

public static void CalculateOptimalWindowUsingAnnot(){
windowSize = 50;
System.out.println("WS\tOverall\tSn\tSp\tCorrectjunk");
while(windowSize<1500){
sbeg=-1;
send=-1;
sAnnotbeg=-1;
sAnnotend=-1;
ProduceAnnotation(false);
System.out.println(windowSize+"\t"+(float)good1/all+"\t"+(float)(goodexon1)/(all-alljunk1)+"\t"+(float)goodexon1/(goodexon1+alljunk1-goodjunk1)+"\t"+(float)(goodjunk1)/(alljunk1));
windowSize+=50;
}
}


//-------------------------------------- CODON USAGE ------------
public static void  CalculateCodonUsage(){

int n = 0,k;

int fr01[] = new int[64];
int fr02[] = new int[64];
int fr11[] = new int[64];
int fr12[] = new int[64];
int fr21[] = new int[64];
int fr22[] = new int[64];

int cfr01[] = new int[64];
int cfr02[] = new int[64];
int cfr11[] = new int[64];
int cfr12[] = new int[64];
int cfr21[] = new int[64];
int cfr22[] = new int[64];


float ffr0[] = new float[64];
float ffr1[] = new float[64];
float ffr2[] = new float[64];

float ffr02[] = new float[64];
float ffr12[] = new float[64];
float ffr22[] = new float[64];

float cfr0[] = new float[64];
float cfr1[] = new float[64];
float cfr2[] = new float[64];

float cfr0_2[] = new float[64];
float cfr1_2[] = new float[64];
float cfr2_2[] = new float[64];


Vector jfr0 = new Vector();
Vector jfr1 = new Vector();
Vector jfr2 = new Vector();

int jfr01[] = new int[64];
int jfr02[] = new int[64];
int jfr11[] = new int[64];
int jfr12[] = new int[64];
int jfr21[] = new int[64];
int jfr22[] = new int[64];

float jffr0[] = new float[64];
float jffr1[] = new float[64];
float jffr2[] = new float[64];
float jffr02[] = new float[64];
float jffr12[] = new float[64];
float jffr22[] = new float[64];

jfr0.addElement(jfr01);
jfr0.addElement(jfr02);
jfr1.addElement(jfr11);
jfr1.addElement(jfr12);
jfr2.addElement(jfr21);
jfr2.addElement(jfr22);



Vector fr0 = new Vector();
Vector fr1 = new Vector();
Vector fr2 = new Vector();

Vector comfr0 = new Vector();
Vector comfr1 = new Vector();
Vector comfr2 = new Vector();

fr0.addElement(fr01);
fr0.addElement(fr02);
fr1.addElement(fr11);
fr1.addElement(fr12);
fr2.addElement(fr21);
fr2.addElement(fr22);

comfr0.addElement(cfr01);
comfr0.addElement(cfr02);
comfr1.addElement(cfr11);
comfr1.addElement(cfr12);
comfr2.addElement(cfr21);
comfr2.addElement(cfr22);


int l=0;
StringBuffer ss = new StringBuffer("");

do{
    for(Iterator i=seq.features(); i.hasNext(); )
      {
      Feature Ft = (Feature) i.next();
      if ((Ft.getType().equals("CDS"))/*&&(Ft.getLocation().isContiguous())*/)
         {
         try{
         k=CalcCodonUsageForSeq(Ft.getSymbols(),fr0,fr1,fr2);
//         ComplementaryCalcCodonUsageForSeq(Ft.getSymbols(),comfr0,comfr1,comfr2);
//         System.out.print("["+Ft.getLocation().getMin()+".."+Ft.getLocation().getMax()+"],");
//         if(Utils.CalcIsGeneComplement(seq,Ft.getLocation().getMin()+2)==-1)
//            {
//            System.out.println();
            //System.out.println(Utils.GetFeatureName(seq,Ft.getLocation().getMin()+2));
            //System.out.println(Ft.getSymbols().seqString());
//            }
         //System.out.print(l+" ");
         n+=k;
         } catch (Exception e)
         {System.out.print("\n"+e+" : ");
          System.out.println(Ft);}
         }
      l++;
      }
ss.append(Utils.getJunkString(seq));
try{
seq = seqI.nextSequence();
} catch(Exception e) {System.out.println(e);}
}while(seqI.hasNext());

// Junk ------------------
// String ss = Utils.getJunkString(seq);
// ss = Utils.FormatSeq(new StringBuffer(ss));
// System.out.println("\n"+ss);
int nj = CalcCodonUsageForSeq(ss.toString(),jfr0,jfr1,jfr2);

System.out.println();
System.out.println("Junk Phase 0:");
CalcCodonUsageForVec(jfr0,nj,jffr0,jffr02);
System.out.println("Junk Phase 1:");
CalcCodonUsageForVec(jfr1,nj,jffr1,jffr12);
System.out.println("Junk Phase 2:");
CalcCodonUsageForVec(jfr2,nj,jffr2,jffr22);

//------------------------

//if(taskModifier!=1)
{
System.out.println();
System.out.println("Phase 0:");
CalcCodonUsageForVec(fr0,n,ffr0,ffr02);
System.out.println("Phase 1:");
CalcCodonUsageForVec(fr1,n,ffr1,ffr12);
System.out.println("Phase 2:");
CalcCodonUsageForVec(fr2,n,ffr2,ffr22);

/*System.out.println();
System.out.println("Complementary Phase 0:");
CalcCodonUsageForVec(comfr0,n,cfr0,cfr0_2);
System.out.println("Complementary Phase 1:");
CalcCodonUsageForVec(comfr1,n,cfr1,cfr1_2);
System.out.println("Complementary Phase 2:");
CalcCodonUsageForVec(comfr2,n,cfr2,cfr2_2);*/
}

/*for(int i=0;i<64;i++){
ffr0[i]=(float)fr01[i]/(windowSize/3);
ffr1[i]=(float)fr11[i]/(windowSize/3);
ffr2[i]=(float)fr21[i]/(windowSize/3);
}*/

float ff0r1[] = Utils.ReconstructPhaseDistribution(ffr0,1);
float ff0r2[] = Utils.ReconstructPhaseDistribution(ffr0,2);

//if(taskModifier!=1)
{
System.out.println();
System.out.println("Distribution phase 0");
printDistribution(ffr0);

System.out.println();
System.out.println("Reconstructed phase 1");
printDistribution(ff0r1);

System.out.println();
System.out.println("Reconstructed phase 2");
printDistribution(ff0r2);

System.out.println();
System.out.println("Reconstruction Quality (phase 1): "+Utils.DistributionDistance(ffr1,ff0r1)/Utils.DistributionNorm(ffr0));
System.out.println("Reconstruction Quality (phase 2): "+Utils.DistributionDistance(ffr2,ff0r2)/Utils.DistributionNorm(ffr0));

System.out.println();
System.out.println("Norm 0: "+Utils.DistributionNorm(ffr0));
System.out.println("Norm 1: "+Utils.DistributionNorm(ff0r1));
System.out.println("Norm 2: "+Utils.DistributionNorm(ff0r2));
System.out.println("01: "+Utils.DistributionDistance(ffr0,ff0r1)+" ("+Utils.DistributionDistance(ffr0,ffr1)+")");
System.out.println("02: "+Utils.DistributionDistance(ffr0,ff0r2)+" ("+Utils.DistributionDistance(ffr0,ffr2)+")");
System.out.println("12: "+Utils.DistributionDistance(ff0r1,ff0r2)+" ("+Utils.DistributionDistance(ffr1,ffr2)+")");
System.out.println("Goodness: "+Utils.DistributionDistance(ffr0,ff0r1)/Utils.DistributionNorm(ffr0));
System.out.println("Dispersion ratio for exons: "+Utils.DistributionNorm(ffr02)/Utils.DistributionDistance(ffr0,ffr1));
System.out.println("Dispersion ratio for junk: "+Utils.DistributionNorm(jffr02)/Utils.DistributionDistance(ffr0,ffr1));
System.out.println("\nFor ViDaExpert:");
}

System.out.print("CodonUsage sphere ");
for(int i=0;i<64;i++)
   System.out.print(ffr0[i]*factor+" ");
System.out.println("0.4 255 0 0 0");

System.out.print("Reconstrph1 sphere ");
for(int i=0;i<64;i++)
if(taskModifier==2)
   System.out.print(ff0r1[i]*factor+" ");
else
   System.out.print(ffr1[i]*factor+" ");
System.out.println("0.3 255 100 0 0");

System.out.print("Reconstrph2 sphere ");
for(int i=0;i<64;i++)
if(taskModifier==2)
   System.out.print(ffr2[i]*factor+" ");
else
   System.out.print(ff0r2[i]*factor+" ");
System.out.println("0.3 255 100 0 0");

if (taskModifier>=1){
cfr0 = Utils.ComplementaryDistribution(ffr0);
if(taskModifier==2){
cfr1 = Utils.ComplementaryDistribution(ffr2);
cfr2 = Utils.ComplementaryDistribution(ffr1);}
else{
cfr1 = Utils.ReconstructPhaseDistribution(cfr0,1);
cfr2 = Utils.ReconstructPhaseDistribution(cfr0,2);}


System.out.println("Assymetry: "+Utils.ComplementaryAsymmetry(ffr0,ff0r1,ff0r2,cfr0,cfr1,cfr2));

System.out.print("ComplCodonUsage sphere ");
for(int i=0;i<64;i++)
   System.out.print(cfr0[i]*factor+" ");
System.out.println("0.4 0 255 0 1");

System.out.print("ComplReconstrph1 sphere ");
for(int i=0;i<64;i++)
   System.out.print(cfr1[i]*factor+" ");
System.out.println("0.3 100 255 100 1");

System.out.print("ComplReconstrph2 sphere ");
for(int i=0;i<64;i++)
   System.out.print(cfr2[i]*factor+" ");
System.out.println("0.3 100 255 100 1");

//float junkf[] = Utils.CalcJunkDistribution(ffr0,ff0r1,ff0r2,cfr0,cfr1,cfr2);
System.out.print("Junk0 sphere ");
for(int i=0;i<64;i++)
if(taskModifier==2)
   System.out.print(jffr0[i]*factor+" ");
else
   System.out.print(jffr0[i]*factor+" ");
System.out.println("0.25 0 0 255");
System.out.print("Junk1 sphere ");
for(int i=0;i<64;i++)
   System.out.print(jffr1[i]*factor+" ");
System.out.println("0.25 0 0 255");
System.out.print("Junk2 sphere ");
for(int i=0;i<64;i++)
   System.out.print(jffr2[i]*factor+" ");
System.out.println("0.25 0 0 255");



System.out.println("Conn01 cylinder connect CodonUsage Reconstrph1 0.1 0.1 255 0 0");
System.out.println("Conn02 cylinder connect CodonUsage Reconstrph2 0.1 0.1 255 0 0");
System.out.println("Conn12 cylinder connect Reconstrph1 Reconstrph2 0.08 0.08 255 100 0");

System.out.println("CConn01 cylinder connect ComplCodonUsage ComplReconstrph1 0.1 0.1 0 255 0");
System.out.println("CConn02 cylinder connect ComplCodonUsage ComplReconstrph2 0.1 0.1 0 255 0");
System.out.println("CConn12 cylinder connect ComplReconstrph1 ComplReconstrph2 0.08 0.08 100 255 100");

}

}

public static void CalcCodonUsageForVec(Vector f,int Number,float fr[],float fr2[]){
int ifr[] = (int[])f.elementAt(0);
int ifr2[] = (int[])f.elementAt(1);;
int m[];

for(int i=0;i<64;i++)
  {
  fr[i] = (float)ifr[i]/Number;
  fr2[i]=(float)Math.sqrt((float)ifr2[i]/Number-fr[i]*fr[i]);
  fr[i]/=windowSize/3;
  fr2[i]/=windowSize/3;
  }

for(int j=0;j<8;j++){
for(int i=0;i<8;i++)
  {
  System.out.print(Utils.TripletName(j*8+i)+" : "+fr[j*8+i]/**(windowSize/3)*/+"("+fr2[j*8+i]+") ");
  } System.out.println();}
}

public static int CalcCodonUsageForSeq(String seq, Vector fr0, Vector fr1, Vector fr2){
int freq[], freq1[], freq2[];
int n=0;
try{
if(seq.length()<windowSize){
freq = Utils.calcTripletFreq(seq,1,seq.length());
freq1 = (int[])fr0.elementAt(0);
freq2 = (int[])fr0.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}

freq = Utils.calcTripletFreq(seq,2,seq.length());
freq1 = (int[])fr1.elementAt(0);
freq2 = (int[])fr1.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}

freq = Utils.calcTripletFreq(seq,3,seq.length());
freq1 = (int[])fr2.elementAt(0);
freq2 = (int[])fr2.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}
n=1;
}
else{

for(int i=1;i<=seq.length()-windowSize-10;i+=everyPosition){
freq = Utils.calcTripletFreq(seq,i,i+windowSize);
freq1 = (int[])fr0.elementAt(0);
freq2 = (int[])fr0.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}

freq = Utils.calcTripletFreq(seq,i+1,i+1+windowSize);
freq1 = (int[])fr1.elementAt(0);
freq2 = (int[])fr1.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}

freq = Utils.calcTripletFreq(seq,i+2,i+2+windowSize);
freq1 = (int[])fr2.elementAt(0);
freq2 = (int[])fr2.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}
n++;
}

}}catch(Exception e){
System.out.println(e);
}
return n;
}

public static int CalcCodonUsageForSeq(SymbolList seq, Vector fr0, Vector fr1, Vector fr2){
int freq[], freq1[], freq2[];
int n=0;
try{
if(seq.length()<windowSize){
freq = Utils.calcTripletFreq(seq,1,seq.length());
freq1 = (int[])fr0.elementAt(0);
freq2 = (int[])fr0.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}

freq = Utils.calcTripletFreq(seq,2,seq.length());
freq1 = (int[])fr1.elementAt(0);
freq2 = (int[])fr1.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}

freq = Utils.calcTripletFreq(seq,3,seq.length());
freq1 = (int[])fr2.elementAt(0);
freq2 = (int[])fr2.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}
n=1;
}
else{

for(int i=1;i<=seq.length()-windowSize-10;i+=everyPosition){
freq = Utils.calcTripletFreq(seq,i,i+windowSize);
freq1 = (int[])fr0.elementAt(0);
freq2 = (int[])fr0.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}

freq = Utils.calcTripletFreq(seq,i+1,i+1+windowSize);
freq1 = (int[])fr1.elementAt(0);
freq2 = (int[])fr1.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}

freq = Utils.calcTripletFreq(seq,i+2,i+2+windowSize);
freq1 = (int[])fr2.elementAt(0);
freq2 = (int[])fr2.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}
n++;
}

}}catch(Exception e){
System.out.println(e);
}
return n;
}


public static int ComplementaryCalcCodonUsageForSeq(SymbolList seq, Vector fr0, Vector fr1, Vector fr2){
int freq[], freq1[], freq2[];
int n=0;
try{
if(seq.length()<windowSize){
freq = Utils.calcTripletFreqComplementary(seq,1,seq.length());
freq1 = (int[])fr0.elementAt(0);
freq2 = (int[])fr0.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}

freq = Utils.calcTripletFreqComplementary(seq,2,seq.length());
freq1 = (int[])fr1.elementAt(0);
freq2 = (int[])fr1.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}

freq = Utils.calcTripletFreqComplementary(seq,3,seq.length());
freq1 = (int[])fr2.elementAt(0);
freq2 = (int[])fr2.elementAt(1);
for(int i=0;i<64;i++){
  freq1[i]+=freq[i]; freq2[i]+=freq[i]*freq[i];}
n=1;
}
else{

for(int i=1;i<=seq.length()-windowSize-10;i+=everyPosition){
freq = Utils.calcTripletFreqComplementary(seq,i,i+windowSize);
freq1 = (int[])fr0.elementAt(0);
freq2 = (int[])fr0.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}

freq = Utils.calcTripletFreqComplementary(seq,i+1,i+1+windowSize);
freq1 = (int[])fr1.elementAt(0);
freq2 = (int[])fr1.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}

freq = Utils.calcTripletFreqComplementary(seq,i+2,i+2+windowSize);
freq1 = (int[])fr2.elementAt(0);
freq2 = (int[])fr2.elementAt(1);
for(int j=0;j<64;j++){
  freq1[j]+=freq[j]; freq2[j]+=freq[j]*freq[j];}
n++;
}

}}catch(Exception e){
System.out.println(e);
}
return n;
}
//-------------------------------------- CODON USAGE ------------
//----------- ANNOTATION WITH KNOWN CODON USAGE -------------------
public static void ProduceAnnotationWithKnownCodonUsage(){
try{

  if(sbeg==-1) sbeg=windowSize/2+2;
  if(send==-1) send=seq.length()-windowSize/2-2;

int fr01[] = new int[64];
int fr02[] = new int[64];
int fr11[] = new int[64];
int fr12[] = new int[64];
int fr21[] = new int[64];
int fr22[] = new int[64];

float ffr0[] = new float[64];
float ffr1[] = new float[64];
float ffr2[] = new float[64];

float ffr02[] = new float[64];
float ffr12[] = new float[64];
float ffr22[] = new float[64];

Vector fr0 = new Vector();
Vector fr1 = new Vector();
Vector fr2 = new Vector();

Vector jfr0 = new Vector();
Vector jfr1 = new Vector();
Vector jfr2 = new Vector();

int jfr01[] = new int[64];
int jfr02[] = new int[64];
int jfr11[] = new int[64];
int jfr12[] = new int[64];
int jfr21[] = new int[64];
int jfr22[] = new int[64];

float jffr0[] = new float[64];
float jffr1[] = new float[64];
float jffr2[] = new float[64];
float jffr02[] = new float[64];
float jffr12[] = new float[64];
float jffr22[] = new float[64];

jfr0.addElement(jfr01);
jfr0.addElement(jfr02);
jfr1.addElement(jfr11);
jfr1.addElement(jfr12);
jfr2.addElement(jfr21);
jfr2.addElement(jfr22);


fr0.addElement(fr01);
fr0.addElement(fr02);
fr1.addElement(fr11);
fr1.addElement(fr12);
fr2.addElement(fr21);
fr2.addElement(fr22);

float cfr0[] = new float[64];
float cfr1[] = new float[64];
float cfr2[] = new float[64];

int l=0,k=0,n=0;

System.out.println("Calculating Codon Usage...");

    for(Iterator i=seq.features(); i.hasNext(); )
      {
      Feature Ft = (Feature) i.next();
      if ((Ft.getType().equals("CDS"))&&(Ft.getLocation().isContiguous()))
         {
         try{
         k=CalcCodonUsageForSeq(Ft.getSymbols(),fr0,fr1,fr2);
//         ComplementaryCalcCodonUsageForSeq(Ft.getSymbols(),comfr0,comfr1,comfr2);
//         System.out.print("["+Ft.getLocation().getMin()+".."+Ft.getLocation().getMax()+"],");
//         if(Utils.CalcIsGeneComplement(seq,Ft.getLocation().getMin()+2)==-1)
//            {
//            System.out.println();
            //System.out.println(Utils.GetFeatureName(seq,Ft.getLocation().getMin()+2));
            //System.out.println(Ft.getSymbols().seqString());
//            }
         System.out.print(l+" ");
         n+=k;
         } catch (Exception e)
         {System.out.print("\n"+e+" : ");
          System.out.println(Ft);}
         }
      l++;
      }

System.out.println("Phase 0:");
CalcCodonUsageForVec(fr0,n,ffr0,ffr02);
System.out.println("Phase 1:");
CalcCodonUsageForVec(fr1,n,ffr1,ffr12);
System.out.println("Phase 2:");
CalcCodonUsageForVec(fr2,n,ffr2,ffr22);

String ss = Utils.getJunkString(seq);
int nj = CalcCodonUsageForSeq(ss,jfr0,jfr1,jfr2);

System.out.println();
System.out.println("Junk Phase 0:");
CalcCodonUsageForVec(jfr0,nj,jffr0,jffr02);
System.out.println("Junk Phase 1:");
CalcCodonUsageForVec(jfr1,nj,jffr1,jffr12);
System.out.println("Junk Phase 2:");
CalcCodonUsageForVec(jfr2,nj,jffr2,jffr22);


System.out.println("Reconstructing phases...");

System.out.println("Phase 1...");
float ff0r1[] = Utils.ReconstructPhaseDistribution(ffr0,1);
System.out.println("Phase 2...");
float ff0r2[] = Utils.ReconstructPhaseDistribution(ffr0,2);

cfr0 = Utils.ComplementaryDistribution(ffr0);
System.out.println("Complementary Phase 1...");
cfr1 = Utils.ReconstructPhaseDistribution(cfr0,1);
System.out.println("Complementary Phase 2...");
cfr2 = Utils.ReconstructPhaseDistribution(cfr0,2);

System.out.println("Creating dataset...");

  VD = new VDataSet();
  n = (send-sbeg)/everyPosition;
  VD.massif = new float[n][64];
  VD.setPointCount(n);
  VD.setCoordCount(64);

  System.out.println("Calculating... "+n+" points");

  int veci[];

  k=0;
  for(int i=sbeg;i<send;i+=everyPosition)
  {
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  if(k<n)
    for(int j=0;j<64;j++) VD.massif[k][j]=(float)veci[j]/(windowSize/3);
  k++;
  }

  VStatistics stat = VD.calcStatistics();
//  System.out.println(stat);
  N = new VNormalization(stat);
  VD.clearPreProcess();
  VD.addToPreProcess(N);
  VD.preProcessData();

VKMeansClusterization KM = new VKMeansClusterization(numberOfClasses);

KM.setClusterizationQuality(0.00001);
KM.doClusterization(VD);
System.out.println("IN: "+KM.iterationsNumber+" Disp: "+KM.calculateAverageIntraClassStdDev(VD));
KM.rearrangeClusterNumbers();

float cfr0r1[] = Utils.ComplementaryDistribution(ffr2);
float cfr0r2[] = Utils.ComplementaryDistribution(ffr1);

//float junkf[] = Utils.CalcJunkDistribution(ffr0,ffr1,ffr2,cfr0,cfr0r1,cfr0r2);
//float junkf[] = Utils.CalcJunkDistribution(ffr0,ff0r1,ff0r2,cfr0,cfr1,cfr2);
//float junkf[] = new float[64];

System.out.println("Setting centroids...");
N.process(jffr0);
N.process(ffr0);
N.process(ffr1);
N.process(ffr2);
N.process(ff0r1);
N.process(ff0r2);
N.process(cfr0);
N.process(cfr1);
N.process(cfr2);
N.process(cfr0r1);
N.process(cfr0r2);


System.out.println("Direct Triangle: "+VVectorCalc.Distance(ffr0,ffr1)+" "+VVectorCalc.Distance(ffr0,ffr2)+" "+VVectorCalc.Distance(ffr2,ffr1));

int cc = KM.calculateClosestCentroid(jffr0);
double rjunkf[] = KM.getCentroidVector(cc);
System.out.println("Closest to junk: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),jffr0));
cc = KM.calculateClosestCentroid(ffr0);
double rffr0[] = KM.getCentroidVector(cc);
System.out.println("Closest to ffr0: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),ffr0));
cc = KM.calculateClosestCentroid(ff0r1);
System.out.println("Closest to ff0r1: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),ff0r1));
cc = KM.calculateClosestCentroid(ff0r2);
System.out.println("Closest to ff0r2: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),ff0r2));

cc = KM.calculateClosestCentroid(ffr1);
System.out.println("Closest to ffr1: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),ffr1));
cc = KM.calculateClosestCentroid(ffr2);
System.out.println("Closest to ffr2: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),ffr2));

cc = KM.calculateClosestCentroid(cfr0);
System.out.println("Closest to cfr0: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),cfr0));
cc = KM.calculateClosestCentroid(cfr1);
System.out.println("Closest to cfr1: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),cfr1));
cc = KM.calculateClosestCentroid(cfr2);
System.out.println("Closest to cfr2: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),cfr2));

cc = KM.calculateClosestCentroid(cfr0r1);
System.out.println("Closest to cfr0r1: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),cfr0r1));
cc = KM.calculateClosestCentroid(cfr0r2);
System.out.println("Closest to cfr0r2: "+cc+" "+VVectorCalc.Distance(KM.getCentroidVector(cc),cfr0r2));


KM.clearCentroids();

KM.addCentroid(jffr0);
KM.addCentroid(ffr0);
//KM.addCentroid(ff0r1);
//KM.addCentroid(ff0r2);
KM.addCentroid(ffr1);
KM.addCentroid(ffr2);

if(numberOfClasses==7){
KM.addCentroid(cfr0);
//KM.addCentroid(cfr1);
//KM.addCentroid(cfr2);
KM.addCentroid(cfr0r1);
KM.addCentroid(cfr0r2);}


float vecf[] = new float[64];
int cn,cn1,cn2,cnf;

all=0;good1=0;good2=0;good3=0;alljunk1=0;alljunk2=0;alljunk3=0;goodjunk1=0;goodjunk2=0;goodjunk3=0;
goodexon1=0;goodexon2=0;goodexon3=0;

System.out.println("Annotation from "+sbeg+" to "+send+" every "+parameter+"....");
for(int i=sbeg;i<send;i+=parameter)
if(i<seq.length()-windowSize/2-2)
  {
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  for(int j=0;j<vecf.length;j++) vecf[j]=(float)veci[j]/(windowSize/3);
  N.process(vecf);
  cn = KM.calculateClosestCentroid(vecf);

  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2+1,i+windowSize/2);
  for(int j=0;j<vecf.length;j++) vecf[j]=(float)veci[j]/(windowSize/3);
  N.process(vecf);
  cn1 = KM.calculateClosestCentroid(vecf);

  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2+2,i+windowSize/2);
  for(int j=0;j<vecf.length;j++) vecf[j]=(float)veci[j]/(windowSize/3);
  N.process(vecf);
  cn2 = KM.calculateClosestCentroid(vecf);

//  if((cn1==0)||(cn2==0)) cn=0;
  cnf=0;
//  if (((cn==1)&&(cn1==2)&&(cn2==3))||((cn==2)&&(cn1==3)&&(cn2==1))||((cn==3)&&(cn1==1)&&(cn2==2))) cnf=1;
//  if (((cn==4)&&(cn1==5)&&(cn2==6))||((cn==5)&&(cn1==6)&&(cn2==4))||((cn==6)&&(cn1==4)&&(cn2==5))) cnf=1;
//  if((cn>0)||(cn1>0)||(cn2>0)) cnf=1;
  cnf=cn;

//  cn = (int)((numberOfClasses+1)*Math.random());
  String fn = new String(bioseqtools.Utils.GetFeatureName(seq,i));
  int cd = bioseqtools.Utils.CalcIsCoding(seq,i);
  int cd1 = bioseqtools.Utils.CalcIsCodingGene(seq,i);
  int cd2 = bioseqtools.Utils.CalcIsCoding2(seq,i);
  int compl = bioseqtools.Utils.CalcIsGeneComplement(seq,i);
  System.out.println(i+"\t"+cn+"\t"+cn1+"\t"+cn2+"\t"+cnf+"\t"+fn+"\t"+(float)cd*6.5+"\t"+(float)cd1*6.5+"\t"+cd2+"\t"+(float)compl*0.5/*+VVectorCalc.Distance(KM.getCentroidVector(0),vecf)*/);
  all++;
  if(((cnf==0)&&(cd==0))||((cnf>0)&&(cd==1))) good1++;
  if(((cnf==0)&&(cd==0))||((cnf>0)&&(cd1==1))) good2++;
  if(((cnf==0)&&(cd==0))||((cnf>0)&&(cd2==1))) good3++;
  if(cd==0) alljunk1++;
  if(cd1==0) alljunk2++;
  if(cd2==0) alljunk3++;
  if((cd==0)&&(cnf==0)) goodjunk1++;
  if((cd1==0)&&(cnf==0)) goodjunk2++;
  if((cd2==0)&&(cnf==0)) goodjunk3++;
  if((cd==1)&&(cnf>0)) goodexon1++;
  if((cd1==1)&&(cnf>0)) goodexon2++;
  if((cd2==1)&&(cnf>0)) goodexon3++;
  }

System.out.println("% of coding: "+(float)(all-alljunk1)/all);
System.out.println("Quality:");
System.out.println(" CDS only: "+(float)good1/all+" in junk: "+(float)goodjunk1/alljunk1+"(from "+alljunk1+") in exons(Sn): "+(float)(goodexon1)/(all-alljunk1)+" (from "+(all-alljunk1)+")");
System.out.println("     Sn: "+(float)goodexon1/(all-alljunk1)+" Sp: "+(float)goodexon1/(goodexon1+alljunk1-goodjunk1));
System.out.println(" Whole genes: "+(float)good2/all+" in junk: "+(float)goodjunk2/alljunk2+"(from "+alljunk2+") in genes: "+(float)(goodexon2)/(all-alljunk2)+" (from "+(all-alljunk2)+")");
System.out.println("     Sn: "+(float)goodexon2/(all-alljunk2)+" Sp: "+(float)goodexon2/(goodexon2+alljunk2-goodjunk2));
System.out.println(" CDS+tRNA+rRNA: "+(float)good3/all+" in junk: "+(float)goodjunk3/alljunk3+"(from "+alljunk3+") in exons: "+(float)(goodexon3)/(all-alljunk3)+" (from "+(all-alljunk3)+")");
System.out.println("     Sn: "+(float)goodexon3/(all-alljunk3)+" Sp: "+(float)goodexon3/(goodexon3+alljunk3-goodjunk3));
} catch (Exception e) {System.out.println(e);};

}
//--------------------------------- ProduceAnnotationMultiply TASK 5
public static void ProduceAnnotationMultiply(boolean doOut){
try{

  if(sbeg==-1) sbeg=windowSize/2+2;
  if(send==-1) send=seq.length()-windowSize/2-2;
  if(sAnnotbeg==-1) sAnnotbeg=windowSize/2+2;
  if(sAnnotend==-1) sAnnotend=seq.length()-windowSize/2-2;

/*  VDataTable VT = new VDataTable();
  VT.LoadFromVDatFile(datFile);
  VD = new VDataSet();
  VTableSelector VS = new VTableSelector();
  VS.selectFirstNumericalColumns(VT,64);
  System.out.print("Loading dataset...");
  VD.loadFromDataTable(VT,VS);
  System.out.println(" Rec: "+VD.getPointCount());
  N = new VNormalization(VD.calcStatistics());
  VD.addToPreProcess(N);
  VD.preProcessData();*/

  if(taskModifier==1)
    VVectorCalc.distanceType=1;

  int NumberOfPoints = 0;

  while(seqI.hasNext()){
  seq = seqI.nextSequence();
  NumberOfPoints+=(int)seq.length()/everyPosition;
  }

  VD = new VDataSet();
  int n = NumberOfPoints;
  VD.massif = new float[n][64];
  VD.setPointCount(n);
  VD.setCoordCount(64);

  if (doOut) System.out.println("Calculating... "+n+" points");

  int veci[];

  LoadSeq();

  int k=0;
  while(seqI.hasNext()){
  seq = seqI.nextSequence();
  sbeg=windowSize/2+2;
  send=seq.length()-windowSize/2-2;
  for(int i=sbeg;i<send;i+=everyPosition)
  {
  //System.out.print(i+" ");
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  if(k<n)
    for(int j=0;j<64;j++) VD.massif[k][j]=veci[j];
  k++;
  }
  }

  VStatistics stat = VD.calcStatistics();
//  System.out.println(stat);
  N = new VNormalization(stat);
  VD.clearPreProcess();
  VD.addToPreProcess(N);
  VD.preProcessData();

VKMeansClusterization KM = new VKMeansClusterization(numberOfClasses);
if (doOut) System.out.println("Clusterization... "+numberOfClasses+" classes");
KM.setClusterizationQuality(0.00001);
KM.doClusterization(VD);
if (doOut) System.out.println("IN: "+KM.iterationsNumber+" Disp: "+KM.calculateAverageIntraClassStdDev(VD));
KM.rearrangeClusterNumbers();

//int veci[];
float vecf[] = new float[64];
int cn;

all=0;good1=0;good2=0;good3=0;alljunk1=0;alljunk2=0;alljunk3=0;goodjunk1=0;goodjunk2=0;goodjunk3=0;
goodexon1=0;goodexon2=0;goodexon3=0;

  LoadSeq();

while(seqI.hasNext()){
seq = seqI.nextSequence();

sAnnotbeg=windowSize/2+2;
sAnnotend=seq.length()-windowSize/2-2;

if (doOut) System.out.println("Annotation from "+sAnnotbeg+" to "+sAnnotbeg+" every "+parameter+"....");
for(int i=sAnnotbeg;i<sAnnotend;i+=parameter)
if(i<seq.length()-windowSize/2-2)
  {
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  for(int j=0;j<vecf.length;j++) vecf[j]=veci[j];
  N.process(vecf);
  cn = KM.calculateClosestCentroid(vecf);
//  cn = (int)((numberOfClasses+1)*Math.random());
  String fn = new String(bioseqtools.Utils.GetFeatureName(seq,i));
  int cd = bioseqtools.Utils.CalcIsCoding(seq,i);
  int cd1 = bioseqtools.Utils.CalcIsCodingGene(seq,i);
  int cd2 = bioseqtools.Utils.CalcIsCoding2(seq,i);
  int compl = bioseqtools.Utils.CalcIsGeneComplement(seq,i);
  if (doOut) System.out.println(i+"\t"+cn+"\t"+fn+"\t"+cd+"\t"+cd1+"\t"+cd2+"\t"+(float)compl*0.5);
  all++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd==1))) good1++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd1==1))) good2++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd2==1))) good3++;
  if(cd==0) alljunk1++;
  if(cd1==0) alljunk2++;
  if(cd2==0) alljunk3++;
  if((cd==0)&&(cn==0)) goodjunk1++;
  if((cd1==0)&&(cn==0)) goodjunk2++;
  if((cd2==0)&&(cn==0)) goodjunk3++;
  if((cd==1)&&(cn>0)) goodexon1++;
  if((cd1==1)&&(cn>0)) goodexon2++;
  if((cd2==1)&&(cn>0)) goodexon3++;
  }
}

if (doOut) System.out.println("% of coding: "+(float)(all-alljunk1)/all);
if (doOut) System.out.println("Quality:");
if (doOut) System.out.println(" CDS only: "+(float)good1/all+" in junk: "+(float)goodjunk1/alljunk1+"(from "+alljunk1+") in exons(Sn): "+(float)(goodexon1)/(all-alljunk1)+" (from "+(all-alljunk1)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon1/(all-alljunk1)+" Sp: "+(float)goodexon1/(goodexon1+alljunk1-goodjunk1));
if (doOut) System.out.println(" Whole genes: "+(float)good2/all+" in junk: "+(float)goodjunk2/alljunk2+"(from "+alljunk2+") in genes: "+(float)(goodexon2)/(all-alljunk2)+" (from "+(all-alljunk2)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon2/(all-alljunk2)+" Sp: "+(float)goodexon2/(goodexon2+alljunk2-goodjunk2));
if (doOut) System.out.println(" CDS+tRNA+rRNA: "+(float)good3/all+" in junk: "+(float)goodjunk3/alljunk3+"(from "+alljunk3+") in exons: "+(float)(goodexon3)/(all-alljunk3)+" (from "+(all-alljunk3)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon3/(all-alljunk3)+" Sp: "+(float)goodexon3/(goodexon3+alljunk3-goodjunk3));

} catch (Exception e) {System.out.println(e);}
}

//--------------------------------------------------------------------------------------

public static void ProduceAnnotationBayesian(boolean doOut){
try{

  if(sbeg==-1) sbeg=windowSize/2+2;
  if(send==-1) send=seq.length()-windowSize/2-2;
  if(sAnnotbeg==-1) sAnnotbeg=windowSize/2+2;
  if(sAnnotend==-1) sAnnotend=seq.length()-windowSize/2-2;

  if(taskModifier==1)
    VVectorCalc.distanceType=1;

  VD = new VDataSet();
  int n = (send-sbeg)/everyPosition;
  VD.massif = new float[n][64];
  VD.setPointCount(n);
  VD.setCoordCount(64);

  if (doOut) System.out.println("Calculating... "+n+" points");

  int veci[];
  int veci1[];
  int veci2[];
  int veci3[];

  int k=0;
  for(int i=sbeg;i<send;i+=everyPosition)
  {
  //System.out.print(i+" ");
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  if(k<n)
    for(int j=0;j<64;j++) VD.massif[k][j]=(float)veci[j]/(windowSize/3);
  k++;
  }

  VStatistics stat = VD.calcStatistics();
//  System.out.println(stat);
  N = new VNormalization(stat);
  VD.clearPreProcess();
  VD.addToPreProcess(N);
  VD.preProcessData();

VKMeansClusterization KM = new VKMeansClusterization(numberOfClasses);
if (doOut) System.out.println("Clusterization... "+numberOfClasses+" classes");
KM.setClusterizationQuality(0.00001);
KM.doClusterization(VD);
if (doOut) System.out.println("IN: "+KM.iterationsNumber+" Disp: "+KM.calculateAverageIntraClassStdDev(VD));
KM.rearrangeClusterNumbers();

//int veci[];
float vecf[] = new float[64];
int cn;

all=0;good1=0;good2=0;good3=0;alljunk1=0;alljunk2=0;alljunk3=0;goodjunk1=0;goodjunk2=0;goodjunk3=0;
goodexon1=0;goodexon2=0;goodexon3=0;

int koefBayes = 6;

if (doOut) System.out.println("Annotation from "+sAnnotbeg+" to "+sAnnotbeg+" every "+parameter+"....");
for(int i=sAnnotbeg;i<sAnnotend;i+=parameter)
if(i<seq.length()-windowSize/2-2)
  {
  veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
  veci1 = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2/koefBayes+2,i+windowSize/2/koefBayes);
//  veci2 = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2/koefBayes+1,i+windowSize/2/koefBayes);
//  veci3 = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2/koefBayes+2,i+windowSize/2/koefBayes);

//  for(int l=0;l<64;l++) System.out.print(veci[l]+" "); System.out.println();
//  for(int l=0;l<64;l++) System.out.print(veci1[l]+" "); System.out.println();
//  for(int l=0;l<64;l++) System.out.print(veci2[l]+" "); System.out.println();
//  for(int l=0;l<64;l++) System.out.print(veci3[l]+" "); System.out.println();

  for(int j=0;j<vecf.length;j++) vecf[j]=(float)veci[j]/(windowSize/3);
  N.process(vecf);
  cn = KM.calculateClosestCentroid(vecf);
//  for(int j=0;j<vecf.length;j++)
//      vecf[j]=(float)(KM.getCentroidVector(cn)[j]+(float)(veci1[j]/(windowSize/3)))/(1+1/koefBayes);
  N.unprocess(vecf);
  for(int j=0;j<vecf.length;j++)
      vecf[j]=(float)(vecf[j]+(float)veci1[j]/(windowSize/3))/(1+1/koefBayes);
  N.process(vecf);

//  cn = (int)((numberOfClasses+1)*Math.random());
  int cnb = KM.calculateClosestCentroid(vecf);
  cn=cnb;
  String fn = new String(bioseqtools.Utils.GetFeatureName(seq,i));
  int cd = bioseqtools.Utils.CalcIsCoding(seq,i);
  int cd1 = bioseqtools.Utils.CalcIsCodingGene(seq,i);
  int cd2 = bioseqtools.Utils.CalcIsCoding2(seq,i);
  int compl = bioseqtools.Utils.CalcIsGeneComplement(seq,i);
  if (doOut) System.out.println(i+"\t"+cn+"b"+cnb+"\t"+fn+"\t"+cd+"\t"+cd1+"\t"+cd2+"\t"+(float)compl*0.5);
  all++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd==1))) good1++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd1==1))) good2++;
  if(((cn==0)&&(cd==0))||((cn>0)&&(cd2==1))) good3++;
  if(cd==0) alljunk1++;
  if(cd1==0) alljunk2++;
  if(cd2==0) alljunk3++;
  if((cd==0)&&(cn==0)) goodjunk1++;
  if((cd1==0)&&(cn==0)) goodjunk2++;
  if((cd2==0)&&(cn==0)) goodjunk3++;
  if((cd==1)&&(cn>0)) goodexon1++;
  if((cd1==1)&&(cn>0)) goodexon2++;
  if((cd2==1)&&(cn>0)) goodexon3++;
  }

if (doOut) System.out.println("% of coding: "+(float)(all-alljunk1)/all);
if (doOut) System.out.println("Quality:");
if (doOut) System.out.println(" CDS only: "+(float)good1/all+" in junk: "+(float)goodjunk1/alljunk1+"(from "+alljunk1+") in exons(Sn): "+(float)(goodexon1)/(all-alljunk1)+" (from "+(all-alljunk1)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon1/(all-alljunk1)+" Sp: "+(float)goodexon1/(goodexon1+alljunk1-goodjunk1));
if (doOut) System.out.println(" Whole genes: "+(float)good2/all+" in junk: "+(float)goodjunk2/alljunk2+"(from "+alljunk2+") in genes: "+(float)(goodexon2)/(all-alljunk2)+" (from "+(all-alljunk2)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon2/(all-alljunk2)+" Sp: "+(float)goodexon2/(goodexon2+alljunk2-goodjunk2));
if (doOut) System.out.println(" CDS+tRNA+rRNA: "+(float)good3/all+" in junk: "+(float)goodjunk3/alljunk3+"(from "+alljunk3+") in exons: "+(float)(goodexon3)/(all-alljunk3)+" (from "+(all-alljunk3)+")");
if (doOut) System.out.println("     Sn: "+(float)goodexon3/(all-alljunk3)+" Sp: "+(float)goodexon3/(goodexon3+alljunk3-goodjunk3));

} catch (Exception e) {System.out.println(e);}
}
//------------------------------------------------------------------

public static void ExtractCodonUsage() throws Exception{
 if(sbeg==-1) sbeg=windowSize/2+2;
 if(send==-1) send=seq.length()-windowSize/2-2;
 if(sAnnotbeg==-1) sAnnotbeg=windowSize/2+2;
 if(sAnnotend==-1) sAnnotend=seq.length()-windowSize/2-2;
 if(taskModifier==1)
   VVectorCalc.distanceType=1;
 VD = new VDataSet();
 int n = (send-sbeg)/everyPosition;
 VD.massif = new float[n][64];
 VD.setPointCount(n);
 VD.setCoordCount(64);
 System.out.println("Calculating... "+n+" points");
 int veci[];
 int k=0;
 for(int i=sbeg;i<send;i+=everyPosition)
 {
 veci = bioseqtools.Utils.calcTripletFreq(seq,i-windowSize/2,i+windowSize/2);
 if(k<n)
   for(int j=0;j<64;j++) VD.massif[k][j]=(float)veci[j]/(windowSize/3.0f);
   //for(int j=0;j<64;j++) VD.massif[k][j]=veci[j];
 k++;
 }

 VStatistics stat = VD.calcStatistics();
 N = new VNormalization(stat);
 VD.clearPreProcess();
 VD.addToPreProcess(N);
 VD.preProcessData();

VKMeansClusterization KM = new VKMeansClusterization(numberOfClasses);
System.out.println("Clusterization... "+numberOfClasses+" classes");
KM.setClusterizationQuality(0.00001);
KM.doClusterization(VD);
System.out.println("IN: "+KM.iterationsNumber+" Disp: "+KM.calculateAverageIntraClassStdDev(VD));
KM.rearrangeClusterNumbers();

int stopcod1 = Utils.TripletNumber("taa");
int stopcod2 = Utils.TripletNumber("tag");
int stopcod3 = Utils.TripletNumber("tga");
int cstopcod1 = Utils.TripletNumber("taa");
int cstopcod2 = Utils.TripletNumber("tag");
int cstopcod3 = Utils.TripletNumber("tga");

//float ;
System.out.println("Junk Phase 0:");
float v[] = KM.getCentroidVectorF(0);
N.unprocess(v);
printDistribution(v);
System.out.println("Junk Phase 1:");
printDistribution(v);
System.out.println("Junk Phase 2:");
printDistribution(v);

float min = 10f;
int p1cl = -1;
for(int i=1;i<numberOfClasses;i++){
  v = KM.getCentroidVectorF(i);
  if((v[stopcod1]+v[stopcod2]+v[stopcod3])<min) { min=v[stopcod1]+v[stopcod2]+v[stopcod3];
  p1cl = i;}
}
min = 10f;
int c1cl = -1;
for(int i=1;i<numberOfClasses;i++){
  v = KM.getCentroidVectorF(i);
  if((v[cstopcod1]+v[cstopcod2]+v[cstopcod3])<min) { min=v[cstopcod1]+v[cstopcod2]+v[cstopcod3]; c1cl = i;}
}

System.out.println("Phase 0:");
v = KM.getCentroidVectorF(p1cl);
N.unprocess(v);
printDistribution(v);
k=1;
for(int i=1;(i<numberOfClasses)&&(k<3);i++){
if((i!=p1cl)&&(i!=c1cl)){
  System.out.println("Phase "+k+":");
  v = KM.getCentroidVectorF(i);
  N.unprocess(v);
  printDistribution(v);
  k++;
}
}

for(int i=1;(i<numberOfClasses);i++){
  System.out.println("All "+i+":");
  v = KM.getCentroidVectorF(i);
  N.unprocess(v);
  printDistribution(v);
}

}

private static void printDistribution(float f[]){
DecimalFormat df = new DecimalFormat();
df.applyPattern("#,#####0.00000");
for(int i=0;i<8;i++)
 {
 for(int j=0;j<8;j++){
   System.out.print(Utils.TripletName(i*8+j)+":"+df.format(f[i*8+j])+"()\t");
   }
System. out.println();
 }
}

private static String addNoise(String s, float insertionLevel, float deletionLevel, float changeLevel){
  StringBuffer sb = new StringBuffer(s);
  int insN = (int)(insertionLevel*s.length());
  int delN = (int)(deletionLevel*s.length());
  int chgN = (int)(changeLevel*s.length());
  int ac=0,cc=0,gc=0,tc=0;
  for(int i=0;i<sb.length();i++){
    switch(sb.charAt(i)){
       case 'a': ac++; break;
       case 'c': cc++; break;
       case 'g': gc++; break;
       case 't': tc++; break;
   }
  }
  System.out.println("Counted a:"+ac+" c:"+cc+" g:"+gc+" t:"+tc);
  Random r = new Random();
  for(int i=0;i<insN;i++){
    double d = r.nextDouble();
    if(i-(int)((float)i/1000f)*1000==0){
      System.out.println(i);
    }
    int pos = (int)(sb.length()*d);
    if((pos>=0)&&(pos<sb.length())){
      char c = chooseLetter(sb.length(),ac,cc,gc,tc);
      sb.insert(pos,c);
    }
  }
  System.out.println("Insertion done: "+sb.length());
  for(int i=0;i<delN;i++){
    if(i-(int)((float)i/1000f)*1000==0){
      System.out.println(i);
    }
    double d = r.nextDouble();
    int pos = (int)(sb.length()*d);
    if((pos>=0)&&(pos<sb.length())){
      char c = chooseLetter(sb.length(),ac,cc,gc,tc);
      sb.delete(pos,pos+1);
    }
  }
  System.out.println("Deletion done: "+sb.length());
  for(int i=0;i<chgN;i++){
    if(i-(int)((float)i/1000f)*1000==0){
      System.out.println(i);
    }
    double d = r.nextDouble();
    int pos = (int)(sb.length()*d);
    if((pos>=0)&&(pos<sb.length())){
      char c = chooseLetter(sb.length(),ac,cc,gc,tc);
      sb.setCharAt(pos,c);
    }
  }
  System.out.println("Change done: "+sb.length());
  return sb.toString();
}

private static char chooseLetter(int len, int ac, int cc, int gc, int tc){
  Random r = new Random();
  char res = 'a';
  double d = r.nextDouble();
  if(d<ac/len) res = 'a';
  if((d>=ac/len)&&(d<(ac+cc)/len)) res = 'c';
  if((d>=(ac+cc)/len)&&(d<(ac+cc+gc)/len)) res = 'g';
  if((d>=(ac+cc+gc)/len)&&(d<(ac+cc+gc+tc)/len)) res = 't';
  return res;
}

}