//
// Copyright (C) 2009-2012 Institute for Computational Biomedicine,
//                         Weill Medical College of Cornell University
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//

#include <string>
#include <iostream>
#define PCRE2_CODE_UNIT_WIDTH 8
#include <pcre2.h>
#include <stdio.h>

#include "Reads.h"
#include "C_Alignments.h"
#include "MessageChunks.h"
#include "SamFlags.h"
#include "C_Gsnap.h"

using namespace std;

#undef C_WRITE_API_WRITE_ALIGNMENT_DEBUG
#ifdef C_WRITE_API_WRITE_ALIGNMENT_DEBUG
#define debug(x) x
#else
#define debug(x)
#endif

#ifdef __GNUC__
#define GOBY_CAPTURE_VIA_OPEN_MEMSTREAM
#endif

// TODO: When reading from COMPACT-READS, one should call
// TODO: addQueryIdentifierWithInt() but ONLY if there is a text
// TODO: identifier for the read in question, otherwise just
// TODO: use the actual queryIndex that comes from CompactReads
// TODO: as the query index and don't bother with addQueryIdentifier*
// TODO: the addQueryIdentifier() should only be used when reading
// TODO: from FA/FQ and what you have is a real string identifier.


/**
 * C API to enable writing Goby compact-alignments in C.
 */
extern "C" {

	void gobyAlignments_openAlignmentsWriterDefaultEntriesPerChunk(char *basename, CAlignmentsWriterHelper **writerHelperpp) {
	    gobyAlignments_openAlignmentsWriter(basename, GOBY_DEFAULT_NUMBER_OF_ENTRIES_PER_CHUNK, writerHelperpp);
	}

	void gobyAlignments_openAlignmentsWriter(char *basename, unsigned int number_of_entries_per_chunk, CAlignmentsWriterHelper **writerHelperpp) {
        debug(fprintf(stderr,"Writing alignment to basename, entries per chunk=%u\n", basename, number_of_entries_per_chunk));
        *writerHelperpp = new CAlignmentsWriterHelper;
        CAlignmentsWriterHelper *writerHelper = *writerHelperpp;
        string basenameStr(basename);
        
        writerHelper->alignmentWriter = new goby::AlignmentWriter(basenameStr, number_of_entries_per_chunk);
        writerHelper->tmhWriter = new goby::TooManyHitsWriter(basenameStr, 1);
        writerHelper->alignmentEntry = NULL;
        writerHelper->sequenceVariation = NULL;
        writerHelper->lastSeqVarRefPosition = -1;
        writerHelper->lastSeqVarReadChar = '\0';
        writerHelper->lastSeqVarRefChar = '\0';
        writerHelper->smallestQueryIndex = -1;
        writerHelper->largestQueryIndex = -1;
        writerHelper->numberOfAlignedReads = 0;
        writerHelper->numberOfMisParsedReads = 0;
        // Adjust this if the incoming quality scores are
        // Illumina (-64) or Sanger (-33), otherwise Phred
        // scores (0-based) will be assumed.
        writerHelper->qualityAdjustment = 0;
        writerHelper->samHelper = NULL;
        writerHelper->captureFile = NULL;
        writerHelper->captureBuffer = NULL;
        writerHelper->captureBufferSize = 0;
        writerHelper->captureIgnoredFile = NULL;
        writerHelper->captureIgnoredBuffer = NULL;
        writerHelper->captureIgnoredBufferSize = 0;
        writerHelper->alignerToGobyTargetIndexMap = NULL;
        writerHelper->alignmentWriter->setQueryLengthsStoredInEntries(true);
        writerHelper->gsnapAlignment = NULL;
    }

    void gobyCapture_open(CAlignmentsWriterHelper *writerHelper, int openIgnoredFile) {
        gobyCapture_close(writerHelper);
#ifdef GOBY_CAPTURE_VIA_OPEN_MEMSTREAM
        debug(fprintf(stderr, "Opening gobyCapture via open_memstream()\n"));
        writerHelper->captureFile = open_memstream(&writerHelper->captureBuffer, &writerHelper->captureBufferSize);
        if (openIgnoredFile) {
            writerHelper->captureIgnoredFile = open_memstream(&writerHelper->captureIgnoredBuffer, &writerHelper->captureIgnoredBufferSize);
        }
#else
        debug(fprintf(stderr, "Opening gobyCapture via temporary files, slower but non-gcc compatible.\n"));
        writerHelper->captureFile = tmpfile();
        if (openIgnoredFile) {
            writerHelper->captureIgnoredFile = tmpfile();
        }
#endif
    }

    FILE *gobyCapture_fileHandle(CAlignmentsWriterHelper *writerHelper) {
        return writerHelper->captureFile;
    }

    FILE *gobyCapture_ignoredFileHandle(CAlignmentsWriterHelper *writerHelper) {
        return writerHelper->captureIgnoredFile;
    }

    /**
     * Start a new section of output.
     */
    void gobyCapture_startNew(CAlignmentsWriterHelper *writerHelper) {
        if (writerHelper->captureFile) {
            rewind(writerHelper->captureFile);
        }
        if (writerHelper->captureIgnoredFile) {
            rewind(writerHelper->captureIgnoredFile);
        }
    }

    char *gobyCapture_capturedData(CAlignmentsWriterHelper *writerHelper) {
        return writerHelper->captureBuffer;
    }

    char *gobyCapture_ignoredData(CAlignmentsWriterHelper *writerHelper) {
        return writerHelper->captureIgnoredBuffer;
    }

    /**
     * A section of output is complete. Populate what was written to the files into
     * writerHelper->captureBuffer / writerHelper->captureIgnoredBuffer.
     * The output should be processed. After you're done processing the output,
     * when you're ready to start writing new output call gobyCapture_startNew().
     */
    void gobyCapture_flush(CAlignmentsWriterHelper *writerHelper) {
	    if (writerHelper->captureFile) {
            fflush(writerHelper->captureFile);
        }
        if (writerHelper->captureIgnoredFile) {
            fflush(writerHelper->captureIgnoredFile);
        }
#ifndef GOBY_CAPTURE_VIA_OPEN_MEMSTREAM
        size_t fileSize;
        if (writerHelper->captureFile) {
            fileSize = ftell(writerHelper->captureFile);
            if (fileSize + 1 > writerHelper->captureBufferSize) {
                writerHelper->captureBufferSize = fileSize + 1;
                writerHelper->captureBuffer = (char *) realloc(writerHelper->captureBuffer, writerHelper->captureBufferSize);
            }
            rewind(writerHelper->captureFile);
            fread(writerHelper->captureBuffer, 1, fileSize, writerHelper->captureFile);
            writerHelper->captureBuffer[fileSize] = '\0';
        }
        if (writerHelper->captureIgnoredFile) {
            fileSize = ftell(writerHelper->captureIgnoredFile);
            if (fileSize + 1 > writerHelper->captureIgnoredBufferSize) {
                writerHelper->captureIgnoredBufferSize = fileSize + 1;
                writerHelper->captureIgnoredBuffer = (char *) realloc(writerHelper->captureIgnoredBuffer, writerHelper->captureIgnoredBufferSize);
            }
            rewind(writerHelper->captureIgnoredFile);
            fread(writerHelper->captureIgnoredBuffer, 1, fileSize, writerHelper->captureIgnoredFile);
            writerHelper->captureIgnoredBuffer[fileSize] = '\0';
        }
#endif
    }

    void gobyCapture_close(CAlignmentsWriterHelper *writerHelper) {
        if (writerHelper->captureFile || writerHelper->captureIgnoredFile) {
#ifdef GOBY_CAPTURE_VIA_OPEN_MEMSTREAM
            debug(fprintf(stderr, "Closing gobyCapture via open_memstream()\n");)
#else
            debug(fprintf(stderr, "Closing gobyCapture via temporary files\n");)
#endif
        }
        if (writerHelper->captureFile) {
            fclose(writerHelper->captureFile);
            writerHelper->captureFile = NULL;
            free(writerHelper->captureBuffer);
            writerHelper->captureBuffer = NULL;
            writerHelper->captureBufferSize = 0;
        }
        if (writerHelper->captureIgnoredFile) {
            fclose(writerHelper->captureIgnoredFile);
            writerHelper->captureIgnoredFile = NULL;
            free(writerHelper->captureIgnoredBuffer);
            writerHelper->captureIgnoredBuffer = NULL;
            writerHelper->captureIgnoredBufferSize = 0;
        }
    }

    int gobyAlignments_getQualityAdjustment(CAlignmentsWriterHelper *writerHelper) {
        return writerHelper->qualityAdjustment;
    }
    void gobyAlignments_setQualityAdjustment(CAlignmentsWriterHelper *writerHelper, int value) {
        debug(fprintf(stderr,"gobyAlEntry_setQualityAdjustment=%d\n", value));
        writerHelper->qualityAdjustment = value;
    }
    void gobyAlignments_setAlignerName(CAlignmentsWriterHelper *writerHelper, char *value) {
        debug(fprintf(stderr,"gobyAlEntry_setAlignerName=%s\n", value));
        string valueStr(value);
        writerHelper->alignmentWriter->setAlignerName(valueStr);
    }
    void gobyAlignments_setAlignerVersion(CAlignmentsWriterHelper *writerHelper, char *value) {
        debug(fprintf(stderr,"gobyAlEntry_setAlignerVersion=%s\n", value));
        string valueStr(value);
        writerHelper->alignmentWriter->setAlignerVersion(valueStr);
    }
    void gobyAlignments_setSorted(CAlignmentsWriterHelper *writerHelper, int sorted /* bool */) {
        debug(fprintf(stderr,"gobyAlignments_setSorted=%d\n", sorted));
        writerHelper->alignmentWriter->setSorted(sorted == 0 ? false : true);
    }
    void gobyAlignments_setIndexed(CAlignmentsWriterHelper *writerHelper, int indexed /* bool */) {
        debug(fprintf(stderr,"gobyAlignments_setIndexed=%d\n", indexed));
        writerHelper->alignmentWriter->setIndexed(indexed == 0 ? false : true);
    }
    void gobyAlignments_addStatisticStr(CAlignmentsWriterHelper *writerHelper, const char *description, const char *value) {
        debug(fprintf(stderr,"gobyAlignments_addStatisticStr %s=%s\n", description, value));
        string descriptionStr(description);
        string valueStr(value);
        writerHelper->alignmentWriter->addStatistic(descriptionStr, valueStr);
    }
    void gobyAlignments_addStatisticInt(CAlignmentsWriterHelper *writerHelper, const char *description, const int value) {
        debug(fprintf(stderr,"gobyAlignments_addStatisticInt %s=%d\n", description, value));
        string descriptionStr(description);
        writerHelper->alignmentWriter->addStatistic(descriptionStr, value);
    }
    void gobyAlignments_addStatisticDouble(CAlignmentsWriterHelper *writerHelper, const char *description, const double value) {
        debug(fprintf(stderr,"gobyAlignments_addStatisticDouble %s=%f\n", description, value));
        string descriptionStr(description);
        writerHelper->alignmentWriter->addStatistic(descriptionStr, value);
    }

    void gobyAlignments_setQueryIndexOccurrencesStoredInEntries(CAlignmentsWriterHelper *writerHelper, int value /* bool */) {
        writerHelper->alignmentWriter->setQueryLengthsStoredInEntries(value == 0 ? false : true);
    }

    /**
     * This should be called once for each target, gobyTargetIndex should start at 0
     * but the alignerTargetIndex can be any value. If you use this version to register
     * targets, when you call gobyAlEntry_setTargetIndex you can provide the
     * alignerTargetIndex and it will be converted to the gobyTargetIndex.
     */
    void gobyAlignments_addTargetWithTranslation(CAlignmentsWriterHelper *writerHelper, const unsigned int gobyTargetIndex, const unsigned int alignerTargetIndex, const char *targetName, const unsigned int targetLength) {
        debug(fprintf(stderr,"gobyAlignments_addTargetIdentifier '%s' = [gi=%u, ai=%u] length=%u\n", targetName, gobyTargetIndex, alignerTargetIndex, targetLength));
        string targetNameStr(targetName);
        writerHelper->alignmentWriter->addTargetIdentifier(targetNameStr, gobyTargetIndex);
        writerHelper->alignmentWriter->addTargetLength(targetLength);
        if (writerHelper->alignerToGobyTargetIndexMap == NULL) {
            writerHelper->alignerToGobyTargetIndexMap = new map<unsigned int, unsigned int>;
        }
        (*(writerHelper->alignerToGobyTargetIndexMap))[alignerTargetIndex] = gobyTargetIndex;
    }

    /**
     * This should be called once for each target, targetIndex should start at 0
     * on the first call and increment by one on each following call.
     */
    void gobyAlignments_addTarget(CAlignmentsWriterHelper *writerHelper, const unsigned int targetIndex, const char *targetName, const unsigned int targetLength) {
        debug(fprintf(stderr,"gobyAlignments_addTargetIdentifier '%s' = [ai=%u] length=%u\n", targetName, targetIndex, targetLength));
        string targetNameStr(targetName);
        writerHelper->alignmentWriter->addTargetIdentifier(targetNameStr, targetIndex);
        writerHelper->alignmentWriter->addTargetLength(targetLength);
    }


    /**
     * Checks if a target identifier has been registered.
     * @return returns 1 if targetName has been registered, otherwise 1.
     */
    int gobyAlignments_isTargetIdentifierRegistered(CAlignmentsWriterHelper *writerHelper, const char *targetName) {
        string targetNameStr(targetName);
        return writerHelper->alignmentWriter->isTargetIdentifierRegistered(targetNameStr);
    }

    /**
     * Look up the target index for the specified identifier. If a targetName is requested
     * that isn't in target_identifiers this will return 0, but but it isn't recommended
     * that this be called for targetName values that haven't been registered by
     * addTargetIdentifier() as 0 is also a valid target index.
     */
    unsigned int gobyAlignments_targetIndexForIdentifier(CAlignmentsWriterHelper *writerHelper, const char *targetName) {
        string targetNameStr(targetName);
        return writerHelper->alignmentWriter->targetIndexForIdentifier(targetNameStr);
    }

    /**
     * If you wish to provide the identifier and have the queryIndex generated
     * automatically. If the identifier has already been registered, you'll get
     * back the same queryIndex as before. This should be called if your INPUT
     * is Fasta or Fastq and your query identifier is a real string and you need
     * to gernate a queryIndex automatically.
     */
    unsigned gobyAlignments_addQueryIdentifier(CAlignmentsWriterHelper *writerHelper, const char *queryIdentifier) {
        string queryIdentifierStr(queryIdentifier);
        unsigned queryIndex = writerHelper->alignmentWriter->addQueryIdentifier(queryIdentifierStr);
        debug(fprintf(stderr,"gobyAlignments_addQueryIdentifier %s=%u\n", queryIdentifier, queryIndex));
        return queryIndex;
    }

    /**
     * If you wish to provide the identifier AND the queryIndex. If the 
     * identifier has already been registered, nothing new will happen.
     * This should be called if your input is compact-reads and you HAVE
     * query-identifier strings for your queries. If you don't have query
     * identifiers, don't bother to call this, you know the queryIndex and
     * that's good enough.
     */
    void gobyAlignments_addQueryIdentifierWithIndex(CAlignmentsWriterHelper *writerHelper, const char *queryIdentifier, unsigned int newQueryIndex) {
        debug(fprintf(stderr,"gobyAlignments_addQueryIdentifierWithIndex %s=%u\n", queryIdentifier, newQueryIndex));
        string queryIdentifierStr(queryIdentifier);
        writerHelper->alignmentWriter->addQueryIdentifierWithIndex(queryIdentifierStr, newQueryIndex);
    }
    
    // get an empty alignment entry to populate
    void gobyAlignments_appendEntry(CAlignmentsWriterHelper *writerHelper) {
        debug(fprintf(stderr,"---------------------------------------\n"));
        debug(fprintf(stderr,"gobyAlignments_appendEntry\n"));
        writerHelper->alignmentEntry = writerHelper->alignmentWriter->appendEntry();
        writerHelper->sequenceVariation = NULL;
	    writerHelper->lastSeqVarRefPosition = -1;
        writerHelper->lastSeqVarReadChar = '\0';
        writerHelper->lastSeqVarRefChar = '\0';
    }
    void gobyAlEntry_setMultiplicity(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setMultiplicity=%u\n", value));
        writerHelper->alignmentEntry->set_multiplicity(value);
    }
    void gobyAlignments_observeQueryIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlignments_observeQueryIndex=%u\n", value));
        if (writerHelper->smallestQueryIndex == -1) {
            writerHelper->smallestQueryIndex = value;
            writerHelper->largestQueryIndex = value;
        } else {
            writerHelper->smallestQueryIndex = min(value, writerHelper->smallestQueryIndex);
            writerHelper->largestQueryIndex = max(value, writerHelper->largestQueryIndex);
        }
    }
    void gobyAlEntry_setQueryIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setQueryIndex=%u\n", value));
        gobyAlignments_observeQueryIndex(writerHelper, value);
        writerHelper->alignmentEntry->set_query_index(value);
    }
    unsigned int gobyAlEntry_getQueryIndex(CAlignmentsWriterHelper *writerHelper) {
        return writerHelper->alignmentEntry->query_index();
    }

    /**
     * The target index. If you registered targets with gobyAlignments_addTargetWithTranslation()
     * value should in the range of the native aligner (can be 1-based or even any-based).
     * If you registered targets with gobyAlignments_addTarget() you should be providing
     * targets that are 0-based.
     */
    void gobyAlEntry_setTargetIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setTargetIndex=%u\n", value));
        unsigned int indexToUse;
        if (writerHelper->alignerToGobyTargetIndexMap == NULL) {
            indexToUse = value;
        } else {
            indexToUse = (*(writerHelper->alignerToGobyTargetIndexMap))[value];
        }
        writerHelper->alignmentEntry->set_target_index(indexToUse);
    }
    void gobyAlEntry_setPosition(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setPosition=%u\n", value));
        writerHelper->alignmentEntry->set_position(value);
    }
    void gobyAlEntry_setMatchingReverseStrand(CAlignmentsWriterHelper *writerHelper, int value /* bool */) {
        debug(fprintf(stderr,"gobyAlEntry_setMatchingReverseStrand=%d\n", value));
        writerHelper->alignmentEntry->set_matching_reverse_strand(value == 0 ? false : true);
    }
    void gobyAlEntry_setQueryPosition(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setQueryPosition=%u\n", value));
        writerHelper->alignmentEntry->set_query_position(value);
    }
    void gobyAlEntry_setScoreInt(CAlignmentsWriterHelper *writerHelper, int value) {
        float fValue = (float) value;
        debug(fprintf(stderr,"gobyAlEntry_setScore=%f\n", fValue));
        writerHelper->alignmentEntry->set_score(fValue);
    }
    void gobyAlEntry_setNumberOfMismatches(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setNumberOfMismatches=%u\n", value));
        writerHelper->alignmentEntry->set_number_of_mismatches(value);
    }
    void gobyAlEntry_setNumberOfIndels(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setNumberOfIndels=%u\n", value));
        writerHelper->alignmentEntry->set_number_of_indels(value);
    }
    void gobyAlEntry_setQueryAlignedLength(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setQueryAlignedLength=%u\n", value));
        writerHelper->alignmentEntry->set_query_aligned_length(value);
    }
    void gobyAlEntry_setTargetAlignedLength(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setTargetAlignedLength=%u\n", value));
        writerHelper->alignmentEntry->set_target_aligned_length(value);
    }
    void gobyAlEntry_setQueryLength(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setQueryLength=%u\n", value));
        writerHelper->alignmentEntry->set_query_length(value);
    }

    void gobyAlEntry_setMappingQuality(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setMappingQuality=%u\n", value));
        writerHelper->alignmentEntry->set_mapping_quality(value);
    }

    /**
     * Define left soft clipped details.
     * @param writerHelper writerHelper
     * @param start the position in query/quality where soft clips start (0-based)
     * @param size the size of the soft clip
     * @param query the query bases
     * @param quality the quality scores, will be adjusted by qualityAdjustment
     */
    void gobyAlEntry_setSoftClippedLeft(CAlignmentsWriterHelper *writerHelper,
            int start, int size, const char *query, const char *quality) {
        if (size > 0) {
            string *basesStr = writerHelper->alignmentEntry->mutable_softclippedbasesleft();
            string *qualityStr;
            if (quality != NULL) {
                qualityStr = writerHelper->alignmentEntry->mutable_softclippedqualityleft();
            }
            int i;
            for (i = start; i < start+size; i++) {
                (*basesStr) += query[i];
                if (quality != NULL) {
                    (*qualityStr) += (quality[i] + writerHelper->qualityAdjustment);
                }
            }
        }
    }

    /**
     * Define right soft clipped details.
     * @param writerHelper writerHelper
     * @param start the position in query/quality where soft clips start (0-based)
     * @param size the size of the soft clip
     * @param query the query bases
     * @param quality the quality scores, will be adjusted by qualityAdjustment
     */
    void gobyAlEntry_setSoftClippedRight(CAlignmentsWriterHelper *writerHelper,
            int start, int size, const char *query, const char *quality) {
        if (size > 0) {
            string *basesStr = writerHelper->alignmentEntry->mutable_softclippedbasesright();
            string *qualityStr;
            if (quality != NULL) {
                qualityStr = writerHelper->alignmentEntry->mutable_softclippedqualityright();
            }
            int i;
            for (i = start; i < start+size; i++) {
                (*basesStr) += query[i];
                if (quality != NULL) {
                    (*qualityStr) += (quality[i] + writerHelper->qualityAdjustment);
                }
            }
        }
    }

    /**
     * Output details about an alignment entries placed unmapped mate.
     * @param writerHelper writerHelper
     * @param length length of the query (and quality if it exists)
     * @param translateQuery if the query needs to be translated for BWA to bases (0-4 to ACGTN).
     * @param reverseStrand if unamped mate is thought to be reverse strand
     * @param query the query bases for the placed unmapped mate
     * @param quality the quality scores for the placed unmapped mate, will be adjusted by qualityAdjustment
     */
    void gobyAlEntry_setPlacedUnmapped(CAlignmentsWriterHelper *writerHelper,
            int length, int translateQuery /*bool*/, int reverseStrand /*bool*/,
            const char *query, const char *quality) {
        if (length > 0) {
            string *basesStr = writerHelper->alignmentEntry->mutable_placedunmappedsequence();
            string *qualityStr;
            if (quality) {
                qualityStr = writerHelper->alignmentEntry->mutable_placedunmappedquality();
            }
            int i;
            for (i = 0; i < length; i++) {
                if (reverseStrand) {
                    if (translateQuery) {
                        (*basesStr) += "TGCAN"[(int)(query[length - 1 - i])];
                    } else {
                        (*basesStr) += GOBY_COMPLEMENT_CODE[query[length - 1 - i]];
                    }
                    if (quality) {
                        (*qualityStr) += (quality[length - 1 - i] + writerHelper->qualityAdjustment);
                    }
                } else {
                    if (translateQuery) {
                        (*basesStr) += "ACGTN"[(int)(query[i])];
                    } else {
                        (*basesStr) += query[i];
                    }
                    if (quality) {
                        (*qualityStr) += (quality[i] + writerHelper->qualityAdjustment);
                    }
                }
            }
        }
    }

    void startNewSequenceVariation(CAlignmentsWriterHelper *writerHelper, unsigned int readIndex, unsigned int refPosition) {
        debug(fprintf(stderr,"... startNewSequenceVariation\n"));
        writerHelper->sequenceVariation = writerHelper->alignmentEntry->add_sequence_variations();
        writerHelper->sequenceVariation->set_read_index(readIndex);
        writerHelper->sequenceVariation->set_position(refPosition);
    }

    /**
     * This method will accumulate sequence variations. If > 1 are registered at readIndex positions next to
     * each other, they will be put into the SAME sequence variation. If a gap of great than one in readIndex
     * is found an additional SeqVar will be created.
     * This will also take "reverse" and "position" into account, which should already have been defined for
     * writerHelper->alignmentEntry.
     * @param readIndex coming in, readIndex is >0< based. When string position and readIndex into the SeqVar
     *        position and readIndex are >1< based.
     * For contiguous sequence variations, this assumes readIndex will increment by one each time.
     */
    void gobyAlEntry_addSequenceVariation(CAlignmentsWriterHelper *writerHelper,
            unsigned int readIndex, unsigned int refPosition, char refChar, char readChar,
            int hasQualCharInt /* bool */, char readQualChar) {
        bool hasQualChar = hasQualCharInt == 0 ? false : true;
        debug(fprintf(stderr,"gobyAlEntry_addSequenceVariation readIndex=%u refPosition=%u ref=%c read=%c hasQualChar=%s readQualChar=%d, adjustedreadQualChar=%d\n",
                readIndex, refPosition, refChar, readChar, hasQualChar ? "true" : "false",
                readQualChar, readQualChar + writerHelper->qualityAdjustment));
        if (writerHelper->sequenceVariation == NULL || writerHelper->lastSeqVarRefPosition == -1) {
            // New sequence variation
            startNewSequenceVariation(writerHelper, readIndex, refPosition);
        } else if (writerHelper->lastSeqVarRefPosition + 1 == refPosition || writerHelper->lastSeqVarRefPosition == refPosition) {
            // Append to prev SequenceVar entry
            debug(fprintf(stderr,"... appending to previous seqVar\n"));
        } else if  ((readChar == '-' && writerHelper->lastSeqVarReadChar != '-') ||
                    (refChar == '-' && writerHelper->lastSeqVarRefChar != '-') ||
                    (writerHelper->lastSeqVarReadChar == '-' && readChar != '-')  ||
                    (writerHelper->lastSeqVarRefChar == '-' && refChar != '-')) {
            // transitioning into or out of an insert or delete, make a new SeqVar
            startNewSequenceVariation(writerHelper, readIndex, refPosition);
        } else {
            // Not contiguous to previous SeqVar
            startNewSequenceVariation(writerHelper, readIndex, refPosition);
        }
        writerHelper->lastSeqVarRefPosition = refPosition;
        writerHelper->lastSeqVarReadChar = readChar;
        writerHelper->lastSeqVarRefChar = refChar;
        string *from = writerHelper->sequenceVariation->mutable_from();
        string *to = writerHelper->sequenceVariation->mutable_to();
        (*from) += refChar;
        (*to) += readChar;
        debug(fprintf(stderr,"... sv->read_index=%u sv->position=%u from=%s to=%s\n",
            writerHelper->sequenceVariation->read_index(), writerHelper->sequenceVariation->position(),
            from->c_str(), to->c_str()));
        if (hasQualChar) {
            string *toQuality = writerHelper->sequenceVariation->mutable_to_quality();
            (*toQuality) += (readQualChar + writerHelper->qualityAdjustment);
        }
    }

    /**
     * !! Before calling this alignments_appendEntry(writerHelper) should
     * !! have already been called to start a new alignment entry.
     * This will output the sequence variations given reference, query,
     * and optional quality string. Upon completion, the sequence variations
     * will be outputted to goby alignment AND outMatches, outSubs, outInserts,
     * and outDeletes will be populated by the number of the same found.
     * Only positions between queryStart and queryEnd (0-based) will be
     * oberseved which leads to a slight deviation in the number of matches
     * found with respect to the number of matches reported by Gsnap,
     * in the Gsnap considers matches even outside of the queryStart/queryEnd
     * boundries (left and right clipping) to still be counted as matches,
     * even though it says to clip them in the alignment.
     * !! IMPORTANT !! The mergedRef, mergedQuery, and mergedQual should all
     * be in the DIRECTION OF THE REFERENCE, NOT the direction of the read.
     * By default, Gsnap outputs data in the direction of the READ, so if
     * the data comes from Gsnap, the reference and query would need be
     * reverse complemented and the quality would need to be reversed before
     * calling this method.
     * @param writerHelper the Goby writer helper
     * @param mergedRef the reference, note that if there are insertions, there
     * should be "-" characters in the appropriate places in mergedRef
     * @param mergedQuery the query/read, note that if there are deletions,
     * there should be "-" characters in the appropriate places in mergedQuery
     * @param mergedQual the merged quality scores or NULL, if there are no
     * quality scores. If there are  deletions, there should be "\0" characters
     * in the appropriate places in mergedQual
     * @param queryStart the inclusive, 0-based start point to start observing
     * for sequence variations. This is also known as the left clip value.
     * @param queryEnd the inclusive, 0-based start point to stop observing
     * for sequence variations. This is also known as the right clip value.
     * @param reverseStrand set to 1 of this query was matched on the
     * reverse strand
     * @param outMatches output value, the nubmer of matches found
     * @param outSubs output value, the number of substitutions found
     * @param outInserts output value, the number of insert bases found
     * @param outDeletes output value, the number of delete bases found
     */
    void gobyAlignments_outputSequenceVariations(
            CAlignmentsWriterHelper *writerHelper,
            const char *reference, const char *query, const char *quality,
            int queryStart, int queryEnd,
            int reverseStrand,
            int *outMatches, int *outSubs, int *outInserts, int *outDeletes) {

        int refPosition = 0;
        int readIndex = 0;
        int hasQualChar;

        int mergedSize = strlen(reference);
        unsigned int refPositions[mergedSize];
        unsigned int readIndexes[mergedSize];
        for (int i = 0; i < mergedSize; i++) {
            char refChar = reference[i];
            if (refChar != '-') {
                refPosition++;
            }
            refPositions[i] = refPosition - queryStart;
            int readIndexI = reverseStrand ? mergedSize - i - 1 : i;
            char queryChar = query[readIndexI];
            if (queryChar != '-') {
                readIndex++;
            }
            readIndexes[readIndexI] = readIndex;
        }

        debug(
            cout << "r=";
            for (int i = 0; i < mergedSize; i++) {
                cout << (refPositions[i] % 10);
            }
            cout << endl;
            cout << "r=" << reference << endl;
            cout << "q=" << query << endl;
            cout << "q=";
            for (int i = 0; i < mergedSize; i++) {
                cout << (readIndexes[i] % 10);
            }
            cout << endl;
        )

        int matches = 0;
        int subs = 0;
        int inserts = 0;
        int deletes = 0;
//        for (int i = merged->queryStart + merged->startClip; 
//                 i <= merged->queryEnd - merged->endClip; i++) {
        for (int i = queryStart;  i <= queryEnd; i++) {
            char refChar = toupper(reference[i]);
            char queryChar = toupper(query[i]);
            char qualChar;
            if (quality != NULL) {
                qualChar = quality[i];
                if (qualChar > GOBY_NO_QUAL) {
                    qualChar += writerHelper->qualityAdjustment;
                }
            } else {
                qualChar = GOBY_NO_QUAL;
            }
            
            if (refChar == queryChar) {
                matches++;
            } else {
                if (refChar == '-') {
                    inserts++;
                } else if (queryChar == '-') {
                    deletes++;
                } else if (queryChar == 'N') {
                    // Done penalize for query=N
                    matches++;
                } else {
                    subs++;
                }
                debug(
                    cout << "Seqvar refPosition=" << refPositions[i] << " " <<
                            "readIndex=" << readIndexes[i] << " " << 
                            "from=" << refChar << " " <<
                            "to=" << queryChar << endl;
                )
                hasQualChar = (qualChar == GOBY_NO_QUAL ? false : true);
                gobyAlEntry_addSequenceVariation(writerHelper,
                    readIndexes[i], refPositions[i], 
                    refChar, queryChar, hasQualChar, qualChar);
            }
        }

        *outMatches = matches;
        *outSubs = subs;
        *outInserts = inserts;
        *outDeletes = deletes;
    }

    /**
     * This method should be called once AFTER you call
     *    gobyAlEntry_setQueryIndex(...)
     *    gobyAlEntry_setQueryAlignedLength(...)
     */
    void gobyAlEntry_appendTooManyHits(CAlignmentsWriterHelper *writerHelper, unsigned int queryIndex, unsigned int alignedLength, int numberOfHits) {
        debug(fprintf(stderr,"gobyAlEntry_appendTooManyHits queryIndex=%u alignedLength=%u numberOfHits=%u\n", queryIndex, alignedLength, numberOfHits));
        writerHelper->tmhWriter->append(queryIndex, numberOfHits, alignedLength);
    }

    void gobyAlEntry_setFragmentIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setFragmentIndex=%u\n", value));
        writerHelper->alignmentEntry->set_fragment_index(value);
    }

    void gobyAlEntry_setInsertSize(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setInsertSize=%u\n", value));
        writerHelper->alignmentEntry->set_insert_size(value);
    }

    void gobyAlEntry_setAmbiguity(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setAmbiguity=%u\n", value));
        writerHelper->alignmentEntry->set_ambiguity(value);
    }
    
    void gobyAlEntry_setQueryIndexOccurrences(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setQueryIndexOccurrences=%u\n", value));
        writerHelper->alignmentEntry->set_query_index_occurrences(value);
    }

    // These are only used when dealing with a Query Pair
    void gobyAlEntry_setPairFlags(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setPairFlags=%u\n", value));
        writerHelper->alignmentEntry->set_pair_flags(value);
    }
    void gobyAlEntry_setPairTargetIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setPairTargetIndex=%u\n", value));
        writerHelper->alignmentEntry->mutable_pair_alignment_link()->set_target_index(value);
    }
    void gobyAlEntry_setPairFragmentIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setPairFragmentIndex=%u\n", value));
        writerHelper->alignmentEntry->mutable_pair_alignment_link()->set_fragment_index(value);
    }
    void gobyAlEntry_setPairPosition(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setPairPosition=%u\n", value));
        writerHelper->alignmentEntry->mutable_pair_alignment_link()->set_position(value);
    }

    // These are only used when dealing with a Splice
    void gobyAlEntry_setSplicedFlags(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedFlags=%u\n", value));
        writerHelper->alignmentEntry->set_spliced_flags(value);
    }
    void gobyAlEntry_setSplicedForwardTargetIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedForwardTargetIndex=%u\n", value));
        writerHelper->alignmentEntry->mutable_spliced_forward_alignment_link()->set_target_index(value);
    }
    void gobyAlEntry_setSplicedForwardFragmentIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedForwardFragmentIndex=%u\n", value));
        writerHelper->alignmentEntry->mutable_spliced_forward_alignment_link()->set_fragment_index(value);
    }
    void gobyAlEntry_setSplicedForwardPosition(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedForwardPosition=%u\n", value));
        writerHelper->alignmentEntry->mutable_spliced_forward_alignment_link()->set_position(value);
    }

    void gobyAlEntry_setSplicedBackwardTargetIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedBackwardTargetIndex=%u\n", value));
        writerHelper->alignmentEntry->mutable_spliced_backward_alignment_link()->set_target_index(value);
    }
    void gobyAlEntry_setSplicedBackwardFragmentIndex(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedBackwardFragmentIndex=%u\n", value));
        writerHelper->alignmentEntry->mutable_spliced_backward_alignment_link()->set_fragment_index(value);
    }
    void gobyAlEntry_setSplicedBackwardPosition(CAlignmentsWriterHelper *writerHelper, unsigned int value) {
        debug(fprintf(stderr,"gobyAlEntry_setSplicedBackwardPosition=%u\n", value));
        writerHelper->alignmentEntry->mutable_spliced_backward_alignment_link()->set_position(value);
    }

    /**
     * Split a string into a vector.
     */
    std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems) {
        std::stringstream ss(s);
        std::string item;
        while(std::getline(ss, item, delim)) {
            elems.push_back(item);
        }
        return elems;
    }

    CSamHelper *samHelper_getResetSamHelper(CAlignmentsWriterHelper *writerHelper) {
        CSamHelper *samHelper = writerHelper->samHelper;
        if (samHelper == NULL) {
            writerHelper->samHelper = new CSamHelper;
            samHelper = writerHelper->samHelper;           
            samHelper->cpp_cigar = new string();
            samHelper->cpp_md = new string();
            samHelper->cpp_sourceQuery = new string();
            samHelper->cpp_sourceQual = new string();
            samHelper->cpp_query = new string();
            samHelper->cpp_qual = new string();
            samHelper->cpp_ref = new string();
            samHelper->minQualValue = writerHelper->qualityAdjustment;
        } else {
            samHelper->cpp_cigar->clear();
            samHelper->cpp_md->clear();
            samHelper->cpp_sourceQuery->clear();
            samHelper->cpp_sourceQual->clear();
            samHelper->cpp_query->clear();
            samHelper->cpp_qual->clear();
            samHelper->cpp_ref->clear();
        }
        samHelper->alignedLength = 0;
        samHelper->numInsertions = 0;
        samHelper->numDeletions = 0;
        samHelper->numMisMatches = 0;
        samHelper->score = 0;
        samHelper->numLeftClipped = 0;
        samHelper->numRightClipped = 0;
        return samHelper;
    }

    void samHelper_setCigar(CSamHelper *samHelper, const char *cigar) {
        if (cigar) {
            (*samHelper->cpp_cigar) += cigar;
        }
    }

    void samHelper_addCigarItem(CSamHelper *samHelper, int length, char op) {
        std:stringstream lengthStream;
        lengthStream << length;
        (*samHelper->cpp_cigar) += lengthStream.str();
        (*samHelper->cpp_cigar) += op;
    }

    void samHelper_setMd(CSamHelper *samHelper, const char *md) {
        if (md) {
            (*samHelper->cpp_md) += md;
        }
    }

    void samHelper_setQueryTranslate(CSamHelper *samHelper, unsigned char *reads, unsigned char *qual, 
            unsigned int length, unsigned int reverseStrand) {
        if (!reverseStrand && qual) {
            (*samHelper->cpp_sourceQual) += (char *) qual;
        }
        int i;
        for (i = 0; i < length; i++) {
            if (reverseStrand) {
                (*samHelper->cpp_sourceQuery) += "TGCAN"[(int)(reads[length - 1 - i])];
                if (qual) {
                    (*samHelper->cpp_sourceQual) += (qual[length - 1 - i]);
                }
            } else {
                (*samHelper->cpp_sourceQuery) += "ACGTN"[(int)(reads[i])];
            }
        }
    }

    void samHelper_setQuery(CSamHelper *samHelper, const char *reads, const char *qual,
            unsigned int length, unsigned int reverseStrand) {
        // TODO: Do we need to reverse this if reverse query?
        if (reads && strlen(reads) > 0) {
            (*samHelper->cpp_sourceQuery) += reads;
        }
        if (qual && strlen(qual) > 0) {
            (*samHelper->cpp_sourceQual) += qual;
        }
    }

    const char *samHelper_getCigarStr(CSamHelper *samHelper) {
        return samHelper->cpp_cigar->c_str();
    }

    void applyCigar(CSamHelper *samHelper) {
        //Setting up pcre2.
        pcre2_code *re;
        string sPattern = string("([0-9]+)([SMID])");
        PCRE2_SPTR pattern = (PCRE2_SPTR)(sPattern.c_str());
        PCRE2_SPTR subject = (PCRE2_SPTR)(samHelper->cpp_cigar->c_str());
        int errornumber;
        int rc;
        PCRE2_SIZE erroroffset;
        PCRE2_SIZE *ovector;
        PCRE2_SIZE subject_length = (PCRE2_SIZE)(sPattern.length());
        pcre2_match_data *match_data;
        re = pcre2_compile(pattern, PCRE2_ZERO_TERMINATED, 0, &errornumber, &erroroffset, NULL);
        match_data = pcre2_match_data_create_from_pattern(re, NULL);
        PCRE2_SIZE start_offset = 0;

        debug (fprintf(stderr, ":: Applying cigar=%s\n", samHelper->cpp_cigar->c_str());)
        int length;
        char op;
        int posInReads = 0;
        int i;
        bool startOfCigar = true;
        samHelper->numInsertions = 0;
        samHelper->numDeletions = 0;
        samHelper->numMisMatches = 0;

        do {
            rc = pcre2_match(re, subject, subject_length, start_offset, 0, match_data, NULL);

            if (rc != PCRE2_ERROR_NOMATCH) {
                //Match succeeded. Converting strings to the needed types.
                ovector = pcre2_get_ovector_pointer(match_data);
                char* sLength = new char[ovector[3] - ovector[2] + 1];
                char* sOp = new char[ovector[5] - ovector[4] + 1];
                sprintf(sLength, "%.*s", (int)(ovector[3] - ovector[2]), (char*)(subject + ovector[2]));
                sprintf(sOp, "%.*s", (int)(ovector[5] - ovector[4]), (char*)(subject + ovector[4]));
                length = stoi(string(sLength));
                op = string(sOp)[0];
                delete[] sLength;
                delete[] sOp;
                //Preparing next match.
                start_offset = ovector[1];
            }

            switch(op) {
                case 'S':
                    // Soft clipping
                    for (i = 0; i < length; i++) {
                        (*samHelper->cpp_ref) += '-';
                        (*samHelper->cpp_query) += '-';
                        (*samHelper->cpp_qual) += ((char) samHelper->minQualValue); // min quality
                    }
                    if (startOfCigar) {
                        samHelper->numLeftClipped = length;
                    } else {
                        samHelper->numRightClipped = length;
                    }
                    posInReads += length;
                    break;
                case 'M':
                    // Account for matches AND mismatches. Any mis-matches will be fixed in applyMd()
                    (*samHelper->cpp_query) += samHelper->cpp_sourceQuery->substr(posInReads, length);
                    if (samHelper->cpp_sourceQual->size() != 0) {
                        (*samHelper->cpp_qual) += samHelper->cpp_sourceQual->substr(posInReads, length);
                    }
                    (*samHelper->cpp_ref) += samHelper->cpp_sourceQuery->substr(posInReads, length);
                    posInReads += length;
                    break;
                case 'I':
                    (*samHelper->cpp_query) += samHelper->cpp_sourceQuery->substr(posInReads, length);
                    if (samHelper->cpp_sourceQual->size() != 0) {
                        (*samHelper->cpp_qual) += samHelper->cpp_sourceQual->substr(posInReads, length);
                    }
                    for (i = 0; i < length; i++) {
                        (*samHelper->cpp_ref) += '-';
                    }
                    samHelper->numInsertions += length;
                    posInReads += length;
                    break;
                case 'D':
                    for (i = 0; i < length; i++) {
                        (*samHelper->cpp_query) += '-';
                        if (samHelper->cpp_sourceQual->size() != 0) {
                            // Minimum qual, placing min value here but it shouldn't get written to
                            // sequence variations
                            (*samHelper->cpp_qual) += ((char) samHelper->minQualValue); // min quality
                        }
                        (*samHelper->cpp_ref) += '?';  // provided in applyMd()
                    }
                    samHelper->numDeletions += length;
                    break;
            }
            startOfCigar = false;
        } while (rc != PCRE2_ERROR_NOMATCH);
    }

    void applyMd(CSamHelper *samHelper) {
        // My RE is simplified from the SAM spec but performs the same task
        // the major difference being mine would allow 5ACG where the current
        // spec would require 5A0C0G0 (which mine will still work with just fine).

        //Setting up pcre2.
        pcre2_code *re;
        string sPattern = string("([0-9]+|[ACGTN]|\\^[ACGTN]+)");
        PCRE2_SPTR pattern = (PCRE2_SPTR)(sPattern.c_str());
        PCRE2_SPTR subject = (PCRE2_SPTR)(samHelper->cpp_md->c_str());
        int errornumber;
        int rc;
        PCRE2_SIZE erroroffset;
        PCRE2_SIZE *ovector;
        PCRE2_SIZE subject_length = (PCRE2_SIZE)(sPattern.length());
        pcre2_match_data *match_data;
        re = pcre2_compile(pattern, PCRE2_ZERO_TERMINATED, 0, &errornumber, &erroroffset, NULL);
        match_data = pcre2_match_data_create_from_pattern(re, NULL);
        PCRE2_SIZE start_offset = 0;

        debug (fprintf(stderr, ":: Applying md=%s\n", samHelper->cpp_md->c_str());)
        string mdPart;
        int position = samHelper->numLeftClipped;
        int i;

        do {
            rc = pcre2_match(re, subject, subject_length, start_offset, 0, match_data, NULL);

            if (rc != PCRE2_ERROR_NOMATCH) {
                //Match succeeded.
                ovector = pcre2_get_ovector_pointer(match_data);
                char* sMdPart = new char[ovector[3] - ovector[2] + 1];
                sprintf(sMdPart, "%.*s", (int)(ovector[3] - ovector[2]), (char*)(subject + ovector[2]));
                mdPart = string(sMdPart);
                delete[] sMdPart;
                //Preparing next match.
                start_offset = ovector[1];
            }

            if (isdigit(mdPart[0])) {
                int length = atoi(mdPart.c_str());
                position += length;
            } else if (mdPart[0] == '^') {
                // Adjust the ref with these characters, ignoring the ^ character so start at 1
                for(i = 1; i < mdPart.size(); i++) {
                    (*samHelper->cpp_ref)[position++] = mdPart[i];
                }
            } else {
                // The regex should only allow a single character here, but we'll accept multiple
                for(i = 0; i < mdPart.size(); i++) {
                    (*samHelper->cpp_ref)[position++] = tolower(mdPart[i]);
                    samHelper->numMisMatches++;
                }
            }
        } while (rc != PCRE2_ERROR_NOMATCH);
    }

    void samHelper_constructRefAndQuery(CSamHelper *samHelper) {
        samHelper->cpp_query->clear();
        samHelper->cpp_qual->clear();
        samHelper->cpp_ref->clear();
        samHelper->alignedLength = 0;
        samHelper->numInsertions = 0;
        samHelper->numDeletions = 0;
        samHelper->numMisMatches = 0;
        samHelper->score = 0;

        debug (
            fprintf(stderr, ":: Reference and query before construction\n");
            fprintf(stderr, ":: read=%s\n", samHelper->cpp_sourceQuery->c_str());
        )

        applyCigar(samHelper);
        applyMd(samHelper);
        samHelper->alignedLength = samHelper->cpp_query->size();
        samHelper->score = samHelper->alignedLength - samHelper->numDeletions - samHelper->numInsertions - samHelper->numMisMatches;
        
        debug (
            fprintf(stderr, ":: Reference and query constructed via SAM\n");
            fprintf(stderr, ":: ref =%s\n", samHelper->cpp_ref->c_str());
            fprintf(stderr, ":: read=%s\n", samHelper->cpp_query->c_str());
        )
        // Figure out start of alignment and alignment length, by observing mismatches at head and tail
        if (samHelper->cpp_query->size() != samHelper->cpp_ref->size()) {
            printf("ERROR! reconstructed reads and refs don't match in size!!\n");
        }
    }

    /**
     * Depends on samHelper_setQueryTranslate() being called.
     */
    const char *samHelper_sourceQuery(CSamHelper *samHelper) {
        return samHelper->cpp_sourceQuery->c_str();
    }
    /**
     * Depends on samHelper_setQueryTranslate() being called.
     */
    const char *samHelper_sourceQual(CSamHelper *samHelper) {
        if (samHelper->cpp_sourceQual->size() == 0) {
            return NULL;
        } else {
            return samHelper->cpp_sourceQual->c_str();
        }
    }
    /**
     * Depends on samHelper_constructRefAndQuery() being called.
     */
    const char *samHelper_constructedQuery(CSamHelper *samHelper) {
        return samHelper->cpp_query->c_str();
    }
    /**
     * Depends on samHelper_constructRefAndQuery() being called.
     */
    const char *samHelper_constructedQual(CSamHelper *samHelper) {
        if (samHelper->cpp_qual->size() == 0) {
            return NULL;
        } else {
            return samHelper->cpp_qual->c_str();
        }
    }

    /**
     * Depends on samHelper_constructRefAndQuery() being called.
     */
    const char *samHelper_constructedRef(CSamHelper *samHelper) {
        return samHelper->cpp_ref->c_str();
    }
    
	void gobyAlignments_finished(CAlignmentsWriterHelper *writerHelper, unsigned int numberOfReads) {
        if (writerHelper->numberOfMisParsedReads > 0) {
            fprintf(stderr, "Number of mis-parsed reads: %d\n",
                    writerHelper->numberOfMisParsedReads);
        }
        debug(fprintf(stderr,"gobyAlignments_finished\n"));
        if (writerHelper != NULL) {
            // Stats for the alignment
            writerHelper->alignmentWriter->setNumberOfQueries(numberOfReads);
            writerHelper->alignmentWriter->setSmallestSplitQueryIndex(writerHelper->smallestQueryIndex);
            writerHelper->alignmentWriter->setLargestSplitQueryIndex(writerHelper->largestQueryIndex);
            writerHelper->alignmentWriter->setNumberOfAlignedReads(writerHelper->numberOfAlignedReads);

            // Close and delete
            writerHelper->alignmentWriter->close();
            delete writerHelper->alignmentWriter;
            writerHelper->tmhWriter->write();
            delete writerHelper->tmhWriter;
            if (writerHelper->samHelper != NULL) {
                delete writerHelper->samHelper->cpp_cigar;
                delete writerHelper->samHelper->cpp_md;
                delete writerHelper->samHelper->cpp_sourceQuery;
                delete writerHelper->samHelper->cpp_sourceQual;
                delete writerHelper->samHelper->cpp_query;
                delete writerHelper->samHelper->cpp_qual;
                delete writerHelper->samHelper->cpp_ref;
                delete writerHelper->samHelper;
            }
            if (writerHelper->alignerToGobyTargetIndexMap != NULL) {
                writerHelper->alignerToGobyTargetIndexMap->clear();
                delete writerHelper->alignerToGobyTargetIndexMap;
            }
            gobyGsnap_destoryAlignment(writerHelper);
            delete writerHelper;

        }
	}

    /** These come from GSNAP. */
    const char* hitTypes[] = {
              "EXACT", "SUB", "INSERTION", "DELETION", "HALFSPLICE_DONOR",
              "HALFSPLICE_ACCEPTOR", "SPLICE", "ONE_THIRD_SHORTEXON",
              "TWO_THIRDS_SHORTEXON", "SHORTEXON", "TERMINAL" };

    /**
     * Debugging to print read and ref sequences including padding.
     * This ALWAYS prints if it is called. Don't call it if you don't want it to print.
     */
    void gobyAlignments_debugSequences(CAlignmentsWriterHelper *writerHelper, int hitType, char *refSequence, char *readSequence, int paddingLeft, int paddingRight) {
        fprintf(stderr,":: type=%s, paddingLeft=%d, paddingRight=%d\n", hitTypes[hitType], paddingLeft, paddingRight);
        fprintf(stderr,":: ref =");
        for (int i = 0; i < paddingLeft; i++) {
            fprintf(stderr, "_");
        }
        fprintf(stderr,"%s", refSequence);
        for (int i = 0; i < paddingRight; i++) {
            fprintf(stderr, "_");
        }
        fprintf(stderr,"\n");

        fprintf(stderr,":: read=");
        for (int i = 0; i < paddingLeft; i++) {
            fprintf(stderr, "_");
        }
        fprintf(stderr,"%s", readSequence);
        for (int i = 0; i < paddingRight; i++) {
            fprintf(stderr, "_");
        }
        fprintf(stderr,"\n");
    }

	/**
	 * Duplicate a string with an optional length to copy. If the incoming
	 * str is NULL, this will return null. If length is -1, the length
	 * that is copied is strlen(str). The size of the returned buffer is
	 * always length + 1 (to include the trailing '\0'). The caller is required
	 * to FREE the string.
	 * @param str the string to copy
	 * @param the maximum length to copy or -1 for the whole string
	 * @return the duplicate string.
	 */

    char *goby_copy_string(char *str, int length) {
        int copy_length = length;
        char *new_str = (char *) NULL;

        if (str != NULL) {
                if (copy_length == -1) {
                        copy_length = strlen(str);
                }
                new_str = (char *) calloc(copy_length + 1, sizeof(char));
                strncpy(new_str, str, copy_length);
                new_str[copy_length] = '\0';
        }

        return new_str;
    }

}