haplo_predict.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <config.h>

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <inttypes.h>

#include <libxml/tree.h>

#ifdef HAPLO_HAVE_PTHREAD
#include <pthread.h>
#endif

#ifdef HAPLO_HAVE_DMALLOC
#include <dmalloc.h>
#endif

#include <jwsc/base/error.h>
#include <jwsc/base/option.h>
#include <jwsc/base/file_io.h>
#include <jwsc/vector/vector.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matblock/matblock.h>

#include "haplo_groups.h"
#include "options.h"
#include "output.h"
#include "input.h"
#include "xml.h"
#include "nb_freq.h"
#include "nb_gauss.h"
#include "nb_gmm.h"
#include "mv_gmm.h"
#ifdef HAPLO_ENABLE_SVM
#include "svm_tree.h"
#endif
#ifdef HAPLO_ENABLE_WEKA
#include "weka.h"
#endif
#include "nearest.h"


#ifdef HAPLO_ENABLE_SVM
#define NUM_SVM_OPTS 2
#else
#define NUM_SVM_OPTS 0
#endif

#ifdef HAPLO_ENABLE_WEKA
#define NUM_WEKA_OPTS 4
#else
#define NUM_WEKA_OPTS 0
#endif

#define  NUM_OPTS_NO_ARG    1 + NUM_SHARED_OPTS_NO_ARG
#define  NUM_OPTS_WITH_ARG  13 + NUM_SVM_OPTS + NUM_WEKA_OPTS + NUM_SHARED_OPTS_WITH_ARG


#define  NUM_ALGOS     8
#define  LABEL_COL     0
#define  OUTPUT_FNAME  "/dev/stdout"


typedef struct 
{ 
    Vector_u32** labels;
    Vector_d**   confs;
    Matrix_i32*  markers; 
    uint32_t     n;
#ifdef HAPLO_HAVE_PTHREAD
    pthread_mutex_t mutex;
#endif
}
Predict_params;


Option_no_arg opts_no_arg[NUM_OPTS_NO_ARG];

Option_with_arg opts_with_arg[NUM_OPTS_WITH_ARG];

static const char* output_fname = OUTPUT_FNAME;


uint32_t get_num_opts_no_arg()
{
    return NUM_OPTS_NO_ARG;
}

uint32_t get_num_opts_with_arg()
{
    return NUM_OPTS_WITH_ARG;
}

void print_usage()
{
    fprintf(stderr, "usage: haplo-predict OPTIONS [data-fname | <stdin>]\n");
    print_options(stderr, 27, NUM_OPTS_NO_ARG, opts_no_arg, NUM_OPTS_WITH_ARG,
            opts_with_arg);
}

Error* process_output_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'output' requires an argument");
    }
    output_fname = arg;
    return NULL;
}

static void init_predict_options(void)
{
    uint32_t i;

    char s_name;
    const char* l_name;
    const char* desc;

    Error* (*fnoarg)();
    Error* (*farg)(const char*);

    init_options(opts_no_arg, opts_with_arg);

    opts.label_col = LABEL_COL;

    i = NUM_SHARED_OPTS_NO_ARG;
    l_name = "exclude-one";
    s_name = 0;
    desc   = "When performing the tandem prediction decision, exclude at most one prediction from the set of classification algorithms. There must be three or more algorithms in play for this to take effect.";
    fnoarg = process_exclude_one_opt;
    init_option_no_arg(&(opts_no_arg[i++]), l_name, s_name, desc, fnoarg);
    assert(i == NUM_OPTS_NO_ARG);

    i = NUM_SHARED_OPTS_WITH_ARG;
    l_name = "output";
    s_name = 0;
    desc   = "File to output the predictions to. The default is stdout.";
    farg   = process_output_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "model-dir";
    s_name = 0;
    desc   = "Directory containing the trained models.";
    farg   = process_model_dir_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nb-freq";
    s_name = 0;
    desc   = "Naive Bayes non-parametric frequency model tree information.";
    farg   = process_nb_freq_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nb-freq-dtd";
    s_name = 0;
    desc   = "Validate the naive Bayes non-parametric frequency model tree information XML file with this DTD.";
    farg   = process_nb_freq_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nb-gauss";
    s_name = 0;
    desc   = "Naive Bayes Gaussian model tree information.";
    farg   = process_nb_gauss_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nb-gauss-dtd";
    s_name = 0;
    desc   = "Validate the naive Bayes Gaussian model tree information XML file with this DTD.";
    farg   = process_nb_gauss_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nb-gmm";
    s_name = 0;
    desc   = "Naive Bayes Gaussian mixture model tree information.";
    farg   = process_nb_gmm_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nb-gmm-dtd";
    s_name = 0;
    desc   = "Validate the naive Bayes Gaussian mixture model tree information XML file with this DTD.";
    farg   = process_nb_gmm_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "mv-gmm";
    s_name = 0;
    desc   = "Multivariate Gaussian mixture model tree information.";
    farg   = process_mv_gmm_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "mv-gmm-dtd";
    s_name = 0;
    desc   = "Validate the multivariate Gaussian mixture model tree information XML file with this DTD.";
    farg   = process_mv_gmm_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

#ifdef HAPLO_ENABLE_SVM
    l_name = "svm";
    s_name = 0;
    desc   = "SVM model tree information.";
    farg   = process_svm_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "svm-dtd";
    s_name = 0;
    desc   = "Validate the SVM model tree information XML file with this DTD.";
    farg   = process_svm_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);
#endif

#ifdef HAPLO_ENABLE_WEKA
    l_name = "weka-j48";
    s_name = 0;
    desc   = "Weka J48 model tree information.";
    farg   = process_weka_j48_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "weka-part";
    s_name = 0;
    desc   = "Weka PART model tree information.";
    farg   = process_weka_part_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "weka-jar";
    s_name = 0;
    desc   = "Weka java archive file. Required for using the Weka algorithms.";
    farg   = process_weka_jar_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "weka-dtd";
    s_name = 0;
    desc   = "Validate the Weka model tree information XML files with this DTD.";
    farg   = process_weka_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);
#endif

    l_name = "nearest";
    s_name = 0;
    desc   = "Nearest neighbor model information.";
    farg   = process_nearest_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nearest-dtd";
    s_name = 0;
    desc   = "Validate the nearest neighbor model information XML file with this DTD.";
    farg   = process_nearest_dtd_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "nearest-max-d";
    s_name = 0;
    desc   = "Maximum distance allowed for a nearest neighbor classification.";
    farg   = process_nearest_max_d_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);
    assert(i == NUM_OPTS_WITH_ARG);
}

static uint8_t num_models_to_predict()
{
    return (opts.nb_freq_fname       != 0) +
           (opts.nb_gauss_fname      != 0) +
           (opts.nb_gmm_fname        != 0) +
           (opts.mv_gmm_fname        != 0) +
#ifdef HAPLO_ENABLE_SVM
           (opts.svm_fname           != 0) +
#endif
#ifdef HAPLO_ENABLE_WEKA
           (opts.weka_j48_fname      != 0) +
           (opts.weka_part_fname     != 0) +
#endif
           (opts.nearest_fname       != 0);
}

static void predict_nb_freq
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
    NB_freq_model_tree* tree = NULL;
    Error* err;

    if (!opts.nb_freq_fname)
        return;

    if ((err = read_nb_freq_model_tree(&tree, opts.nb_freq_fname,
                    opts.nb_freq_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_nb_freq_model_tree(labels_out, 
                    confs_out, markers, tree, 0)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_nb_freq_model_tree(tree);
}

static void predict_nb_gauss
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
    NB_gauss_model_tree* tree = NULL;
    Error* err;

    if (!opts.nb_gauss_fname)
        return;

    if ((err = read_nb_gauss_model_tree(&tree, opts.nb_gauss_fname, 
                    opts.nb_gauss_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_nb_gauss_model_tree(labels_out, 
                    confs_out, markers, tree, 0)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_nb_gauss_model_tree(tree);
}

static void predict_nb_gmm
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
    NB_gmm_model_tree* tree = NULL;
    Error* err;

    if (!opts.nb_gmm_fname)
        return;

    if ((err = read_nb_gmm_model_tree(&tree, opts.nb_gmm_fname, 
                    opts.nb_gmm_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_nb_gmm_model_tree(labels_out, 
                    confs_out, markers, tree, 0)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_nb_gmm_model_tree(tree);
}

static void predict_mv_gmm
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
    MV_gmm_model_tree* tree = NULL;
    Error* err;

    if (!opts.mv_gmm_fname)
        return;

    if ((err = read_mv_gmm_model_tree(&tree, opts.mv_gmm_fname, 
                    opts.mv_gmm_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_mv_gmm_model_tree(labels_out, 
                    confs_out, markers, tree, 0)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_mv_gmm_model_tree(tree);
}

static void predict_svm
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
#ifdef HAPLO_ENABLE_SVM
    SVM_model_tree* tree = NULL;
    Error* err;

    if (!opts.svm_fname)
        return;

    if ((err = read_svm_model_tree(&tree, opts.svm_fname, 
                    opts.svm_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_svm_model_tree(labels_out, confs_out,
                    markers, tree)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_svm_model_tree(tree);
#else
    return;
#endif
}

static void predict_j48
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
#ifdef HAPLO_ENABLE_WEKA
    Weka_model_tree* tree = NULL;
    Error* err;

    if (!opts.weka_j48_fname)
        return;

    if ((err = read_weka_model_tree(&tree, opts.weka_j48_fname, 
                    opts.weka_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_weka_j48_model_tree(labels_out, 
                    confs_out, markers, tree, opts.weka_jar_fname)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_weka_model_tree(tree);
#else
    return;
#endif
}

static void predict_part
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
#ifdef HAPLO_ENABLE_WEKA
    Weka_model_tree* tree = NULL;
    Error* err;

    if (!opts.weka_part_fname)
        return;

    if ((err = read_weka_model_tree(&tree, opts.weka_part_fname, 
                    opts.weka_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_weka_part_model_tree(labels_out, 
                    confs_out, markers, tree, opts.weka_jar_fname)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_weka_model_tree(tree);
#else
    return;
#endif
}

static void predict_nearest
(
    Vector_u32**      labels_out,
    Vector_d**        confs_out,
    const Matrix_i32* markers
)
{
    Nearest_model* model = NULL;
    Error*         err;

    if (!opts.nearest_fname)
        return;

    if ((err = read_nearest_model(&model, opts.nearest_fname,
                    opts.nearest_dtd_fname, opts.model_dirname)) ||
        (err = predict_labels_with_nearest_model(labels_out, confs_out,
                    markers, model)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    free_nearest_model(model);
}

static void write_results_header
(
    const Matblock_u8* ids,
    const Vector_u32*  labels,
    const Vector_u32*  nb_freq_labels,
    const Vector_u32*  nb_gauss_labels,
    const Vector_u32*  nb_gmm_labels,
    const Vector_u32*  mv_gmm_labels,
    const Vector_u32*  svm_labels,
    const Vector_u32*  j48_labels,
    const Vector_u32*  part_labels,
    const Vector_u32*  nearest_labels,
    FILE*              fp
)
{
    uint32_t i;

    if (!opts.header_out || opts.output_format == HAPLO_OUTPUT_XML)
        return;

    if (ids)
    {
        for (i = 0; i < ids->num_rows; i++)
        {
            switch (opts.output_format)
            {
                case HAPLO_OUTPUT_TXT:
                    fprintf(fp, "ID %-7d  ", i+1);
                    break;
                case HAPLO_OUTPUT_CSV:
                    fprintf(fp, "ID %d,", i+1);
                    break;
                case HAPLO_OUTPUT_XML:
                    break;
            }
        }
    }

    if (labels)
    {
        switch (opts.output_format)
        {
            case HAPLO_OUTPUT_TXT:
                fprintf(fp, "%-10s  ", "Label");
                break;
            case HAPLO_OUTPUT_CSV:
                fprintf(fp, "%s,", "Label");
                break;
            case HAPLO_OUTPUT_XML:
                break;
        }
    }

    switch (opts.output_format)
    {
        case HAPLO_OUTPUT_TXT:
            fprintf(fp, "%-10s  %-4s", "Ancestor", "Type");
            if (nb_freq_labels) 
                fprintf(fp, "  %-10s  %-5s", "NB-Freq", "Conf");
            if (nb_gauss_labels) 
                fprintf(fp, "  %-10s  %-5s", "NB-Gauss", "Conf");
            if (nb_gmm_labels) 
                fprintf(fp, "  %-10s  %-5s", "NB-Gmm", "Conf");
            if (mv_gmm_labels) 
                fprintf(fp, "  %-10s  %-5s", "MV-Gmm", "Conf");
            if (svm_labels) 
                fprintf(fp, "  %-10s  %-5s", "SVM", "Conf");
            if (j48_labels) 
                fprintf(fp, "  %-10s  %-5s", "J48", "Conf");
            if (part_labels) 
                fprintf(fp, "  %-10s  %-5s", "PART", "Conf");
            if (nearest_labels) 
                fprintf(fp, "  %-10s  %-5s", "Nearest", "Dist");
            break;
        case HAPLO_OUTPUT_CSV:
            fprintf(fp, "%s,%s", "Ancestor", "Type");
            if (nb_freq_labels) 
                fprintf(fp, ",%s,%s", "NB-Freq", "Conf");
            if (nb_gauss_labels) 
                fprintf(fp, ",%s,%s", "NB-Gauss", "Conf");
            if (nb_gmm_labels) 
                fprintf(fp, ",%s,%s", "NB-Gmm", "Conf");
            if (mv_gmm_labels) 
                fprintf(fp, ",%s,%s", "MV-Gmm", "Conf");
            if (svm_labels) 
                fprintf(fp, ",%s,%s", "SVM", "Conf");
            if (j48_labels) 
                fprintf(fp, ",%s,%s", "J48", "Conf");
            if (part_labels) 
                fprintf(fp, ",%s,%s", "PART", "Conf");
            if (nearest_labels) 
                fprintf(fp, ",%s,%s", "Nearest", "Dist");
            break;
        case HAPLO_OUTPUT_XML:
            break;
    }
    fprintf(fp, "\n");
}

static void find_ancestors
(
    Vector_u32**      ancestor_types_out,
    Vector_u32**      ancestor_labels_out,
    const Vector_u32* labels_1,
    const Vector_u32* labels_2,
    const Vector_u32* labels_3,
    const Vector_u32* labels_4,
    const Vector_u32* labels_5,
    const Vector_u32* labels_6,
    const Vector_u32* labels_7,
    const Vector_u32* labels_8
)
{
    uint32_t            n, nn, N;
    uint32_t            i;
    uint32_t            num_labels;
    uint32_t            ancestor_label;
    Haplo_ancestor_type ancestor_type;
    Vector_u32*         labels;
    Vector_u32*         labelss;

    N = 0;

    if (labels_1) {N++; num_labels = labels_1->num_elts;}
    if (labels_2) {N++; num_labels = labels_2->num_elts;}
    if (labels_3) {N++; num_labels = labels_3->num_elts;}
    if (labels_4) {N++; num_labels = labels_4->num_elts;}
    if (labels_5) {N++; num_labels = labels_5->num_elts;}
    if (labels_6) {N++; num_labels = labels_6->num_elts;}
    if (labels_7) {N++; num_labels = labels_7->num_elts;}
    if (labels_8) {N++; num_labels = labels_8->num_elts;}

    assert(N > 0);

    labels = NULL;
    create_vector_u32(&labels, N);

    create_vector_u32(ancestor_types_out, num_labels);
    create_vector_u32(ancestor_labels_out, num_labels);

    for (i = 0; i < num_labels; i++)
    {
        N = 0;

        if (labels_1) labels->elts[ N++ ] = labels_1->elts[ i ];
        if (labels_2) labels->elts[ N++ ] = labels_2->elts[ i ];
        if (labels_3) labels->elts[ N++ ] = labels_3->elts[ i ];
        if (labels_4) labels->elts[ N++ ] = labels_4->elts[ i ];
        if (labels_5) labels->elts[ N++ ] = labels_5->elts[ i ];
        if (labels_6) labels->elts[ N++ ] = labels_6->elts[ i ];
        if (labels_7) labels->elts[ N++ ] = labels_7->elts[ i ];
        if (labels_8) labels->elts[ N++ ] = labels_8->elts[ i ];

        ancestor_type = find_ancestor_index_of_set(&ancestor_label, labels);

        if (opts.exclude_one && ancestor_type == HAPLO_ANCESTOR_NONE && N > 3)
        {
            labelss = NULL;
            create_vector_u32(&labelss, N-1);
            for (n = 0; ancestor_type == HAPLO_ANCESTOR_NONE && n < N; n++)
            {
                for (nn = 0; nn < N-1; nn++)
                {
                    if (nn < n)
                    {
                        labelss->elts[ nn ] = labels->elts[ nn ];
                    }
                    else
                    {
                        labelss->elts[ nn ] = labels->elts[ nn+1 ];
                    }
                }

                ancestor_type = find_ancestor_index_of_set(&ancestor_label,
                        labelss);
            }
            free_vector_u32(labelss);
        }

        (*ancestor_types_out)->elts[ i ] = ancestor_type;
        (*ancestor_labels_out)->elts[ i ] = ancestor_label;
    }
}

static void write_results
(
    const Matblock_u8* ids,
    const Vector_u32*  labels,
    const Matrix_i32*  markers,
    const Vector_u32*  nb_freq_labels,
    const Vector_d*    nb_freq_confs,
    const Vector_u32*  nb_gauss_labels,
    const Vector_d*    nb_gauss_confs,
    const Vector_u32*  nb_gmm_labels,
    const Vector_d*    nb_gmm_confs,
    const Vector_u32*  mv_gmm_labels,
    const Vector_d*    mv_gmm_confs,
    const Vector_u32*  svm_labels,
    const Vector_d*    svm_confs,
    const Vector_u32*  j48_labels,
    const Vector_d*    j48_confs,
    const Vector_u32*  part_labels,
    const Vector_d*    part_confs,
    const Vector_u32*  nearest_labels,
    const Vector_d*    nearest_dists,
    const Vector_u32*  ancestor_types,
    const Vector_u32*  ancestor_labels
)
{
    uint32_t    i;
    FILE*       fp        = NULL;
    xmlDoc*     xml_doc   = NULL;
    xmlNode*    xml_root  = NULL;
    xmlNode*    xml_node  = NULL;
    char*       xml_buf;
    uint32_t    xml_len;
    Error*      err;

    if ((err = open_output(&fp, &xml_doc, "haplo-predict-out",
                    "haplo-predict-out.dtd", output_fname)))
    {
        print_error_msg("haplo-predict", err->msg);
    }

    write_results_header(ids, labels, nb_freq_labels, nb_gauss_labels,
            nb_gmm_labels, mv_gmm_labels, svm_labels, j48_labels, part_labels,
            nearest_labels, fp);

    if (opts.output_format == HAPLO_OUTPUT_XML)
    {
        xml_root = xmlDocGetRootElement(xml_doc);
    }

    for (i = 0; i < markers->num_rows; i++)
    {
        if (opts.output_format == HAPLO_OUTPUT_XML)
        {
            xml_node = XMLNewChild(xml_root, "sample", NULL);
            xml_len = 256;
            xml_buf = malloc(xml_len*sizeof(char));
            snprintf(xml_buf, xml_len, "%d", i+1);
            XMLNewProp(xml_node, "number", xml_buf);
            free(xml_buf);
        }

        write_ids(ids, i, HAPLO_SEP_SUFFIX, fp, xml_node);
        write_label(labels, i, HAPLO_SEP_SUFFIX, fp, xml_node);

        write_ancestor_label(ancestor_types, ancestor_labels, i, 
                HAPLO_SEP_NONE, fp, xml_node);

        write_prediction("nb-freq", nb_freq_labels, nb_freq_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("nb-gauss", nb_gauss_labels, nb_gauss_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("nb-gmm", nb_gmm_labels, nb_gmm_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("mv-gmm", mv_gmm_labels, mv_gmm_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("svm", svm_labels, svm_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("j48", j48_labels, j48_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("part", part_labels, part_confs, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        write_prediction("nearest", nearest_labels, nearest_dists, i, 
                HAPLO_SEP_PREFIX, fp, xml_node);

        if (opts.output_format != HAPLO_OUTPUT_XML)
        {
            fprintf(fp, "\n");
        }
    }

    if ((err = close_output(fp, xml_doc, output_fname)))
    {
        print_error_msg("haplo-predict", err->msg);
    }
}


static void predict(Predict_params* p)
#ifdef HAPLO_HAVE_PTHREAD
{
    uint32_t n;
    typedef void (*f)(Vector_u32**, Vector_d**, const Matrix_i32*);

    f predict_f[NUM_ALGOS] = { 
                        predict_nb_freq, predict_nb_gauss, 
                        predict_nb_gmm, predict_mv_gmm, 
                        predict_svm, predict_j48,
                        predict_part, predict_nearest };

    pthread_mutex_lock(&(p->mutex));
    n = p->n;

    while (n < NUM_ALGOS)
    {
        pthread_mutex_unlock(&(p->mutex));

        p->labels[ n ] = NULL;
        p->confs[ n ] = NULL;
        predict_f[ n ](&(p->labels[ n ]), &(p->confs[ n ]), p->markers);

        pthread_mutex_lock(&(p->mutex));
        n = ++(p->n);
    }

    pthread_mutex_unlock(&(p->mutex));
}
#else
{
    uint32_t n;
    typedef void (*f)(Vector_u32**, Vector_d**, const Matrix_i32*);

    f predict_f[NUM_ALGOS] = { 
                        predict_nb_freq, predict_nb_gauss, 
                        predict_nb_gmm, predict_mv_gmm, 
                        predict_svm, predict_j48,
                        predict_part, predict_nearest };

    for (n = p->n; n < NUM_ALGOS; n++)
    {
        p->labels[ n ] = NULL;
        p->confs[ n ] = NULL;
        predict_f[ n ](&(p->labels[ n ]), &(p->confs[ n ]), p->markers);
    }
}
#endif


int main(int argc, const char** argv)
{
    int         argi;
    const char* data_fname = "/dev/stdin";
    Error*      err;

    uint32_t n;
    typedef void* (*f)(void*);
    Predict_params p;

#ifdef HAPLO_HAVE_PTHREAD
    uint32_t t;
    pthread_t* threads;
#endif

    Matblock_u8* ids             = NULL;
    Vector_u32*  labels          = NULL;
    Matrix_i32*  markers         = NULL;
    Vector_u32*  ancestor_types  = NULL;
    Vector_u32*  ancestor_labels = NULL;

    Vector_u32*  pred_labels[NUM_ALGOS] = {NULL};
    Vector_d*    pred_confs[NUM_ALGOS]  = {NULL};

    init_predict_options();

    if ((err = process_options(argc, argv, &argi, NUM_OPTS_NO_ARG, opts_no_arg,
                    NUM_OPTS_WITH_ARG, opts_with_arg)) != NULL)
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    if ((argc - argi) == 1)
    {
        data_fname = argv[ argi ];
    }

    if (num_models_to_predict() == 0)
    {
        print_error_msg_exit("haplo-predict", "No models to predict with");
    }

    if ((err = read_haplo_groups(opts.labels_fname)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    if ((err = read_input(&ids, &labels, &markers, data_fname)))
    {
        print_error_msg_exit("haplo-predict", err->msg);
    }

    p.labels  = pred_labels;
    p.confs   = pred_confs;
    p.markers = markers;
    p.n       = 0;

#ifdef HAPLO_HAVE_PTHREAD
    pthread_mutex_init(&(p.mutex), NULL);
    assert(threads = malloc(opts.num_threads*sizeof(pthread_t)));
    for (t = 0; t < opts.num_threads; t++)
    {
        pthread_create(&(threads[ t ]), NULL, (f)predict, &p);
    }
    for (t = 0; t < opts.num_threads; t++)
    {
        pthread_join(threads[ t ], NULL);
    }
    free(threads);
#else
    predict(&p);
#endif

    find_ancestors(&ancestor_types, &ancestor_labels, pred_labels[0],
            pred_labels[1], pred_labels[2], pred_labels[3],
            pred_labels[4], pred_labels[5], pred_labels[6],
            pred_labels[7]);

    write_results(ids, labels, markers, pred_labels[0], pred_confs[0],
            pred_labels[1], pred_confs[1], pred_labels[2], pred_confs[2],
            pred_labels[3], pred_confs[3], pred_labels[4], pred_confs[4],
            pred_labels[5], pred_confs[5], pred_labels[6], pred_confs[6],
            pred_labels[7], pred_confs[7], ancestor_types, ancestor_labels);

    free_matblock_u8(ids);
    free_vector_u32(labels);
    free_matrix_i32(markers);
    free_vector_u32(ancestor_types);
    free_vector_u32(ancestor_labels);

    for (n = 0; n < NUM_ALGOS; n++)
    {
        free_vector_u32(pred_labels[ n ]);
        free_vector_d(pred_confs[ n ]);
    }

    if (get_num_unhandled_errors() > 0)
    {
        print_error_msg_exit("haplo-predict", "Unhandled errors");
    }

    return EXIT_SUCCESS;
}