Add GLprobs and MSAprobs to binaries
[jabaws.git] / binaries / src / MSAProbs-0.9.7 / MSAProbs / MSAReadMatrix.cpp
diff --git a/binaries/src/MSAProbs-0.9.7/MSAProbs/MSAReadMatrix.cpp b/binaries/src/MSAProbs-0.9.7/MSAProbs/MSAReadMatrix.cpp
new file mode 100644 (file)
index 0000000..6ff0643
--- /dev/null
@@ -0,0 +1,174 @@
+/***********************************************
+ * # Copyright 2009-2010. Liu Yongchao
+ * # Contact: Liu Yongchao, School of Computer Engineering,
+ * #                    Nanyang Technological University.
+ * # Emails:    liuy0039@ntu.edu.sg; nkcslyc@hotmail.com
+ * #
+ * # GPL version 3.0 applies.
+ * #
+ * ************************************************/
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "MSAReadMatrix.h"
+
+#define TRACE 0
+
+////////////////////////////////////////////////////////////
+// extern variables for scoring matrix data
+////////////////////////////////////////////////////////////
+extern float g_gap_open1, g_gap_open2, g_gap_ext1, g_gap_ext2;
+extern char *aminos, *bases, matrixtype[20];
+extern int subst_index[26];
+
+extern double sub_matrix[26][26];
+
+extern float TEMPERATURE;
+extern int MATRIXTYPE;
+
+extern float GAPOPEN;
+extern float GAPEXT;
+
+typedef struct {
+       char input[30];
+       int matrix;
+       int N;
+       float T;
+       float beta;
+       char opt;                       //can be 'P' or 'M'
+       float gapopen;
+       float gapext;
+} argument_decl;
+
+//argument support
+extern argument_decl argument;
+
+/////////////////////////////////////////////////////////
+//sets substitution matrix type
+////////////////////////////////////////////////////////
+void setmatrixtype(int le) {
+       switch (le) {
+       case 160:
+               strcpy(matrixtype, "gonnet_160");
+               break;
+       case 4:
+               strcpy(matrixtype, "nuc_simple");
+               break;
+       default:
+               strcpy(matrixtype, "CUSTOM");
+               break;
+
+       };
+
+}
+
+///////////////////////////////////////////////////////////////////
+//sets matrix flag
+///////////////////////////////////////////////////////////////////
+inline int matrixtype_to_int() {
+
+       if (!strcmp(matrixtype, "nuc_simple"))
+               return 4;
+       else if (!strcmp(matrixtype, "gonnet_160"))
+               return 160;
+       else
+               return 1000;
+
+}
+
+/////////////////////////////////////////////////////////////////
+//
+// Can read any scoring matrix as long as it is defined in Matrix.h
+// AND it is a lower triangular 
+// AND the order of amino acids/bases is mentioned
+/////////////////////////////////////////////////////////////////
+
+inline void read_matrix(score_matrix matrx) {
+       int i, j, basecount, position = 0;
+
+       bases = (char *) matrx.monomers;
+
+       basecount = strlen(bases);
+
+       for (i = 0; i < basecount; i++)
+               subst_index[i] = -1;
+
+       for (i = 0; i < basecount; i++)
+               subst_index[bases[i] - 'A'] = i;
+
+       if (TRACE == 1)
+               printf("\nbases read: %d\n", basecount);
+
+       for (i = 0; i < basecount; i++)
+               for (j = 0; j <= i; j++) {
+
+                       double value = exp(argument.beta * matrx.matrix[position++]);
+                       sub_matrix[i][j] = value;
+                       sub_matrix[j][i] = value;
+               }
+
+       if (TRACE)
+               for (i = 0; i < basecount; i++) {
+                       for (j = 0; j < basecount; j++)
+                               printf(" %g ", sub_matrix[i][j]);
+                       printf("\n");
+               }
+
+}
+
+////////////////////////////////////////////////////////////////////////////////// 
+//intialize the arguments (default values)
+//////////////////////////////////////////////////////////////////////////////////
+void init_arguments() {
+       float gap_open = 0, gap_ext = 0;
+       int le;
+
+       le = matrixtype_to_int();
+
+       argument.N = 1;
+       strcpy(argument.input, "tempin");
+       argument.matrix = le;
+       argument.gapopen = GAPOPEN;
+       argument.gapext = GAPEXT;
+       argument.T = TEMPERATURE;
+       argument.beta = 1.0 / TEMPERATURE;
+       argument.opt = 'P';
+
+       if (le == 4)            //NUC OPTION :default is nuc_simple
+                       {
+               read_matrix(nuc_simple);
+               gap_open = -4;
+               gap_ext = -0.25;
+       }
+
+       else if (le == 160)  //PROT option: default is gonnet_160
+                       {
+               if (TRACE)
+                       printf("read matrix\n");
+               read_matrix(gonnet_160);
+               gap_open = -22;
+               gap_ext = -1;
+       } else if (le == 1000) {  //Error handling
+               printf("Error: enter a valid matrix type\n");
+               exit(1);
+               //additional matrices can only be lower triangular
+       }
+
+       //now override the gapopen and gapext
+       if (argument.gapopen != 0.0 || argument.gapext != 0.00)
+
+       {
+               gap_open = -argument.gapopen;
+               gap_ext = -argument.gapext;
+       }
+
+       if (TRACE)
+               printf("%f %f %f %d\n", argument.T, gap_open, gap_ext, le);
+
+       argument.gapopen = gap_open;
+       argument.gapext = gap_ext;
+       argument.opt = 'P';
+
+}