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