nb_freq.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>

#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 "xml.h"
#include "haplo_groups.h"
#include "nb_freq.h"


void train_nb_freq_model
(
    NB_freq_model**   model_out,
    const Vector_u32* labels, 
    const Matrix_i32* markers,
    const Vector_d*   label_priors
)
{
    uint32_t  num_labels;
    uint32_t  num_samples;
    uint32_t  sample;
    uint32_t  marker;
    uint32_t  label;
    uint32_t  i;
    uint32_t  num_markers;
    uint32_t  num_marker_vals;
    uint32_t  min_marker_val;
    uint32_t  max_marker_val;
    int32_t   marker_val;
    double    sum;

    Matrix_d* marker_given_label;
    Vector_d* priors;

    assert(markers->num_rows == labels->num_elts);

    num_labels = get_num_haplo_groups();

    if (*model_out != NULL)
    {
        marker_given_label = (*model_out)->marker_given_label;
        priors = (*model_out)->priors;
    }
    else
    {
        marker_given_label = NULL;
        priors = NULL;
    }

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

    max_marker_val = 1;
    min_marker_val = 1;

    for (sample = 0; sample < num_samples; sample++)
    {
        for (marker = 0; marker < num_markers; marker++)
        {
            if (markers->elts[ sample ][ marker ] > (int32_t)min_marker_val)
            {
                if (markers->elts[ sample ][ marker ] > (int32_t)max_marker_val)
                {
                    max_marker_val = markers->elts[ sample ][ marker ];
                }
            }
        }
    }

    num_marker_vals = max_marker_val - min_marker_val + 1;

    create_zero_matrix_d(&marker_given_label, num_labels, 
            num_marker_vals*num_markers);

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

        for (marker = 0; marker < num_markers; marker++)
        {
            marker_val = markers->elts[ sample ][ marker ];

            /* Ignore non-positive marker values. */
            if (marker_val > min_marker_val)
            {
                i = num_marker_vals*marker + 
                    ((uint32_t)marker_val - min_marker_val);

                marker_given_label->elts[ label ][ i ]++;
            }
        }
    }

    /* If a label doesn't have any representation in the data, leave the
     * entire row of marker_given_label zero for sparse reading and writing.
     * Otherwise, if a label does appear, set all it's "zero" counts to a
     * small value so that there is no zero probability, just a small value. */
    for (label = 0; label < num_labels; label++)
    {
        sum = 0;

        for (i = 0; i < marker_given_label->num_cols; i++)
        {
            sum += marker_given_label->elts[ label ][ i ];
        }

        if (sum)
        {
            for (i = 0; i < marker_given_label->num_cols; i++)
            {
                if (marker_given_label->elts[ label ][ i ] == 0)
                {
                    marker_given_label->elts[ label ][ i ] = 0.001;
                }
            }
        }
    }

    for (label = 0; label < num_labels; label++)
    {
        for (marker = 0; marker < num_markers; marker++)
        {
            sum = 0;

            for (marker_val = min_marker_val; marker_val < num_marker_vals; 
                    marker_val++)
            {
                i = num_marker_vals*marker + 
                    ((uint32_t)marker_val - min_marker_val); 

                sum += marker_given_label->elts[ label ][ i ];
            }

            if (sum == 0)
                break;

            sum = 1.0 / sum;

            for (marker_val = min_marker_val; marker_val < num_marker_vals; 
                    marker_val++)
            {
                i = num_marker_vals*marker + 
                    ((uint32_t)marker_val - min_marker_val); 

                marker_given_label->elts[ label ][ i ] *= sum;
            }
        }
    }

    if (label_priors == NULL)
    {
        create_zero_vector_d(&priors, num_labels);

        for (sample = 0; sample < num_samples; sample++)
        {
            label = labels->elts[ sample ];
            priors->elts[ label ]++;
        }
        normalize_vector_sum_d(&priors, priors);
    }
    else
    {
        normalize_vector_sum_d(&priors, label_priors);
    }


    if (*model_out == NULL)
    {
        assert(*model_out = malloc(sizeof(NB_freq_model)));
    }
    (*model_out)->num_markers        = num_markers;
    (*model_out)->min_marker_val     = min_marker_val;
    (*model_out)->num_marker_vals    = num_marker_vals;
    (*model_out)->marker_given_label = marker_given_label;
    (*model_out)->priors             = priors;
}


Error* predict_label_with_nb_freq_model
(
    uint32_t*            label_out,
    double*              confidence_out,
    const Vector_i32*    markers,
    const NB_freq_model* model,
    uint32_t             order
)
{
    uint32_t num_markers;
    uint32_t marker;
    uint32_t num_labels;
    uint32_t label;
    uint32_t best_label;
    uint32_t i;
    int32_t  marker_val;

    double  map;
    double  likelihood;
    double  prior;
    double  posterior;
    double  posterior_sum;

    Vector_u32* best_labels = NULL;
    Vector_d* best_labels_conf = NULL;

    num_labels = get_num_haplo_groups();
    num_markers = markers->num_elts;

    map = 0;
    posterior_sum = 0;

    create_init_vector_u32(&best_labels, num_labels, 0);
    create_init_vector_d(&best_labels_conf, num_labels, 0);

    for (label = 0; label < num_labels; label++)
    {
        likelihood = 1.0;

        for (marker = 0; marker < num_markers; marker++)
        {
            marker_val = markers->elts[ marker ];

            /* Ignore non-positive or invalid marker values. */
            if (marker_val > 0 && marker_val <= 
                    (model->num_marker_vals + model->min_marker_val))
            {
                i = model->num_marker_vals*marker + 
                    ((uint32_t)marker_val - model->min_marker_val);

                likelihood *= model->marker_given_label->elts[ label ][ i ];
            }
        }

        prior = model->priors->elts[ label ];

        assert(finite(posterior = likelihood * prior));
        posterior_sum += posterior;

        if (posterior > map)
        {
            map = posterior;
            best_label = label;
            for (i = num_labels-1; i > 0; i--)
            {
                best_labels->elts[ i ] = best_labels->elts[ i - 1];
                best_labels_conf->elts[ i ] = best_labels_conf->elts[ i - 1];
            }
            best_labels->elts[ 0 ] = best_label;
            best_labels_conf->elts[ 0 ] = map;
        }
    }

    assert(order < best_labels->num_elts);

    *label_out = best_labels->elts[ order ];

    if (posterior_sum > 0)
    {
        assert(finite(*confidence_out = best_labels_conf->elts[order] / 
            (double)posterior_sum));
    }
    else
    {
        *confidence_out = 0;
    }

    free_vector_u32(best_labels);
    free_vector_d(best_labels_conf);

    return NULL;
}


Error* predict_labels_with_nb_freq_model
(
    Vector_u32**         labels_out,
    Vector_d**           confidence_out,
    const Matrix_i32*    markers,
    const NB_freq_model* model
)
{
    uint32_t num_samples;
    uint32_t sample;
    uint32_t num_markers;
    uint32_t marker;

    Vector_u32* labels;
    Vector_d*   confidence;
    Vector_i32* 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_vector_d(confidence_out, num_samples);
    confidence = *confidence_out;

    markers_v = NULL;
    create_vector_i32(&markers_v, num_markers);

    for (sample = 0; sample < num_samples; sample++)
    {
        for (marker = 0; marker < num_markers; marker++)
        {
            markers_v->elts[ marker ] = markers->elts[ sample ][marker];
        }

        if ((e = predict_label_with_nb_freq_model(&(labels->elts[ sample ]), 
                &(confidence->elts[ sample ]), markers_v, model, 0)) != NULL)
        {
            return e;
        }
    }

    free_vector_i32(markers_v);

    return NULL;
}


Error* read_nb_freq_model(NB_freq_model** model_out, const char* fname)
{
    FILE*    fp;
    uint32_t label;
    uint32_t num_labels;
    uint32_t num_cols;
    uint32_t s;

    Vector_u32*    sparse = NULL;
    Matrix_d*      sparse_marker_given_label = NULL;
    NB_freq_model* model;
    Error*         e;
    int            slen = 256;
    char           str[slen];

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

    if (*model_out == NULL)
    {
        assert((*model_out = malloc(sizeof(NB_freq_model))) != NULL);
        model = *model_out;
        model->marker_given_label = NULL;
        model->priors = NULL;
    }

    if (fscanf(fp, "%u %u %u\n", &(model->num_markers), 
                    &(model->min_marker_val), &(model->num_marker_vals)) != 3)
    {
        snprintf(str, slen, "%s: %s", fname, "Improperly formatted model file");
        return JWSC_EARG(str);
    }

    if ((e = read_vector_with_header_fp_u32(&(sparse), fp)) != NULL)
    {
        snprintf(str, slen, "%s: %s", fname, "Improperly formatted model file");
        return JWSC_EARG(str);
    }

    if ((e = read_matrix_with_header_fp_d(&(sparse_marker_given_label), fp))
            != NULL)
    {
        snprintf(str, slen, "%s: %s", fname, "Improperly formatted model file");
        return JWSC_EARG(str);
    }

    if ((e = read_vector_with_header_fp_d(&(model->priors), fp)) != NULL)
    {
        snprintf(str, slen, "%s: %s", fname, "Improperly formatted model file");
        return JWSC_EARG(str);
    }

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

    num_labels = sparse->num_elts;
    num_cols   = sparse_marker_given_label->num_cols;

    create_zero_matrix_d(&(model->marker_given_label), num_labels, num_cols);

    s = 0;
    for (label = 0; label < num_labels; label++)
    {
        if (sparse->elts[ label ])
        {
            copy_matrix_block_into_matrix_d(model->marker_given_label,
                    label, 0, sparse_marker_given_label, s++, 0, 1, num_cols);
        }
    }

    free_vector_u32(sparse);
    free_matrix_d(sparse_marker_given_label);

    return NULL;
}


Error* write_nb_freq_model(NB_freq_model* model, const char* fname)
{
    FILE* fp;
    uint32_t label;
    uint32_t num_labels;
    uint32_t i;
    uint32_t num_cols;
    uint32_t s;
    int  slen = 256;
    char str[slen];

    Vector_u32* sparse = NULL;
    Matrix_d*   sparse_marker_given_label = NULL;

    num_labels = model->marker_given_label->num_rows;
    num_cols   = model->marker_given_label->num_cols;

    create_zero_vector_u32(&sparse, num_labels);

    s = 0;
    for (label = 0; label < num_labels; label++)
    {
        for (i = 0; i < num_cols; i++)
        {
            if (model->marker_given_label->elts[ label ][ i ] > 0)
            {
                sparse->elts[ label ] = 1;
                s++;
                break;
            }
        }
    }

    create_matrix_d(&sparse_marker_given_label, s, num_cols);

    s = 0;
    for (label = 0; label < num_labels; label++)
    {
        if (sparse->elts[ label ])
        {
            copy_matrix_block_into_matrix_d(sparse_marker_given_label,
                    s++, 0, model->marker_given_label, label, 0, 1, num_cols);
        }
    }

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

    fprintf(fp, "%u %u %u\n", model->num_markers, model->min_marker_val,
            model->num_marker_vals);

    write_vector_with_header_fp_u32(sparse, fp);
    write_matrix_with_header_fp_d(sparse_marker_given_label, fp);
    write_vector_with_header_fp_d(model->priors, fp);

    free_vector_u32(sparse);
    free_matrix_d(sparse_marker_given_label);

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

    return NULL;
}


void free_nb_freq_model(NB_freq_model* model)
{
    if (!model)
        return;

    free_matrix_d(model->marker_given_label);
    free_vector_d(model->priors);
    free(model);
}


static Error* read_nb_freq_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_nb_freq_model_tree_from_xml_node
(
     NB_freq_model_tree**      tree_out, 
     const NB_freq_model_tree* parent,
     uint32_t                  parent_label,
     xmlNode*                  xml_node
)
{
    NB_freq_model_tree* tree;
    uint32_t i, j;
    uint32_t label;
    uint32_t altlabel;
    uint32_t num_altlabels;
    double   p;
    xmlAttr* attr;
    xmlNode* it;
    xmlNode* itt;
    Error*   err;
    const char* fname;

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

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

    tree->parent       = parent;
    tree->parent_label = parent_label;
    tree->num_groups   = 0;
    tree->subtrees     = NULL;
    tree->labels       = NULL;
    tree->altlabels    = NULL;
    tree->priors       = NULL;
    tree->model        = NULL;
    tree->model_fname  = NULL;

    assert(XMLStrEqual(xml_node->name, "model"));

    for (it = xml_node->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "group"))
        {
            tree->num_groups++;
        }
    }

    assert(tree->subtrees = calloc(tree->num_groups, sizeof(void*)));
    assert(tree->altlabels = calloc(tree->num_groups, sizeof(void*)));
    create_vector_u32(&(tree->labels), tree->num_groups);

    for (attr = xml_node->properties; attr; attr = attr->next)
    {
        if (XMLStrEqual(attr->name, "priors") && 
                XMLStrEqual(attr->children->content, "true"))
        {
            create_init_vector_d(&(tree->priors), get_num_haplo_groups(), 1.0);
            break;
        }
    }

    i = 0;
    for (it = xml_node->children; it; it = it->next)
    {
        j = 0;
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "file"))
        {
            fname = (const char*) it->children->content;
            tree->model_fname = malloc((strlen(fname)+1) * sizeof(char));
            strcpy(tree->model_fname, fname);
        }
        else if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "group"))
        {
            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(&(tree->altlabels[ 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)))
                    {
                        return err;
                    }
                    tree->labels->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)))
                    {
                        return err;
                    }
                    tree->altlabels[ i ]->elts[ j++ ] = altlabel;
                }
                else if (tree->priors && itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "prior"))
                {
                    if (sscanf((char*)itt->children->content, "%lf", &p) != 1 ||
                            p < 0 || p > 1)
                    {
                        return JWSC_EARG("Invalid prior");
                    }
                    tree->priors->elts[ i ] = p;
                }
                else if (itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "model"))
                {
                    if ((err = create_nb_freq_model_tree_from_xml_node(
                                 &(tree->subtrees[ i ]), tree, label, itt)))
                    {
                        return err;
                    }
                }
            }
            i++;
        }
    }
    assert(i == tree->num_groups);

    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,
    const Vector_u32* model_labels,
    Vector_u32*const* model_altlabels
)
{
    uint8_t     b;
    uint32_t    i, j, k;
    uint32_t    n;
    Vector_u32* train_labels = NULL;
    Matrix_i32* train_markers = NULL;

    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 < model_labels->num_elts; j++)
        {
            if (is_ancestor(data_labels->elts[ i ], model_labels->elts[ j ]))
            {
                train_labels->elts[ n ] = model_labels->elts[ j ];
                copy_matrix_block_into_matrix_i32(train_markers, n, 0,
                        data_markers, i, 0, 1, data_markers->num_cols);
                n++;
                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->elts[ j ];
                        copy_matrix_block_into_matrix_i32(train_markers, n, 0,
                                data_markers, i, 0, 1, data_markers->num_cols);
                        n++;
                        break;
                    }
                }
                if (b)
                {
                    break;
                }
            }
        }
    }

    if (!n)
    {
        return JWSC_EARG("No data for model");
    }

    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_nb_freq_model_node
(
    NB_freq_model_node* node,
    const Vector_u32*   labels, 
    const Matrix_i32*   markers
)
{
    uint32_t i;
    Vector_u32* node_labels  = NULL;
    Matrix_i32* node_markers = NULL;
    Error* err = NULL;
    int  slen = 256;
    char str[slen];

    if ((err = create_model_training_data(&node_labels, &node_markers, labels,
                    markers, node->labels, node->altlabels)))
    {
        snprintf(str, slen, "%s: %s", node->model_fname, err->msg);
        return JWSC_EARG(str);
    }

    train_nb_freq_model(&(node->model), node_labels, node_markers, 
            node->priors);

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

    free_vector_u32(node_labels);
    free_matrix_i32(node_markers);

    return NULL;
}


Error* train_nb_freq_model_tree
(
    NB_freq_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;
    xmlNode* it;
    Error*   err;
    int      slen = 256;
    char     str[slen];

    assert(tree_xml_fname);

    if ((err = read_nb_freq_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, "model"))
        {
            xml_node = it;
            break;
        }
    }

    if ((err = create_nb_freq_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_nb_freq_model_node(*tree_out, labels, markers)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    return NULL;
}


static Error* recursively_predict_label_in_model_tree
(
    uint32_t*                 label_out,
    double*                   confidence_out,
    const NB_freq_model_tree* tree,
    const Vector_i32*         markers_v,
    uint32_t                  order
)
{
    uint32_t subtree_label;
    uint32_t i;
    double   confidence;
    Error*   err;

    if ((err = predict_label_with_nb_freq_model(label_out, &confidence, 
                    markers_v, tree->model, order)) != NULL)
    {
        return err;
    }
    *confidence_out *= confidence;

    for (i = 0; i < tree->num_groups; i++)
    {
        if (tree->subtrees[ i ])
        {
            subtree_label = tree->subtrees[ i ]->parent_label;

            if (subtree_label == *label_out)
            {
                if ((err = recursively_predict_label_in_model_tree(label_out,
                                confidence_out, tree->subtrees[ i ], markers_v,
                                order)) != NULL)
                {
                    return err;
                }
            }
        }
    }

    return NULL;
}


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

    Vector_u32* labels;
    Vector_d*   confidence;
    Vector_i32* markers_v;
    Error*      err;

    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;
    create_vector_i32(&markers_v, num_markers);

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

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

    free_vector_i32(markers_v);

    return NULL;
}


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

    snprintf(buf, 1024, "%s/%s", model_dirname, node->model_fname);
    if ((err = read_nb_freq_model(&(node->model), buf)))
    {
        return err;
    }

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

    return NULL;
}


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

    assert(tree_xml_fname && model_dirname);

    if ((err = read_nb_freq_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, "model"))
        {
            xml_node = it;
            break;
        }
    }

    if ((err = create_nb_freq_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_nb_freq_model_node(*tree_out, model_dirname)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    return NULL;
}


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

    snprintf(buf, 1024, "%s/%s", model_dirname, tree->model_fname);
    if ((err = write_nb_freq_model(tree->model, buf)))
    {
        return err;
    }

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

    return NULL;
}


void free_nb_freq_model_tree(NB_freq_model_tree* tree)
{
    uint32_t i;

    if (!tree)
        return;

    for (i = 0; i < tree->num_groups; i++)
    {
        free_nb_freq_model_tree(tree->subtrees[ i ]);
        free_vector_u32(tree->altlabels[ i ]);
    }

    free(tree->subtrees);
    free(tree->altlabels);
    free_vector_u32(tree->labels);
    free_vector_d(tree->priors);
    free_nb_freq_model(tree->model);
    free(tree->model_fname);
    free(tree);
}