sources/oc/oc.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <float.h>
   4 #include <string.h>
   5 #include <time.h>
   6 #include <math.h>
   7 #include <errno.h>
   8 #include <limits.h>
   9
  10 /* for regex's */
  11 #include <sys/types.h>
  12 #include <regex.h>
  13
  14 #include "gjutil.h"
  15 #include "oc.h"
  16
  17 /*#define VERY_SMALL DBL_MIN*/
  18 #define VERY_SMALL -DBL_MAX
  19 #define VERY_BIG DBL_MAX
  20 #define MAX_ID_LEN 30
  21
  22 /* Postscript plotting constants */
  23 #define POINT 72
  24 #define X_OFFSET 40
  25 #define Y_OFFSET 100
  26
  27 #define DEBUG
  28
  29 float MAXside = 11.0 * POINT;   /* Max dimensions of paper */
  30 float MINside = 8.0 * POINT;
  31 float WTreeFraction = 0.9;      /* Fraction of max taken up by Tree */
  32 float HTreeFraction = 0.9;      /* (Calculated after subtracting offsets) */
  33 float TreeExtendFac = 0.05;     /* Extension of tree beyond min/max value factor */
  34 double base_val = 1000.0;
  35 float IdOffsetFactor = 0.01;
  36 float NumberOffsetFactor = 0.05;
  37 float NumberFontSize = 12;
  38 float IdFontSize = 12;
  39 float TitleFontSize = 12;
  40
  41 /****************************************************************************
  42    Program oc:
  43
  44    General purpose cluster analysis program.  Written with flexibility
  45    in mind rather than all out efficiency.
  46
  47    Copyright(c) G. J. Barton (1993,1997,2002)
  48
  49    Author: G. J. Barton
  50    School of Life Sciences
  51    University of Dundee
  52    Dow St.
  53    Dundee DD1 5EH
  54    Scotland
  55
  56    email:geoff@compbio.dundee.ac.uk
  57
  58    Example input file:
  59
  60 7
  61 First
  62 Second
  63 Third
  64 Fourth
  65 Fifth
  66 Sixth
  67 Seventh
  68 100.0
  69 100.0
  70 50.0
  71 33.0
  72 25.0
  73 20.0
  74 100.0
  75 50.0
  76 50.0
  77 33.0
  78 33.0
  79 100.0
  80 33.0
  81 20.0
  82 25.0
  83 33.0
  84 20.0
  85 25.0
  86 100.0
  87 100.0
  88 100.0
  89
  90 11/8/93: Add single order option.  This forces the clustering to only add
  91 a single sequence at a time to the growing cluster. All other options are
  92 equivalent to the full tree method. GJB.
  93
  94 15/12/93:  Correct bug in definition of VERY_SMALL ie now = -DBL_MAX. This
  95 allows the complete and means linkage to work with negative numbers.
  96
  97 04/08/2002:Update affiliation and also add gjnoc2 routine to the distribution.
  98 ----------------------------------------------------------------------------*/
  99
 100 int
 101 main (int argc, char *argv[])
 102 {
 103   extern FILE *std_in, *std_out, *std_err;
 104   extern double base_val;
 105   int error;
 106   int n;
 107   int nclust;
 108   int i, j, k, l;
 109   int i1loc, j1loc;
 110   int i2loc, j2loc;
 111   int i3loc, j3loc;
 112   double U1val, U2val, U3val, rval;
 113   struct sc *clust;
 114   double **arr;                 /* upper diagonal array structure */
 115   int *unclust;                 /* list of unclustered entities */
 116   int nunclust;                 /* number of unclustered entities */
 117   char **idents;                /* array of idents for the entities */
 118   int *notparent;               /* list of clusters that have no children */
 119   int nnotparent;               /* number of clusters that have no children */
 120   int *inx;                     /* array of indexes - inx[i] is the re-ordered position */
 121   char *postfile;               /* name of PostScript output file */
 122
 123   /* TIME INFORMATION */
 124   clock_t start_time, end_time, initial_time, final_time;
 125
 126   /* options */
 127   int sim;                      /* similarity scores =1, distances = 2 */
 128   int type;                     /* single linkage = 1, complete linkage = 2, means = 3 */
 129   int option;                   /* use to feed the Grp_Grp function */
 130   int ps;
 131   int dolog;                    /* flag to take logs on output to .ps file */
 132   double cut;
 133   int do_cut;
 134   int showid;
 135   int amps_out;                 /* output tree file and order file for amps */
 136   int timeclus;                 /* set to 1 if timing output to std_err required */
 137   int single_order;
 138
 139   /* amps files */
 140   char *amps_order_file;
 141   char *amps_tree_file;
 142   FILE *ford;
 143   FILE *ftree;
 144   FILE *fp;
 145
 146   GJinitfile ();
 147
 148   sim = 1;
 149   type = 1;
 150   ps = 0;                       /* default is no PostScript output */
 151   showid = 0;
 152   do_cut = 0;
 153   timeclus = 0;
 154   amps_out = 0;
 155   dolog = 0;
 156   single_order = 0;
 157
 158   start_time = clock ();
 159   initial_time = start_time;
 160
 161   if (argc > 1)
 162     {
 163       for (i = 1; i < argc; ++i)
 164         {
 165           if (strcmp (argv[i], "sim") == 0)
 166             {
 167               sim = 1;
 168               base_val = VERY_BIG;
 169               cut = VERY_SMALL;
 170               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 171             }
 172           else if (strcmp (argv[i], "dis") == 0)
 173             {
 174               sim = 2;
 175               base_val = VERY_SMALL;
 176               cut = VERY_BIG;
 177               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 178             }
 179           else if (strcmp (argv[i], "single") == 0)
 180             {
 181               type = 1;
 182               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 183             }
 184           else if (strcmp (argv[i], "complete") == 0)
 185             {
 186               type = 2;
 187               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 188             }
 189           else if (strcmp (argv[i], "means") == 0)
 190             {
 191               type = 3;
 192               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 193             }
 194           else if (strcmp (argv[i], "order") == 0)
 195             {
 196               /* do a single order clustering rather than a full tree */
 197               single_order = 1;
 198               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 199             }
 200           else if (strcmp (argv[i], "ps") == 0
 201                    || strcmp (argv[i], "eps") == 0)
 202             {
 203               ++i;
 204               ps = 1;
 205               postfile = GJcat (2, argv[i], ".eps");
 206               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 207             }
 208           else if (strcmp (argv[i], "log") == 0)
 209             {
 210               dolog = 1;
 211               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 212             }
 213           else if (strcmp (argv[i], "cut") == 0)
 214             {
 215               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 216               ++i;
 217               cut = atof (argv[i]);
 218               do_cut = 1;
 219               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 220             }
 221           else if (strcmp (argv[i], "amps") == 0)
 222             {
 223               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 224               amps_out = 1;
 225               ++i;
 226               if (single_order)
 227                 {
 228                   amps_order_file = GJcat (2, argv[i], ".ord");
 229                 }
 230               else
 231                 {
 232                   amps_tree_file = GJcat (2, argv[i], ".tree");
 233                   amps_order_file = GJcat (2, argv[i], ".tord");
 234                 }
 235               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 236             }
 237           else if (strcmp (argv[i], "id") == 0)
 238             {
 239               showid = 1;       /* output idents rather than index */
 240               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 241             }
 242           else if (strcmp (argv[i], "timeclus") == 0)
 243             {
 244               timeclus = 1;     /* output idents rather than index */
 245               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 246             }
 247           else if (strcmp (argv[i], "maxside") == 0)
 248             {
 249               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 250               ++i;
 251               MAXside = atof (argv[i]) * POINT;
 252               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 253             }
 254           else if (strcmp (argv[i], "minside") == 0)
 255             {
 256               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 257               ++i;
 258               MINside = atof (argv[i]) * POINT;
 259               /*              fprintf(std_err,"Option: argv[%d]=%s\n",i,argv[i]); */
 260             }
 261           else
 262             {
 263               fprintf (std_err, "Unrecognised Option: argv[%d]=%s\n", i, argv[i]);
 264             }
 265         }
 266     }
 267   else
 268     {
 269       fprintf (std_err, "Cluster analysis program\n\n");
 270       fprintf (std_err, "Author:  G. J. Barton, University of Dundee, UK\n");
 271       fprintf (std_err, "         http://www.compbio.dundee.ac.uk\n");
 272       fprintf (std_err, "Copyright:  G. J. Barton (1993,1997,2002,2004)\n\n");
 273       fprintf (std_err, "Please cite: OC a cluster analysis program, Barton, G. J.  1993\n\n");
 274       fprintf (std_err, "Usage: oc <sim/dis> <single/complete/means> <eps file> <cut N>\n\n");
 275       fprintf (std_err, "Version 2.1a (Apr. 2005) - Requires a file to be piped to standard input\n");
 276       fprintf (std_err, "Format:  Line   1:  Number (N) of entities to cluster (e.g. 10)\n");
 277       fprintf (std_err, "Format:  Lines 2 to 2+N-1:  Identifier codes for the entities (e.g. Entity1)\n");
 278       fprintf (std_err, "Format:  N*(N-1)/2:  Distances, or similarities - ie the upper diagonal\n\n");
 279       fprintf (std_err, "Options:\n");
 280       fprintf (std_err, "sim = similarity /  dis = distances\n");
 281       fprintf (std_err, "method = single/complete/means\n");
 282       fprintf (std_err, "ps <file> = plot out portrait dendrogram to <file.eps> \n");
 283       fprintf (std_err, "eps <file> = plot out portrait dendrogram to <file.eps> \n");
 284       fprintf (std_err, "minside <dim> = define length of short (x) side of plot in inches \n");
 285       fprintf (std_err, "maxside <dim> = define length of long  (y) side of plot in inches \n");
 286       fprintf (std_err, "log = take logs before calculation \n");
 287       fprintf (std_err, "cut = only show clusters above/below the cutoff\n");
 288       fprintf (std_err, "id = output identifier codes rather than indexes for entities\n");
 289       fprintf (std_err, "timeclus = output times to generate each cluster\n");
 290       fprintf (std_err, "amps <file> = produce amps <file>.tree and <file>.tord files\n\n");
 291       fprintf (std_err, "See oc_manual.txt file for brief instructions\n");
 292
 293       exit (EXIT_SUCCESS);
 294     }
 295
 296   /* Read in the OC file */
 297   /* allocate and read in the upper diagonal array */
 298   fprintf (std_err, "Reading Upper Diagonal\n");
 299   n = read_num (std_in);
 300   idents = read_idents (std_in, n);
 301   arr = read_up_diag (std_in, n);
 302   if (dolog)
 303     {
 304       up_diag_log (arr, n);
 305     }
 306
 307 /*        write_up_diag(std_out,arr,n);*/
 308   fprintf (std_err, "Read: %d Entries\n", n);
 309   /* File has been read in */
 310
 311   end_time = clock ();
 312   fprintf (std_err, "CPU time: %f seconds\n",
 313            ((float) (end_time - start_time)) / CLOCKS_PER_SEC);
 314   start_time = end_time;
 315
 316   /* holds index of unclustered entities */
 317
 318   fprintf (std_err, "Setting up unclust\n");
 319   unclust = (int *) GJmalloc (sizeof (int) * n);
 320   for (i = 0; i < n; ++i)
 321     {
 322       unclust[i] = i;
 323     }
 324   nunclust = n;
 325
 326   fprintf (std_err, "Setting up notparent\n");
 327   notparent = (int *) GJmalloc (sizeof (int) * n);
 328   for (i = 0; i < n; ++i)
 329     {
 330       notparent[i] = 0;
 331     }
 332   nnotparent = 0;
 333
 334
 335
 336   /* set up the clust array to hold results */
 337   /* there are always n-1 clusters formed - we will also store the nth (all entities) */
 338   fprintf (std_err, "Setting up clust\n");
 339   clust = (struct sc *) GJmalloc (sizeof (struct sc) * n);
 340   end_time = clock ();
 341   fprintf (std_err, "CPU time: %f seconds\n",
 342            ((float) (end_time - start_time)) / CLOCKS_PER_SEC);
 343   start_time = end_time;
 344
 345   nclust = 0;
 346   U1val = VERY_SMALL;
 347   U2val = VERY_SMALL;
 348   U3val = VERY_SMALL;
 349   if (sim == 1)
 350     {
 351       if (type == 1)
 352         {
 353           fprintf (std_err, "Single linkage on similarity\n");
 354           option = 1;
 355         }
 356       else if (type == 2)
 357         {
 358           fprintf (std_err, "Complete linkage on similarity\n");
 359           option = 0;
 360         }
 361       else if (type == 3)
 362         {
 363           fprintf (std_err, "Means linkage on similarity\n");
 364           option = 2;
 365         }
 366     }
 367   else
 368     {
 369       if (type == 1)
 370         {
 371           fprintf (std_err, "Single linkage on distance\n");
 372           option = 0;
 373         }
 374       else if (type == 2)
 375         {
 376           fprintf (std_err, "Complete linkage on distance\n");
 377           option = 1;
 378         }
 379       else if (type == 3)
 380         {
 381           fprintf (std_err, "Means linkage on distance\n");
 382           option = 2;
 383         }
 384     }
 385
 386 /* implement max similarity single linkage first */
 387
 388   fprintf (std_err, "Doing Cluster Analysis...\n");
 389   while (nclust < (n - 1))
 390     {
 391       if (sim == 1)
 392         {
 393           U1val = VERY_SMALL;
 394           U2val = VERY_SMALL;
 395           U3val = VERY_SMALL;
 396         }
 397       else
 398         {
 399           U1val = VERY_BIG;
 400           U2val = VERY_BIG;
 401           U3val = VERY_BIG;
 402         }
 403       if (nunclust > 0)
 404         {
 405 /*                      fprintf(std_err,"Option 1\n");*/
 406           if (!single_order || (single_order && nclust == 0))
 407             {
 408               /* we've not clustered all entities into groups */
 409               /* find max of what is left */
 410               for (i = 0; i < nunclust - 1; ++i)
 411                 {
 412                   for (j = (i + 1); j < nunclust; ++j)
 413                     {
 414                       rval = val_A_B (arr, n, unclust[i], unclust[j]);
 415                       if (sim == 1)
 416                         {
 417                           if (rval > U1val)
 418                             {
 419                               U1val = rval;
 420                               i1loc = unclust[i];
 421                               j1loc = unclust[j];
 422                             }
 423                         }
 424                       else
 425                         {
 426                           if (rval < U1val)
 427                             {
 428                               U1val = rval;
 429                               i1loc = unclust[i];
 430                               j1loc = unclust[j];
 431                             }
 432                         }
 433                     }
 434                 }               /* ... for(i=0... */
 435             }                   /* ...if(!single_order... */
 436         }
 437       if (nnotparent > 0 && nunclust > 0)
 438         {
 439 /*                      fprintf(std_err,"Option 2\n");*/
 440           /* we have some clusters */
 441           /* so compare them to the unclustered entitites */
 442           for (i = 0; i < nunclust; ++i)
 443             {
 444               for (j = 0; j < nnotparent; ++j)
 445                 {
 446                   k = notparent[j];
 447                   rval =
 448                     val_Grp_B (arr, n, clust[k].members, clust[k].n,
 449                                unclust[i], option);
 450                   if (sim == 1)
 451                     {
 452                       if (rval >= U2val)
 453                         {
 454                           /* need >= because want biggest previous cluster */
 455                           U2val = rval;
 456                           i2loc = unclust[i];
 457                           j2loc = k;
 458                         }
 459                     }
 460                   else
 461                     {
 462                       if (rval <= U2val)
 463                         {
 464                           /* need <= because want biggest previous cluster */
 465                           U2val = rval;
 466                           i2loc = unclust[i];
 467                           j2loc = k;
 468                         }
 469                     }
 470                 }
 471             }
 472         }
 473       if (nnotparent > 1)
 474         {
 475 /*                      fprintf(std_err,"Option 3\n");*/
 476           /* compare clusters to clusters if there are more than 1 */
 477           for (i = 0; i < nnotparent - 1; ++i)
 478             {
 479               k = notparent[i];
 480               for (j = (i + 1); j < nnotparent; ++j)
 481                 {
 482                   l = notparent[j];
 483                   rval = val_Grp_Grp (arr, n, clust[k].members, clust[k].n,
 484                                       clust[l].members, clust[l].n, option);
 485                   if (sim == 1)
 486                     {
 487                       if (rval >= U3val)
 488                         {
 489                           U3val = rval;
 490                           i3loc = k;
 491                           j3loc = l;
 492                         }       /* if(rval >= U3val... */
 493                     }
 494                   else
 495                     {
 496                       if (rval <= U3val)
 497                         {
 498                           U3val = rval;
 499                           i3loc = k;
 500                           j3loc = l;
 501                         }       /* if(rval >= U3val... */
 502                     }
 503                 }               /* if j=(i+1) ... */
 504             }                   /* for(i=0... */
 505         }
 506       /* find which search gave the biggest (or smallest) value */
 507       if ((sim == 1 && U3val >= U1val && U3val >= U2val) ||
 508           (sim == 2 && U3val <= U1val && U3val <= U2val))
 509         {
 510           /* joining of  two clusters */
 511           clust[nclust].score = U3val;
 512           clust[nclust].n = clust[i3loc].n + clust[j3loc].n;
 513           clust[nclust].members =
 514             (int *) GJmalloc (sizeof (int) * clust[nclust].n);
 515           for (i = 0; i < clust[i3loc].n; ++i)
 516             {
 517               clust[nclust].members[i] = clust[i3loc].members[i];
 518             }
 519           for (i = 0; i < clust[j3loc].n; ++i)
 520             {
 521               clust[nclust].members[i + clust[i3loc].n] =
 522                 clust[j3loc].members[i];
 523             }
 524           clust[nclust].parentA = &clust[i3loc];
 525           clust[nclust].parentB = &clust[j3loc];
 526           clust[nclust].lab = 1;
 527           remove_notparent (notparent, &nnotparent, i3loc);
 528           sub_notparent (notparent, &nnotparent, j3loc, nclust);
 529         }
 530       else if ((sim == 1 && U2val >= U1val && U2val >= U3val) ||
 531                (sim == 2 && U2val <= U1val && U2val <= U3val))
 532         {
 533           /* joining of  a single entity to an existing cluster */
 534           clust[nclust].score = U2val;
 535           clust[nclust].n = 1 + clust[j2loc].n;
 536           clust[nclust].members =
 537             (int *) GJmalloc (sizeof (int) * clust[nclust].n);
 538           /* copy previous cluster list into this one - not strictly neccessary */
 539           /* but done to simplify cluster comparisons above */
 540           for (i = 1; i < clust[nclust].n; ++i)
 541             {
 542               clust[nclust].members[i] = clust[j2loc].members[i - 1];
 543             }
 544           /* stick solitary cluster member on the beginning */
 545           clust[nclust].members[0] = i2loc;
 546           /* put in parent pointers */
 547           clust[nclust].parentA = make_entity (i2loc, base_val);
 548           clust[nclust].parentB = &clust[j2loc];
 549           clust[nclust].lab = 1;
 550           sub_notparent (notparent, &nnotparent, j2loc, nclust);
 551           remove_unclust (unclust, &nunclust, i2loc);
 552         }
 553       else if ((sim == 1 && U1val >= U2val && U1val >= U3val) ||
 554                (sim == 2 && U1val <= U2val && U1val <= U3val))
 555         {
 556           /* unclustered pair were biggest */
 557           clust[nclust].score = U1val;
 558           clust[nclust].n = 2;
 559           clust[nclust].members =
 560             (int *) GJmalloc (sizeof (int) * clust[nclust].n);
 561           clust[nclust].members[0] = i1loc;
 562           clust[nclust].members[1] = j1loc;
 563           clust[nclust].parentA = make_entity (i1loc, base_val);
 564           clust[nclust].parentB = make_entity (j1loc, base_val);
 565           clust[nclust].lab = 1;
 566           add_notparent (notparent, &nnotparent, nclust);
 567           remove_unclust (unclust, &nunclust, i1loc);
 568           remove_unclust (unclust, &nunclust, j1loc);
 569         }
 570
 571 /*              fprintf(std_err,"Current cluster\n");*/
 572 /*              show_entity(&clust[nclust],std_err);*/
 573
 574       if (timeclus)
 575         {
 576           fprintf (std_err, "Cluster: %d DONE!  ", nclust);
 577           end_time = clock ();
 578           fprintf (std_err, "CPU time: %f seconds\n",
 579                    ((float) (end_time - start_time)) / CLOCKS_PER_SEC);
 580           start_time = end_time;
 581         }
 582
 583       ++nclust;
 584     }
 585   /* now we should have all clusters stored in the clust array */
 586   /* set up the inx array */
 587   inx = (int *) GJmalloc (sizeof (int) * n);
 588   for (i = 0; i < clust[nclust - 1].n; ++i)
 589     {
 590       inx[clust[nclust - 1].members[i]] = i;
 591     }
 592
 593   if (do_cut)
 594     {
 595       /* only output the biggest clusters that are above the cutoff */
 596       if (sim == 1)
 597         {
 598           for (i = (nclust - 1); i >= 0; --i)
 599             {
 600               if (clust[i].lab == 1)
 601                 {
 602                   /* only output clusters that aren't already output as part of something bigger */
 603                   if (clust[i].score >= cut)
 604                     {
 605                       fprintf (std_out, "## %d %g %d\n", i, clust[i].score,
 606                                clust[i].n);
 607                       show_id_entity (&clust[i], idents, std_out);
 608                       mark_parents (&clust[i]);
 609                     }
 610                 }
 611             }
 612         }
 613       else
 614         {
 615           for (i = (nclust - 1); i >= 0; --i)
 616             {
 617               if (clust[i].lab == 1)
 618                 {
 619                   /* only output clusters that aren't already output as part of something bigger */
 620                   if (clust[i].score <= cut)
 621                     {
 622                       fprintf (std_out, "## %d %g %d\n", i, clust[i].score,
 623                                clust[i].n);
 624                       show_id_entity (&clust[i], idents, std_out);
 625                       mark_parents (&clust[i]);
 626                     }
 627                 }
 628             }
 629         }
 630       /* now output the unclustered entities */
 631       fprintf (std_out, "\nUNCLUSTERED ENTITIES\n");
 632       write_unclustered (&clust[nclust - 1], idents, std_out);
 633
 634     }
 635   else
 636     {
 637       /* output all clusters from smallest to largest */
 638       if (amps_out)
 639         {
 640 /*          fprintf(std_err,"amps_out not implemented yet\n");*/
 641           if (!single_order)
 642             {
 643               fprintf (std_err,
 644                        "Opening AMPS order and tree files for writing\n");
 645               ford = GJfopen (amps_order_file, "w", 1);
 646               ftree = GJfopen (amps_tree_file, "w", 1);
 647             }
 648           else
 649             {
 650               fprintf (std_err, "Opening AMPS order file for writing\n");
 651               ford = GJfopen (amps_order_file, "w", 1);
 652             }
 653           /* output order file */
 654           if (single_order)
 655             {
 656               /* order file must be reversed */
 657               for (i = nclust; i > -1; --i)
 658                 {
 659                   fprintf (ford, "%11d%20.2f%10d\n",
 660                            clust[nclust - 1].members[i] + 1, 0.0, 0);
 661                 }
 662             }
 663           else
 664             {
 665               /* order file for use with the tree */
 666               for (i = 0; i < (nclust + 1); ++i)
 667                 {
 668                   fprintf (ford, "%11d%20.2f%10d\n",
 669                            clust[nclust - 1].members[i] + 1, 0.0, 0);
 670                 }
 671             }
 672           GJfclose (ford, 1);
 673           if (!single_order)
 674             {
 675               /* output the tree_file */
 676               for (i = 0; i < nclust; ++i)
 677                 {
 678                   print_amps_cluster (clust[i].parentA, inx, ftree);
 679                   print_amps_cluster (clust[i].parentB, inx, ftree);
 680                 }
 681               GJfclose (ftree, 1);
 682             }
 683         }
 684       else
 685         {
 686           for (i = 0; i < nclust; ++i)
 687             {
 688               fprintf (std_out, "## %d %g %d\n", i, clust[i].score,
 689                        clust[i].n);
 690               if (showid)
 691                 {
 692                   show_id_entity (&clust[i], idents, std_out);
 693                 }
 694               else
 695                 {
 696                   show_entity (&clust[i], std_out);
 697                 }
 698
 699               /*            show_inx_entity(&clust[i],inx,std_out); */
 700
 701               /*            fprintf(std_out,"Parents:\n");
 702                  show_entity(clust[i].parentA,std_out);
 703                  show_entity(clust[i].parentB,std_out);
 704                */
 705             }
 706         }
 707     }
 708
 709   if (ps == 1)
 710     {
 711       fp = GJfopen (postfile, "w", 1);
 712       draw_dendrogram (fp, clust, idents, inx, nclust, sim);
 713     }
 714
 715   final_time = clock ();
 716   fprintf (std_err, "Total CPU time: %f seconds\n",
 717            (float) (final_time - initial_time) / CLOCKS_PER_SEC);
 718   exit (EXIT_SUCCESS);
 719
 720 }
 721
 722 /*********************************************************************
 723 up_diag_log:  take logs for the upper diagonal structure
 724 -------------------------------------------------------------------*/
 725
 726 void
 727 up_diag_log (double **arr, int n)
 728 {
 729   int i, j;
 730
 731   for (i = 0; i < (n - 1); ++i)
 732     {
 733       for (j = 0; j < (n - i - 1); ++j)
 734         {
 735           arr[i][j] = log (arr[i][j]);
 736         }
 737     }
 738 }
 739
 740 /**********************************************************************
 741 write_up_diag:
 742 writes out n, the number of entities, followed by
 743 n*(n-1)/2 numbers that are the pair distances or similarities between
 744 the n entities.
 745
 746 -----------------------------------------------------------------------*/
 747 void
 748 write_up_diag (FILE * infile, double **ret_val, int n)
 749 {
 750   int i, j;
 751
 752   fprintf (infile, "Number of Entities: %d\n", n);
 753   for (i = 0; i < (n - 1); ++i)
 754     {
 755       for (j = 0; j < (n - i - 1); ++j)
 756         {
 757           fprintf (infile, " %f", ret_val[i][j]);
 758         }
 759       fprintf (infile, "\n");
 760     }
 761 }
 762
 763 /**********************************************************************
 764 show_entity:
 765
 766 print out the score, number of members and the members for
 767 the entity
 768 ----------------------------------------------------------------------*/
 769 void
 770 show_entity (struct sc *entity, FILE * out)
 771 {
 772   int i;
 773
 774   if (entity == NULL)
 775     {
 776       fprintf (out, "No entity\n");
 777       return;
 778     }
 779
 780   fprintf (out, "Entity Score: %g  Number of members: %d\n",
 781            entity->score, entity->n);
 782
 783   for (i = 0; i < entity->n; ++i)
 784     {
 785       fprintf (out, " %d", entity->members[i]);
 786     }
 787   fprintf (out, "\n");
 788 }
 789
 790 /**********************************************************************
 791 show_id_entity:
 792
 793 print out the score, number of members and the members for
 794 the entity using the id codes rather than numerical positions
 795 ----------------------------------------------------------------------*/
 796 void
 797 show_id_entity (struct sc *entity, char **id, FILE * out)
 798 {
 799   int i;
 800
 801   if (entity == NULL)
 802     {
 803       fprintf (out, "No entity\n");
 804       return;
 805     }
 806
 807 /*    fprintf(out,"Entity Score: %f  Number of members: %d\n",
 808             entity->score,entity->n);
 809 */
 810
 811   for (i = 0; i < entity->n; ++i)
 812     {
 813       fprintf (out, " %s", id[entity->members[i]]);
 814     }
 815   fprintf (out, "\n");
 816 }
 817
 818 /**********************************************************************
 819 show_inx_entity:
 820
 821 print out the score, number of members and the members for
 822 the entity - followed by the index for the entity
 823 ----------------------------------------------------------------------*/
 824 void
 825 show_inx_entity (struct sc *entity, int *inx, FILE * out)
 826 {
 827   int i;
 828
 829   if (entity == NULL)
 830     {
 831       fprintf (out, "No entity\n");
 832       return;
 833     }
 834
 835   fprintf (out, "Entity Score: %g  Number of members: %d\n",
 836            entity->score, entity->n);
 837   for (i = 0; i < entity->n; ++i)
 838     {
 839       fprintf (out, " %d", entity->members[i]);
 840     }
 841   fprintf (out, "\n");
 842   for (i = 0; i < entity->n; ++i)
 843     {
 844       fprintf (out, " %d", inx[entity->members[i]]);
 845     }
 846   fprintf (out, "\n");
 847 }
 848
 849 /**********************************************************************
 850 print_amps_cluster:
 851
 852 print out the cluster in AMPS tree file format.
 853 ie 1x,20i5  Fortran format.
 854
 855 inx is necessary to so that the identifying codes are re-ordered
 856 appropriately for AMPS
 857
 858 Need +1 cos AMPS uses sequence numbering from 1...N not 0...N-1
 859
 860 ----------------------------------------------------------------------*/
 861 void
 862 print_amps_cluster (struct sc *entity, int *inx, FILE * fp)
 863 {
 864   int i;
 865   int j;
 866
 867   fprintf (fp, " %5d\n", entity->n);
 868
 869   j = 0;
 870   for (i = 0; i < entity->n; ++i)
 871     {
 872       if (j == 20)
 873         {
 874           j = 0;
 875         }
 876       if (j == 0)
 877         {
 878           fprintf (fp, " %5d", inx[entity->members[i]] + 1);
 879         }
 880       else
 881         {
 882           fprintf (fp, "%5d", inx[entity->members[i]] + 1);
 883         }
 884       if (j == 19)
 885         {
 886           fprintf (fp, "\n");
 887         }
 888       ++j;
 889     }
 890   /* change from < to <= to avoid blunders at end of line */
 891   if (j <= 19)
 892     fprintf (fp, "\n");
 893 }
 894
 895 /******************************************************************
 896 make_entity: return pointer to malloc'd space for a struct sc
 897         structure for a single entity. i.
 898         base_val: The value form minimum distance, or maximum similarity.
 899 -------------------------------------------------------------------*/
 900 struct sc *
 901 make_entity (int i, double base_val)
 902 {
 903   struct sc *ret_val;
 904   ret_val = (struct sc *) GJmalloc (sizeof (struct sc));
 905   ret_val->score = base_val;
 906   ret_val->n = 1;
 907   ret_val->members = (int *) GJmalloc (sizeof (int));
 908   ret_val->members[0] = i;
 909   ret_val->parentA = NULL;
 910   ret_val->parentB = NULL;
 911   ret_val->lab = 1;
 912   return ret_val;
 913 }
 914
 915 /******************************************************************
 916 val_A_B: Given the upper diagonal array, return the value stored
 917          in the notional arr[i][j].
 918          Where 0 <= i < n-1
 919          and   i < j < n
 920
 921 Note:  This could/should be turned into a macro - without the bounds
 922        check.
 923 ------------------------------------------------------------------
 924
 925 double val_A_B(
 926         double **arr,
 927         int n,
 928         int i,
 929         int j)
 930
 931 {
 932     extern FILE *std_err;
 933     if(!(i >= 0 && i < n-1 && i < j && j < n)){
 934         GJerror("Invalid request in val_A_B");
 935         fprintf(std_err,"i,j,n %d %d %d \n",i,j,n);
 936         exit(-1);
 937     }
 938
 939     return arr[i][j-i-1];
 940 }
 941 */
 942 /*******************************************************************
 943 val_Grp_B: Given upper diagonal array,
 944         find min, max, or mean of comparison between the
 945         entities in grp and j
 946         choice = 0 =  find min value of grp-grp comparison
 947                  1 =  find max ...
 948                  2 =  find mean ...
 949
 950 Note:  there is some duplication in the code to avoid testing choice
 951        on every step in the loop
 952 --------------------------------------------------------------------*/
 953 double
 954 val_Grp_B (double **arr,        /* upper diagonal array */
 955            int n,               /* side of array */
 956            int *grpA,           /* array containing entities of group */
 957            int ng,              /* number of members in grpA */
 958            int j,               /* entity to compare to grpA */
 959            int choice           /* see above */
 960   )
 961 {
 962   int k;
 963   double val;
 964   double rval;
 965
 966   if (choice == 0)
 967     {
 968       /* min case */
 969       val = VERY_BIG;
 970       for (k = 0; k < ng; ++k)
 971         {
 972           if (grpA[k] < j)
 973             {
 974               rval = val_A_B (arr, n, grpA[k], j);
 975             }
 976           else
 977             {
 978               rval = val_A_B (arr, n, j, grpA[k]);
 979             }
 980           if (rval < val)
 981             val = rval;
 982         }
 983     }
 984   else if (choice == 1)
 985     {
 986       val = VERY_SMALL;
 987       for (k = 0; k < ng; ++k)
 988         {
 989           if (grpA[k] < j)
 990             {
 991               rval = val_A_B (arr, n, grpA[k], j);
 992             }
 993           else
 994             {
 995               rval = val_A_B (arr, n, j, grpA[k]);
 996             }
 997           if (rval > val)
 998             val = rval;
 999         }
1000     }
1001   else if (choice == 2)
1002     {
1003       val = 0.0;
1004       for (k = 0; k < ng; ++k)
1005         {
1006           if (grpA[k] < j)
1007             {
1008               rval = val_A_B (arr, n, grpA[k], j);
1009             }
1010           else
1011             {
1012               rval = val_A_B (arr, n, j, grpA[k]);
1013             }
1014           val += rval;
1015         }
1016       val /= ng;
1017     }
1018   return val;
1019 }
1020
1021 /*******************************************************************
1022 val_Grp_Grp: Given upper diagonal array,
1023         find min, max, or mean of comparison between the
1024         entities in grpA and grpB
1025         choice = 0 =  find min value of grp-grp comparison
1026                  1 =  find max ...
1027                  2 =  find mean ...
1028
1029 Note:  there is some duplication in the code to avoid testing choice
1030        on every step in the loop
1031 --------------------------------------------------------------------*/
1032 double
1033 val_Grp_Grp (double **arr,      /* upper diagonal array */
1034              int n,             /* side of array */
1035              int *grpA,         /* array containing entities of group */
1036              int ngA,           /* number of members in grpA */
1037              int *grpB,         /* array containing entities of group B */
1038              int ngB,           /* number of members in grpB */
1039              int choice         /* see above */
1040   )
1041 {
1042   int k, l;
1043   double val;
1044   double rval;
1045
1046   if (choice == 0)
1047     {
1048       /* min case */
1049       val = VERY_BIG;
1050       for (k = 0; k < ngA; ++k)
1051         {
1052           for (l = 0; l < ngB; ++l)
1053             {
1054               if (grpA[k] < grpB[l])
1055                 {
1056                   rval = val_A_B (arr, n, grpA[k], grpB[l]);
1057                 }
1058               else
1059                 {
1060                   rval = val_A_B (arr, n, grpB[l], grpA[k]);
1061                 }
1062               if (rval < val)
1063                 val = rval;
1064             }
1065         }
1066     }
1067   else if (choice == 1)
1068     {
1069       val = VERY_SMALL;
1070       for (k = 0; k < ngA; ++k)
1071         {
1072           for (l = 0; l < ngB; ++l)
1073             {
1074               if (grpA[k] < grpB[l])
1075                 {
1076                   rval = val_A_B (arr, n, grpA[k], grpB[l]);
1077                 }
1078               else
1079                 {
1080                   rval = val_A_B (arr, n, grpB[l], grpA[k]);
1081                 }
1082               if (rval > val)
1083                 val = rval;
1084             }
1085         }
1086     }
1087   else if (choice == 2)
1088     {
1089       val = 0.0;
1090       for (k = 0; k < ngA; ++k)
1091         {
1092           for (l = 0; l < ngB; ++l)
1093             {
1094               if (grpA[k] < grpB[l])
1095                 {
1096                   rval = val_A_B (arr, n, grpA[k], grpB[l]);
1097                 }
1098               else
1099                 {
1100                   rval = val_A_B (arr, n, grpB[l], grpA[k]);
1101                 }
1102               val += rval;
1103             }
1104         }
1105       val /= (ngA * ngB);
1106     }
1107   return val;
1108 }
1109
1110 /*****************************************************************
1111 remove_unclust:
1112
1113 Take the array containing integers showing members that have
1114 yet to be clustered and remove an element.  nunclust is returned
1115 reduced by one.
1116 -----------------------------------------------------------------*/
1117 void
1118 remove_unclust (int *unclust, int *nunclust, int val)
1119 {
1120   int i, j, found;
1121   i = 0;
1122   found = 0;
1123
1124   while (i < *nunclust)
1125     {
1126       if (unclust[i] == val)
1127         {
1128           ++i;
1129           found = 1;
1130         }
1131       if (found)
1132         {
1133           j = i - 1;
1134         }
1135       else
1136         {
1137           j = i;
1138         }
1139       unclust[j] = unclust[i];
1140       ++i;
1141     }
1142   --(*nunclust);
1143 }
1144
1145 void
1146 iprintarr (int *arr, int n, FILE * out)
1147 {
1148   int i;
1149   for (i = 0; i < n; ++i)
1150     {
1151       fprintf (out, " %d", arr[i]);
1152     }
1153   fprintf (out, "\n");
1154 }
1155
1156 /******************************************************************
1157 no_share:  return 1 if grpA and grpB have no common members
1158                  0 if there are common members
1159 ----------------------------------------------------------------*/
1160
1161 int
1162 no_share (int *grpA, int ngA, int *grpB, int ngB)
1163 {
1164   int i, j;
1165
1166   for (i = 0; i < ngA; ++i)
1167     {
1168       for (j = 0; j < ngB; ++j)
1169         {
1170           if (grpA[i] == grpB[j])
1171             return 0;
1172         }
1173     }
1174   return 1;
1175 }
1176
1177 /*****************************************************************
1178 sub_notparent:  given the notparent list substitutes A by B
1179 -----------------------------------------------------------------*/
1180 void
1181 sub_notparent (int *np, int *n, int A, int B)
1182 {
1183   int i;
1184   for (i = 0; i < *n; ++i)
1185     {
1186       if (np[i] == A)
1187         {
1188           np[i] = B;
1189           return;
1190         }
1191     }
1192 }
1193
1194 /*****************************************************************
1195 remove_notparent:  given the notparent list removes A from the list
1196
1197 -----------------------------------------------------------------*/
1198 void
1199 remove_notparent (int *np, int *n, int A)
1200 {
1201   remove_unclust (np, n, A);
1202 }
1203
1204 /*****************************************************************
1205 add_notparent:  given the notparent list adds A to the list
1206
1207 -----------------------------------------------------------------*/
1208 void
1209 add_notparent (int *np, int *n, int A)
1210 {
1211   np[*n] = A;
1212   ++(*n);
1213 }
1214
1215 /**************************************************************************
1216 Return a pointer to space for an upper diagonal square array stored rowwise
1217 G. J. B. 16 May 1993
1218
1219 double precision array.
1220 side = n;
1221 --------------------------------------------------------------------------*/
1222 double **
1223 GJDudarr (int n)
1224 {
1225   double **ret_val;
1226   int i, nm1;
1227
1228   nm1 = n - 1;
1229   ret_val = (double **) GJmalloc (sizeof (double *) * nm1);
1230
1231
1232   for (i = 0; i < nm1; ++i)
1233     {
1234       ret_val[i] = (double *) GJmalloc (sizeof (double) * (nm1 - i));
1235     }
1236
1237   return ret_val;
1238 }
1239
1240 /***************************************************************************
1241 draw_dendrogram:
1242
1243 Output the dendrogram to fout using PostScript commands
1244
1245 ---------------------------------------------------------------------------*/
1246 void
1247 draw_dendrogram (FILE * fout,
1248                  struct sc *clust, char **idents, int *inx, int n, int sim)
1249 {
1250   extern double base_val;
1251   extern float MINside, MAXside;
1252   extern float WTreeFraction, HTreeFraction, TreeExtendFac;
1253
1254   float Twidth, Theight;        /* x and y limits for the tree part of the plot */
1255   float Xoffset, Yoffset;       /* x and y offset (points) */
1256
1257   int i;
1258   double x1, x2, x3, y1, y2;
1259   double Xrange, Yrange, Xmin, Xmax, Ymin, Ymax;        /*,TXmax; */
1260 /*    double Xscale,Yscale;*/
1261   double Xid, Yid;
1262
1263   char *buff;
1264
1265   buff = GJstrcreate (MAX_ID_LEN, NULL);
1266
1267   Xoffset = X_OFFSET;
1268   Yoffset = Y_OFFSET;
1269
1270   Twidth = (MINside * WTreeFraction) - Xoffset;
1271   Theight = (MAXside * HTreeFraction) - Yoffset;
1272
1273   Xmin = VERY_BIG;
1274   Xmax = VERY_SMALL;
1275   Ymin = VERY_BIG;
1276   Ymax = VERY_SMALL;
1277
1278   /* find the min and max in X excluding the base_val */
1279   for (i = 0; i < n; ++i)
1280     {
1281       if (sim == 1)
1282         {
1283           if (clust[i].score < Xmin)
1284             Xmin = clust[i].score;
1285           if (clust[i].score != base_val && clust[i].score > Xmax)
1286             Xmax = clust[i].score;
1287         }
1288       else
1289         {
1290           if (clust[i].score != base_val && clust[i].score < Xmin)
1291             Xmin = clust[i].score;
1292           if (clust[i].score > Xmax)
1293             Xmax = clust[i].score;
1294         }
1295     }
1296   /* set the max/min value slightly bigger/smaller than the actual max/min */
1297
1298   if (sim == 1)
1299     Xmax *= (1.0 + TreeExtendFac);
1300   if (sim == 2)
1301     Xmin *= (1.0 - TreeExtendFac);
1302   Xrange = Xmax - Xmin;
1303   Ymin = 0.0;
1304   Ymax = n;
1305   Yrange = (double) Ymax - Ymin;
1306
1307   /* write out the preamble */
1308
1309   fprintf (fout, "%%!PS-Adobe-2.0 EPSF-2.0\n");
1310   fprintf (fout, "%%%%Title: oc output file\n");
1311   fprintf (fout,
1312            "%%%%Creator: oc by G. J. Barton: http://www.compbio.dundee.ac.uk\n");
1313   /*    fprintf(fout,"%%%%Pages: 1\n"); */
1314   fprintf (fout, "%%%%DocumentFonts: Times-Roman\n");
1315   fprintf (fout, "%%%%BoundingBox: %d %d %d %d\n", (int) Xoffset,
1316            (int) Yoffset, (int) MINside, (int) MAXside);
1317   fprintf (fout, "%%%%EndComments\n");
1318
1319   fprintf (fout, "/Times-Roman findfont 12 scalefont setfont\n");
1320
1321   /* get the X position for id codes then write them out */
1322   if (sim == 1)
1323     {
1324       Xid = Xmax + (Xrange * IdOffsetFactor);
1325     }
1326   else
1327     {
1328       Xid = Xmin - (Xrange * IdOffsetFactor);
1329     }
1330   Xid = trans_x (Xid, sim, Xmin, Xmax, Xrange, Twidth, Xoffset);
1331
1332   for (i = 0; i < (n + 1); ++i)
1333     {
1334       Yid = trans_y ((double) i, Ymin, Ymax, Yrange, Theight, Yoffset);
1335       PSPutText ((float) Xid, (float) Yid, idents[clust[n - 1].members[i]],
1336                  fout);
1337     }
1338
1339   /* write min and max scores below the plot */
1340   x1 = Xmin;
1341   x2 = Xmax;
1342   y1 = Ymin - (Yrange * NumberOffsetFactor);
1343   x1 = trans_x (x1, sim, Xmin, Xmax, Xrange, Twidth, Xoffset);
1344   x2 = trans_x (x2, sim, Xmin, Xmax, Xrange, Twidth, Xoffset);
1345   y1 = trans_y (y1, Ymin, Ymax, Yrange, Theight, Yoffset);
1346   sprintf (buff, "%g", Xmin);
1347   PSPutText ((float) x1, (float) y1, buff, fout);
1348   sprintf (buff, "%g", Xmax);
1349   PSPutText ((float) x2, (float) y1, buff, fout);
1350
1351   for (i = 0; i < n; ++i)
1352     {
1353       float bob = 1. / 0.;
1354
1355       x1 = clust[i].parentA->score;
1356       x2 = clust[i].score;
1357       x3 = clust[i].parentB->score;
1358       if (x1 == base_val && sim == 1)
1359         x1 = Xmax;
1360       if (x1 == base_val && sim == 2)
1361         x1 = Xmin;
1362       if (x3 == base_val && sim == 1)
1363         x3 = Xmax;
1364       if (x3 == base_val && sim == 2)
1365         x3 = Xmin;
1366       y1 = get_mean (clust[i].parentA[0].members, inx, clust[i].parentA[0].n);
1367       y2 = get_mean (clust[i].parentB[0].members, inx, clust[i].parentB[0].n);
1368       /* transform x and y to fit on the page, then */
1369       x1 = trans_x (x1, sim, Xmin, Xmax, Xrange, Twidth, Xoffset);
1370       x2 = trans_x (x2, sim, Xmin, Xmax, Xrange, Twidth, Xoffset);
1371
1372       /* test for inf */
1373       if (bob == x2)
1374         {
1375           fprintf (stderr, "%e, %e, %e, %e, %e, %e, %e\n", x1, (double)sim, Xmin,
1376                    Xmax, Xrange, Twidth, Xoffset);
1377           exit (EXIT_FAILURE);
1378         }
1379
1380       x3 = trans_x (x3, sim, Xmin, Xmax, Xrange, Twidth, Xoffset);
1381       y1 = trans_y (y1, Ymin, Ymax, Yrange, Theight, Yoffset);
1382       y2 = trans_y (y2, Ymin, Ymax, Yrange, Theight, Yoffset);
1383
1384       draw_arch (fout, (float) x1, (float) x2, (float) x3, (float) y1,
1385                  (float) y2);
1386     }
1387   fprintf (fout, "showpage\n");
1388 }
1389
1390 double
1391 get_mean (int *arr, int *inx, int n)
1392 {
1393   double sum;
1394   int i;
1395   sum = 0.0;
1396   for (i = 0; i < n; ++i)
1397     {
1398       sum += (double) inx[arr[i]];
1399     }
1400   return sum / n;
1401 }
1402
1403 void
1404 draw_arch (FILE * fout, float x1, float x2, float x3, float y1, float y2)
1405 {
1406   fprintf (fout, "%f %f moveto\n", x1, y1);
1407   fprintf (fout, "%f %f lineto\n", x2, y1);
1408   fprintf (fout, "%f %f lineto\n", x2, y2);
1409   fprintf (fout, "%f %f lineto\n", x3, y2);
1410   fprintf (fout, "stroke newpath\n");
1411 }
1412
1413 double
1414 trans_x (double X,
1415          int sim,
1416          double Xmin, double Xmax, double Xrange, float Twidth, float Xoffset)
1417 {
1418   if (sim == 1)
1419     {
1420       return (((X - Xmin) / Xrange) * (double) Twidth) + Xoffset;
1421     }
1422   else
1423     {
1424       return (((Xmax - X) / Xrange) * (double) Twidth) + Xoffset;
1425     }
1426 }
1427
1428 double
1429 trans_y (double X,
1430          double Xmin, double Xmax, double Xrange, float Twidth, float Xoffset)
1431 {
1432   return (((X - Xmin) / Xrange) * (double) Twidth) + Xoffset;
1433 }
1434
1435 void
1436 PSPutText (float x, float y, char *text, FILE * outf)
1437 {
1438   fprintf (outf, " %f %f moveto (%s) show \n", x, y, text);
1439 }
1440
1441 /**********************************************************************
1442 mark_parents:  This works back up the tree setting the entity->lab
1443 values to 0 for all parents.
1444 Hopefully the recursion won't run out of space for the size of problem
1445 we typically deal with...
1446 ----------------------------------------------------------------------*/
1447
1448 void
1449 mark_parents (struct sc *entity)
1450 {
1451   entity->lab = 0;
1452   if (entity->parentA != NULL)
1453     {
1454       mark_parents (entity->parentA);
1455     }
1456   if (entity->parentB != NULL)
1457     {
1458       mark_parents (entity->parentB);
1459     }
1460 }
1461
1462 /*********************************************************************
1463 write_unclustered:
1464
1465 Work back down the tree and print out the identifiers of entitites that
1466 are not clustered into any group.  If mark_parents has not been called,
1467 then all entities will be output.  Otherwise, only those that have
1468 not been marked will be output.
1469 --------------------------------------------------------------------*/
1470
1471 void
1472 write_unclustered (struct sc *entity, char **idents, FILE * fout)
1473 {
1474   if (entity->lab == 1)
1475     {
1476       /* this could be an unclustered entity or a cluster
1477          containing unclustered entities  */
1478       if (entity->parentA == NULL && entity->parentB == NULL)
1479         {
1480           fprintf (fout, "%s\n", idents[entity->members[0]]);
1481           return;
1482         }
1483       /* we don't need both if's since any entitity with one NULL parent
1484          must also have the other parent NULL
1485        */
1486       if (entity->parentA != NULL)
1487         {
1488           write_unclustered (entity->parentA, idents, fout);
1489         }
1490       if (entity->parentB != NULL)
1491         {
1492           write_unclustered (entity->parentB, idents, fout);
1493         }
1494     }
1495 }
1496
1497 /*********************************************************************
1498  * OC FILE READING FUNCTIONS
1499 --------------------------------------------------------------------*/
1500
1501
1502 /*********************************************************************
1503 sys_fatal_error:
1504
1505 Pass a string to the function, prints it and the system error, then exits with
1506 failure.
1507 --------------------------------------------------------------------*/
1508
1509 void
1510 sys_fatal_error (const char *error)
1511 {
1512   perror (error);
1513   exit (EXIT_FAILURE);
1514 }
1515
1516 /*********************************************************************
1517 getline:
1518
1519 Returns a null terminated string representing a line, from a file pointer. Or,
1520 if the file is empty, NULL.
1521 --------------------------------------------------------------------*/
1522
1523 #define BUFFER 80               /* Size of buffer for reading in a line */
1524
1525 char *
1526 getline_loc (FILE * file)
1527 {
1528   int c, i, j;
1529   char *line;
1530   i = j = 0;
1531
1532   if ((line = malloc (sizeof (char) * BUFFER)) == NULL)
1533     sys_fatal_error ("malloc() failed\n");
1534
1535   while ((c = getc (file)) != EOF)
1536     {
1537
1538       if (i + 1 % BUFFER == 0)
1539         {
1540           if ((line = realloc (line, sizeof (char) * BUFFER * ++j)) == NULL)
1541             sys_fatal_error ("realloc() failed\n");
1542         }
1543
1544       if (c == (int) '\n')
1545         break;
1546
1547       line[i++] = (char) c;
1548     }
1549
1550   if (i == 0)
1551     return NULL;
1552
1553   if ((line = realloc (line, sizeof (char) * i)) == NULL)
1554     sys_fatal_error ("realloc() failed\n");
1555
1556   line[i] = '\0';
1557
1558   return line;
1559 }
1560
1561 /**********************************************************************
1562 read_num:
1563 Reads in the number of sequences that are found in the OC file.
1564
1565 Returns an int
1566
1567 -----------------------------------------------------------------------*/
1568 int
1569 read_num (FILE * infile)
1570 {
1571   int n;
1572   char *line;
1573
1574
1575   /*      regex_t cmpregex;
1576      const char pattern[] = "^ *[0-9]\\+ *$";
1577    */
1578
1579   if ((line = getline_loc (infile)) == NULL)
1580     {
1581       fprintf (stderr, "File empty.\n");
1582       exit (EXIT_FAILURE);
1583     }
1584
1585   /* comment this stuff out as it fails on SunOS
1586
1587      if (regcomp(&cmpregex, pattern, REG_NOSUB) != 0) {
1588      fprintf(stderr, "Regex compilation failed! Check source code for error.\n");
1589      exit(EXIT_FAILURE);
1590      }
1591
1592      if (regexec(&cmpregex, line, 0, NULL, 0) == REG_NOMATCH) {
1593      fprintf(stderr, "File doesn't contain integer number of entries on first line. Line is:\n%s\n", line);
1594      exit(EXIT_FAILURE);
1595      }
1596    */
1597
1598   n = atoi (line);
1599
1600   /*      regfree(&cmpregex); */
1601   free (line);
1602   return n;
1603 }
1604
1605 /**********************************************************************
1606 read_idents:
1607 reads n identifiers from infile
1608
1609 returns a pointer to an array of identifiers
1610 -----------------------------------------------------------------------*/
1611
1612 char **
1613 read_idents (FILE * infile, int n)
1614 {
1615   char **ret_val;
1616   int i;
1617   int error;
1618   char *buff;
1619
1620   buff = GJstrcreate (MAX_ID_LEN, NULL);
1621   ret_val = (char **) GJmalloc (sizeof (char *) * n);
1622
1623   for (i = 0; i < n; ++i)
1624     {
1625       error = fscanf (infile, "%s", buff);
1626       if (error == EOF)
1627         {
1628           fprintf (stderr, "Premature end of file reached in identifiers\n");
1629           exit (EXIT_FAILURE);
1630         }
1631       ret_val[i] = GJstrdup (buff);
1632     }
1633   return ret_val;
1634 }
1635
1636 /**********************************************************************
1637 read_up_diag:
1638 reads a file containing n, the number of entities, followed by
1639 n*(n-1)/2 numbers that are the pair distances or similarities between
1640 the n entities.
1641
1642 Returns an upper diagonal type pointer to pointer array
1643
1644 -----------------------------------------------------------------------*/
1645 double **
1646 read_up_diag (FILE * infile, int n)
1647 {
1648   double **ret_val;
1649   int error;
1650   int i, j;
1651
1652   ret_val = (double **) GJDudarr (n);
1653   for (i = 0; i < (n - 1); ++i)
1654     {
1655       for (j = 0; j < (n - i - 1); ++j)
1656         {
1657           error = fscanf (infile, "%lf", &ret_val[i][j]);
1658           if (error == EOF)
1659             {
1660               fprintf (stderr,
1661                        "Premature end of file reached in upper diagonal\n");
1662               exit (EXIT_FAILURE);
1663             }
1664         }
1665     }
1666   return ret_val;
1667 }