svm_tree.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 <inttypes.h>
#include <assert.h>
#include <string.h>
#include <errno.h>
#include <math.h>

#include <libxml/tree.h>
#include <libxml/parser.h>
#include <libxml/valid.h>

#if defined HAPLO_HAVE_LIBSVM_H
#include <libsvm.h>
#elif defined HAPLO_HAVE_SVM_H
#include <svm.h>
#endif

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

#include <jwsc/base/error.h>
#include <jwsc/base/limits.h>
#include <jwsc/base/file_io.h>
#include <jwsc/vector/vector.h>
#include <jwsc/vector/vector_io.h>
#include <jwsc/vector/vector_math.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matrix/matrix_io.h>
#include <jwsc/matblock/matblock.h>
#include <jwsc/matblock/matblock_io.h>
#include <jwsc/stat/gmm.h>

#include "xml.h"
#include "haplo_groups.h"
#include "svm_tree.h"


void train_svm_model
(
    SVM_model**        model_out,
    const Vector_u32*  labels,
    const Matrix_i32*  markers,
    double             cost,
    double             gamma
)
{
    uint32_t s, num_samples;
    uint32_t m, num_markers;
    int32_t  marker_val;

    SVM_model*            model;
    struct svm_parameter* param;

    if (*model_out != NULL)
    {
        free_svm_model(*model_out);
    }

    assert(*model_out = malloc(sizeof(SVM_model)));
    model = *model_out;

    num_samples = markers->num_rows;
    num_markers = markers->num_cols;

    assert(model->prob = malloc(sizeof(struct svm_problem)));
    assert(model->prob->y = malloc(num_samples*sizeof(double)));
    assert(model->prob->x = malloc(num_samples*sizeof(struct svm_node*)));

    model->prob->l = num_samples;
    for (s = 0; s < num_samples; s++)
    {
        model->prob->y[ s ] = labels->elts[ s ];
        model->prob->x[ s ] = malloc((num_markers+1)*sizeof(struct svm_node));
        assert(model->prob->x[ s ]);

        for (m = 0; m < num_markers; m++)
        {
            model->prob->x[ s ][ m ].index = m+1;
            marker_val = markers->elts[ s ][ m ];
            model->prob->x[ s ][ m ].value = marker_val;
        }
        model->prob->x[ s ][ m ].index = -1;
    }

    assert(param = malloc(sizeof(struct svm_parameter)));

    param->svm_type = C_SVC;
    param->kernel_type = RBF;
    param->gamma = gamma;
    param->C = cost;
    param->nr_weight = 0;
    param->probability = 1;
    param->shrinking = 1;       /* svm-train default */
    param->eps = 0.001;         /* svm-train default */
    param->cache_size = 40;     /* svm-train default */

    assert(model->svm = svm_train(model->prob, param));

    free(param);
}


Error* predict_label_with_svm_model
(
    uint32_t*               label_out,
    double**                confidence_out,
    const struct svm_node*  markers,
    const SVM_model*        model
)
{
    assert(*confidence_out);
    *label_out = svm_predict_probability(model->svm, markers, *confidence_out);

    return NULL;
}


static void get_predicted_labels_confidence
(
    Vector_d**        confidence_out,
    const Vector_u32* labels_v,
    const Matrix_d*   confidence_matrix, 
    const SVM_model*  model
)
{
    uint32_t  s, num_samples;
    uint32_t  l, num_labels;
    uint32_t  label;
    int*      labels;
    Vector_d* confidence;

    num_samples = labels_v->num_elts;

    create_zero_vector_d(confidence_out, num_samples);
    confidence = *confidence_out;

    num_labels = (uint32_t)svm_get_nr_class(model->svm);
    assert(labels = malloc(num_labels*sizeof(uint32_t)));
    svm_get_labels(model->svm, labels);

    for (s = 0; s < num_samples; s++)
    {
        label = labels_v->elts[ s ];

        for (l = 0; l < num_labels; l++)
        {
            if ((uint32_t)labels[ l ] == label)
            {
                confidence->elts[ s ] = confidence_matrix->elts[ s ][ l ];
            }
        }
    }

    free(labels);
}


Error* predict_labels_with_svm_model
(
    Vector_u32**      labels_out,
    Vector_d**        confidence_out,
    const Matrix_i32* markers,
    const SVM_model*  model
)
{
    uint32_t s, num_samples;
    uint32_t m, num_markers;

    Vector_u32* labels;
    Matrix_d*   confidence_matrix;

    struct svm_node* markers_v;

    num_samples = markers->num_rows;
    num_markers = markers->num_cols;

    assert(markers_v = malloc((num_markers+1)*sizeof(struct svm_node)));

    create_zero_vector_u32(labels_out, num_samples);
    labels = *labels_out;

    confidence_matrix = NULL;
    create_zero_matrix_d(&confidence_matrix, num_samples, 
            svm_get_nr_class(model->svm));

    for (s = 0; s < num_samples; s++)
    {
        for (m = 0; m < num_markers; m++)
        {
            markers_v[ m ].index = m+1;
            markers_v[ m ].value = markers->elts[ s ][ m ];
        }
        markers_v[ m].index = -1;

        predict_label_with_svm_model(&(labels->elts[ s ]),
                &(confidence_matrix->elts[ s ]), markers_v, model);
    }

    get_predicted_labels_confidence(confidence_out, labels, confidence_matrix, 
            model);

    free(markers_v);
    free_matrix_d(confidence_matrix);

    return NULL;
}


Error* read_svm_model(SVM_model** model_out, const char* fname)
{
    int  slen = 256;
    char str[slen];

    if (*model_out)
    {
        free_svm_model(*model_out);
    }

    assert(*model_out = malloc(sizeof(SVM_model)));
    (*model_out)->prob = NULL;

    if (!((*model_out)->svm = svm_load_model(fname)))
    {
        snprintf(str, slen, "%s: %s", fname, "Could not read model");
        return JWSC_EARG(str);
    }

    return NULL;
}


Error* write_svm_model(SVM_model* model, const char* fname)
{
    int  slen = 256;
    char str[slen];

    if (svm_save_model(fname, model->svm) < 0)
    {
        snprintf(str, slen, "%s: %s", fname, "Could not write model");
        return JWSC_EARG(str);
    }

    return NULL;
}


Error* write_svm_model_training_data
(
    const Vector_u32* labels,
    const Matrix_i32* markers,
    const char*       fname
)
{
    FILE*    fp;
    uint32_t s, num_samples;
    uint32_t m, num_markers;
    int  slen = 256;
    char str[slen];

    num_samples = markers->num_rows;
    num_markers = markers->num_cols;

    if ((fp = fopen(fname, "w")) == NULL)
    {
        snprintf(str, slen, "%s: %s", fname, strerror(errno));
        return JWSC_EIO(str);
    }

    for (s = 0; s < num_samples; s++)
    {
        fprintf(fp, "%-4d", labels->elts[ s ]);
        for (m = 0; m < num_markers; m++)
        {
            fprintf(fp, " %2d:%-4d", m+1, markers->elts[ s ][ m ]);
        }
        fprintf(fp, "\n");
    }

    if (fclose(fp) != 0)
    {
        snprintf(str, slen, "%s: %s", fname, strerror(errno));
        return JWSC_EIO(str);
    }

    return NULL;
}


void free_svm_model(SVM_model* model)
{
    uint32_t i;

    if (!model)
        return;

    svm_free_and_destroy_model(&(model->svm));

    if (model->prob)
    {
        for (i = 0; i < model->prob->l; i++ )
        {
            free(model->prob->x[ i ]);
        }
        free(model->prob->x);
        free(model->prob->y);
        free(model->prob);
    }

    free(model);
}


static Error* read_svm_xml_doc
(
    xmlDoc**    xml_doc_out,
    const char* xml_fname,
    const char* dtd_fname
)
{
    xmlParserCtxt* xml_parse_ctxt;
    xmlValidCtxt*  xml_valid_ctxt;
    xmlDtd*        xml_dtd;
    int            slen = 256;
    char           str[slen];

    assert(xml_parse_ctxt = xmlNewParserCtxt());

    if (!(*xml_doc_out = xmlCtxtReadFile(xml_parse_ctxt, xml_fname, NULL, 0)))
    {
        snprintf(str, slen, "%s: %s", xml_fname, "Could not parse file");
        return JWSC_EARG(str);
    } 

    xmlFreeParserCtxt(xml_parse_ctxt);

    if (dtd_fname)
    {
        assert(xml_valid_ctxt = xmlNewValidCtxt());

        if (!(xml_dtd = xmlParseDTD(NULL, (xmlChar*)dtd_fname)))
        {
            snprintf(str, slen, "%s: %s", dtd_fname, "Could not parse DTD");
            return JWSC_EARG(str);
        }

        if (!xmlValidateDtd(xml_valid_ctxt, *xml_doc_out, xml_dtd))
        {
            snprintf(str, slen, "%s: %s", xml_fname, "XML file not valid");
            return JWSC_EARG(str);
        }

        xmlFreeValidCtxt(xml_valid_ctxt);
        xmlFreeDtd(xml_dtd);
    }

    return NULL;
}

static Error* create_svm_model_tree_from_xml_node
(
     SVM_model_tree**      tree_out, 
     const SVM_model_tree* parent,
     uint32_t              parent_label,
     xmlNode*              xml_node
);

static Error* create_svm_model_node_from_xml_node
(
     xmlNode*        xml_node,
     SVM_model_node* svm_node,
     uint32_t        i
)
{
    uint32_t    label;
    uint32_t    altlabel;
    uint32_t    num_altlabels;
    uint32_t    g;
    uint32_t    a;
    const char* fname;
    xmlNode*    it;
    xmlNode*    itt;
    Error*      err;
    int         slen = 256;
    char        str[slen];

    g = 0;
    for (it = xml_node; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "file"))
        {
            fname = (const char*) it->children->content;
            svm_node->model_fnames[ i ] = malloc((strlen(fname)+1) * 
                    sizeof(char));
            strcpy(svm_node->model_fnames[ i ], fname);
        }
        else if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "cost"))
        {
            if (sscanf((char*)it->children->content, "%lf", 
                        &(svm_node->cost->elts[ i ])) != 1)
            {
                snprintf(str, slen, "%s: %s", fname, "Invalid cost");
                return JWSC_EARG(str);
            }
        }
        else if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "gamma"))
        {
            if (sscanf((char*)it->children->content, "%lf", 
                        &(svm_node->gamma->elts[ i ])) != 1)
            {
                snprintf(str, slen, "%s: %s", fname, "Invalid gamma");
                return JWSC_EARG(str);
            }
        }
        else if (it->type == XML_ELEMENT_NODE && 
                (XMLStrEqual(it->name, "group-one") ||
                 XMLStrEqual(it->name, "group-all")))
        {
            a = 0;
            num_altlabels = 0;
            for (itt = it->children; itt; itt = itt->next)
            {
                if (itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "altlabel"))
                {
                    num_altlabels++;
                }
            }
            if (num_altlabels)
            {
                create_vector_u32(&(svm_node->altlabels[ g ][ i ]), 
                        num_altlabels);
            }

            for (itt = it->children; itt; itt = itt->next)
            {
                if (itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "label"))
                {
                    if ((err = lookup_haplo_group_index_from_label(&label, 
                                    (const char*) itt->children->content)))
                    {
                        snprintf(str, slen, "%s: %s", fname, err->msg);
                        return JWSC_EARG(str);
                    }
                    svm_node->labels[ g ]->elts[ i ] = label;
                }
                else if (itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "altlabel"))
                {
                    if ((err = lookup_haplo_group_index_from_label(&altlabel, 
                                    (const char*) itt->children->content)))
                    {
                        snprintf(str, slen, "%s: %s", fname, err->msg);
                        return JWSC_EARG(str);
                    }
                    svm_node->altlabels[ g ][ i ]->elts[ a++ ] = altlabel;
                }
                else if (itt->type == XML_ELEMENT_NODE && 
                         (XMLStrEqual(itt->name, "binary-model") ||
                          XMLStrEqual(itt->name, "one-vs-all-model") ||
                          XMLStrEqual(itt->name, "one-vs-one-model")))
                {
                    if ((err = create_svm_model_tree_from_xml_node(
                                    &(svm_node->subtrees[ g ][ i ]), svm_node,
                                    label, itt)))
                    {
                        snprintf(str, slen, "%s: %s", fname, err->msg);
                        return JWSC_EARG(str);
                    }
                    break;
                }
            }
            g++;
        }
    }

    return NULL;
}


static Error* create_svm_model_tree_from_xml_node
(
     SVM_model_tree**      tree_out, 
     const SVM_model_tree* parent,
     uint32_t              parent_label,
     xmlNode*              xml_node
)
{
    SVM_model_tree* tree;
    uint32_t i;
    xmlNode* it;
    Error*   err;

    if (*tree_out)
    {
        free_svm_model_tree(*tree_out);
    }

    assert(*tree_out = malloc(sizeof(SVM_model_tree)));
    tree = *tree_out;

    tree->parent       = parent;
    tree->parent_label = parent_label;
    tree->subtrees[0]  = NULL;
    tree->subtrees[1]  = NULL;
    tree->labels[0]    = NULL;
    tree->labels[1]    = NULL;
    tree->altlabels[0] = NULL;
    tree->altlabels[1] = NULL;
    tree->cost         = NULL;
    tree->gamma        = NULL;
    tree->num_models   = 0;
    tree->models       = NULL;
    tree->model_fnames = NULL;

    assert(XMLStrEqual(xml_node->name, "binary-model") ||
           XMLStrEqual(xml_node->name, "one-vs-all-model") ||
           XMLStrEqual(xml_node->name, "one-vs-one-model"));

    for (it = xml_node; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && 
                (XMLStrEqual(it->name, "binary-model") ||
                 XMLStrEqual(it->name, "one-vs-all-model") ||
                 XMLStrEqual(it->name, "one-vs-one-model")))
        {
            tree->num_models++;
        }
    }

    assert(tree->subtrees[0] = calloc(tree->num_models, sizeof(void*)));
    assert(tree->subtrees[1] = calloc(tree->num_models, sizeof(void*)));
    assert(tree->altlabels[0] = calloc(tree->num_models, sizeof(void*)));
    assert(tree->altlabels[1] = calloc(tree->num_models, sizeof(void*)));
    assert(tree->model_fnames = calloc(tree->num_models, sizeof(void*)));
    assert(tree->models = calloc(tree->num_models, sizeof(void*)));

    create_vector_u32(&(tree->labels[0]), tree->num_models);
    create_vector_u32(&(tree->labels[1]), tree->num_models);
    create_vector_d(&(tree->cost), tree->num_models);
    create_vector_d(&(tree->gamma), tree->num_models);

    i = 0;
    for (it = xml_node; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && 
                (XMLStrEqual(it->name, "binary-model") ||
                 XMLStrEqual(it->name, "one-vs-all-model") ||
                 XMLStrEqual(it->name, "one-vs-one-model")))
        {
            if ((err = create_svm_model_node_from_xml_node(it->children, tree, 
                        i++)))
            {
                return err;
            }
        }
    }
    assert(i == tree->num_models);

    return NULL;
}


static Error* create_model_training_data
(
    Vector_u32**       train_labels_out,
    Matrix_i32**       train_markers_out,
    const Vector_u32*  data_labels,
    const Matrix_i32*  data_markers,
    uint32_t*          model_labels,
    const Vector_u32** model_altlabels
)
{
    uint8_t     b;
    uint32_t    i, j, k;
    uint32_t    n;
    uint32_t    m[2] = {0};
    const char* label_str;
    Vector_u32* train_labels = NULL;
    Matrix_i32* train_markers = NULL;
    int         slen = 256;
    char        str[slen];

    copy_vector_u32(&train_labels, data_labels);
    copy_matrix_i32(&train_markers, data_markers);

    n = 0;
    for (i = 0; i < data_labels->num_elts; i++)
    {
        for (j = 0; j < 2; j++)
        {
            if (is_ancestor(data_labels->elts[ i ], model_labels[ j ]))
            {
                train_labels->elts[ n ] = model_labels[ j ];
                copy_matrix_block_into_matrix_i32(train_markers, n, 0,
                        data_markers, i, 0, 1, data_markers->num_cols);
                n++;
                m[j]++;
                break;
            }
            else if (model_altlabels[ j ])
            {
                for (k = 0; k < model_altlabels[ j ]->num_elts; k++)
                {
                    if ((b = is_ancestor(data_labels->elts[ i ], 
                                    model_altlabels[ j ]->elts[ k ])))
                    {
                        train_labels->elts[ n ] = model_labels[ j ];
                        copy_matrix_block_into_matrix_i32(train_markers, n, 0,
                                data_markers, i, 0, 1, data_markers->num_cols);
                        n++;
                        m[j]++;
                        break;
                    }
                }
                if (b)
                {
                    break;
                }
            }
        }
    }

    assert(n == (m[0] + m[1]));

    if (!(m[j=0]) || !(m[j=1]))
    {
        lookup_haplo_group_label_from_index(&label_str, model_labels[j]);
        snprintf(str, slen, "No data for model label %s", label_str);
        return JWSC_EARG(str);
    }

    copy_vector_section_u32(train_labels_out, train_labels, 0, n);
    copy_matrix_block_i32(train_markers_out, train_markers, 0, 0, n, 
            train_markers->num_cols);

    free_vector_u32(train_labels);
    free_matrix_i32(train_markers);

    return NULL;
}


static Error* train_svm_model_node
(
    SVM_model_node*   node,
    const Vector_u32* labels, 
    const Matrix_i32* markers
)
{
    uint32_t          i, j;
    uint32_t          label[2];
    const Vector_u32* altlabel[2];
    Vector_u32*       node_labels  = NULL;
    Matrix_i32*       node_markers = NULL;
    Error*            err;
    int               slen = 256;
    char              str[slen];

    for (i = 0; i < node->num_models; i++)
    {
        label[0] = node->labels[0]->elts[ i ];
        label[1] = node->labels[1]->elts[ i ];

        altlabel[0] = node->altlabels[0][ i ];
        altlabel[1] = node->altlabels[1][ i ];

        if ((err = create_model_training_data(&node_labels, &node_markers, 
                        labels, markers, label, altlabel)))
        {
            snprintf(str, slen, "%s: %s", node->model_fnames[i], err->msg);
            return JWSC_EARG(str);
        }

        train_svm_model(&(node->models[ i ]), node_labels, node_markers,
                node->cost->elts[ i ], node->gamma->elts[ i ]);

        for (j = 0; j < 2; j++)
        {
            if (node->subtrees[ j ][ i ])
            {
                if ((err = train_svm_model_node(node->subtrees[ j ][ i ],
                                labels, markers)))
                {
                    return err;
                }
            }
        }
    }

    free_vector_u32(node_labels);
    free_matrix_i32(node_markers);

    return NULL;
}


Error* train_svm_model_tree
(
    SVM_model_tree**  tree_out,
    const Vector_u32* labels, 
    const Matrix_i32* markers,
    const char*       tree_xml_fname,
    const char*       tree_dtd_fname
)
{
    xmlDoc*  xml_doc;
    xmlNode* xml_node = NULL;
    xmlNode* it;
    Error*   err;
    int      slen = 256;
    char     str[slen];

    assert(tree_xml_fname);

    if ((err = read_svm_xml_doc(&xml_doc, tree_xml_fname, tree_dtd_fname)))
    {
        return err;
    }

    for (it = xmlDocGetRootElement(xml_doc)->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && 
                (XMLStrEqual(it->name, "binary-model") ||
                 XMLStrEqual(it->name, "one-vs-all-model") ||
                 XMLStrEqual(it->name, "one-vs-one-model")))
        {
            xml_node = it;
            break;
        }
    }

    assert(xml_node);

    if ((err = create_svm_model_tree_from_xml_node(tree_out, NULL, 0, 
                    xml_node)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    if ((err = train_svm_model_node(*tree_out, labels, markers)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    return NULL;
}


static double get_predicted_label_confidence
(
    uint32_t         label,
    double*          confs,
    const SVM_model* model
)   
{   
    uint32_t l, num_labels;
    int*     labels;
    double   conf;

    num_labels = (uint32_t)svm_get_nr_class(model->svm);
    assert(labels = malloc(num_labels*sizeof(int)));
    svm_get_labels(model->svm, labels);

    for (l = 0; l < num_labels; l++)
    {
        if (labels[ l ] == label)
        {
           conf = confs[ l ];
        }
    }

    free(labels);

    return conf;
}


static uint32_t get_best_one_against_all_model
(
    const SVM_model_tree* tree,
    Vector_u32*           labels,
    Vector_d*             confs
)
{
    uint32_t    label, other_label;
    uint32_t    num_models;
    uint32_t    m;
    uint32_t    num_labels;
    int32_t     best_model           = -1;
    double      max_conf             = 0;
    double      min_conf             = 1;
    double      label_conf;
    double      normalizing_c        = 0;

    Vector_d* label_conf_sum = NULL;

    num_models = tree->num_models;
    num_labels = get_num_haplo_groups();

    assert(num_models > 1);

    create_init_vector_d(&label_conf_sum, num_labels, 0);

    /* Use the sum the confidence values for each predicted label to decide 
     * which is best. If there is a tie, the first occuring label is choosen. */
    for (m = 0; m < num_models; m++)
    {
        label = labels->elts[ m ];
        label_conf = confs->elts[ m ];

        if (tree->labels[0]->elts[ m ] == label)
        {
            label_conf_sum->elts[ label ] += label_conf;
            other_label = tree->labels[1]->elts[ m ];
        }
        else
        {
            label_conf_sum->elts[ tree->labels[0]->elts[ m ] ] += 1-label_conf;
            other_label = tree->labels[0]->elts[ m ];
        }
    }

    for (m = 0; m < num_models; m++)
    {
        label = labels->elts[ m ];
        label_conf = label_conf_sum->elts[ label ];

        if ((label == tree->labels[0]->elts[ m ]) && (max_conf < label_conf))
        {
            max_conf = label_conf;
            best_model = m;
        }
    }

    /* It's possible that nothing was choosen because the 'other_label' was  
     * selected in all groups. In this case, choose the predicted label whose
     * model had the least confidence in its 'other_label' prediction. */
    if (best_model < 0)
    {
        for (m = 0; m < num_models; m++)
        {
            label = labels->elts[ m ];
            assert(label == tree->labels[1]->elts[ m ]);
            label_conf = confs->elts[ m ];

            if (min_conf > label_conf)
            {
                min_conf = label_conf;
                best_model = m;
            }
        }

        /* Hack the predicted_labels and confidence data structures to reflect
         * that 'other_label' was not choosen. */
        labels->elts[ best_model ] = tree->labels[0]->elts[ best_model ];

        label_conf = confs->elts[ best_model ];
        confs->elts[ best_model ] = 1.0 - label_conf;
    }

    /* Need to normalize the confidence value for the best label prediction so 
     * that it is a probability. At this level in the tree, one of the labels 
     * will be chosen, so the confidence value for these labels needs to sum
     * to one. */
    for (label = 0; label < num_labels; label++)
    {
        normalizing_c += label_conf_sum->elts[ label ];
    }
    assert(normalizing_c > 1.0e-16);
    confs->elts[ best_model ] /= normalizing_c;

    free_vector_d(label_conf_sum);

    return best_model;
}


static uint32_t get_best_one_against_one_model
(
    const SVM_model_tree* tree,
    Vector_u32*           labels,
    Vector_d*             confs
)
{
    uint32_t    m, num_models;
    uint32_t    l, num_labels;
    uint32_t    best_model;
    uint32_t    max;
    uint32_t    the_winner;
    uint32_t   num_winners;
    double max_conf;

    Vector_u32* winners = NULL;

    num_models = tree->num_models;
    num_labels = get_num_haplo_groups();

    assert(num_models > 1);

    create_init_vector_u32(&winners, num_labels, 0);

    for (m = 0; m < num_models; m++)
    {
        winners->elts[ labels->elts[ m ] ]++;
    }

    max = 0;
    the_winner = 0;
    num_winners = 0;

    for (l = 0; l < num_labels; l++)
    {
        if (winners->elts[ l ] > max)
        {
            the_winner = l;
            max = winners->elts[ l ];
        }
    }
    for (l = 0; l < num_labels; l++)
    {
        if (max == winners->elts[ l ])
        {
            num_winners++;
        }
    }

    assert(max > 0);

    max_conf = 0;

    for (m = 0; m < num_models; m++)
    {
        if (labels->elts[ m ] == the_winner && confs->elts[ m ] > max_conf)
        {
            best_model = m;
        }
    }

    free_vector_u32(winners);

    return best_model;
}


static Error* recursively_predict_label_in_model_tree
(
    uint32_t*                label_out,
    double*                  confidence_out,
    const SVM_model_tree*    tree,
    const struct svm_node*   markers
)
{
    uint32_t i;
    uint32_t m, num_models;
    uint32_t subtree_label;
    Error*   err;

    Vector_u32* labels = NULL;
    Vector_d*   confs  = NULL;
    Vector_d*   v      = NULL;

    num_models = tree->num_models;

    create_vector_u32(&labels, num_models);
    create_vector_d(&confs, num_models);

    for (m = 0; m < num_models; m++)
    {
        create_zero_vector_d(&v, svm_get_nr_class(tree->models[ m ]->svm));

        if ((err = predict_label_with_svm_model(&(labels->elts[ m ]),
                        &(v->elts), markers, tree->models[ m ])))
        {
            return err;
        }

        confs->elts[ m ] = get_predicted_label_confidence(labels->elts[ m ], 
                v->elts, tree->models[ m ]);
    }
    free_vector_d(v);

    // If there are multiple models in the node, get the prediction from
    // the best one.
    m = (num_models > 1) ? get_best_one_against_all_model(tree, labels, confs) 
        : 0;

    *label_out = labels->elts[ m ];
    *confidence_out *= confs->elts[ m ];

    free_vector_u32(labels);
    free_vector_d(confs);

    // If the best model has a subtree, recursively predict in the subtree.
    for (m = 0; m < num_models; m++)
    {
        for (i = 0; i < 2; i++)
        {
            if (!(tree->subtrees[ i ][ m ]))
                continue;

            subtree_label = tree->subtrees[ i ][ m ]->parent_label;

            if (subtree_label == *label_out)
            {
                if ((err = recursively_predict_label_in_model_tree(label_out,
                                confidence_out, tree->subtrees[ i ][ m ],
                                markers)))
                {
                    return err;
                }
            }
        }
    }

    return NULL;
}


Error* predict_labels_with_svm_model_tree
(
    Vector_u32**          labels_out,
    Vector_d**            confidence_out,
    const Matrix_i32*     markers,
    const SVM_model_tree* tree
)
{
    uint32_t s, num_samples;
    uint32_t m, num_markers;

    Vector_u32*      labels;
    Vector_d*        confidence;
    struct svm_node* markers_v;
    Error*           e;

    num_samples = markers->num_rows;
    num_markers = markers->num_cols;

    create_vector_u32(labels_out, num_samples);
    labels = *labels_out;

    create_init_vector_d(confidence_out, num_samples, 1.0);
    confidence = *confidence_out;

    markers_v = NULL;
    assert(markers_v = malloc((num_markers+1)*sizeof(struct svm_node)));

    for (s = 0; s < num_samples; s++)
    {
        for (m  = 0; m  < num_markers; m ++)
        {
            markers_v[ m  ].index = m +1;
            markers_v[ m  ].value = markers->elts[ s ][ m  ];
        }
        markers_v[ m ].index = -1;

        if ((e = recursively_predict_label_in_model_tree(
                        &(labels->elts[ s ]), 
                        &(confidence->elts[ s ]), tree, markers_v)) 
                != NULL)
        {
            return e;
        }
    }

    free(markers_v);

    return NULL;
}


static Error* read_svm_model_node
(
     SVM_model_node* node, 
     const char*     model_dirname
)
{
    uint32_t i, j;
    char     buf[1024] = {0};
    Error*   err;

    for (i = 0; i < node->num_models; i++)
    {
        snprintf(buf, 1024, "%s/%s", model_dirname, node->model_fnames[ i ]);
        if ((err = read_svm_model(&(node->models[ i ]), buf)))
        {
            return err;
        }

        for (j = 0; j < 2; j++)
        {
            if (node->subtrees[ j ][ i ])
            {
                if ((err = read_svm_model_node(node->subtrees[ j ][ i ],
                                model_dirname)))
                {
                    return err;
                }
            }
        }
    }

    return NULL;
}


Error* read_svm_model_tree
(
    SVM_model_tree** tree_out,
    const char*      tree_xml_fname,
    const char*      tree_dtd_fname,
    const char*      model_dirname
)
{
    xmlDoc*  xml_doc;
    xmlNode* xml_node = NULL;
    xmlNode* it;
    Error*   err;
    int      slen = 256;
    char     str[slen];

    assert(tree_xml_fname && model_dirname);

    if ((err = read_svm_xml_doc(&xml_doc, tree_xml_fname, tree_dtd_fname)))
    {
        return err;
    }

    for (it = xmlDocGetRootElement(xml_doc)->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && 
                (XMLStrEqual(it->name, "binary-model") ||
                 XMLStrEqual(it->name, "one-vs-all-model") ||
                 XMLStrEqual(it->name, "one-vs-one-model")))
        {
            xml_node = it;
            break;
        }
    }

    assert(xml_node);

    if ((err = create_svm_model_tree_from_xml_node(tree_out, NULL, 0,
                    xml_node)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    if ((err = read_svm_model_node(*tree_out, model_dirname)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    return NULL;
}


Error* write_svm_model_tree
(
    const SVM_model_tree* tree,
    const char*           model_dirname
)
{
    uint32_t i, j;
    char     buf[1024] = {0};
    Error*   err;
    int      slen = 256;
    char     str[slen];

    for (i = 0; i < tree->num_models; i++)
    {
        snprintf(buf, 1024, "%s/%s", model_dirname, tree->model_fnames[ i ]);
        if ((err = write_svm_model(tree->models[ i ], buf)))
        {
            return err;
        }

        for (j = 0; j < 2; j++)
        {
            if (tree->subtrees[ j ][ i ])
            {
                if ((err = write_svm_model_tree(tree->subtrees[ j ][ i ], 
                                model_dirname)))
                {
                    snprintf(str, slen, "%s: %s", tree->model_fnames[ i ], 
                            err->msg);
                    return JWSC_EARG(str);
                }
            }
        }
    }

    return NULL;
}


static Error* write_svm_model_node_training_data
(
    SVM_model_node*   node,
    const Vector_u32* labels, 
    const Matrix_i32* markers,
    const char*       data_dirname
)
{
    uint32_t          i, j;
    uint32_t          label[2];
    char              buf[1024] = {0};
    const Vector_u32* altlabel[2];
    Vector_u32*       node_labels  = NULL;
    Matrix_i32*       node_markers = NULL;
    Error*            err;
    int               slen = 256;
    char              str[slen];

    for (i = 0; i < node->num_models; i++)
    {
        label[0] = node->labels[0]->elts[ i ];
        label[1] = node->labels[1]->elts[ i ];

        altlabel[0] = node->altlabels[0][ i ];
        altlabel[1] = node->altlabels[1][ i ];

        if ((err = create_model_training_data(&node_labels, &node_markers, 
                        labels, markers, label, altlabel)))
        {
            snprintf(str, slen, "%s: %s", node->model_fnames[i], err->msg);
            return JWSC_EARG(str);
        }

        snprintf(buf, 1024, "%s/%s", data_dirname, node->model_fnames[ i ]);

        if ((err = write_svm_model_training_data(node_labels, node_markers,
                        buf)))
        {
            return err;
        }

        for (j = 0; j < 2; j++)
        {
            if (node->subtrees[ j ][ i ])
            {
                if ((err = write_svm_model_node_training_data(
                                node->subtrees[ j ][ i ], labels, markers,
                                data_dirname)))
                {
                    return err;
                }
            }
        }
    }

    free_vector_u32(node_labels);
    free_matrix_i32(node_markers);

    return NULL;
}


Error* write_svm_model_tree_training_data
(
    const Vector_u32* labels, 
    const Matrix_i32* markers,
    const char*       tree_xml_fname,
    const char*       tree_dtd_fname,
    const char*       data_dirname
)
{
    SVM_model_tree* tree = NULL;
    xmlDoc*  xml_doc;
    xmlNode* xml_node = NULL;
    xmlNode* it;
    Error*   err;
    int      slen = 256;
    char     str[slen];

    assert(tree_xml_fname);

    if ((err = read_svm_xml_doc(&xml_doc, tree_xml_fname, tree_dtd_fname)))
    {
        return err;
    }

    for (it = xmlDocGetRootElement(xml_doc)->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && 
                (XMLStrEqual(it->name, "binary-model") ||
                 XMLStrEqual(it->name, "one-vs-one-model") ||
                 XMLStrEqual(it->name, "one-vs-all-model")))
        {
            xml_node = it;
            break;
        }
    }

    assert(xml_node);

    if ((err = create_svm_model_tree_from_xml_node(&tree, NULL, 0, xml_node)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    if ((err = write_svm_model_node_training_data(tree, labels, markers,
                    data_dirname)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    free_svm_model_tree(tree);

    return NULL;
}


void free_svm_model_tree(SVM_model_tree* tree)
{
    uint32_t i;

    if (!tree)
        return;

    for (i = 0; i < tree->num_models; i++) 
    {
        free_svm_model_tree(tree->subtrees[0][ i ]);
        free_svm_model_tree(tree->subtrees[1][ i ]);
        free_vector_u32(tree->altlabels[0][ i ]);
        free_vector_u32(tree->altlabels[1][ i ]);
        free_svm_model(tree->models[ i ]);
        free(tree->model_fnames[ i ]);
    }

    free(tree->subtrees[0]);
    free(tree->subtrees[1]);
    free_vector_u32(tree->labels[0]);
    free_vector_u32(tree->labels[1]);
    free(tree->altlabels[0]);
    free(tree->altlabels[1]);
    free_vector_d(tree->cost);
    free_vector_d(tree->gamma);
    free(tree->models);
    free(tree->model_fnames);
    free(tree);
}