92f77f7bef111331d3c88d14cf98ef2e64c765f8
kent
  Tue Dec 21 12:06:48 2021 -0800
Improved a structure name and removed some debugging.  Added some comments.

diff --git src/utils/matrixClusterColumns/matrixClusterColumns.c src/utils/matrixClusterColumns/matrixClusterColumns.c
index 2ba453a..d0f47bd 100644
--- src/utils/matrixClusterColumns/matrixClusterColumns.c
+++ src/utils/matrixClusterColumns/matrixClusterColumns.c
@@ -1,585 +1,592 @@
 /* matrixClusterColumns - Group the columns of a matrix into clusters, and output a matrix 
  * the with same number of rows and generally much fewer columns.. */
 #include "common.h"
 #include "linefile.h"
 #include "hash.h"
 #include "options.h"
 #include "obscure.h"
 #include "fieldedTable.h"
 #include "sqlNum.h"
 #include "pthreadDoList.h"
 
 #define clThreadCount 50
 #define chunkItemsPerThread 5
 #define chunkMaxSize (clThreadCount * chunkItemsPerThread)
 
 void usage()
 /* Explain usage and exit. */
 {
 errAbort(
   "matrixClusterColumns - Group the columns of a matrix into clusters, and output a matrix with\n"
   "the same number of rows and generally much fewer columns. Combines columns by taking mean.\n"
   "usage:\n"
   "   matrixClusterColumns inMatrix.tsv meta.tsv cluster outMatrix.tsv outStats.tsv [cluster2 outMatrix2.tsv outStats2.tsv ... ]\n"
   "where:\n"
   "   inMatrix.tsv is a file in tsv format with cell labels in first row and gene labels in first column\n"
   "   meta.tsv is a table where the first row is field labels and the first column is sample ids\n"
   "   cluster is the name of the field with the cluster names\n"
   "You can produce multiple clusterings in the same pass through the input matrix by specifying\n"
   "additional cluster/outMatrix/outStats triples in the command line.\n"
   "options:\n"
   "   -makeIndex=index.tsv - output index tsv file with <matrix-col1><input-file-pos><line-len>\n"
   "   -median if set ouput median rather than mean cluster value\n"
   );
 }
 
 /* Command line validation table. */
 static struct optionSpec options[] = {
    {"makeIndex", OPTION_STRING},
    {"median", OPTION_BOOLEAN},
    {NULL, 0},
 };
 
 
 void readLineArray(char *fileName, int *retCount, char ***retLines)
 /* Return an array of strings, one for each line of file.  Return # of lines in file too */
 {
 /* This is sloppy with memory but it doesn't matter since we won't free it. */
 struct slName *el, *list = readAllLines(fileName);
 if (list == NULL)
     errAbort("%s is empty", fileName);
 int count = slCount(list);
 char **lines;
 AllocArray(lines, count);
 int i;
 for (i=0, el=list; i<count; ++i, el = el->next)
     {
     lines[i] = el->name;
     }
 *retCount = count;
 *retLines = lines;
 }
 
 int countNonzero(double *a, int size)
 /* Return number of nonzero items in array a */
 {
 int count = 0;
 while (--size >= 0)
     if (*a++ != 0.0)
         ++count;
 return count;
 }
 
 double sumArray(double *a, int size)
 /* Return sum of all items in array */
 {
 double sum = 0.0;
 while (--size >= 0)
     sum += *a++;
 return sum;
 }
 
 struct ccMatrix
 /* Local matrix structure - little wrapper around a fielded table/lineFile combination */
     {
     struct ccMatrix *next;
 
         /* kind of private fields */
     struct lineFile *lf;	    // Line file for tab-sep case
     struct fieldedTable *ft;	    // fielded table if a tab-sep file
 
 	/* From below are fields that yu can read but not change */
     int colCount;		    // Number of columns in a row
     char **colLabels;		    // A label for each column 
     int curRow;			    // Current row we are processing
     };
 
 struct ccMatrix *ccMatrixOpen(char *matrixFile)
 /* Local helper matrix structure.  Could be simplified */
 /* Figure out if it's a mtx or tgz file and open it */
 {
 /* Read in labels if there are any */
 struct ccMatrix *v = needMem(sizeof(*v));
 struct lineFile *lf = v->lf = lineFileOpen(matrixFile, TRUE);
 struct fieldedTable *ft = v->ft = fieldedTableReadTabHeader(lf, NULL, 0);
 v->colCount = ft->fieldCount-1;	    // Don't include row label field 
 v->colLabels = ft->fields+1;	// +1 to skip over row label field
 return v;
 }
 
 void ccMatrixClose(struct ccMatrix **pV)
 /* Close up ccMatrix */
 {
 struct ccMatrix *v = *pV;
 if (v != NULL)
     {
     fieldedTableFree(&v->ft);
     freez(pV);
     }
 }
 
 struct clustering
 /* Stuff we need to cluster something.  This is something we might do
  * repeatedly to same input matrix */
     {
     struct clustering *next;
     char *clusterField;	    /* Field to cluster on */
     char *outMatrixFile;    /* Where to put matrix result */
     char *outStatsFile;	    /* Where to put stats result */
     int clusterMetaIx;	    /* Index of cluster field in meta table */
     int *colToCluster;	    /* Maps input matrix column to clustered column */
     int clusterCount;   /* Number of different values in column we are clustering */
     double *clusterTotal;  /* A place to hold totals for each cluster */
     double *clusterGrandTotal;	/* Total over all rows */
     int *clusterElements;  /* A place to hold counts for each cluster */
     double *clusterResults; /* Results after clustering */
     struct hash *clusterSizeHash;   /* Keyed by cluster, int value is elements in cluster */
     char **clusterNames;	    /* Holds name of each cluster */
     int *clusterSizes;	    /* An array that holds size of each cluster */
 
     /* Things needed by median handling */
     boolean doMedian;	/* If true we calculate median */
     double **clusterSamples; /* An array that holds an array big enough for all vals in cluster. */
 
     struct dyString **chunkLinesOut; /* parallel output */
     double **chunkSubtotals;	     /* Totals */
     FILE *matrixFile;		     /* serial output */
     };
 
 
 struct clustering *clusteringNew(char *clusterField, char *outMatrixFile, char *outStatsFile,
     struct fieldedTable *metaTable, struct ccMatrix *v, boolean doMedian)
 /* Make up a new clustering structure */
 {
 /* Check that all column names in matrix are unique */
 int colCount = v->colCount;
 char **colLabels = v->colLabels;
 struct hash *uniqColHash = hashNew(0);
 int colIx;
 for (colIx=0; colIx < colCount; colIx = colIx+1)
     {
     char *label = colLabels[colIx];
     if (hashLookup(uniqColHash, label) == NULL)
 	hashAdd(uniqColHash, label, NULL);
     else
         errAbort("Duplicated column label %s in input matrix", label);
     }
 
 struct clustering *clustering;
 AllocVar(clustering);
 clustering->clusterField = clusterField;
 clustering->outMatrixFile = outMatrixFile;
 clustering->outStatsFile = outStatsFile;
 int clusterFieldIx = clustering->clusterMetaIx = fieldedTableMustFindFieldIx(metaTable, clusterField);
 
 /* Make up hash of sample names with cluster name values 
  * and also hash of cluster names with size values */
 struct hash *sampleHash = hashNew(0);	/* Keyed by sample value is cluster */
 struct hash *clusterSizeHash = clustering->clusterSizeHash = hashNew(0);
 struct fieldedRow *fr;
 for (fr = metaTable->rowList; fr != NULL; fr = fr->next)
     {
     char **row = fr->row;
     char *sample = row[0];
     if (!hashLookup(uniqColHash, sample))
         errAbort("%s is in %s but not input matrix", sample, metaTable->name);
 
     hashAdd(sampleHash, sample, row[clusterFieldIx]);
     hashIncInt(clusterSizeHash, row[clusterFieldIx]);
     }
 
 /* Find all uniq cluster names */
 struct slName *nameList = NULL;
 struct hash *uniqHash = hashNew(0);
 for (fr = metaTable->rowList; fr != NULL; fr = fr->next)
     {
     char *cluster = fr->row[clusterFieldIx];
     if (hashLookup(uniqHash, cluster) == NULL)
         {
 	slNameAddHead(&nameList, cluster);
 	hashAdd(uniqHash, cluster, NULL);
 	}
     }
 hashFree(&uniqHash);
 
 /* Just alphabetize names for now */
 slNameSort(&nameList);
 slSort(&nameList, slNameCmpWordsWithEmbeddedNumbers);
 
 /* Make up hash that maps cluster names to cluster ids */
 struct hash *clusterIxHash = hashNew(0);	/* Keyed by cluster, no value */
 int i;
 struct slName *name;
 for (name = nameList, i=0; name != NULL; name = name->next, ++i)
     hashAddInt(clusterIxHash, name->name, i);
 int clusterCount = clustering->clusterCount = clusterIxHash->elCount;
 
 /* Make up array that holds size of each cluster */
 AllocArray(clustering->clusterSizes, clusterCount);
 AllocArray(clustering->clusterNames, clusterCount);
 for (i = 0, name = nameList; i < clusterCount; ++i, name = name->next)
     {
     clustering->clusterSizes[i] = hashIntVal(clustering->clusterSizeHash, name->name);
     clustering->clusterNames[i] = name->name;
     verbose(2, "clusterSizes[%d] = %d\n", i, clustering->clusterSizes[i]);
     }
 
 if (doMedian)
     {	
     /* Allocate arrays to hold number of samples and all sample vals for each cluster */
     clustering->doMedian = doMedian;
     AllocArray(clustering->clusterSamples, clusterCount);
     int clusterIx;
     for (clusterIx = 0; clusterIx < clusterCount; ++clusterIx)
 	{
 	double *samples;
 	AllocArray(samples, clustering->clusterSizes[clusterIx]);
 	clustering->clusterSamples[clusterIx] = samples;
 	}
     }
 
 
 /* Make up array that has -1 where no cluster available, otherwise output index, also
  * hash up all column labels. */
 int *colToCluster = clustering->colToCluster = needHugeMem(colCount * sizeof(colToCluster[0]));
 int unclusteredColumns = 0, missCount = 0;
 for (colIx=0; colIx < colCount; colIx = colIx+1)
     {
     char *colName = colLabels[colIx];
     char *clusterName = hashFindVal(sampleHash, colName);
     colToCluster[colIx] = -1;
     if (clusterName != NULL)
         {
 	int clusterId = hashIntValDefault(clusterIxHash, clusterName, -1);
 	colToCluster[colIx] = clusterId;
 	if (clusterId == -1)
 	    {
 	    verbose(3, "%s is in expression matrix but not in sample cluster file", clusterName);
 	    ++missCount;
 	    }
 	}
     else
 	unclusteredColumns += 1;
     }
 verbose(1, "%d total columns, %d unclustered, %d misses\n", 
     colCount, unclusteredColumns, missCount);
 
 /* Allocate space for results for clustering one line */
 clustering->clusterResults = needHugeMem(clusterCount * sizeof(clustering->clusterResults[0]));
 
 /* Allocate a few more things */
 clustering->clusterTotal = needMem(clusterCount*sizeof(clustering->clusterTotal[0]));
 clustering->clusterGrandTotal = needMem(clusterCount*sizeof(clustering->clusterGrandTotal[0]));
 clustering->clusterElements = needMem(clusterCount*sizeof(clustering->clusterElements[0]));
 
 /* Open file - and write out header */
 FILE *f = clustering->matrixFile = mustOpen(clustering->outMatrixFile, "w");
 if (v->ft->startsSharp)
     fputc('#', f);
 
 /* First field name agrees with first column of matrix */
 fputs( v->ft->fields[0],f);
 
 /* Use our clusters for the rest of the names */
 for (name = nameList; name != NULL; name = name->next) 
     {
     fputc('\t', f);
     fputs(name->name, f);
     }
 fputc('\n', f);
 
 /* Allocate parallel output buffers */
 AllocArray(clustering->chunkSubtotals, chunkMaxSize);
 for (i=0; i<chunkMaxSize; ++i)
     {
     AllocArray(clustering->chunkSubtotals[i], clusterCount);
     }
 AllocArray(clustering->chunkLinesOut, chunkMaxSize);
 for (i=0; i<chunkMaxSize; ++i)
     clustering->chunkLinesOut[i] = dyStringNew(colCount * 3);  
 
 /* Clean up and return result */
 hashFree(&sampleHash);
 hashFree(&clusterIxHash);
 return clustering;
 }
 
 void outputClusterStats(struct clustering *clustering)
 /* Output statistics on each cluster in this clustering. */
 {
 FILE *f = mustOpen(clustering->outStatsFile, "w");
 fprintf(f, "#cluster\tcount\ttotal\n");
 int i;
 for (i=0; i<clustering->clusterCount; ++i)
     {
     fprintf(f, "%s\t%d\t%g\n", clustering->clusterNames[i], 
 	clustering->clusterSizes[i], clustering->clusterGrandTotal[i]);
     }
 carefulClose(&f);
 }
 
 void clusterRow(struct clustering *clustering, int colCount, char *rowLabel, 
     double *a, double *totalingTemp, int *elementsTemp, struct dyString *out,
     double *subtotals)
 /* Process a row in a, outputting in clustering->file */
 {
 /* Zero out cluster histogram */
 double *clusterTotal = totalingTemp;
 int *clusterElements = elementsTemp;
 int clusterCount = clustering->clusterCount;
 int i;
 for (i=0; i<clusterCount; ++i)
     {
     clusterTotal[i] = 0.0;
     clusterElements[i] = 0;
     }
 
 /* Loop through rest of row filling in histogram */
 int *colToCluster = clustering->colToCluster;
 boolean doMedian = clustering->doMedian;
 for (i=0; i<colCount; ++i)
     {
     int clusterIx = colToCluster[i];
     if (clusterIx >= 0)
 	{
 	double val = a[i];
 	int valCount = clusterElements[clusterIx];
 	clusterElements[clusterIx] = valCount+1;
 	clusterTotal[clusterIx] += val;
 	if (doMedian)
 	    {
 	    if (valCount >= clustering->clusterSizes[clusterIx])
 		internalErr();
 	    clustering->clusterSamples[clusterIx][valCount] = val;
 	    }
 	}
     }
 
 /* Do output to outstrng and do grand totalling */
 dyStringClear(out);
 dyStringPrintf(out, "%s", rowLabel);
 for (i=0; i<clusterCount; ++i)
     {
     dyStringAppendC(out, '\t');
     double total = clusterTotal[i];
     subtotals[i] += total;
     int elements = clusterElements[i];
     double val;
     if (doMedian)
 	{
 	val = doubleMedian(elements, clustering->clusterSamples[i]);
 	dyStringPrintf(out, "%g", val);
 	}
     else
 	{
 	if (total > 0)
 	    {
 	    val = total/elements;
 	    dyStringPrintf(out, "%g", val);
 	    }
 	else
 	    dyStringAppendC(out, '0');
 	}
     }
 dyStringAppendC(out, '\n');
 }
 
 
 static void addRowToIndex(FILE *fIndex, char *rowLabel, struct lineFile *lf)
 /* Write out info to index file about where this row begins */
 {
 if (fIndex)
   {
   fprintf(fIndex, "%s", rowLabel);
   fprintf(fIndex, "\t%lld\t%lld\n",  (long long)lineFileTell(lf), (long long)lineFileTellSize(lf));
   }
 }
 
-struct lineIoInfo
-/* Enough rows to do things in parallel we hope? */
+struct lineIoItem
+/* This is an item fed to a parallel worker.  It corresponds to a single line of
+ * input matrix */
     {
-    struct lineIoInfo *next;
+    struct lineIoItem *next;	/* Pointer to next in list */
+
+    /* Information about input file and where we are in it. */
     char *fileName;
-    struct clustering *clusteringList;
-    int lineIx;	    /* Index of line in file */
-    int chunkIx;  /* Index of line in chunk */
+    int lineIx;	    /* Index of line in input file */
     long long lineStartOffset;	/* Start offset within file */
     long long lineSize;	/* Size of line */
     long long lineEndOffset;	
+    int chunkIx;  /* Index of line in chunk */
+
+    /* Slightly parsed input. */
     struct dyString *rowLabel;	/* Just the row label of input */
     struct dyString *lineIn;     /* Unparsed rest of input line */
-    double *vals;		 /* Array of values parsed from string  */
-    double *totalingTemp;        /* Buffer for parllel computation */
+
+    struct clustering *clusteringList;  /* Instructions on how to cluster and output */
+    
+    /* Temporary values used for calculating output */
+    double *vals;		 /* Parse out input matrix values for this line */
+    double *totalingTemp;        /* Buffer for parallel computation of totals */
     int *elementsTemp;           /* buffers for parallel computation */
     };
 
 void lineWorker(void *item, void *context)
 /* A worker to execute a single column clustering  */
 {
-struct lineIoInfo *lii = item;
+struct lineIoItem *lii = item;
 struct ccMatrix *v = context;
 int xSize = v->colCount;
 char *s = lii->lineIn->string;
 char *rowLabel = lii->rowLabel->string;;
 
 /* Convert ascii to floating point, with little optimization for the many zeroes we usually see */
 int i;
 double *vals = lii->vals;
 for (i=0; i<xSize; ++i)
     {
     char *str = nextTabWord(&s);
     if (str == NULL)
         errAbort("not enough fields in input matrix line %d", lii->lineIx);
     double val = ((str[0] == '0' && str[1] == 0) ? 0.0 : sqlDouble(str));
     vals[i] = val;
     }
 
 /* Then do the clustering */
 struct clustering *clustering;
 for (clustering = lii->clusteringList; clustering != NULL; clustering = clustering->next)
       {
       int chunkIx = lii->chunkIx;
       clusterRow(clustering, xSize, rowLabel, 
 	    lii->vals, lii->totalingTemp, lii->elementsTemp,
 	    clustering->chunkLinesOut[chunkIx], clustering->chunkSubtotals[chunkIx]);
       }
 }
 
 void matrixClusterColumns(char *matrixFile, char *metaFile, char *sampleField,
     int outputCount, char **clusterFields, char **outMatrixFiles, char **outStatsFiles,
     char *outputIndex, boolean doMedian)
 /* matrixClusterColumns - Group the columns of a matrix into clusters, and output a matrix 
  * the with same number of rows and generally much fewer columns.. */
 {
 FILE *fIndex = NULL;
 if (outputIndex)
     fIndex = mustOpen(outputIndex, "w");
 
 /* Load up metadata and make sure we have all of the cluster fields we need 
  * and fill out array of clusterIx corresponding to clusterFields in metaFile. */
 struct fieldedTable *metaTable = fieldedTableFromTabFile(metaFile, metaFile, NULL, 0);
 struct hash *metaHash = fieldedTableIndex(metaTable, sampleField);
 verbose(1, "Read %d rows from %s\n", metaHash->elCount, metaFile);
 
 /* Load up input matrix first line at least */
 struct ccMatrix *v = ccMatrixOpen(matrixFile);
 verbose(1, "matrix %s has %d fields\n", matrixFile, v->colCount);
 
 /* Create a clustering for each output and find index in metaTable for each. */
 struct clustering *clusteringList = NULL, *clustering;
 int i;
 for (i=0; i<outputCount; ++i)
     {
     clustering = clusteringNew(clusterFields[i], outMatrixFiles[i], outStatsFiles[i], 
 			metaTable, v, doMedian);
     slAddTail(&clusteringList, clustering);
     }
 
 /* Set up buffers for pthread workers */
-struct lineIoInfo chunks[chunkMaxSize];
+struct lineIoItem chunks[chunkMaxSize];
 for (i=0; i<chunkMaxSize; ++i)
     {
-    struct lineIoInfo *chunk = &chunks[i];
+    struct lineIoItem *chunk = &chunks[i];
     chunk->fileName = matrixFile;
     chunk->clusteringList = clusteringList;
     chunk->chunkIx = i;
     chunk->lineIn = dyStringNew(0);
     chunk->rowLabel = dyStringNew(0);
     AllocArray(chunk->vals, v->colCount);
     AllocArray(chunk->totalingTemp, v->colCount);
     AllocArray(chunk->elementsTemp, v->colCount);
     }
 
 /* Chug through big matrix a row at a time clustering */
 dotForUserInit(1);
 
 boolean atEof = FALSE;
 struct lineFile *lf = v->lf;
-uglyf("Starting main loop on %d columns\n", v->colCount);
 while (!atEof)
     {
     /* Read a chunk of lines of the file */
-    struct lineIoInfo *chunkList = NULL, *chunk;
+    struct lineIoItem *chunkList = NULL, *chunk;
     int chunkSize;
     for (chunkSize = 0; chunkSize < chunkMaxSize; chunkSize += 1)
 	{
 	char *line;
 	if (!lineFileNextReal(lf, &line))
 	   {
 	   atEof = TRUE;
 	   break;
 	   }
 	chunk = &chunks[chunkSize];
 	chunk->lineIx = lf->lineIx;
 	char *rowLabel = nextTabWord(&line);
 	addRowToIndex(fIndex, rowLabel, lf);
 	dyStringClear(chunk->rowLabel);
 	dyStringAppend(chunk->rowLabel, rowLabel);
 	dyStringClear(chunk->lineIn);
 	dyStringAppend(chunk->lineIn, line);
 	slAddHead(&chunkList, chunk);
 	}
     slReverse(&chunkList);
 
     /* Calculate strings to print in parallel */
     pthreadDoList(clThreadCount, chunkList, lineWorker, v);
 
     /* Do the actual file writes in serial and add subtotals to grand total */
     for (chunk = chunkList; chunk != NULL; chunk = chunk->next)
 	{
 	struct clustering *clustering;
 	for (clustering = clusteringList; clustering != NULL; clustering = clustering->next)
 	    {
 	    fputs(clustering->chunkLinesOut[chunk->chunkIx]->string, clustering->matrixFile);
 
 	    /* Add subtotals calculated in parallel to grandTotal. */
 	    int i;
 	    double *grandTotal = clustering->clusterGrandTotal;
 	    double *subtotal = clustering->chunkSubtotals[chunk->chunkIx];
 	    for (i=0; i<clustering->clusterCount; ++i)
 		{
 	        grandTotal[i] += subtotal[i];
 		subtotal[i] = 0;
 		}
 	    }
 	}
     dotForUser();
     }
 
 
 dotForUserEnd();
 
 /* Do stats and close files */
 for (clustering = clusteringList; clustering != NULL; clustering = clustering->next)
     {
     outputClusterStats(clustering);
     carefulClose(&clustering->matrixFile);
     }
 
 ccMatrixClose(&v);
 carefulClose(&fIndex);
 }
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 optionInit(&argc, argv, options);
 char *makeIndex = optionVal("makeIndex", NULL);
 int minArgc = 6;
 if (argc < minArgc || ((argc-minArgc)%3)!=0)  // Force minimum, even number
     usage();
 int outputCount = 1 + (argc-minArgc)/3;	      // Add one since at minimum have 1
 char *clusterFields[outputCount], *outMatrixFiles[outputCount], *outStatsFiles[outputCount];
 int i;
 char **triples = argv + minArgc - 3;
 for (i=0; i<outputCount; ++i)
     {
     clusterFields[i] = triples[0];
     outMatrixFiles[i] = triples[1];
     outStatsFiles[i] = triples[2];
     triples += 3;
     }
 matrixClusterColumns(argv[1], argv[2], argv[3], 
     outputCount, clusterFields, outMatrixFiles, outStatsFiles, makeIndex, optionExists("median"));
 return 0;
 }