nb_gauss.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 <jwsc/prob/pdf.h>

#include "xml.h"
#include "haplo_groups.h"
#include "nb_gauss.h"


void train_nb_gauss_model
(
    NB_gauss_model**  model_out,
    const Vector_u32* labels, 
    const Matrix_i32* markers,
    const Vector_d*   label_priors
)
{
    uint32_t  num_labels;
    uint32_t  label;
    uint32_t  num_samples;
    uint32_t  sample;
    uint32_t  num_markers;
    uint32_t  marker;
    double    marker_val;

    Matrix_d*   mu;
    Matrix_d*   sigma;
    Vector_d*   priors;
    Matrix_u32* marker_counts;

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

    num_labels = get_num_haplo_groups();

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

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

    create_zero_matrix_d(&mu, num_labels, num_markers);
    create_init_matrix_d(&sigma, num_labels, num_markers, 0.1);

    marker_counts = NULL;
    create_zero_matrix_u32(&marker_counts, num_labels, 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 ];

            if (marker_val > 0)
            {
                marker_counts->elts[ label ][ marker ]++;
                mu->elts[ label ][ marker ] += marker_val;
                sigma->elts[ label ][ marker ] += marker_val*marker_val;
            }
        }
    }
    for (label = 0; label < num_labels; label++)
    {
        for (marker = 0; marker < num_markers; marker++)
        {
            if (marker_counts->elts[ label ][ marker ] > 0)
            {
                mu->elts[ label ][ marker ] /= 
                    marker_counts->elts[ label ][ marker ];

                sigma->elts[ label ][ marker ] /= 
                    marker_counts->elts[ label ][ marker ];

                sigma->elts[ label ][ marker ] -= 
                    mu->elts[ label ][ marker ] * mu->elts[ label ][ marker ];

                sigma->elts[ label ][ marker ] = 
                    sqrt(sigma->elts[ label ][ marker ]);
            }
        }
    }

    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_gauss_model)));
    }
    (*model_out)->num_markers  = num_markers;
    (*model_out)->mu           = mu;
    (*model_out)->sigma        = sigma;
    (*model_out)->priors       = priors;

    free_matrix_u32(marker_counts);
}


Error* predict_label_with_nb_gauss_model
(
    uint32_t*             label_out,
    double*               confidence_out,
    const Vector_i32*     markers,
    const NB_gauss_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;
    double   marker_val;
    double   mu, sigma;
    double   log_map;
    double   log_ll;
    double   log_prior;
    double   log_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;

    log_map = log(JWSC_MIN_LOG_ARG);
    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++)
    {
        log_ll = 0;

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

            /* Ignore non-positive marker values. */
            if (marker_val > 0)
            {
                mu = model->mu->elts[ label ][ marker ];
                sigma = model->sigma->elts[ label ][ marker ];

                log_ll += log_gaussian_pdf_d(mu, sigma, marker_val);
            }
        }

        if (model->priors->elts[ label ] < JWSC_MIN_LOG_ARG)
        {
            log_prior = log(JWSC_MIN_LOG_ARG);
        }
        else
        {
            log_prior = log(model->priors->elts[ label ]);
        }

        log_posterior = log_ll + log_prior;
        posterior_sum += exp(log_posterior);

        if (log_posterior > log_map)
        {
            log_map = log_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 ] = log_map;
        }
    }

    assert(order < best_labels->num_elts);

    *label_out = best_labels->elts[ order ];

    if (posterior_sum > 0)
    {
        assert(finite(*confidence_out = exp(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_gauss_model
(
    Vector_u32**          labels_out,
    Vector_d**            confidence_out,
    const Matrix_i32*     markers,
    const NB_gauss_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_gauss_model(&(labels->elts[ sample ]), 
                &(confidence->elts[ sample ]), markers_v, model, 0)) != NULL)
        {
            return e;
        }
    }

    free_vector_i32(markers_v);

    return NULL;
}


Error* read_nb_gauss_model(NB_gauss_model** model_out, const char* fname)
{
    FILE*           fp;
    NB_gauss_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_EIO(str);
    }

    if (*model_out == NULL)
    {
        assert((*model_out = malloc(sizeof(NB_gauss_model))) != NULL);
        model = *model_out;
        model->mu = NULL;
        model->sigma = NULL;
        model->priors = NULL;
    }

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

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

    if ((e = read_matrix_with_header_fp_d(&(model->sigma), fp)) != NULL)
    {
        snprintf(str, slen, "%s: %s", fname, "Improperly formatted model file");
        return JWSC_EIO(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_EIO(str);
    }

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

    return NULL;
}


Error* write_nb_gauss_model(NB_gauss_model* model, const char* fname)
{
    FILE* fp;
    int   slen = 256;
    char  str[slen];

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

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

    write_matrix_with_header_fp_d(model->mu, fp);
    write_matrix_with_header_fp_d(model->sigma, fp);
    write_vector_with_header_fp_d(model->priors, fp);

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

    return NULL;
}


void free_nb_gauss_model(NB_gauss_model* model)
{
    if (!model)
        return;

    free_matrix_d(model->mu);
    free_matrix_d(model->sigma);
    free_vector_d(model->priors);
    free(model);
}


static Error* read_nb_gauss_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_EIO(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_EIO(str);
        }

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

        xmlFreeValidCtxt(xml_valid_ctxt);
        xmlFreeDtd(xml_dtd);
    }

    return NULL;
}


static Error* create_nb_gauss_model_tree_from_xml_node
(
     NB_gauss_model_tree**      tree_out, 
     const NB_gauss_model_tree* parent,
     uint32_t                   parent_label,
     xmlNode*                   xml_node
)
{
    NB_gauss_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_gauss_model_tree(*tree_out);
    }

    assert(*tree_out = malloc(sizeof(NB_gauss_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_gauss_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_gauss_model_node
(
    NB_gauss_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;
    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_gauss_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_gauss_model_node(node->subtrees[ i ], labels,
                        markers)))
            {
                return err;
            }
        }
    }

    free_vector_u32(node_labels);
    free_matrix_i32(node_markers);

    return NULL;
}


Error* train_nb_gauss_model_tree
(
    NB_gauss_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_gauss_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_gauss_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_gauss_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_gauss_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_gauss_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_gauss_model_tree
(
    Vector_u32**              labels_out,
    Vector_d**                confidence_out,
    const Matrix_i32*         markers,
    const NB_gauss_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_gauss_model_node
(
     NB_gauss_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_gauss_model(&(node->model), buf)))
    {
        return err;
    }

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

    return NULL;
}


Error* read_nb_gauss_model_tree
(
    NB_gauss_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_gauss_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_gauss_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_gauss_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_gauss_model_tree
(
    const NB_gauss_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_gauss_model(tree->model, buf)))
    {
        return err;
    }

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

    return NULL;
}


void free_nb_gauss_model_tree(NB_gauss_model_tree* tree)
{
    uint32_t i;

    if (!tree)
        return;

    for (i = 0; i < tree->num_groups; i++)
    {
        free_nb_gauss_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_gauss_model(tree->model);
    free(tree->model_fname);
    free(tree);
}