/* PROSPERO V.1.3

Copyright 

Richard Mott 2000

Wellcome Trust Centre For Human Genetics
University of Oxford
Roosevelt Drive
Oxford OX3 7AD
UK

The software package ARIADNE is distributed in the hope that it will be
useful, but in order that the University as a charitable foundation
protects its assets for the benefit of its educational and research
purposes, the University makes clear that no condition is made or to
be implied, nor is any warranty given or to be implied, as to the
accuracy of ARIADNE, or that it will be suitable for
any particular purpose or for use under any specific conditions, or that
the content or use of ARIADNE will not constitute
or result in infringement of third-party rights.  Furthermore, the
University disclaims all responsibility for the use which is made of 
ARIADNE.


*/

#include<stdio.h>
#include<math.h>
#include"cl.h"
#include"seq_util.h"
#include"gapstat.h"
#include"ariadne.h"




int main( int argc, char **argv ) {


  SEQUENCE *seq1=NULL, *seq2=NULL;
  PROFILE *pro=NULL;
  FILE *seqfp1, *seqfp2, *profp;
  char seqfilename1[256];
  char seqfilename2[256];
  char profilename[256];
  char matname[256];
  int **matrix;
  int top = 5;
  int starts, gaps, len;
  int score;
  int stats = 1;
  double *freq1, *freq2;
  double lambda0, K0, Kplus, s, r, pval, H, alpha, theta, logkappa;
  int A=11, B=1, pseudocount=100;
  int self=0;
  float ethresh=0;

  seqfp1 = ArgFile("-seq1", "r", argc, argv, seqfilename1 );
  seqfp2 = ArgFile("-seq2", "r", argc, argv, seqfilename2 );
  profp = ArgFile("-pro", "r", argc, argv, profilename );

  strcpy(matname, "blosum62");
  GetArg("-matrix", matname, argc, argv);
  GetInt("-A", &A, argc, argv);
  GetInt("-B", &B, argc, argv);
  GetInt("-top",&top,argc,argv);
  GetInt("-pseudocount", &pseudocount,argc,argv);
  GetFloat("-ethresh",&ethresh,argc,argv);

  GetHelp(argc,argv);

  printf( "***** PROSPERO v%s  %s *****\n\n", ARIADNE_VERSION, Now() );

  printf( "%s", ARIADNE_PREAMBLE );

  printf("using gap penalty %d+%dk\n", A, B );
  
  if ( profp ) {
    if (! (pro = ReadProfile( profp ))) {
      fprintf(stderr, "ERROR Could not read profile %s\n", profilename );
      exit(1);
    }
  }
  else {
    if ( ! *matname ) {
      fprintf(stderr, "ERROR: No matrix\n");
      exit(1);
    }
    else if ( ! (matrix = ReadBlastMatrix(matname)) ) {
      fprintf(stderr,"ERROR: Could not find matrix %s\n", matname );
      exit(1);
    }
    else {
      printf( "using matrix %s\n", matname );
    }
  }
  if ( ethresh > 0 ) 
    printf("printing all alignments with eval < %f\n", ethresh );
  else
    printf("showing top %d alignments\n", top );

  if ( seqfp1 ) {
    seq1 = GetFastaSeq( seqfp1 );
    UpCase( seq1->s);
    printf ("using sequence1 %s\n", seq1->name );
  }
  else {
    fprintf(stderr, "Error no sequence\n");
    exit(1);
  }

  if ( seqfp2 ) {
    if (! ( seq2 = GetFastaSeq( seqfp2 ) ) ) {
      fprintf(stderr, "Error no sequence\n");
      exit(1);
    }
    UpCase(seq2->s);
    printf ("using sequence2 %s\n", seq2->name );
  }
  
  if ( seq1 && (! seq2) && (! pro ) ) {
    self = 1;
    seq2 = seq1;
    printf( "using self-comparison\n");
  }
  if ( pro && seq1 ) {
    freq1 = PseudoResidueFrequencies( seq1, pseudocount );
    printf( "using profile %s\n", pro->name );
    if ( KarlinAltschulProfileStatistics( pro, freq1, &lambda0, &K0, &Kplus, &H, &r, &s ) ) {
      alpha = 2*s*exp(-lambda0*(A+B))/(1-exp(-lambda0*B));
      if ( alpha > 0.25 ) {
	printf( "WARNING: scoring scheme too lax (alpha = %e) for reliable statistics. \nTry increasing the gap penalty\n", alpha);
	stats = 0;
      }
    }
    else {
      printf( "WARNING: scoring scheme too lax for statistics.\n");
      stats = 0;
    } 

    TopAlignProteinToProfile( seq1, pro, A, B, top, ethresh, self, lambda0, K0, H, alpha );
    return 0;
  }
  else if ( seq1 && seq2 ) {
    pro = SeqToProfile(seq2, matrix); /* turn the query into a profile */
    freq1 = PseudoResidueFrequencies( seq1, pseudocount );
    freq2 = PseudoResidueFrequencies( seq2, pseudocount );

    if ( KarlinAltschulStatistics( matrix, freq1, freq2, &lambda0, &K0, &Kplus, &H, &r, &s ) ) {
      alpha = 2*s*exp(-lambda0*(A+B))/(1-exp(-lambda0*B));
      if ( alpha > 0.25 ) {
	printf( "WARNING: scoring scheme too lax (alpha = %e) for reliable statistics. \nTry increasing the gap penalty\n", alpha);
	stats = 0;
      }
    }
    else {
      printf( "WARNING: scoring scheme too lax for statistics.\n");
      stats = 0;
    } 

    TopAlignProteinToProfile( seq1, pro, A, B, top, ethresh, self, lambda0, K0, H, alpha );
    return 0;
  }
}


