binaries/src/clustalo/src/clustal/symmatrix.c

   1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4;  indent-tabs-mode: nil -*- */
   2
   3 /*********************************************************************
   4  * Clustal Omega - Multiple sequence alignment
   5  *
   6  * Copyright (C) 2010 University College Dublin
   7  *
   8  * Clustal-Omega is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This file is part of Clustal-Omega.
  14  *
  15  ********************************************************************/
  16
  17 /*
  18  *  RCS $Id: symmatrix.c 230 2011-04-09 15:37:50Z andreas $
  19  *
  20  *
  21  * Functions for symmetric (square) matrices including diagonal.
  22  * supports the notion of non-square sub-matrices of a symmetric
  23  * matrix, i.e. where |rows|<|cols|.
  24  *
  25  * FIXME Allocating one big chunk of memory is probably
  26  * much faster and also easier to maintain.
  27  */
  28
  29 #ifdef HAVE_CONFIG_H
  30 #include "config.h"
  31 #endif
  32
  33 #include <stdlib.h>
  34 #include <stdio.h>
  35 #include <ctype.h>
  36 #include <string.h>
  37 #include <assert.h>
  38 #include "symmatrix.h"
  39
  40
  41 #if 0
  42 #define DEBUG
  43 #endif
  44
  45 #define MAX_BUF_SIZE 65536
  46
  47 /**
  48  * @brief Allocates symmat and its members and initialises them. Data
  49  * will be calloced, i.e. initialised with zeros.
  50  *
  51  * @param[out] symmat
  52  * newly allocated and initialised symmatrix instance
  53  * @param[in] nrows
  54  * number of rows
  55  * @param[in]
  56  * ncols number of columns
  57  *
  58  * @return: non-zero on error
  59  *
  60  * @see FreeSymMatrix()
  61  *
  62  * @note: symmat data will be of fake shape nrows x ncols
  63  *
  64  */
  65 int
  66 NewSymMatrix(symmatrix_t **symmat, int nrows, int ncols)
  67 {
  68     int i; /* aux */
  69
  70     assert(nrows>0 && ncols>0 && ncols>=nrows);
  71     assert(ncols>0 && ncols>=nrows);
  72
  73     (*symmat) = (symmatrix_t *) malloc(1*sizeof(symmatrix_t));
  74     if (NULL == (*symmat)) {
  75         fprintf(stderr, "Couldn't allocate memory (%s|%s)\n",
  76                 __FILE__, __FUNCTION__);
  77         return -1;
  78     }
  79
  80     (*symmat)->data = (double **) malloc(nrows * sizeof(double *));
  81     if (NULL == (*symmat)->data) {
  82         fprintf(stderr, "Couldn't allocate memory (%s|%s)\n",
  83                 __FILE__, __FUNCTION__);
  84         free(*symmat);
  85         *symmat = NULL;
  86         return -1;
  87     }
  88     for (i=0; i<nrows; i++) {
  89         (*symmat)->data[i] = (double *) calloc((ncols-i), sizeof(double));
  90         if (NULL == (*symmat)->data[i]) {
  91             fprintf(stderr, "Couldn't allocate memory (%s|%s)\n",
  92                     __FILE__, __FUNCTION__);
  93             while (0!=--i) {
  94                 free((*symmat)->data[i]);
  95             }
  96             free((*symmat)->data);
  97             free(*symmat);
  98             *symmat = NULL;
  99             return -1;
 100         }
 101 #ifdef TRACE
 102         fprintf(stderr, "DEBUG(%s|%s():%d) initialising symmat with the number of the beast\n",
 103                 __FILE__, __FUNCTION__, __LINE__);
 104         {
 105             int j;
 106             for (j=0; j<ncols-i; j++) {
 107                 (*symmat)->data[i][j] = -666.0;
 108             }
 109
 110         }
 111 #endif
 112     }
 113
 114     (*symmat)->nrows = nrows;
 115     (*symmat)->ncols = ncols;
 116
 117     return 0;
 118 }
 119 /***   end: new_symmatrix   ***/
 120
 121
 122 /**
 123  * @brief Sets symmat data of given index to given value
 124  *
 125  * @param[in] symmat
 126  * symmatrix_t whose data is to be set
 127  * @param[in] i
 128  * first index
 129  * @param[in] j
 130  * second index
 131  * @param[in] value
 132  * value used to set data point
 133  *
 134  * @see SymMatrixGetValue()
 135  *
 136  * @note This is a convenience function that checks index order.
 137  *
 138  */
 139 void
 140 SymMatrixSetValue(symmatrix_t *symmat, const int i, const int j, const double value)
 141 {
 142     assert(NULL != symmat);
 143
 144     if (i<=j) {
 145         assert(i < symmat->nrows && j < symmat->ncols);
 146         symmat->data[i][j-i] = value;
 147     } else {
 148         assert(j < symmat->nrows && i < symmat->ncols);
 149         symmat->data[j][i-j] = value;
 150     }
 151 }
 152 /***   end: symmatrix_setvalue   ***/
 153
 154
 155
 156 /**
 157  * @brief Returns element of symmat corresponding to given indices
 158  *
 159  * @param[in] symmat
 160  * symmatrix_t of question
 161  * @param[in] i
 162  * index i
 163  * @param[in] j
 164  * index j
 165  *
 166  * @return requested value
 167  *
 168  * @see SymMatrixSetValue()
 169  *
 170  * @note This is a convenience function that checks index order.
 171  */
 172 double
 173 SymMatrixGetValue(symmatrix_t *symmat, const int i, const int j)
 174 {
 175     assert(NULL != symmat);
 176
 177     if (i<=j) {
 178         assert(i < symmat->nrows && j < symmat->ncols);
 179         return symmat->data[i][j-i];
 180     } else {
 181         assert(j < symmat->nrows && i < symmat->ncols);
 182         return symmat->data[j][i-j];
 183     }
 184 }
 185 /***   end: symmatrix_getvalue   ***/
 186
 187
 188
 189
 190
 191 /**
 192  * @brief Returns a pointer to an element of symmat corresponding to
 193  * given indices
 194  *
 195  * @param[out] val
 196  * Value to be set
 197  * @param[in] symmat
 198  * symmatrix_t of question
 199  * @param[in] i
 200  * index i
 201  * @param[in] j
 202  * index j
 203  *
 204  * @return pointer to value
 205  *
 206  * @see SymMatrixGetValue()
 207  *
 208  * @note This is a convenience function that checks index order.
 209  *
 210  */
 211 void
 212 SymMatrixGetValueP(double **val, symmatrix_t *symmat,
 213                    const int i, const int j)
 214 {
 215     assert(NULL != symmat);
 216
 217     if (i<=j) {
 218         assert(i < symmat->nrows && j < symmat->ncols);
 219         *val = &(symmat->data[i][j-i]);
 220     } else {
 221         assert(j < symmat->nrows && i < symmat->ncols);
 222         *val = &(symmat->data[j][i-j]);
 223     }
 224 }
 225 /***   end: symmatrix_getvaluep   ***/
 226
 227
 228 /**
 229  * @brief Frees memory allocated by data members of symmat and symmat
 230  * itself.
 231  *
 232  * @param[in] symmat
 233  * symmatrix_t to be freed
 234  *
 235  * @note Use in conjunction with NewSymMatrix()
 236  *
 237  * @see NewSymMatrix()
 238  *
 239  */
 240 void
 241 FreeSymMatrix(symmatrix_t **symmat)
 242 {
 243     int i;
 244     if (NULL != (*symmat)) {
 245         if (NULL != (*symmat)->data) {
 246             for (i=0; i<(*symmat)->nrows; i++) {
 247                 free((*symmat)->data[i]);
 248             }
 249             free((*symmat)->data);
 250         }
 251     }
 252     free(*symmat);
 253     *symmat = NULL;
 254 }
 255 /***   end: free_symmatrix   ***/
 256
 257
 258
 259 /**
 260  * @brief Print out a symmat in phylip style. Distances are printed on
 261  * one line per sequence/object. Since we also support matrices with
 262  * rows\<cols, the first line can also be nrows by ncolumns instead
 263  * of just nrows
 264  *
 265  * @param[in] symmat
 266  * the symmatrix_t to print
 267  * @param[in] labels
 268  * sequence/objects labels/names. must be at least of length symmat nrows
 269  * @param[in] path
 270  * filename or NULL. If NULL stdout will be used.
 271  *
 272  */
 273 void
 274 SymMatrixPrint(symmatrix_t *symmat, char **labels,  const char *path)
 275 {
 276     FILE *fp = NULL;
 277     int i, j; /* aux */
 278     int max_label_len = 0;
 279
 280     if (NULL==symmat || NULL==labels) {
 281         fprintf(stderr,
 282                 "One of the provided arguments is empty or NULL (print_matrix)\n");
 283         return;
 284     }
 285     if (NULL==path) {
 286         fp = stdout;
 287     } else {
 288         if (NULL==(fp=fopen(path, "w"))) {
 289             fprintf(stderr, "Couldn't open %s for writing.", path);
 290             return;
 291         }
 292     }
 293
 294     /* find maximum label length for pretty printing
 295      */
 296     for (i=0; i<symmat->nrows; i++) {
 297         int this_len;
 298         assert(NULL != labels[i]);
 299         this_len = strlen(labels[i]);
 300         if (this_len>max_label_len)
 301             max_label_len = this_len;
 302     }
 303
 304     if (symmat->ncols==symmat->nrows) {
 305         fprintf(fp, "%u\n", symmat->ncols);
 306     } else {
 307         /* this is not standard phylip, but necessary to indicate
 308            seed matrices */
 309         fprintf(fp, "%u x %u\n", symmat->nrows, symmat->ncols);
 310     }
 311     for (i=0; i<symmat->nrows; i++) {
 312         /* phylip restriction is 10 characters. we don't care here
 313          */
 314         fprintf(fp, "%-*s", max_label_len, labels[i]);
 315         /* we are lazy and do it like fastdist: write all distances
 316          *  for one seq to one line
 317          */
 318         for (j=0; j<symmat->ncols; j++) {
 319             fprintf(fp, " %f", SymMatrixGetValue(symmat, i, j));
 320         }
 321         fprintf(fp, "%s", "\n");
 322     }
 323     if (NULL != path) {
 324         (void) fclose(fp);
 325     } else {
 326         (void) fflush(fp);
 327     }
 328 }
 329 /***   end: SymMatrixPrint   ***/
 330
 331
 332
 333 /**
 334  *
 335  * @brief Read a distance matrix in phylip format
 336  *
 337  * @param[in] pcFileIn
 338  * distance matrix filename
 339  * @param[out] prSymMat_p
 340  * the symmatrix_t. will be allocated here.
 341  * @return:
 342  * non-zero on error
 343  *
 344  * @note:
 345  * FIXME untested code
 346  */
 347 int
 348 SymMatrixRead(char *pcFileIn, symmatrix_t **prSymMat_p)
 349 {
 350     FILE *prFilePointer;
 351     char *buf;
 352     /* number of parsed sequences */
 353     int iNParsedSeq = 0;
 354     /* number of parsed distances per sequence */
 355     int iNParsedDists = 0;
 356     /* total number of sequences/objects */
 357     int iTotalNSeq = 0;
 358     int iRetCode = 0;
 359
 360     assert(NULL!=pcFileIn);
 361
 362     /* FIXME: test correctness and implement label checking */
 363     fprintf(stderr, "WARNING: Reading of distance matrix from file not thoroughly tested!\n");
 364     fprintf(stderr, "WARNING: Assuming same order of sequences in sequence file and distance matrix file (matching of labels not implemented)\n");
 365
 366     if (NULL == (buf = (char *) malloc(MAX_BUF_SIZE * sizeof(char)))) {
 367         fprintf(stderr, "ERROR: couldn't allocate memory at %s:%s:%d\n",
 368                 __FILE__, __FUNCTION__, __LINE__);
 369         return -1;
 370     }
 371
 372     if (NULL == (prFilePointer = fopen(pcFileIn, "r"))) {
 373         fprintf(stderr, "ERROR: Couldn't open %s for reading\n", pcFileIn);
 374         free(buf);
 375         return -1;
 376     }
 377
 378     /* get number of sequences from first line and allocate memory for
 379      * distance matrix
 380      *
 381      */
 382     if (NULL == fgets(buf, MAX_BUF_SIZE, prFilePointer) ) {
 383         fprintf(stderr, "Couldn't read first line from %s\n", pcFileIn);
 384         iRetCode = -1;
 385         goto closefile_and_freebuf;
 386     }
 387     if (strlen(buf)==MAX_BUF_SIZE-1) {
 388         fprintf(stderr, "%s\n", "Looks like I couldn't read complete line. Wrong format (or too small MAX_BUF_SIZE)");
 389         iRetCode = -1;
 390         goto closefile_and_freebuf;
 391     }
 392     if (sscanf(buf, "%d", &iTotalNSeq)!=1) {
 393         fprintf(stderr, "ERROR: couldn't parse number of sequences from first line of %s\n", pcFileIn);
 394         iRetCode = -1;
 395         goto closefile_and_freebuf;
 396     }
 397
 398 #if TRACE
 399     Log(&rLog, LOG_FORCED_DEBUG, "iTotalNSeq parsed from %s is %d\n", pcFileIn, iTotalNSeq);
 400 #endif
 401
 402     if (NewSymMatrix(prSymMat_p, iTotalNSeq, iTotalNSeq)) {
 403         fprintf(stderr, "FATAL %s", "Memory allocation for distance matrix failed");
 404         iRetCode = -1;
 405         goto closefile_and_freebuf;
 406     }
 407
 408
 409     /* parse file line by line
 410      *
 411      */
 412     while (NULL != fgets(buf, MAX_BUF_SIZE, prFilePointer)) {
 413         char *szToken;
 414         int is_newseq;
 415
 416         if (MAX_BUF_SIZE-1 == strlen(buf)) {
 417             fprintf(stderr, "%s\n", "Looks like I couldn't read complete line. Wrong format (or too small MAX_BUF_SIZE)");
 418             iRetCode = -1;
 419             goto closefile_and_freebuf;
 420         }
 421
 422 #ifdef DEBUG
 423         Log(&rLog, LOG_FORCED_DEBUG, "Got line: %s\n", buf);
 424 #endif
 425
 426         /* new sequence label at beginning of line?
 427          */
 428         if (isblank(buf[0])) {
 429             is_newseq = 0;
 430         } else {
 431             is_newseq = 1;
 432         }
 433
 434         /* tokenise line and treat new sequence specially
 435          */
 436         szToken = strtok(buf, " \t");
 437         if (is_newseq==1) {
 438             iNParsedSeq++;
 439             iNParsedDists=0;
 440
 441             /* if format is lower dimensional matrix then first
 442              * sequence has no distances but might have newline
 443              * character at it's end.
 444              */
 445             while (isspace(szToken[strlen(szToken)-1])) {
 446                 szToken[strlen(szToken)-1]='\0';
 447             }
 448             /* FIXME save label? */
 449             szToken = strtok(NULL, " \t");
 450         }
 451         /* from here on it's only parsing of distances */
 452         while (szToken != NULL) {
 453             double dist;
 454             iNParsedDists++;
 455
 456             /* only parse what's needed */
 457             if (iNParsedDists!=iNParsedSeq) {
 458                 /* parse and store distance
 459                  */
 460                 if (sscanf(szToken, "%lf", &dist)!=1) {
 461                     fprintf(stderr, "Couldn't parse float from entry '%s'\n", szToken);
 462                     iRetCode = -1;
 463                     goto closefile_and_freebuf;
 464                 }
 465 #if TRACE
 466                 Log(&rLog, LOG_FORCED_DEBUG, "Parsed distance %d for seq %d = %f\n", iNParsedDists-1, iNParsedSeq-1, dist);
 467 #endif
 468                 SymMatrixSetValue(*prSymMat_p, iNParsedSeq-1, iNParsedDists-1, dist);
 469                 SymMatrixSetValue(*prSymMat_p, iNParsedDists-1, iNParsedSeq-1, dist);
 470             }
 471             szToken = strtok(NULL, " \t");
 472         }
 473     }
 474
 475     if (iTotalNSeq!=iNParsedSeq) {
 476         fprintf(stderr, "expected %d seqs, but only parsed %d\n", iTotalNSeq, iNParsedSeq);
 477         iRetCode = -1;
 478         goto closefile_and_freebuf;
 479     }
 480 #if TRACE
 481     for (i=0; i<iNParsedSeq; i++) {
 482         int j;
 483         for (j=0; j<iNParsedSeq; j++) {
 484             Log(&rLog, LOG_FORCED_DEBUG, "prSymMat_p[%d][%d]=%f\n", i, j, (*prSymMat_p)[i][j]);
 485         }
 486     }
 487 #endif
 488
 489 closefile_and_freebuf:
 490     fclose(prFilePointer);
 491     free(buf);
 492
 493     return iRetCode;
 494 }
 495 /***   end: SymMatrixRead   ***/
 496