@@ -37,7 +37,7 @@

 #include <stdio.h>    // printf

 #include <stdlib.h>   // exit

-#ifdef WINMATHH

+#ifdef MATH_H_CLIB

 #include <math.h>     // sqrt, pow

 #else

 #include <cmath>     // sqrt, pow

@@ -36,8 +36,13 @@

 #include <string.h>     // strcmp, strstr

 #include <stdio.h>    // printf

 #include <stdlib.h>   // exit

-//#include <math.h>     // sqrt, pow

+

+#ifdef WINMATHH

+#include <math.h>     // sqrt, pow

+#else

 #include <cmath>     // sqrt, pow

+#endif

+

 #include <limits.h>   // INT_MIN

 #include <float.h>    // FLT_MIN

 #include <time.h>     // clock

@@ -36,7 +36,8 @@

 #include <string.h>     // strcmp, strstr

 #include <stdio.h>    // printf

 #include <stdlib.h>   // exit

-#include <math.h>     // sqrt, pow

+//#include <math.h>     // sqrt, pow

+#include <cmath>     // sqrt, pow

 #include <limits.h>   // INT_MIN

 #include <float.h>    // FLT_MIN

 #include <time.h>     // clock

@@ -918,11 +918,10 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL,

         par.gapg = par.gapgV;

         par.gaph = par.gaphV;

         par.gapi = par.gapiV;

-

-        par.gapExtension = rHhalignPara.gapExtension;

-        par.gapOpening   = rHhalignPara.gapOpening;

-

     }

+

+    par.gapExtension = rHhalignPara.gapExtension;

+    par.gapOpening   = rHhalignPara.gapOpening;

     // Read input file (HMM, HHM, or alignment format), and add pseudocounts etc.

     q.cQT = 'q';

new file mode 100644

@@ -0,0 +1,1424 @@

+/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */

+

+/*********************************************************************

+ * Clustal Omega - Multiple sequence alignment

+ *

+ * Copyright (C) 2010 University College Dublin

+ *

+ * Clustal-Omega 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 file is part of Clustal-Omega.

+ *

+ ********************************************************************/

+

+/*

+ *  RCS $Id: hhalign.cpp 284 2013-06-12 10:10:11Z fabian $

+ */

+

+/* hhalign.C: 

+ * Align a multiple alignment to an alignment or HMM 

+ * Print out aligned input sequences in a3m format

+ * Compile:              g++ hhalign.C -o hhalign -I/usr/include/ -L/usr/lib -lpng -lz -O3 -fno-strict-aliasing 

+ * Compile with efence:  g++ hhalign.C -o hhalign -I/usr/include/ -lefence -L/usr/lib -lpng -lz -O -g  

+ * 

+ * Error codes: 0: ok  1: file format error  2: file access error  

+ *              3: memory error  4: internal numeric error  

+ *              5: command line error

+ */

+#undef PNG           /* include options for making png files? 

+			(will need the png library) */

+#define MAIN

+#include <iostream>   // cin, cout, cerr

+#include <fstream>    // ofstream, ifstream 

+#include <string.h>     // strcmp, strstr

+#include <stdio.h>    // printf

+#include <stdlib.h>   // exit

+#include <math.h>     // sqrt, pow

+#include <limits.h>   // INT_MIN

+#include <float.h>    // FLT_MIN

+#include <time.h>     // clock

+#include <ctype.h>    // islower, isdigit etc

+

+using std::cout;

+using std::cerr;

+using std::endl;

+using std::ios;

+using std::ifstream;

+using std::ofstream;

+

+int iAux_GLOBAL;

+

+#include "general.h"

+#include "util-C.h"     /* imax, fmax, iround, iceil, ifloor, strint, strscn, 

+			 strcut, substr, uprstr, uprchr, Basename etc. */

+#include "list-C.h"     // list data structure

+#include "hash-C.h"     // hash data structure

+

+#include "hhdecl-C.h"      // Constants, global variables, struct Parameters

+#include "hhutil-C.h"      /* MatchChr, InsertChr, aa2i, i2aa, log2, 

+			    fast_log2, WriteToScreen, */

+#include "hhmatrices-C.h"  // BLOSUM50, GONNET, HSDM

+#include "hhhmm.h"       // class HMM

+#include "hhhit.h"       // class Hit

+#include "hhalignment.h" // class Alignment

+#include "hhhalfalignment.h" // class HalfAlignment

+#include "hhfullalignment.h" // class FullAlignment

+#include "hhhitlist.h"   // class Hit

+#include "hhhmm-C.h"       // class HMM

+#include "hhalignment-C.h" // class Alignment

+#include "hhhit-C.h"       // class Hit

+#include "hhhalfalignment-C.h" // class HalfAlignment

+#include "hhfullalignment-C.h" // class FullAlignment

+#include "hhhitlist-C.h"   // class HitList 

+#include "hhfunc-C.h"      // some functions common to hh programs

+

+#ifdef PNG

+#include "pngwriter.h"   //PNGWriter (http://pngwriter.sourceforge.net/)

+#include "pngwriter.cc"  //PNGWriter (http://pngwriter.sourceforge.net/)

+#endif	    

+

+////////////////////////////////////////////////////////////////////////

+// Global variables 

+////////////////////////////////////////////////////////////////////////

+HMM *q;                 // Create query HMM with maximum of MAXRES match states

+HMM *t;                 /* Create template HMM with maximum of MAXRES 

+			  match states  */

+Alignment *qali;        /* (query alignment might be needed outside of hhfunc.C 

+			  for -a option) */

+Hit hit;               // Ceate new hit object pointed at by hit

+HitList hitlist;       /* list of hits with one Hit object for each 

+			  pairwise comparison done */

+char aliindices[256];  /* hash containing indices of all alignments 

+			  which to show in dot plot */

+char* dmapfile=NULL;   /* where to write the coordinates for the HTML map file 

+			  (to click the alignments) */

+char* alitabfile=NULL; // where to write pairs of aligned residues

+char* strucfile=NULL;  // where to read structure scores

+char* pngfile=NULL;    // pointer to pngfile

+char* tcfile=NULL;     // TCoffee output file name

+float probmin_tc=0.05; /* 5% minimum posterior probability for printing 

+			  pairs of residues for TCoffee */

+

+int dotW=10;           // average score of dot plot over window [i-W..i+W]

+float dotthr=0.5;      // score threshold for dot plot

+int dotscale=600;      // size scale of dotplot

+char dotali=0;         // show no alignments in dotplot

+float dotsat=0.3;      // saturation of grid and alignments in dot plot

+float pself=0.001;     // maximum p-value of 2nd and following self-alignments

+int Nstochali=0;       // # of stochastically traced alignments in dot plot

+float** Pstruc=NULL;   /* structure matrix which can be multiplied to prob 

+			  ratios from aa column comparisons in Forward() */

+float** Sstruc=NULL;   /* structure matrix which can be added to log odds 

+			  from aa column comparisons in Forward() */

+bool nucleomode = false;

+

+////////////////////////////////////////////////////////////////////////////

+// Help functions

+////////////////////////////////////////////////////////////////////////////

+/**

+ * @brief  general help function, not accessible in Clustal-Omega

+ */

+#if 0

+void 

+help()

+{

+  printf("\n");

+  printf("HHalign %s\n",VERSION_AND_DATE);

+  printf("Align a query alignment/HMM to a template alignment/HMM by HMM-HMM alignment\n");

+  printf("If only one alignment/HMM is given it is compared to itself and the best\n");

+  printf("off-diagonal alignment plus all further non-overlapping alignments above \n");

+  printf("significance threshold are shown.\n");

+  printf("%s",REFERENCE);

+  printf("%s",COPYRIGHT);

+  printf("\n");

+  printf("Usage: %s -i query [-t template] [options]  \n",program_name);

+  printf(" -i <file>     input query alignment  (fasta/a2m/a3m) or HMM file (.hhm)\n");

+  printf(" -t <file>     input template alignment (fasta/a2m/a3m) or HMM file (.hhm)\n");

+#ifdef PNG

+  printf(" -png <file>   write dotplot into PNG-file (default=none)           \n");

+#endif

+  printf("\n");         

+  printf("Output options:                                                           \n");

+  printf(" -o <file>     write output alignment to file\n"); 

+  printf(" -ofas <file>  write alignments in FASTA, A2M (-oa2m) or A3M (-oa3m) format   \n"); 

+  printf(" -a <file>     write query alignment in a3m format to file (default=none)\n");

+  printf(" -aa <file>    append query alignment in a3m format to file (default=none)\n");

+  printf(" -atab <file>  write alignment as a table (with posteriors) to file (default=none)\n");

+  printf(" -v <int>      verbose mode: 0:no screen output  1:only warings  2: verbose\n");

+  printf(" -seq  [1,inf[ max. number of query/template sequences displayed  (def=%i)  \n",par.nseqdis);

+  printf(" -nocons       don't show consensus sequence in alignments (default=show) \n");

+  printf(" -nopred       don't show predicted 2ndary structure in alignments (default=show) \n");

+  printf(" -nodssp       don't show DSSP 2ndary structure in alignments (default=show) \n");

+  printf(" -aliw int     number of columns per line in alignment list (def=%i)\n",par.aliwidth);

+  printf(" -P <float>    for self-comparison: max p-value of alignments (def=%.2g\n",pself);

+  printf(" -p <float>    minimum probability in summary and alignment list (def=%G) \n",par.p);

+  printf(" -E <float>    maximum E-value in summary and alignment list (def=%G)     \n",par.E);

+  printf(" -Z <int>      maximum number of lines in summary hit list (def=%i)       \n",par.Z);

+  printf(" -z <int>      minimum number of lines in summary hit list (def=%i)       \n",par.z);

+  printf(" -B <int>      maximum number of alignments in alignment list (def=%i)    \n",par.B);

+  printf(" -b <int>      minimum number of alignments in alignment list (def=%i)    \n",par.b);

+  printf(" -rank int     specify rank of alignment to write with -a or -aa option (default=1)\n");

+  printf("\n");         

+#ifdef PNG

+  printf("Dotplot options:\n");

+  printf(" -dwin <int>   average score in dotplot over window [i-W..i+W] (def=%i)      \n",dotW);

+  printf(" -dthr <float> score threshold for dotplot (default=%.2f)                    \n",dotthr);

+  printf(" -dsca <int>   if value <= 20: size of dot plot unit box in pixels           \n");

+  printf("               if value > 20: maximum dot plot size in pixels (default=%i)   \n",dotscale);

+  printf(" -dali <list>  show alignments with indices in <list> in dot plot            \n");

+  printf("               <list> = <index1> ... <indexN>  or  <list> = all              \n");

+  printf("\n");         

+#endif

+  printf("Filter input alignment (options can be combined):                         \n");

+  printf(" -id   [0,100] maximum pairwise sequence identity (%%) (def=%i)   \n",par.max_seqid);

+  printf(" -diff [0,inf[ filter most diverse set of sequences, keeping at least this    \n");

+  printf("               many sequences in each block of >50 columns (def=%i)\n",par.Ndiff);

+  printf(" -cov  [0,100] minimum coverage with query (%%) (def=%i) \n",par.coverage);

+  printf(" -qid  [0,100] minimum sequence identity with query (%%) (def=%i) \n",par.qid);

+  printf(" -qsc  [0,100] minimum score per column with query  (def=%.1f)\n",par.qsc);

+//   printf(" -csc  [0,100] minimum score per column with core alignment (def=%-.2f)\n",par.coresc);

+//   printf(" -qscc [0,100] minimum score per column of core sequence with query (def=%-.2f)\n",par.qsc_core);

+  printf("\n");         

+  printf("Input alignment format:                                                     \n");

+  printf(" -M a2m        use A2M/A3M (default): upper case = Match; lower case = Insert;\n");         

+  printf("               '-' = Delete; '.' = gaps aligned to inserts (may be omitted)   \n");

+  printf(" -M first      use FASTA: columns with residue in 1st sequence are match states\n");

+  printf(" -M [0,100]    use FASTA: columns with fewer than X%% gaps are match states   \n");

+  printf("\n");    

+  printf("HMM-HMM alignment options:                                                  \n");

+  printf(" -glob/-loc    global or local alignment mode (def=local)         \n");

+  printf(" -alt <int>    show up to this number of alternative alignments (def=%i)    \n",par.altali);

+  printf(" -vit          use Viterbi algorithm for alignment instead of MAC algorithm \n");

+  printf(" -mac          use Maximum Accuracy (MAC) alignment (default)  \n");

+  printf(" -mact [0,1[   posterior probability threshold for MAC alignment (def=%.3f) \n",par.mact);

+  printf("               A threshold value of 0.0 yields global alignments.\n");

+  printf(" -sto <int>    use global stochastic sampling algorithm to sample this many alignments\n");

+  printf(" -excl <range> exclude query positions from the alignment, e.g. '1-33,97-168'\n");

+  printf(" -shift [-1,1] score offset (def=%-.3f)                                      \n",par.shift);

+  printf(" -corr [0,1]   weight of term for pair correlations (def=%.2f)               \n",par.corr);

+  printf(" -ssm  0-4     0:no ss scoring [default=%i]               \n",par.ssm);

+  printf("               1:ss scoring after alignment                                  \n");

+  printf("               2:ss scoring during alignment                                 \n");

+  printf(" -ssw  [0,1]   weight of ss score  (def=%-.2f)                               \n",par.ssw);

+  printf("\n");

+  printf(" -def          read default options from ./.hhdefaults or <home>/.hhdefault. \n");

+  printf("\n");

+  printf("Example: %s -i T0187.a3m -t d1hz4a_.hhm -png T0187pdb.png \n",program_name);

+  cout<<endl;

+//   printf("More help:                                                         \n");

+//   printf(" -h out        output options                                      \n");

+//   printf(" -h hmm        options for building HMM from multiple alignment    \n");

+//   printf(" -h gap        options for setting gap penalties                   \n");

+//   printf(" -h ali        options for HMM-HMM alignment                       \n");

+//   printf(" -h all        all options \n");

+}

+#endif

+

+/**

+ * @brief  helpt for output, not accessible in Clustal-Omega

+ */

+#if 0

+void 

+help_out()

+{

+  printf("\n");

+  printf("Output options: \n");

+  printf(" -v            verbose mode (default: show only warnings)                 \n");

+  printf(" -v 0          suppress all screen output                                 \n");

+  printf(" -nocons       don't show consensus sequence in alignments (default=show) \n");

+  printf(" -nopred       don't show predicted 2ndary structure in alignments (default=show) \n");

+  printf(" -nodssp       don't show DSSP SS 2ndary structure in alignments (default=show) \n");

+  printf(" -seq  [1,inf[ max. number of query/template sequences displayed  (def=%i)  \n",par.nseqdis);

+  printf(" -aliw [40,..[ number of columns per line in alignment list (def=%i)\n",par.aliwidth);

+  printf(" -P <float>    for self-comparison: max p-value of alignments (def=%.2g\n",pself);

+  printf(" -p <float>    minimum probability in summary and alignment list (def=%G) \n",par.p);

+  printf(" -E <float>    maximum E-value in summary and alignment list (def=%G)     \n",par.E);

+  printf(" -Z <int>      maximum number of lines in summary hit list (def=%i)       \n",par.Z);

+  printf(" -z <int>      minimum number of lines in summary hit list (def=%i)       \n",par.z);

+  printf(" -B <int>      maximum number of alignments in alignment list (def=%i)    \n",par.B);

+  printf(" -b <int>      minimum number of alignments in alignment list (def=%i)    \n",par.b);

+  printf(" -rank int     specify rank of alignment to write with -a or -aa option (def=1)\n");

+  printf(" -tc <file>    write a TCoffee library file for the pairwise comparison   \n");         

+  printf(" -tct [0,100]  min. probobability of residue pairs for TCoffee (def=%i%%)\n",iround(100*probmin_tc));         

+  printf("\n");         

+#ifdef PNG

+  printf("Dotplot options:\n");

+  printf(" -dwin int     average score in dotplot over window [i-W..i+W] (def=%i)   \n",dotW);

+  printf(" -dthr float   score threshold for dotplot (default=%.2f)                 \n",dotthr);

+  printf(" -dsca int     size of dot plot box in pixels  (default=%i)               \n",dotscale);

+  printf(" -dali <list>  show alignments with indices in <list> in dot plot\n");

+  printf("               <list> = <index1> ... <indexN>  or  <list> = all              \n");

+  printf(" -dmap <file>  print list of coordinates in png plot  \n");

+#endif

+}

+#endif

+

+/**

+ * @brief  help hit HMM options, not accessible in Clustal-Omega

+ */

+#if 0

+void 

+help_hmm()

+{

+  printf("\n");

+  printf("Options to filter input alignment (options can be combined):              \n");

+  printf(" -id   [0,100] maximum pairwise sequence identity (%%) (def=%i)   \n",par.max_seqid);

+  printf(" -diff [0,inf[ filter most diverse set of sequences, keeping at least this    \n");

+  printf("               many sequences in each block of >50 columns (def=%i)\n",par.Ndiff);

+  printf(" -cov  [0,100] minimum coverage with query (%%) (def=%i) \n",par.coverage);

+  printf(" -qid  [0,100] minimum sequence identity with query (%%) (def=%i) \n",par.qid);

+  printf(" -qsc  [0,100] minimum score per column with query  (def=%.1f)\n",par.qsc);

+//   printf(" -csc  [0,100] minimum score per column with core alignment (def=%-.2f)\n",par.coresc);

+//   printf(" -qscc [0,100] minimum score per column of core sequence with query (def=%-.2f)\n",par.qsc_core);

+  printf("                                                                          \n");

+  printf("HMM-building options:                                                     \n");

+  printf(" -M a2m        use A2M/A3M (default): upper case = Match; lower case = Insert;\n");         

+  printf("               '-' = Delete; '.' = gaps aligned to inserts (may be omitted)   \n");

+  printf(" -M first      use FASTA: columns with residue in 1st sequence are match states\n");

+  printf(" -M [0,100]    use FASTA: columns with fewer than X%% gaps are match states   \n");

+  printf(" -tags         do NOT neutralize His-, C-myc-, FLAG-tags, and \n");

+  printf("               trypsin recognition sequence to background distribution    \n");

+  printf("                                                                          \n");

+  printf("Pseudocount options:                                                      \n");

+  printf(" -Gonnet       use the Gonnet substitution matrix (default)               \n");

+  printf(" -Blosum50     use the Blosum50 substitution matrix                       \n");

+  printf(" -Blosum62     use the Blosum62 substitution matrix                       \n");

+  printf(" -HSDM         use the structure-derived HSDM substitution matrix         \n");

+  printf(" -pcm  0-3     Pseudocount mode (default=%-i)                             \n",par.pcm);

+  printf("               tau = substitution matrix pseudocount admixture            \n");

+  printf("               0: no pseudo counts:     tau = 0                           \n");

+  printf("               1: constant              tau = a                           \n");

+  printf("               2: divergence-dependent: tau = a/(1 + ((Neff-1)/b)^c)       \n");

+  printf("                  Neff=( (Neff_q^d+Neff_t^d)/2 )^(1/d)                       \n");

+  printf("                  Neff_q = av number of different AAs per column in query  \n");

+  printf("               3: column-specific:      tau = \'2\' * (Neff(i)/Neff)^(w*Neff/20)\n");

+  printf(" -pca  [0,1]   set a (overall admixture) (def=%-.1f)                      \n",par.pca);

+  printf(" -pcb  [1,inf[ set b (threshold for Neff) (def=%-.1f)                      \n",par.pcb);

+  printf(" -pcc  [0,3]   set c (extinction exponent for tau(Neff))  (def=%-.1f)      \n",par.pcc);

+  printf(" -pcw  [0,3]   set w (weight of pos-specificity for pcs) (def=%-.1f)      \n",par.pcw);

+}

+#endif

+

+/**

+ * @brief  help with gap handling, not accessible in Clustal-Omega

+ */

+#if 0

+void 

+help_gap()

+{

+  printf("\n");

+  printf("Gap cost options:                                                         \n");

+  printf(" -gapb [0,inf[ transition pseudocount admixture (def=%-.2f)               \n",par.gapb);

+  printf(" -gapd [0,inf[ Transition pseudocount admixture for opening gap (default=%-.2f)\n",par.gapd);

+  printf(" -gape [0,1.5] Transition pseudocount admixture for extending gap (def=%-.1f)\n",par.gape);

+  printf(" -gapf ]0,inf] factor for increasing/reducing the gap open penalty for deletes (def=%-.2f)\n",par.gapf);

+  printf(" -gapg ]0,inf] factor for increasing/reducing the gap open penalty for deletes (def=%-.2f)\n",par.gapg);

+  printf(" -gaph ]0,inf] factor for increasing/reducing the gap extension penalty for deletes(def=%-.2f)\n",par.gaph);

+  printf(" -gapi ]0,inf] factor for increasing/reducing the gap extension penalty for inserts(def=%-.2f)\n",par.gapi);

+  printf(" -egq  [0,inf[ penalty (bits) for end gaps aligned to query residues (def=%-.2f)\n",par.egq);

+  printf(" -egt  [0,inf[ penalty (bits) for end gaps aligned to template residues (def=%-.2f)\n",par.egt);

+}

+#endif

+

+/**

+ * @brief  help with alignment options, not accessible in Clustal-Omega

+ */

+#if 0

+void 

+help_ali()

+{

+  printf("\n");

+  printf("Alignment options:  \n");

+  printf(" -glob/-loc    global or local alignment mode (def=global)                \n");

+  printf(" -mac          use Maximum Accuracy (MAC) alignment instead of Viterbi\n");

+  printf(" -mact [0,1]   posterior prob threshold for MAC alignment (def=%.3f)      \n",par.mact);

+  printf(" -sto <int>    use global stochastic sampling algorithm to sample this many alignments\n");

+  printf(" -sc   <int>   amino acid score         (tja: template HMM at column j) (def=%i)\n",par.columnscore);

+  printf("        0      = log2 Sum(tja*qia/pa)   (pa: aa background frequencies)    \n");

+  printf("        1      = log2 Sum(tja*qia/pqa)  (pqa = 1/2*(pa+ta) )               \n");

+  printf("        2      = log2 Sum(tja*qia/ta)   (ta: av. aa freqs in template)     \n");

+  printf("        3      = log2 Sum(tja*qia/qa)   (qa: av. aa freqs in query)        \n");

+   printf(" -corr [0,1]   weight of term for pair correlations (def=%.2f)            \n",par.corr);

+  printf(" -shift [-1,1] score offset (def=%-.3f)                                   \n",par.shift);

+  printf(" -r            repeat identification: multiple hits not treated as independent\n");

+  printf(" -ssm  0-2     0:no ss scoring [default=%i]               \n",par.ssm);

+  printf("               1:ss scoring after alignment                               \n");

+  printf("               2:ss scoring during alignment                              \n");

+  printf(" -ssw  [0,1]   weight of ss score compared to column score (def=%-.2f)    \n",par.ssw);

+  printf(" -ssa  [0,1]   ss confusion matrix = (1-ssa)*I + ssa*psipred-confusion-matrix [def=%-.2f)\n",par.ssa);

+}

+#endif

+

+/**

+ * @brief  general help menu, not accessible in Clustal-Omega

+ */

+#if 0

+void 

+help_all()

+{

+  help();

+  help_out();

+  help_hmm();

+  help_gap();

+  help_ali();

+  printf("\n");

+  printf("Default options can be specified in './.hhdefaults' or '~/.hhdefaults'\n");

+}

+#endif

+

+/////////////////////////////////////////////////////////////////////////////

+//// Processing input options from command line and .hhdefaults file

+/////////////////////////////////////////////////////////////////////////////

+/**

+ * @brief  process arguments from commandline, not accessible from Clustal-Omega

+ */

+#if 0

+void 

+ProcessArguments(int argc, char** argv)

+{

+  //Processing command line input

+  for (int i=1; i<argc; i++)

+    { 

+      if (v>=4) cout<<i<<"  "<<argv[i]<<endl; //PRINT

+      if (!strcmp(argv[i],"-i"))

+	{

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no query file following -i\n"; throw 4;}

+	  else strcpy(par.infile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-t"))

+	{

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no template file following -d\n"; throw 4;}

+	  else strcpy(par.tfile,argv[i]); /* FS, ProcessArguments */

+	}

+      else if (!strcmp(argv[i],"-o"))

+	{

+	  if (++i>=argc) 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no filename following -o\n"; throw 4;}

+	  else strcpy(par.outfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-ofas"))

+	{

+	  par.outformat=1;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -o\n"; throw 4;}

+	  else strcpy(par.pairwisealisfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-oa2m"))

+	{

+	  par.outformat=2;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -o\n"; throw 4;}

+	  else strcpy(par.pairwisealisfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-oa3m"))

+	{

+	  par.outformat=3;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -o\n"; throw 4;}

+	  else strcpy(par.pairwisealisfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-rank") && (i<argc-1)) par.hitrank=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-Oa3m"))

+	{

+	  par.append=0;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Oa3m\n"; throw 4;}

+	  else strcpy(par.alnfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-Aa3m"))

+	{

+	  par.append=1;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Aa3m\n"; throw 4;}

+	  else strcpy(par.alnfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-Ohhm"))

+	{

+	  par.append=0;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Ohhm\n"; throw 4;}

+	  else strcpy(par.hhmfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-Ahhm"))

+	{

+	  par.append=1;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Ahhm\n"; throw 4;}

+	  else strcpy(par.hhmfile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-Opsi"))

+	{

+	  par.append=0;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Opsi\n"; throw 4;}

+	  else strcpy(par.psifile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-Apsi"))

+	{

+	  par.append=1;

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Apsi\n"; throw 4;}

+	  else strcpy(par.psifile,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-png"))

+	{

+	  if (++i>=argc) 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no filename following -png\n"; throw 4;}

+	  else 

+	    {

+	      pngfile = new(char[strlen(argv[i])+1]);

+	      strcpy(pngfile,argv[i]);

+	    }

+	}

+      else if (!strcmp(argv[i],"-Struc"))

+	{

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no query file following -Struc\n"; throw 4;}

+	  else 

+	    {

+	      strucfile = new(char[strlen(argv[i])+1]);

+	      strcpy(strucfile,argv[i]);

+	    }

+	}

+      else if (!strcmp(argv[i],"-atab") || !strcmp(argv[i],"-Aliout"))

+	{

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no query file following -Struc\n"; throw 4;}

+	  else 

+	    {

+	      alitabfile = new(char[strlen(argv[i])+1]);

+	      strcpy(alitabfile,argv[i]);

+	    }

+	}

+      else if (!strcmp(argv[i],"-tc"))

+	{

+	  if (++i>=argc || argv[i][0]=='-') 

+	    {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Opsi\n"; throw 4;}

+	  else 

+	    {

+	      tcfile = new(char[strlen(argv[i])+1]);

+	      strcpy(tcfile,argv[i]);

+	    }

+	}

+      else if (!strcmp(argv[i],"-h")|| !strcmp(argv[i],"--help"))

+	{

+	  if (++i>=argc) {help(); throw 0;}

+	  if (!strcmp(argv[i],"out")) {help_out(); throw 0;}

+	  if (!strcmp(argv[i],"hmm")) {help_hmm(); throw 0;}

+	  if (!strcmp(argv[i],"gap")) {help_gap(); throw 0;}

+	  if (!strcmp(argv[i],"ali")) {help_ali(); throw 0;}

+	  if (!strcmp(argv[i],"all")) {help_all(); throw 0;}

+	  else {help(); throw 0;}

+	}

+      else if (!strcmp(argv[i],"-excl"))

+	{

+	  if (++i>=argc) {help(); throw 4;}

+	  par.exclstr = new(char[strlen(argv[i])+1]);

+	  strcpy(par.exclstr,argv[i]);

+	}

+      else if (!strcmp(argv[i],"-v") && (i<argc-1) && argv[i+1][0]!='-' ) v=atoi(argv[++i]);

+      else if (!strcmp(argv[i],"-v"))  v=2;

+      else if (!strcmp(argv[i],"-v0")) v=0;

+      else if (!strcmp(argv[i],"-v1")) v=1;

+      else if (!strcmp(argv[i],"-v2")) v=2;

+      else if (!strcmp(argv[i],"-v3")) v=3;

+      else if (!strcmp(argv[i],"-v4")) v=4;

+      else if (!strcmp(argv[i],"-v5")) v=5;

+      else if (!strcmp(argv[i],"-P") && (i<argc-1)) pself=atof(argv[++i]);

+      else if (!strcmp(argv[i],"-p") && (i<argc-1)) par.p = atof(argv[++i]);

+      else if (!strcmp(argv[i],"-e") && (i<argc-1)) par.E = atof(argv[++i]);

+      else if (!strcmp(argv[i],"-E") && (i<argc-1)) par.E = atof(argv[++i]);

+      else if (!strcmp(argv[i],"-b") && (i<argc-1)) par.b = atoi(argv[++i]);

+      else if (!strcmp(argv[i],"-B") && (i<argc-1)) par.B = atoi(argv[++i]);

+      else if (!strcmp(argv[i],"-z") && (i<argc-1)) par.z = atoi(argv[++i]);

+      else if (!strcmp(argv[i],"-Z") && (i<argc-1)) par.Z = atoi(argv[++i]);

+      else if (!strncmp(argv[i],"-nocons",7)) par.showcons=0;

+      else if (!strncmp(argv[i],"-nopred",7)) par.showpred=0;

+      else if (!strncmp(argv[i],"-nodssp",7)) par.showdssp=0;

+      else if (!strncmp(argv[i],"-mark",7)) par.mark=1;

+      else if (!strcmp(argv[i],"-seq") && (i<argc-1))  par.nseqdis=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-aliw") && (i<argc-1)) par.aliwidth=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-id") && (i<argc-1))   par.max_seqid=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-tct") && (i<argc-1))  probmin_tc=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-dwin") && (i<argc-1)) dotW=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-dsca") && (i<argc-1)) dotscale=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-dthr") && (i<argc-1)) dotthr=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-dali") && (i<argc-1))  

+	{

+	  dotali=1; 

+	  for (int index=0; index<256; index++) aliindices[index]=0;

+	  while (i+1<argc && argv[i+1][0]!='-') // adds index to hash aliindices

+	    {

+	      i++;

+	      if (strcmp(argv[i],"all")) aliindices[atoi(argv[i])]=1;	      

+	      else dotali=2;

+	    }

+	}

+      else if (!strcmp(argv[i],"-dmap"))  

+	{

+	  if (++i>=argc) 

+	    {help(); cerr<<endl<<"Error in "<<program_name<<": no filename following -o\n"; throw 4;}

+	  else 

+	    {

+	      dmapfile = new(char[strlen(argv[i])+1]);

+	      strcpy(dmapfile,argv[i]);

+	    }

+	}

+      else if (!strcmp(argv[i],"-dsat") && (i<argc-1)) dotsat=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-qid") && (i<argc-1))  par.qid=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-qsc") && (i<argc-1))  par.qsc=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-cov") && (i<argc-1))  par.coverage=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-diff") && (i<argc-1)) par.Ndiff=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-qscc") && (i<argc-1))    par.qsc_core=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-csc") && (i<argc-1))     par.coresc=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-Gonnet")) par.matrix=0; 

+      else if (!strcmp(argv[i],"-HSDM")) par.matrix=1; 

+      else if (!strcmp(argv[i],"-BLOSUM50")) par.matrix=2; 

+      else if (!strcmp(argv[i],"-Blosum50")) par.matrix=2; 

+      else if (!strcmp(argv[i],"-B50")) par.matrix=2; 

+      else if (!strcmp(argv[i],"-BLOSUM62")) par.matrix=3; 

+      else if (!strcmp(argv[i],"-Blosum62")) par.matrix=3; 

+      else if (!strcmp(argv[i],"-B62")) par.matrix=3; 

+      else if (!strcmp(argv[i],"-pcm") && (i<argc-1)) par.pcm=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-pca") && (i<argc-1)) par.pca=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-pcb") && (i<argc-1)) par.pcb=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-pcc") && (i<argc-1)) par.pcc=atof(argv[++i]);  

+      else if (!strcmp(argv[i],"-pcw") && (i<argc-1)) par.pcw=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-gapb") && (i<argc-1)) { par.gapb=atof(argv[++i]); if (par.gapb<=0.01) par.gapb=0.01;} 

+      else if (!strcmp(argv[i],"-gapd") && (i<argc-1)) par.gapd=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-gape") && (i<argc-1)) par.gape=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-gapf") && (i<argc-1)) par.gapf=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-gapg") && (i<argc-1)) par.gapg=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-gaph") && (i<argc-1)) par.gaph=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-gapi") && (i<argc-1)) par.gapi=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-egq") && (i<argc-1)) par.egq=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-egt") && (i<argc-1)) par.egt=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-ssgap")) par.ssgap=1;

+      else if (!strcmp(argv[i],"-ssgapd") && (i<argc-1)) par.ssgapd=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-ssgape") && (i<argc-1)) par.ssgape=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-ssgapi") && (i<argc-1)) par.ssgapi=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-ssm") && (i<argc-1)) par.ssm=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-ssw") && (i<argc-1)) par.ssw=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-ssa") && (i<argc-1)) par.ssa=atof(argv[++i]); 

+      else if (!strncmp(argv[i],"-glo",3)) {par.loc=0; if (par.mact>0.3 && par.mact<0.301) {par.mact=0;} }

+      else if (!strncmp(argv[i],"-loc",3)) par.loc=1;

+      else if (!strncmp(argv[i],"-alt",4) && (i<argc-1)) par.altali=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-map") || !strcmp(argv[i],"-MAP")) par.forward=2; 

+      else if (!strcmp(argv[i],"-mac") || !strcmp(argv[i],"-MAC")) par.forward=2; 

+      else if (!strcmp(argv[i],"-vit")) par.forward=0; 

+      else if (!strcmp(argv[i],"-sto") && (i<argc-1))  {Nstochali=atoi(argv[++i]); par.forward=1;}

+      else if (!strcmp(argv[i],"-r")) par.repmode=1; 

+      else if (!strcmp(argv[i],"-M") && (i<argc-1)) 

+	if (!strcmp(argv[++i],"a2m") || !strcmp(argv[i],"a3m"))  par.M=1; 

+	else if(!strcmp(argv[i],"first"))  par.M=3; 

+	else if (argv[i][0]>='0' && argv[i][0]<='9') {par.Mgaps=atoi(argv[i]); par.M=2;}

+	else cerr<<endl<<"WARNING: Ignoring unknown argument: -M "<<argv[i]<<"\n";

+      else if (!strcmp(argv[i],"-shift") && (i<argc-1)) par.shift=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-mact") && (i<argc-1)) {par.mact=atof(argv[++i]); par.forward=2;}

+      else if (!strcmp(argv[i],"-wstruc") && (i<argc-1)) par.wstruc=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-opt") && (i<argc-1)) par.opt=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-sc") && (i<argc-1)) par.columnscore=atoi(argv[++i]); 

+      else if (!strcmp(argv[i],"-corr") && (i<argc-1)) par.corr=atof(argv[++i]); 

+      else if (!strcmp(argv[i],"-def")) par.readdefaultsfile=1; 

+      else if (!strcmp(argv[i],"-ovlp") && (i<argc-1)) par.min_overlap=atoi(argv[++i]);

+      else if (!strcmp(argv[i],"-tags")) par.notags=0;

+      else if (!strcmp(argv[i],"-notags")) par.notags=1;

+      else cerr<<endl<<"WARNING: Ignoring unknown option "<<argv[i]<<" ...\n";

+      if (v>=4) cout<<i<<"  "<<argv[i]<<endl; //PRINT

+    } // end of for-loop for command line input

+}

+#endif

+

+

+

+

+/////////////////////////////////////////////////////////////////////////////

+//// MAIN PROGRAM

+/////////////////////////////////////////////////////////////////////////////

+/**

+ *

+ * hhalign()

+ *

+ * @param[in,out] ppcFirstProf

+ * first profile to be aligned

+ * @param[in] iFirstCnt

+ * number of sequences in 1st profile

+ * @param[in,out] ppcSecndProf

+ * second profile to be aligned

+ * @param[in] iSecndCnt

+ * number of sequences in 2nd profile

+ * @param[out] dScore_p

+ * score of alignment

+ * @param[in] prHMM

+ * HMM info of external HMM (background)

+ * @param[in] pcPrealigned1

+ * (gapped) 1st sequence aligned to background HMM

+ * @param[in] pcRepresent1

+ * (gapped) sequence representing background HMM aligned to 1st sequence

+ * @param[in] pcPrealigned2

+ * (gapped) 2nd sequence aligned to background HMM

+ * @param[in] pcRepresent2

+ * (gapped) sequence representing background HMM aligned to 2nd sequence

+ * @param[in] rHhalignPara,

+ * various parameters passed down from commandline

+ * e.g., iMaxRamMB

+ * @param[out] rHHscores

+ * qualify goodness of alignment

+ *

+ * iFlag,zcAux,zcError are debugging arguments

+ *

+ * @return Non-zero on error

+ */

+extern "C" int 

+hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL,

+        char **ppcSecndProf, int iSecndCnt, double *pdWeightsR,

+        double *dScore_p, hmm_light *prHMM,

+        char *pcPrealigned1, char *pcRepresent1, 

+        char *pcPrealigned2, char *pcRepresent2, 

+        hhalign_para rHhalignPara, hhalign_scores *rHHscores, 

+        int iFlag, int iVerbosity, 

+        char zcAux[], char zcError[]) {

+

+	v = 3;

+

+#ifdef CLUSTALO

+    int iRetVal = RETURN_OK;

+    iAux_GLOBAL = iFlag;

+#ifndef CLUSTALO_NOFILE

+    int argc = 0;

+    char **argv = NULL;

+    argv = (char **)malloc(argc*sizeof(char *));

+    for (int i = 0; i < argc; i++){

+        argv[i] = (char *)malloc(100);

+    }

+    strcpy(argv[0], "./hhalign");

+    strcpy(argv[1], "-t");

+    strcpy(argv[2], "hhalign.C");

+    strcpy(argv[3], "-i");

+    strcpy(argv[4], "hhalign.C");

+    strcpy(argv[5], "-ofas");

+    strcpy(argv[6], "out");

+#endif

+#endif

+    /*char* argv_conf[MAXOPT];*/ /* Input arguments from .hhdefaults file 

+      (first=1: argv_conf[0] is not used) */

+    /*int argc_conf;*/           // Number of arguments in argv_conf 

+    /*char inext[IDLEN]="";*/    // Extension of query input file (hhm or a3m) 

+    /*char text[IDLEN]="";*/     // Extension of template input file (hhm or a3m)

+    /*int** ali=NULL;*/          // ali[i][j]=1 if (i,j) is part of an alignment

+    /*int** alisto=NULL;*/       // ali[i][j]=1 if (i,j) is part of an alignment

+    int Nali;                    /* number of normally backtraced alignments 

+                                    in dot plot */

+    

+    SetDefaults(); 

+    strcpy(par.tfile,""); /* FS, Initialise Argument */

+    strcpy(par.alnfile,""); 

+    par.p=0.0 ;                  /* minimum threshold for inclusion in hit list 

+                                    and alignment listing */

+    par.E=1e6;                   /* maximum threshold for inclusion in hit list 

+                                    and alignment listing */

+    par.b=1;                     // min number of alignments

+    par.B=100;                   // max number of alignments

+    par.z=1;                     // min number of lines in hit list

+    par.Z=100;                   // max number of lines in hit list

+    par.append=0;                /* append alignment to output file 

+                                    with -a option */

+    par.altali=1;                /* find only ONE (possibly overlapping) 

+                                    subalignment  */

+    par.hitrank=0;               /* rank of hit to be printed as a3m alignment 

+                                    (default=0) */

+    par.outformat=3;             // default output format for alignment is a3m

+    hit.self=0;                  // no self-alignment

+    par.forward=2;               /* 0: Viterbi algorithm; 

+                                    1: Viterbi+stochastic sampling; 

+                                    2:Maximum Accuracy (MAC) algorithm */

+    

+    // Make command line input globally available

+#ifndef CLUSTALO_NOFILE

+    par.argv=argv; 

+    par.argc=argc;

+    RemovePathAndExtension(program_name,argv[0]); 

+#endif

+    

+#ifndef CLUSTALO_NOFILE

+    /* Enable changing verbose mode before defaults file 

+       and command line are processed */

+    for (int i=1; i<argc; i++)

+        { 

+            if (!strcmp(argv[i],"-def")) par.readdefaultsfile=1;

+            else if (argc>1 && !strcmp(argv[i],"-v0")) v=0;

+            else if (argc>1 && !strcmp(argv[i],"-v1")) v=1;

+            else if (argc>2 && !strcmp(argv[i],"-v")) v=atoi(argv[i+1]);

+        }

+    

+    // Read .hhdefaults file?

+    if (par.readdefaultsfile) 

+        {

+            // Process default otpions from .hhconfig file

+            ReadDefaultsFile(argc_conf,argv_conf);

+            ProcessArguments(argc_conf,argv_conf);

+        }

+#endif

+    /* Process command line options 

+       (they override defaults from .hhdefaults file) */

+#ifndef CLUSTALO_NOFILE

+    ProcessArguments(argc,argv);

+#endif

+    

+    

+#ifdef CLUSTALO

+    int iAuxLen1 = strlen(ppcFirstProf[0]);

+    int iAuxLen2 = strlen(ppcSecndProf[0]); 

+    if ( (0 == iAuxLen1) || (0 == iAuxLen2) ){ /* problem with empty profiles, FS, r249 -> r250 */

+        sprintf(zcError, "%s:%s:%d: strlen(prof1)=%d, strlen(prof2)=%d -- Nothing to align!\n",

+                __FUNCTION__, __FILE__, __LINE__, iAuxLen1, iAuxLen2);

+        iRetVal = RETURN_UNKNOWN;

+        /* Note: at this stage cannot do   'goto this_is_the_end;' 

+           because would cross initialisation of several variables   */

+        return iRetVal;

+    }

+    par.maxResLen  = iAuxLen1 > iAuxLen2 ? iAuxLen1 : iAuxLen2;

+    const int ciGoodMeasureSeq = 10;

+    par.maxResLen += ciGoodMeasureSeq;

+    par.maxColCnt  = iAuxLen1 + iAuxLen2 + ciGoodMeasureSeq;

+    //par.max_seqid=iFirstCnt+iSecndCnt+3; /* -id     */

+    par.max_seqid=DEFAULT_FILTER;        /* -id     */

+    par.loc=0; par.mact=0;              /* -glob   */

+    par.nseqdis=iFirstCnt+iSecndCnt;   /* -seq    */

+    par.showcons=0;                   /* -nocons */

+    par.showdssp=0;                  /* -nodssp */

+    par.Mgaps=100;                  /* -M      */

+    par.M=2;                       /* -M      */

+    par.pdWg1=pdWeightsL;         /* tree wg */

+    par.pdWg2=pdWeightsR;        /* tree wg */

+    v = 0;                      /* -v0     */

+    /* NOTE: *qali declared globally but only created here,

+       pass in number of sequences to get rid of statically 

+       defined MAXSEQ (FS)

+    */

+    Alignment qali(iFirstCnt+iSecndCnt);

+    HMM q(iFirstCnt+iSecndCnt);

+    HMM t(iFirstCnt+iSecndCnt);

+    //initialize

+    static char const empty[] = "";

+    strcpy(q.file, empty);

+    strcpy(t.file, empty);

+    strcpy(t.sfam, empty);

+    strcpy(q.sfam, empty);

+    strcpy(q.fold, empty);

+    strcpy(t.fold, empty);

+    strcpy(t.cl, empty);

+    strcpy(q.cl, empty);

+#endif

+    

+    

+#ifndef CLUSTALO_NOFILE

+    // Check command line input and default values

+    if (!*par.infile)