Code reorganization to be less horrible.

LICENSE.txt

@@ -0,0 +1,19 @@
+Copyright (C) 2017 by Benedict Paten (benedictpaten@gmail.com) & Arthur Rand (arand@soe.ucsc.edu)
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.

README

@@ -11,6 +11,9 @@ Features:
 - add code to break up haplotypes where less certain
 - add interface code
 
+Organization:
+- Documentation / readme
+
 Optimization:
 - investigate more effective pruning in the cross product function using the probabilities
     DO PRUNING IN CROSS PRODUCT FUNCTION?

impl/column.c

@@ -0,0 +1,150 @@
+#include "stRPHmm.h"
+
+/*
+ * Read partitioning hmm column (stRPColumn) functions
+ */
+
+stRPColumn *stRPColumn_construct(int64_t refStart, int64_t length, int64_t depth,
+        stProfileSeq **seqHeaders, uint8_t **seqs) {
+    stRPColumn *column = st_calloc(1, sizeof(stRPColumn));
+
+    // Reference coordinates
+    column->refStart = refStart;
+    column->length = length;
+
+    // Sequences
+    column->depth = depth;
+    column->seqHeaders = seqHeaders;
+    column->seqs = seqs;
+
+    // Initially contains not states
+    column->head = NULL;
+
+    return column;
+}
+
+void stRPColumn_destruct(stRPColumn *column) {
+    // Clean up the contained cells
+    stRPCell *cell = column->head;
+    while(cell != NULL) {
+        stRPCell *pCell = cell;
+        cell = cell->nCell;
+        stRPCell_destruct(pCell);
+    }
+
+    free(column->seqHeaders);
+    free(column->seqs);
+    free(column);
+}
+
+void stRPColumn_print(stRPColumn *column, FILE *fileHandle, bool includeCells) {
+    /*
+     * Print a description of the column. If includeCells is true then print the
+     * state of the cells too.
+     */
+    fprintf(fileHandle, "\tCOLUMN: REF_START: %" PRIi64
+            " REF_LENGTH: %" PRIi64 " DEPTH: %" PRIi64
+            " TOTAL_PROB: %f\n",
+            column->refStart, column->length, column->depth,
+            (float)column->totalLogProb);
+    for(int64_t i=0; i<column->depth; i++) {
+        stProfileSeq_print(column->seqHeaders[i], fileHandle, 0);
+    }
+    if(includeCells) {
+        stRPCell *cell = column->head;
+        do {
+            fprintf(fileHandle, "\t\t");
+            stRPCell_print(cell, fileHandle);
+        } while((cell = cell->nCell) != NULL);
+    }
+}
+
+void stRPColumn_split(stRPColumn *column, int64_t firstHalfLength, stRPHmm *hmm) {
+    /*
+     * Split the column into two to divide into two smaller reference intervals
+     */
+
+    // Create column
+    stProfileSeq **seqHeaders = st_malloc(sizeof(stProfileSeq *) * column->depth);
+    memcpy(seqHeaders, column->seqHeaders, sizeof(stProfileSeq *) * column->depth);
+    uint8_t **seqs = st_malloc(sizeof(uint8_t *) * column->depth);
+    // Update the pointers to the seqs
+    for(int64_t i=0; i<column->depth; i++) {
+        seqs[i] = &(column->seqs[i][firstHalfLength * ALPHABET_SIZE]);
+    }
+    assert(firstHalfLength > 0); // Non-zero length for first half
+    assert(column->length-firstHalfLength > 0); // Non-zero length for second half
+    stRPColumn *rColumn = stRPColumn_construct(column->refStart+firstHalfLength,
+            column->length-firstHalfLength, column->depth, seqHeaders, seqs);
+
+    // Create merge column
+    uint64_t acceptMask = makeAcceptMask(column->depth);
+    stRPMergeColumn *mColumn = stRPMergeColumn_construct(acceptMask, acceptMask);
+
+    // Copy cells
+    stRPCell *cell = column->head;
+    stRPCell **pCell = &(rColumn->head);
+    do {
+        *pCell = stRPCell_construct(cell->partition);
+        stRPMergeCell_construct(cell->partition, cell->partition, mColumn);
+        pCell = &((*pCell)->nCell);
+    } while((cell = cell->nCell) != NULL);
+
+    // Create links
+    rColumn->pColumn = mColumn;
+    mColumn->nColumn = rColumn;
+
+    // If is the last column
+    if(column->nColumn == NULL) {
+       assert(hmm->lastColumn == column);
+       hmm->lastColumn = rColumn;
+       assert(rColumn->nColumn == NULL);
+    }
+    else {
+        column->nColumn->pColumn = rColumn;
+        rColumn->nColumn = column->nColumn;
+    }
+    column->nColumn = mColumn;
+    mColumn->pColumn = column;
+
+    // Increase column number
+    hmm->columnNumber++;
+
+    // Adjust length of previous column
+    column->length = firstHalfLength;
+}
+
+/*
+ * Read partitioning hmm state (stRPCell) functions
+ */
+
+stRPCell *stRPCell_construct(int64_t partition) {
+    stRPCell *cell = st_calloc(1, sizeof(stRPCell));
+    cell->partition = partition;
+    return cell;
+}
+
+void stRPCell_destruct(stRPCell *cell) {
+    free(cell);
+}
+
+void stRPCell_print(stRPCell *cell, FILE *fileHandle) {
+    /*
+     * Prints a debug representation of the cell.
+     */
+    char *partitionString = intToBinaryString(cell->partition);
+    fprintf(fileHandle, "CELL PARTITION: %s FORWARD_PROB: %f BACKWARD_PROB: %f\n",
+            partitionString, (float)cell->forwardLogProb, (float)cell->backwardLogProb);
+    free(partitionString);
+}
+
+double stRPCell_posteriorProb(stRPCell *cell, stRPColumn *column) {
+    /*
+     * Calculate the posterior probability of visiting the given cell. Requires that the
+     * forward and backward algorithms have been run.
+     */
+    double p = exp(cell->forwardLogProb + cell->backwardLogProb - column->totalLogProb);
+    assert(p <= 1.1);
+    assert(p >= 0.0);
+    return p > 1.0 ? 1.0 : p;
+}