+++ /dev/null
-/*
- * dmat.c
- *
- * $Id$
- *
- *****************************************************************************
- *
- * Copyright (c) 2004, Luke Sheneman
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * + Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * + Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in
- * the documentation and/or other materials provided with the
- * distribution.
- * + The names of its contributors may not be used to endorse or promote
- * products derived from this software without specific prior
- * written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- *
- *****************************************************************************
- *
- * Distance matrix parser
- *
- *****************************************************************************
- *
- * AUTHOR:
- *
- * Luke Sheneman
- * sheneman@cs.uidaho.edu
- *
- */
-
-
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <float.h>
-#include <errno.h>
-#include <string.h>
-
-
-
-#include "common.h"
-#include "clearcut.h"
-
-#include "dmat.h"
-
-
-
-
-
-/*
- *
- * NJ_is_alpha() - determine if character is an alphabetic character
- *
- * INPUT:
- * ------
- * c -- character to test
- *
- * RETURN:
- * -------
- * int -- 1 if character is alphabetic (A-Z || a-z)
- * 0 if character is NOT alphabetic
- *
- */
-static inline
-int
-NJ_is_alpha(char c) {
-
- if( (c >= 'A' && c <= 'Z') ||
- (c >= 'a' && c <= 'z') ) {
- return(1);
- } else {
- return(0);
- }
-}
-
-
-
-/*
- *
- * NJ_is_whitespace() - determine if character is a whitespace character
- *
- * INPUT:
- * ------
- * c -- character to test
- *
- * RETURN:
- * -------
- * int -- 1 if character is whitespace (space, tab, CR, LF)
- * 0 if character is NOT whitespace
- *
- */
-static inline
-int
-NJ_is_whitespace(char c) {
-
- if( c == ' ' || /* space */
- c == '\n' || /* newline */
- c == '\r' || /* carriage-return */
- c == '\v' || /* vertical tab */
- c == '\f' || /* form feed */
- c == '\t' ) { /* horizontal tab */
- return(1);
- } else {
- return(0);
- }
-
-}
-
-
-
-
-/*
- *
- * NJ_is_number() - determine if character is a number
- *
- * INPUT:
- * ------
- * c -- character to test
- *
- * RETURN:
- * -------
- * int -- 1 if character is a number (0-9)
- * 0 if character is NOT a number
- *
- */
-static inline
-int
-NJ_is_number(char c) {
- if(c >= '0' && c <= '9') {
- return(1);
- } else {
- return(0);
- }
-}
-
-
-
-/*
- * NJ_is_distance() - check if string is a properly formatted distance value
- *
- */
-static inline
-int
-NJ_is_distance(char *token) {
-
- int i;
- char c;
- int exponent_state;
- int expsign_state;
- int dpoint_state;
-
- /* if token is NULL return failure */
- if(!token) {
- return(0);
- }
-
- exponent_state = 0;
- expsign_state = 0;
- dpoint_state = 0;
-
- /* The first character must be a number, a decimal point or a sign */
- c = token[0];
- if(!NJ_is_number(c) &&
- c != '.' &&
- c != '-' &&
- c != '+' ) {
-
- goto BAD;
- }
-
- /*
- * if the first character is not a number, and string is only one
- * character long, then we return failure.
- */
- if(strlen(token) == 1) {
- if(!NJ_is_number(c)) {
- goto BAD;
- }
- }
-
- for(i=0;i<strlen(token);i++) {
-
- c = token[i];
-
- /* make sure that all chars in dist string are in list of valid chars */
- if(!NJ_is_number(c) &&
- c != '.' &&
- c != '-' &&
- c != '+' &&
- c != 'e' &&
- c != 'E' ) {
-
- goto BAD;
- }
-
- /* not the first char and we are not in exponent state but read (+,-) */
- if(i>0 && !exponent_state) {
- if(c == '-' ||
- c == '+') {
- goto BAD;
- }
- }
-
- /* if we are in the exponent state, and we've already seen a sign */
- if(exponent_state && expsign_state) {
- if(c == '-' ||
- c == '+') {
- goto BAD;
- }
- }
-
- /* if we are in the exponent state and we see a decimal point */
- if(exponent_state) {
- if(c == '.') {
- goto BAD;
- }
- }
-
- /* if we are in the exponent state and see another e or E */
- if(exponent_state) {
- if(c == 'e' ||
- c == 'E') {
- goto BAD;
- }
- }
-
- /* if we are dpoint_state and see another decimal point */
- if(dpoint_state) {
- if(c == '.') {
- goto BAD;
- }
- }
-
-
- /* enter the exponent state if we need to */
- if(!exponent_state) {
- if(c == 'e' ||
- c == 'E') {
- exponent_state = 1;
- }
- }
-
- /* enter the expsign_state if we need to */
- if(exponent_state && !expsign_state) {
- if(c == '-' ||
- c == '+') {
- expsign_state = 1;
- }
- }
-
- /* if not in dpoint state and we see a dpoint */
- if(!dpoint_state) {
- if(c == '.') {
- dpoint_state = 1;
- }
- }
-
- }
-
- /* the token must end in a number char */
- if(!NJ_is_number(token[strlen(token)-1])) {
- goto BAD;
- }
-
- /* token is a valid numerical distance */
- return(1);
-
- BAD:
-
- /* token is invalid distance format */
- return(0);
-}
-
-
-
-
-/*
- * NJ_is_label() -
- *
- * Simply, if token is not a valid number, then it is a name
- *
- */
-static inline
-int
-NJ_is_label(char *token) {
- if(NJ_is_distance(token)) {
- return(0);
- } else {
- return(1);
- }
-}
-
-
-
-/*
- * NJ_get_token() - get a token from an input stream
- *
- */
-static inline
-int
-NJ_get_token(FILE *fp,
- NJ_DIST_TOKEN *token) {
-
- char c;
- int index;
-
- c = fgetc(fp);
- if(feof(fp)) {
- token->type = NJ_EOF_STATE;
- return(token->type);
- }
-
- if(NJ_is_whitespace(c)) {
- token->buf[0] = c;
- token->buf[1] = '\0';
- token->type = NJ_WS_STATE;
-
- return NJ_WS_STATE;
- }
-
- index = 0;
- while(!NJ_is_whitespace(c)) {
-
- /* reallocate our buffer if necessary */
- if(index >= token->bufsize) {
- token->bufsize *= 2;
- token->buf = (char *)realloc(token->buf, token->bufsize*sizeof(char));
- if(!token->buf) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_get_token()\n");
- exit(-1);
- }
- }
-
- token->buf[index++] = c;
-
- c = fgetc(fp);
- if(feof(fp)) {
- token->type = NJ_EOF_STATE;
- break;
- }
- }
-
- token->buf[index] = '\0';
-
- if(token->type != NJ_EOF_STATE) {
-
- if(NJ_is_distance(token->buf)) {
- token->type = NJ_FLOAT_STATE;
- } else {
- token->type = NJ_NAME_STATE;
- }
-
- }
-
- return(token->type);
-}
-
-
-
-
-
-
-/*
- * NJ_parse_distance_matrix() -- Takes a filename and returns a distance matrix
- *
- *
- * INPUT:
- * ------
- * nj_args -- a pointer to a structure containing the command-line arguments
- *
- * OUTPUT:
- * -------
- * <DMAT *> -- NULL (failure)
- * -- A pointer to a populated distance matrix
- *
- * DESCRIPTION:
- * ------------
- * This function implements a simple state machine to parse a distance matrix
- * in approximate PHYLIP format. This function auto-detects whether the
- * distance matrix is in upper, lower, or fully-symmetric format and handles
- * it accordingly. For full/symmetric matrices, values must be symmetric
- * around the diagonal, which is required to be zeroes. Names and values must
- * be separated by whitespace (space, tab, newlines, etc.). Taxon labels can
- * include numbers, but must start with non-numerical symbols.
- *
- *
- * *** UPPER FORMAT EXAMPLE ***
- *
- * 4
- * seq1 0.2 0.3 0.1
- * seq2 0.2 0.3
- * seq3 0.1
- * seq4
- *
- * *** LOWER FORMAT EXAMPLE ***
- *
- * 4
- * seq1
- * seq2 0.3
- * seq3 0.2 0.4
- * seq4 0.3 0.1 0.3
- *
- * *** SYMMETRIC (FULL) EXAMPLE ***
- *
- * 4
- * seq1 0.0 0.3 0.5 0.3
- * seq2 0.3 0.0 0.1 0.2
- * seq3 0.5 0.1 0.0 0.9
- * seq4 0.3 0.2 0.9 0.0
- *
- * Values in the distance matrix can be positive or negative, integers or
- * real values. Values can also be parsed in exponential notation form.
- *
- */
-DMAT *
-NJ_parse_distance_matrix(NJ_ARGS *nj_args) {
-
- DMAT *dmat = NULL;
- FILE *fp = NULL;
- NJ_DIST_TOKEN *token = NULL;
-
- int state, dmat_type;
- int row;
- int fltcnt;
- int x, y, i;
- int numvalread;
- int expectedvalues = -1;
- float val;
- int first_state = 0;
-
-
- /* allocate our distance matrix and token structure */
- dmat = (DMAT *)calloc(1, sizeof(DMAT));
- token = (NJ_DIST_TOKEN *)calloc(1, sizeof(NJ_DIST_TOKEN));
- if(token) {
- token->bufsize = NJ_INITIAL_BUFSIZE;
- token->buf = (char *)calloc(token->bufsize, sizeof(char));
- }
- if(!dmat || !token || !token->buf) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
- goto XIT_BAD;
- }
-
- /* open distance matrix file here */
- if(nj_args->stdin_flag) {
- fp = stdin;
- } else {
- fp = fopen(nj_args->infilename, "r");
- if(fp==NULL) {
- fprintf(stderr, "Clearcut: Could not open distance matrix: %s\n", nj_args->infilename);
- perror("Clearcut");
- goto XIT_BAD;
- }
- }
-
- /* get the number of taxa in this file */
- fscanf(fp, "%ld\n", &dmat->ntaxa);
- if(dmat->ntaxa < 2) {
- fprintf(stderr, "Clearcut: Invalid number of taxa in distance matrix\n");
-
- goto XIT_BAD;
- }
-
- /* set our initial working size according to the # of taxa */
- dmat->size = dmat->ntaxa;
-
- /* allocate space for the distance matrix values here */
- dmat->val =
- (float *)calloc(NJ_NCELLS(dmat->ntaxa), sizeof(float));
- if(!dmat->val) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
- goto XIT_BAD;
- }
-
- /* taxa names */
- dmat->taxaname = (char **)calloc(dmat->ntaxa, sizeof(char *));
- if(!dmat->taxaname) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
- goto XIT_BAD;
- }
-
- /* set the initial state of our state machine */
- dmat_type = NJ_PARSE_UNKNOWN;
- row = -1;
- fltcnt = 0;
- numvalread = 0;
-
-
- /* read the input one character at a time to drive simple state machine */
- state = NJ_get_token(fp, token);
- while(state != NJ_EOF_STATE) {
-
- switch(state) {
-
- case NJ_NAME_STATE:
-
- if(first_state == 0) {
- first_state = 1;
- }
-
- row++;
-
- if(row > 0 && dmat_type == NJ_PARSE_UNKNOWN) {
-
- if(fltcnt == dmat->ntaxa) {
-
- dmat_type = NJ_PARSE_SYMMETRIC;
- expectedvalues = dmat->ntaxa * dmat->ntaxa;
-
- } else if (fltcnt == dmat->ntaxa-1) {
-
- dmat_type = NJ_PARSE_UPPER;
- expectedvalues = ((dmat->ntaxa) * (dmat->ntaxa-1)) / 2;
-
- /* shift everything in first row by one char */
- for(i=dmat->ntaxa-2;i>=0;i--) {
- dmat->val[i+1] = dmat->val[i];
- }
-
- } else if (fltcnt == 0) {
-
- dmat_type = NJ_PARSE_LOWER;
- expectedvalues = ((dmat->ntaxa) * (dmat->ntaxa-1)) / 2;
-
- } else {
- goto XIT_BAD;
- }
- }
-
- if(row >= dmat->ntaxa) {
- goto XIT_BAD;
- }
-
- /* allocate space for this taxon label */
- dmat->taxaname[row] = (char *)calloc(strlen(token->buf)+1, sizeof(char));
- if(!dmat->taxaname[row]) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
- goto XIT_BAD;
- }
-
- strcpy(dmat->taxaname[row], token->buf);
-
- fltcnt = 0;
-
- break;
-
-
- case NJ_FLOAT_STATE:
-
- if(first_state == 0) {
- goto XIT_BAD;
- }
-
- val = atof(token->buf);
- if(errno) {
- fprintf(stderr, "Clearcut: Distance value out-of-range.\n");
- goto XIT_BAD;
- }
-
- x = row;
- y = fltcnt;
-
- switch(dmat_type) {
-
- case NJ_PARSE_UNKNOWN:
-
- dmat->val[NJ_MAP(x, y, dmat->size)] = val;
-
- break;
-
- case NJ_PARSE_SYMMETRIC:
-
- if(fltcnt >= dmat->ntaxa) {
- fprintf(stderr, "Clearcut: Incorrect number of distance values on row.\n");
- goto XIT_BAD;
- }
-
- if(x < y) {
- dmat->val[NJ_MAP(x, y, dmat->size)] = val;
- } else if(x > y) {
- if(!NJ_FLT_EQ(val, dmat->val[NJ_MAP(y, x, dmat->size)])) {
- fprintf(stderr, "Clearcut: Full matrices must be symmetric.\n");
- goto XIT_BAD;
- }
- } else {
- if(!NJ_FLT_EQ(val, 0.0)) {
- fprintf(stderr, "Clearcut: Values along the diagonal in a symmetric matrix must be zero.\n");
- goto XIT_BAD;
-
- }
- }
-
- break;
-
- case NJ_PARSE_UPPER:
-
- if(fltcnt > dmat->ntaxa-row) {
- fprintf(stderr, "Clearcut: Incorrect number of distance values on row.\n");
- goto XIT_BAD;
- }
-
- dmat->val[NJ_MAP(x, x+y+1, dmat->size)] = val;
-
- break;
-
- case NJ_PARSE_LOWER:
-
- if(fltcnt > row-1) {
- fprintf(stderr, "Clearcut: Incorrect number of distance values on row.\n");
- goto XIT_BAD;
- }
-
- dmat->val[NJ_MAP(y, x, dmat->size)] = val;
-
- break;
-
- default:
-
- goto XIT_BAD;
-
- break;
- }
-
- fltcnt++;
- numvalread++;
-
- break;
-
- case NJ_WS_STATE:
-
- break;
-
- case NJ_EOF_STATE:
-
- if(first_state == 0) {
- goto XIT_BAD;
- }
-
- break;
-
- default:
-
- fprintf(stderr, "Clearcut: Unknown state in distance matrix parser.\n");
- break;
-
- }
-
- /* get next token from stream */
- state = NJ_get_token(fp, token);
- }
-
-
- /*
- * At the end, if we have not read the number of values that we predicted
- * we would need, then there was a problem and we need to punt.
- */
- if(numvalread != expectedvalues) {
- fprintf(stderr, "Clearcut: Incorrect number of values in the distance matrix.\n");
- goto XIT_BAD;
- }
-
- /* special check to make sure first value read is 0.0 */
- if(dmat_type == NJ_PARSE_SYMMETRIC) {
- if(!NJ_FLT_EQ(dmat->val[NJ_MAP(0, 0, dmat->size)], 0.0)) {
- fprintf(stderr, "Clearcut: Values along the diagonal in a symmetric matrix must be zero.\n");
- goto XIT_BAD;
- }
- }
-
-
- /* now lets allocate space for the r and r2 columns */
- dmat->r = (float *)calloc(dmat->ntaxa, sizeof(float));
- dmat->r2 = (float *)calloc(dmat->ntaxa, sizeof(float));
- if(!dmat->r || !dmat->r2) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_parse_distance_matrix()\n");
- goto XIT_BAD;
- }
-
- /* track some memory addresses */
- dmat->rhandle = dmat->r;
- dmat->r2handle = dmat->r2;
- dmat->valhandle = dmat->val;
-
- /* close matrix file here */
- if(!nj_args->stdin_flag) {
- fclose(fp);
- }
-
- if(token) {
- if(token->buf) {
- free(token->buf);
- }
- free(token);
- }
-
- return(dmat);
-
-
-
- /* clean up our partial progress */
- XIT_BAD:
-
- if(fp) {
- fprintf(stderr, "Clearcut: Syntax error in distance matrix at offset %ld.\n", ftell(fp));
- }
-
- /* close matrix file here */
- if(!nj_args->stdin_flag) {
- if(fp) {
- fclose(fp);
- }
- }
-
- /* if we have a valid dmat (partial or complete), we need to free it */
- if(dmat) {
- NJ_free_dmat(dmat);
- }
-
- if(token) {
- if(token->buf) {
- free(token->buf);
- }
- free(token);
- }
-
- return(NULL);
-}
-
-
-
-
-
-
-
-/*
- * NJ_output_matrix() - Output a distance matrix to the specified file
- *
- *
- * INPUTS:
- * -------
- * nj_args -- a pointer to a data structure holding the command-line args
- * dmat -- a pointer to a distance matrix
- *
- *
- * RETURNS:
- * --------
- * NOTHING
- *
- *
- * DESCRIPTION:
- * ------------
- * If the appropriate flag was specified in the command-line, this function
- * now outputs the parsed or computed distance matrix to a file. This
- * can be useful if generating a distance matrix was the primary goal of
- * running the program, or if one wanted to debug and/or verify the
- * correctness of the program.
- *
- * Currently this function outputs full/symmetric matrices only.
- *
- */
-void
-NJ_output_matrix(NJ_ARGS *nj_args,
- DMAT *dmat) {
-
- FILE *fp = NULL;
- long int i, j;
-
-
-
- /* if we haven't specieid matrixout, return immediately */
- if(!nj_args->matrixout) {
- return;
- }
-
- /* open the specified matrix file for writing */
- fp = fopen(nj_args->matrixout, "w");
- if(!fp) {
- fprintf(stderr, "Clearcut: Could not open matrix file %s for output.\n", nj_args->matrixout);
- return;
- }
-
- /* output the number of taxa in the matrix */
- fprintf(fp, " %ld\n", dmat->size);
-
- fprintf(fp, "%s\n", dmat->taxaname[0]); // print the first taxon name outside of the main loop
-
- for(i=1;i<dmat->size;i++) {
-
- /* output taxaname */
- fprintf(fp, "%s\t", dmat->taxaname[i]);
-
- for(j=0;j<i;j++) {
- if(nj_args->expdist) { /* exponential notation (or not) */
- fprintf(fp, "%e ", dmat->val[NJ_MAP(j,i,dmat->size)]);
- } else {
- fprintf(fp, "%f ", dmat->val[NJ_MAP(j,i,dmat->size)]);
- }
- }
-
- fprintf(fp, "\n");
- }
-
-#ifdef FULL_SYMMETRIC_MATRIX
-
- /* output the number of taxa in the matrix */
- fprintf(fp, " %ld\n", dmat->size);
- for(i=0;i<dmat->size;i++) {
-
- /* output taxaname */
- fprintf(fp, "%s\t", dmat->taxaname[i]);
-
- for(j=0;j<dmat->size;j++) {
- if(i>j) {
- if(nj_args->expdist) { /* exponential notation (or not) */
- fprintf(fp, "%e ", dmat->val[NJ_MAP(j,i,dmat->size)]);
- } else {
- fprintf(fp, "%f ", dmat->val[NJ_MAP(j,i,dmat->size)]);
- }
- } else if(i<j) {
- if(nj_args->expdist) { /* exponential notation (or not) */
- fprintf(fp, "%e ", dmat->val[NJ_MAP(i,j,dmat->size)]);
- } else {
- fprintf(fp, "%f ", dmat->val[NJ_MAP(i,j,dmat->size)]);
- }
- } else {
- if(nj_args->expdist) { /* exponential notation (or not) */
- fprintf(fp, "%e ", 0.0);
- } else {
- fprintf(fp, "%f ", 0.0);
- }
- }
- }
-
- fprintf(fp, "\n");
- }
-
-#endif // FULL_SYMMETRIC_MATRIX
-
- /* close the file here */
- if(fp) {
- fclose(fp);
- }
-
- return;
-}
-
-
-
-
-
projects / mothur.git / blobdiff