]> git.donarmstrong.com Git - mothur.git/blob - dmat.cpp
working on pam
[mothur.git] / dmat.cpp
1 /*
2  * dmat.c
3  *
4  * $Id$
5  *
6  *****************************************************************************
7  *
8  * Copyright (c) 2004,  Luke Sheneman
9  * All rights reserved.
10  *
11  * Redistribution and use in source and binary forms, with or without 
12  * modification, are permitted provided that the following conditions 
13  * are met:
14  * 
15  *  + Redistributions of source code must retain the above copyright 
16  *    notice, this list of conditions and the following disclaimer. 
17  *  + Redistributions in binary form must reproduce the above copyright 
18  *    notice, this list of conditions and the following disclaimer in 
19  *    the documentation and/or other materials provided with the 
20  *    distribution. 
21  *  + The names of its contributors may not be used to endorse or promote 
22  *    products derived  from this software without specific prior 
23  *    written permission. 
24  *
25  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
26  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
27  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
28  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
29  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
30  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
31  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
32  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
33  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
34  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
35  * POSSIBILITY OF SUCH DAMAGE.  
36  *
37  *****************************************************************************
38  *
39  * Distance matrix parser
40  *
41  *****************************************************************************
42  *
43  * AUTHOR:
44  * 
45  *   Luke Sheneman
46  *   sheneman@cs.uidaho.edu
47  *
48  */
49
50
51
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <float.h>
55 #include <errno.h>
56 #include <string.h>
57
58
59
60 #include "common.h"
61 #include "clearcut.h"
62
63 #include "dmat.h"
64
65
66
67
68
69 /*
70  *
71  * NJ_is_alpha() - determine if character is an alphabetic character
72  *
73  * INPUT:
74  * ------
75  *  c -- character to test
76  *
77  * RETURN:
78  * -------
79  *   int -- 1 if character is alphabetic (A-Z || a-z)
80  *          0 if character is NOT alphabetic
81  *
82  */
83 static inline
84 int 
85 NJ_is_alpha(char c) {
86
87   if( (c >= 'A' && c <= 'Z') ||
88       (c >= 'a' && c <= 'z') ) {
89     return(1);
90   } else {
91     return(0);
92   }
93 }
94
95
96
97 /*
98  *
99  * NJ_is_whitespace() - determine if character is a whitespace character
100  *
101  * INPUT:
102  * ------
103  *  c -- character to test
104  *
105  * RETURN:
106  * -------
107  *   int -- 1 if character is whitespace (space, tab, CR, LF)
108  *          0 if character is NOT whitespace
109  *
110  */
111 static inline
112 int
113 NJ_is_whitespace(char c) {
114
115   if( c == ' '  ||   /* space           */
116       c == '\n' ||   /* newline         */
117       c == '\r' ||   /* carriage-return */
118       c == '\v' ||   /* vertical tab    */
119       c == '\f' ||   /* form feed       */
120       c == '\t' ) {  /* horizontal tab  */
121     return(1);
122   } else {
123     return(0);
124   }
125
126 }
127
128
129
130
131 /*
132  *
133  * NJ_is_number() - determine if character is a number
134  *
135  * INPUT:
136  * ------
137  *  c -- character to test
138  *
139  * RETURN:
140  * -------
141  *   int -- 1 if character is a number (0-9)
142  *          0 if character is NOT a number
143  *
144  */
145 static inline
146 int
147 NJ_is_number(char c) {
148   if(c >= '0' && c <= '9') {
149     return(1);
150   } else {
151     return(0);
152   }
153 }
154
155
156
157 /*
158  * NJ_is_distance() - check if string is a properly formatted distance value
159  *
160  */
161 static inline
162 int
163 NJ_is_distance(char *token) {
164   
165   int i;
166   char c;
167   int exponent_state;
168   int expsign_state;
169   int dpoint_state;
170
171   /* if token is NULL return failure */
172   if(!token) {
173     return(0);
174   }
175   
176   exponent_state = 0;
177   expsign_state  = 0;
178   dpoint_state   = 0;
179   
180   /* The first character must be a number, a decimal point or a sign */
181   c = token[0];
182   if(!NJ_is_number(c) &&
183      c != '.'         &&
184      c != '-'         &&
185      c != '+' )  {
186     
187     goto BAD;
188   } 
189   
190   /* 
191    * if the first character is not a number, and string is only one 
192    * character long, then we return failure.
193    */
194   if(strlen(token) == 1) {
195     if(!NJ_is_number(c)) {
196       goto BAD;
197     }
198   }
199   
200   for(i=0;i<strlen(token);i++) {
201
202     c = token[i];
203     
204     /* make sure that all chars in dist string are in list of valid chars */
205     if(!NJ_is_number(c) &&
206        c != '.'         &&
207        c != '-'         &&
208        c != '+'         &&
209        c != 'e'         &&
210        c != 'E'         ) {
211
212       goto BAD;
213     }
214
215     /* not the first char and we are not in exponent state but read (+,-) */
216     if(i>0 && !exponent_state) {
217       if(c == '-' || 
218          c == '+') {
219         goto BAD;
220       }
221     }
222
223     /* if we are in the exponent state, and we've already seen a sign */
224     if(exponent_state && expsign_state) {
225       if(c == '-' ||
226          c == '+') {
227         goto BAD;
228       }
229     }
230     
231     /* if we are in the exponent state and we see a decimal point */
232     if(exponent_state) {
233       if(c == '.') {
234         goto BAD;
235       }
236     }
237     
238     /* if we are in the exponent state and see another e or E */
239     if(exponent_state) {
240       if(c == 'e' ||
241          c == 'E') {
242         goto BAD;
243       }
244     }
245     
246     /* if we are dpoint_state and see another decimal point */
247     if(dpoint_state) {
248       if(c == '.') {
249         goto BAD;
250       }
251     }
252     
253     
254     /* enter the exponent state if we need to */
255     if(!exponent_state) {
256       if(c == 'e' ||
257          c == 'E') {
258         exponent_state = 1;
259       }
260     }
261
262     /* enter the expsign_state if we need to */
263     if(exponent_state && !expsign_state) {
264       if(c == '-' ||
265          c == '+') {
266         expsign_state = 1;
267       }
268     }
269
270     /* if not in dpoint state and we see a dpoint */
271     if(!dpoint_state) {
272       if(c == '.') {
273         dpoint_state = 1;
274       }
275     }
276
277   }
278   
279   /* the token must end in a number char */
280   if(!NJ_is_number(token[strlen(token)-1])) {
281     goto BAD;
282   }
283   
284   /* token is a valid numerical distance */
285   return(1);
286   
287  BAD:
288
289   /* token is invalid distance format */
290   return(0);
291 }
292
293
294
295
296 /*
297  * NJ_is_label() - 
298  *
299  * Simply, if token is not a valid number, then it is a name
300  *
301  */
302 static inline
303 int
304 NJ_is_label(char *token) {
305   if(NJ_is_distance(token)) {
306     return(0);
307   } else {
308     return(1);
309   }
310 }
311
312
313
314 /*
315  * NJ_get_token() - get a token from an input stream 
316  *
317  */
318 static inline
319 int
320 NJ_get_token(FILE *fp,
321              NJ_DIST_TOKEN *token) {
322
323   char c;
324   int index;
325
326   c = fgetc(fp);
327   if(feof(fp)) {
328     token->type = NJ_EOF_STATE;
329     return(token->type);
330   }
331
332   if(NJ_is_whitespace(c)) {
333     token->buf[0] = c;
334     token->buf[1] = '\0';
335     token->type   = NJ_WS_STATE;
336
337     return NJ_WS_STATE;
338   } 
339
340   index = 0;
341   while(!NJ_is_whitespace(c)) {
342
343     /* reallocate our buffer if necessary */
344     if(index >= token->bufsize) {
345       token->bufsize *= 2;
346       token->buf = (char *)realloc(token->buf, token->bufsize*sizeof(char));
347       if(!token->buf) {
348         fprintf(stderr, "Clearcut: Memory allocation error in NJ_get_token()\n");
349         exit(-1);
350       }
351     }
352
353     token->buf[index++] = c;
354     
355     c = fgetc(fp);
356     if(feof(fp)) {
357       token->type = NJ_EOF_STATE;
358       break;
359     }
360   }
361   
362   token->buf[index] = '\0';
363   
364   if(token->type != NJ_EOF_STATE) {
365
366     if(NJ_is_distance(token->buf)) {
367       token->type = NJ_FLOAT_STATE;
368     } else {
369       token->type = NJ_NAME_STATE;
370     }
371
372   }
373   
374   return(token->type);
375 }
376
377
378
379
380  
381
382 /* 
383  * NJ_parse_distance_matrix() -- Takes a filename and returns a distance matrix
384  *
385  *
386  * INPUT:
387  * ------
388  *   nj_args -- a pointer to a structure containing the command-line arguments
389  *
390  * OUTPUT:
391  * -------
392  *   <DMAT *> -- NULL  (failure)
393  *            -- A pointer to a populated distance matrix
394  *
395  * DESCRIPTION:
396  * ------------
397  *   This function implements a simple state machine to parse a distance matrix
398  *   in approximate PHYLIP format.  This function auto-detects whether the 
399  *   distance matrix is in upper, lower, or fully-symmetric format and handles
400  *   it accordingly.  For full/symmetric matrices, values must be symmetric 
401  *   around the diagonal, which is required to be zeroes.  Names and values must
402  *   be separated by whitespace (space, tab, newlines, etc.).  Taxon labels can
403  *   include numbers, but must start with non-numerical symbols.
404  *
405  *
406  *   *** UPPER FORMAT EXAMPLE ***
407  * 
408  *   4
409  *   seq1 0.2 0.3 0.1
410  *   seq2     0.2 0.3
411  *   seq3         0.1
412  *   seq4
413  * 
414  *   *** LOWER FORMAT EXAMPLE ***
415  *
416  *   4
417  *   seq1
418  *   seq2 0.3
419  *   seq3 0.2 0.4
420  *   seq4 0.3 0.1 0.3
421  *
422  *   *** SYMMETRIC (FULL) EXAMPLE ***
423  *  
424  *   4
425  *   seq1 0.0 0.3 0.5 0.3
426  *   seq2 0.3 0.0 0.1 0.2
427  *   seq3 0.5 0.1 0.0 0.9
428  *   seq4 0.3 0.2 0.9 0.0
429  *
430  *  Values in the distance matrix can be positive or negative, integers or
431  *  real values.  Values can also be parsed in exponential notation form.
432  * 
433  */
434 DMAT *
435 NJ_parse_distance_matrix(NJ_ARGS *nj_args) {
436   
437   DMAT *dmat           = NULL;
438   FILE *fp            = NULL;
439   NJ_DIST_TOKEN *token = NULL;
440
441   int state, dmat_type;
442   int row;
443   int fltcnt;
444   int x, y, i;
445   int numvalread;
446   int expectedvalues = -1;
447   float val;
448   int first_state = 0;
449
450
451   /* allocate our distance matrix and token structure */
452   dmat = (DMAT *)calloc(1, sizeof(DMAT));
453   token = (NJ_DIST_TOKEN *)calloc(1, sizeof(NJ_DIST_TOKEN));
454   if(token) {
455     token->bufsize = NJ_INITIAL_BUFSIZE;
456     token->buf     = (char *)calloc(token->bufsize, sizeof(char));
457   }
458   if(!dmat || !token || !token->buf) {
459     fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
460     goto XIT_BAD;
461   }
462   
463   /* open distance matrix file here */
464   if(nj_args->stdin_flag) {
465     fp = stdin;
466   } else {
467     fp = fopen(nj_args->infilename, "r");
468     if(fp==NULL) {
469       fprintf(stderr, "Clearcut: Could not open distance matrix: %s\n", nj_args->infilename);
470       perror("Clearcut");
471       goto XIT_BAD;
472     }
473   }
474
475   /* get the number of taxa in this file */
476   fscanf(fp, "%ld\n", &dmat->ntaxa);
477   if(dmat->ntaxa < 2) {
478     fprintf(stderr, "Clearcut: Invalid number of taxa in distance matrix\n");
479
480     goto XIT_BAD;
481   }
482
483   /* set our initial working size according to the # of taxa */
484   dmat->size = dmat->ntaxa;
485
486   /* allocate space for the distance matrix values here */
487   dmat->val = 
488     (float *)calloc(NJ_NCELLS(dmat->ntaxa), sizeof(float));
489   if(!dmat->val) {
490     fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
491     goto XIT_BAD;
492   }
493
494   /*  taxa names */
495   dmat->taxaname = (char **)calloc(dmat->ntaxa, sizeof(char *));
496   if(!dmat->taxaname) {
497     fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
498     goto XIT_BAD;
499   }
500
501   /* set the initial state of our state machine */
502   dmat_type   = NJ_PARSE_UNKNOWN;
503   row         = -1;
504   fltcnt      = 0;
505   numvalread  = 0;
506
507
508   /* read the input one character at a time to drive simple state machine */
509   state = NJ_get_token(fp, token);
510   while(state != NJ_EOF_STATE) {
511     
512     switch(state) {
513
514     case NJ_NAME_STATE:
515       
516       if(first_state == 0) {
517         first_state = 1;
518       }
519
520       row++;
521
522       if(row > 0 && dmat_type == NJ_PARSE_UNKNOWN) {
523
524         if(fltcnt == dmat->ntaxa) {
525
526           dmat_type = NJ_PARSE_SYMMETRIC;
527           expectedvalues = dmat->ntaxa * dmat->ntaxa;
528
529         } else if (fltcnt == dmat->ntaxa-1) {
530
531           dmat_type = NJ_PARSE_UPPER;
532           expectedvalues = ((dmat->ntaxa) * (dmat->ntaxa-1)) / 2;
533
534           /* shift everything in first row by one char */
535           for(i=dmat->ntaxa-2;i>=0;i--) {
536             dmat->val[i+1] = dmat->val[i];
537           }
538
539         } else if (fltcnt == 0) {
540
541           dmat_type = NJ_PARSE_LOWER;
542           expectedvalues = ((dmat->ntaxa) * (dmat->ntaxa-1)) / 2;
543
544         } else {
545           goto XIT_BAD;
546         }
547       }
548       
549       if(row >= dmat->ntaxa) {
550         goto XIT_BAD;
551       }
552       
553       /* allocate space for this taxon label */
554       dmat->taxaname[row] = (char *)calloc(strlen(token->buf)+1, sizeof(char));
555       if(!dmat->taxaname[row]) {
556         fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
557         goto XIT_BAD;
558       }
559       
560       strcpy(dmat->taxaname[row], token->buf); 
561       
562       fltcnt = 0;
563
564       break;
565
566
567     case NJ_FLOAT_STATE:
568
569       if(first_state == 0) {
570         goto XIT_BAD;
571       }
572
573       val = atof(token->buf);
574       if(errno) {
575         fprintf(stderr, "Clearcut: Distance value out-of-range.\n");
576         goto XIT_BAD;
577       }
578       
579       x = row;
580       y = fltcnt;
581
582       switch(dmat_type) {
583
584       case NJ_PARSE_UNKNOWN:
585
586         dmat->val[NJ_MAP(x, y, dmat->size)] = val;
587
588         break;
589
590       case NJ_PARSE_SYMMETRIC:
591         
592         if(fltcnt >= dmat->ntaxa) {
593           fprintf(stderr, "Clearcut: Incorrect number of distance values on row.\n");
594           goto XIT_BAD;
595         }
596
597         if(x < y) {
598           dmat->val[NJ_MAP(x, y, dmat->size)] = val;
599         } else if(x > y) {
600           if(!NJ_FLT_EQ(val, dmat->val[NJ_MAP(y, x, dmat->size)])) {
601             fprintf(stderr, "Clearcut: Full matrices must be symmetric.\n");
602             goto XIT_BAD;
603           }
604         } else {
605           if(!NJ_FLT_EQ(val, 0.0)) {
606             fprintf(stderr, "Clearcut: Values along the diagonal in a symmetric matrix must be zero.\n");
607             goto XIT_BAD;
608
609             }
610         }
611
612         break;
613         
614       case NJ_PARSE_UPPER:
615
616         if(fltcnt > dmat->ntaxa-row) {
617           fprintf(stderr, "Clearcut: Incorrect number of distance values on row.\n");
618           goto XIT_BAD;
619         }
620         
621         dmat->val[NJ_MAP(x, x+y+1, dmat->size)] = val;
622         
623         break;
624
625       case NJ_PARSE_LOWER:
626         
627         if(fltcnt > row-1) {
628           fprintf(stderr, "Clearcut: Incorrect number of distance values on row.\n");
629           goto XIT_BAD;
630         }
631
632         dmat->val[NJ_MAP(y, x, dmat->size)] = val;
633         
634         break;
635         
636       default:
637         
638         goto XIT_BAD;
639         
640         break;
641       }
642
643       fltcnt++;
644       numvalread++;
645       
646       break;
647
648     case NJ_WS_STATE:
649
650       break;
651
652     case NJ_EOF_STATE:
653
654       if(first_state == 0) {
655         goto XIT_BAD;
656       }
657
658       break;
659
660     default:
661
662       fprintf(stderr, "Clearcut: Unknown state in distance matrix parser.\n");
663       break;
664
665     }
666
667     /* get next token from stream */
668     state = NJ_get_token(fp, token);
669   }
670
671
672   /* 
673    * At the end, if we have not read the number of values that we predicted
674    * we would need, then there was a problem and we need to punt.
675    */
676   if(numvalread != expectedvalues) {
677     fprintf(stderr, "Clearcut: Incorrect number of values in the distance matrix.\n");
678     goto XIT_BAD;
679   }
680   
681   /* special check to make sure first value read is 0.0 */
682   if(dmat_type == NJ_PARSE_SYMMETRIC) {
683     if(!NJ_FLT_EQ(dmat->val[NJ_MAP(0, 0, dmat->size)], 0.0)) {
684       fprintf(stderr, "Clearcut: Values along the diagonal in a symmetric matrix must be zero.\n");
685       goto XIT_BAD;
686     }
687   }
688
689   
690   /* now lets allocate space for the r and r2 columns */
691   dmat->r  = (float *)calloc(dmat->ntaxa, sizeof(float));
692   dmat->r2 = (float *)calloc(dmat->ntaxa, sizeof(float));
693   if(!dmat->r || !dmat->r2) {
694     fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
695     goto XIT_BAD;
696   }
697   
698   /* track some memory addresses */
699   dmat->rhandle   = dmat->r;
700   dmat->r2handle  = dmat->r2;
701   dmat->valhandle = dmat->val;
702
703   /* close matrix file here */
704   if(!nj_args->stdin_flag) {
705     fclose(fp);
706   }
707   
708   if(token) {
709     if(token->buf) {
710       free(token->buf);
711     }
712     free(token);
713   }
714
715   return(dmat);
716
717
718
719   /* clean up our partial progress */
720  XIT_BAD:
721
722   if(fp) {
723     fprintf(stderr, "Clearcut: Syntax error in distance matrix at offset %ld.\n", ftell(fp));
724   }
725
726   /* close matrix file here */
727   if(!nj_args->stdin_flag) {
728     if(fp) {
729       fclose(fp);
730     }
731   }
732
733   /* if we have a valid dmat (partial or complete), we need to free it */
734   if(dmat) {
735     NJ_free_dmat(dmat);
736   }
737   
738   if(token) {
739     if(token->buf) {
740       free(token->buf);
741     }
742     free(token);
743   }
744   
745   return(NULL);
746 }
747
748
749
750
751
752
753
754 /*
755  * NJ_output_matrix() - Output a distance matrix to the specified file 
756  *
757  *
758  * INPUTS:
759  * -------
760  *  nj_args -- a pointer to a data structure holding the command-line args
761  *     dmat -- a pointer to a distance matrix
762  *
763  *
764  * RETURNS:
765  * --------
766  *   NOTHING
767  *
768  *
769  * DESCRIPTION:
770  * ------------
771  *   If the appropriate flag was specified in the command-line, this function
772  *   now outputs the parsed or computed distance matrix to a file.  This 
773  *   can be useful if generating a distance matrix was the primary goal of 
774  *   running the program, or if one wanted to debug and/or verify the
775  *   correctness of the program.
776  *
777  *   Currently this function outputs full/symmetric matrices only.
778  *
779  */
780 void
781 NJ_output_matrix(NJ_ARGS *nj_args,
782                  DMAT *dmat) {
783   
784   FILE *fp = NULL;
785   long int i, j;
786
787
788   
789   /* if we haven't specieid matrixout, return immediately */
790   if(!nj_args->matrixout) {
791     return;
792   }
793   
794   /* open the specified matrix file for writing */
795   fp = fopen(nj_args->matrixout, "w");
796   if(!fp) {
797     fprintf(stderr, "Clearcut: Could not open matrix file %s for output.\n", nj_args->matrixout);
798     return;
799   }
800
801   /* output the number of taxa in the matrix */
802   fprintf(fp, "   %ld\n", dmat->size);
803
804   fprintf(fp, "%s\n", dmat->taxaname[0]); // print the first taxon name outside of the main loop
805
806   for(i=1;i<dmat->size;i++) {
807     
808     /* output taxaname */
809     fprintf(fp, "%s\t", dmat->taxaname[i]);
810
811     for(j=0;j<i;j++) {
812       if(nj_args->expdist) {  /* exponential notation (or not) */
813         fprintf(fp, "%e ", dmat->val[NJ_MAP(j,i,dmat->size)]);  
814       } else {
815         fprintf(fp, "%f ", dmat->val[NJ_MAP(j,i,dmat->size)]);
816       }
817     }
818     
819     fprintf(fp, "\n");
820   }
821
822 #ifdef FULL_SYMMETRIC_MATRIX 
823
824   /* output the number of taxa in the matrix */
825   fprintf(fp, "   %ld\n", dmat->size);
826   for(i=0;i<dmat->size;i++) {
827     
828     /* output taxaname */
829     fprintf(fp, "%s\t", dmat->taxaname[i]);
830
831     for(j=0;j<dmat->size;j++) {
832       if(i>j) {
833         if(nj_args->expdist) {  /* exponential notation (or not) */
834           fprintf(fp, "%e ", dmat->val[NJ_MAP(j,i,dmat->size)]);  
835         } else {
836           fprintf(fp, "%f ", dmat->val[NJ_MAP(j,i,dmat->size)]);
837         }
838       } else if(i<j) {
839         if(nj_args->expdist) {  /* exponential notation (or not) */
840           fprintf(fp, "%e ", dmat->val[NJ_MAP(i,j,dmat->size)]);
841         } else {
842           fprintf(fp, "%f ", dmat->val[NJ_MAP(i,j,dmat->size)]);
843         }
844       } else {
845         if(nj_args->expdist) {  /* exponential notation (or not) */
846           fprintf(fp, "%e ", 0.0);
847         } else {
848           fprintf(fp, "%f ", 0.0);
849         }
850       }
851     }
852     
853     fprintf(fp, "\n");
854   }
855
856 #endif // FULL_SYMMETRIC_MATRIX
857   
858   /* close the file here */
859   if(fp) {
860     fclose(fp);
861   }
862   
863   return;
864 }
865
866
867
868
869