Reading QUASAR scoring matrices in C++

/**
 * QUASARReader.cpp
 * An efficient reader for the QUASAR substitution matrix format.
 * See http://www.bio.ifi.lmu.de/QUASAR/
 * 
 * Requires Boost String algorithms (header-only)
 * 
 * Compile like this: clang++ -o quasar quasar.cpp
 * 
 * Version 1.0.3
 * 
 * Copyright (c) 2013 Uli Koehler
 * This file is released as public domain.
 * Please mention the original author if used in your software.
 */
#include <fstream>
#include <iostream>
#include <cassert>
#include <string>
#include <cstring>
#include <map>
#include <sstream>
#include <algorithm>
#include <vector>
#include <boost/algorithm/string.hpp>

/**
 * A quadratic substitution matrix.
 * This class only supports symmetric substitution
 * matrices that map equal character sets.
 */
class SubstitutionMatrix {
public:
    /**
     * Construct a new substitution matrix by character index and default score.
     */
    SubstitutionMatrix(const std::string& index, double defaultScore = 0) : index(index) {
        size = index.size();
        matrix = new double[size * size];
        for(size_t i = 0; i < size * size; i++) {
            matrix[i] = defaultScore;
        }
    }
    ~SubstitutionMatrix() {
        if(matrix != nullptr) {
            delete[] matrix;
        }
    }
    /**
     * For a given character, determine the index
     * in the matrix.
     * @return The index, or string::npos if not found
     */
    int getIndexForChar(char c) {
        return index.find(c);
    }
    /**
     * Set the score of replacing character a with character b.
     */
    void setScore(char a, char b, double score) {
        setScore(getIndexForChar(a), getIndexForChar(b), score);
    }
    /**
     * Set a score by index.
     * @param a The index of character a
     * @param b The index of character b
     * @param score The score of replacing a by b (or b by a)
     */
    void setScore(int a, int b, double score) {
        assert(a >= 0);
        assert(b >= 0);
        assert(a < size);
        assert(b < size);
        //Symmetric matrix - a->b == b-> a
        matrix[a + size*b] = score;
        matrix[b + size*a] = score;
    }
    /**
     * Get the score for replacing
     * the character a with the character b.
     * @return The score
     */
    double getScore(char a, char b) {
        //Replacing a and b does not matter (symmetric matrix)
        return getScore(getIndexForChar(a), getIndexForChar(b));
    }
    /**
     * Get the score for replacing
     * the character at position a with
     * the character at position b
     * @return The score
     */
    double getScore(int a, int b) {
        assert(a >= 0);
        assert(b >= 0);
        assert(a < size);
        assert(b < size);
        //Replacing a and b does not matter (symmetric matrix)
        return matrix[b + size * a];
    }
    /**
     * Get the matrix size.
     * The matrix has getMatrixSize()**2 entries.
     */
    size_t getMatrixSize() {
        return size * size;
    }
    /**
     * Get a string of all characters in the matrix alphabet
     */
    std::string getAlphabet() {
        return index;
    }
    /**
     * Get the raw matrix buffer.
     * The valid length is sizeof(double) * size * size
     */
    double* getMatrix() {
        return matrix;
    }
private:
    /**
     * The two dimensional matrix that contains the scores.
     * Row-first or Column-first indexing does not matter here
     * because the matrix is symmetric.
     */
    double* matrix;
    size_t size;
    /**
     * The list of alphabet characters,
     * assigned to the array indices of the matrix.
     * The first character is assigned to index 0
     */
    std::string index;
};

SubstitutionMatrix* readSubstitutionMatrix(const char* filename) {
    static const char rowIndexPrefix[] = "ROWINDEX";
    static const char colIndexPrefix[] = "COLINDEX";
    static const char matrixPrefix[] = "MATRIX";
    std::ifstream infile(filename);
    std::string line = "", rowIndex = "", colIndex = "";
    int matrixRowIndex = 0; //The number of MATRIX rows already processed
    SubstitutionMatrix* matrix = nullptr;
    while(std::getline(infile, line)) {
        if(line.compare(0, sizeof(rowIndexPrefix) - 1, rowIndexPrefix) == 0) {
            rowIndex = line.substr(sizeof(rowIndexPrefix));
            boost::algorithm::trim(rowIndex);
        } else if(line.compare(0, sizeof(colIndexPrefix) - 1, colIndexPrefix) == 0) {
            colIndex = line.substr(sizeof(colIndexPrefix));
            boost::algorithm::trim(colIndex);
        } else if(line.compare(0, sizeof(matrixPrefix) - 1, matrixPrefix) == 0) {
            //Initialize matrix if not done yet
            if(matrix == nullptr) {
                //Matrices with different row and column indices are currently unsupported
                assert(!rowIndex.empty());
                assert(!colIndex.empty());
                assert(rowIndex == colIndex);
                matrix = new SubstitutionMatrix(rowIndex);
            }
            //Split elements
            int matrixColIndex = 0; //The number of columns in the current row already processed
            std::vector<std::string> elements;
            std::string lineContent = boost::trim_copy(line.substr(sizeof(matrixPrefix)));
            boost::split(elements, lineContent, boost::is_any_of("\t "), boost::token_compress_on);
            for(std::string elem : elements) {
                boost::trim(elem);
                double score = atof(elem.c_str());
                //Row and column index are interchangable
                matrix->setScore(matrixRowIndex, matrixColIndex, score);
                matrixColIndex++;
            }
            matrixRowIndex++;
        }
    }
    infile.close();
    return matrix;
}

using std::cout;
using std::cerr;
using std::endl;

/**
 * The main function just tests if the matrix reader does not fail for a given matrix.
 */
int main(int argc, char** argv) {
    if(argc < 2) {
        cerr << "Usage: " << argv[0] << " <QUASAR input file>" << endl;
        return 1;
    }
    SubstitutionMatrix* matrix = readSubstitutionMatrix(argv[1]);
    delete matrix;
    return 0;
}