mv_gmm.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/matblock/matblock.h>
#include <jwsc/matblock/matblock_io.h>
#include <jwsc/stat/gmm.h>

#include "xml.h"
#include "haplo_groups.h"
#include "input.h"
#include "mv_gmm.h"


void train_mv_gmm_model
(
    MV_gmm_model**    model_out,
    const Vector_u32* labels, 
    const Matrix_i32* markers,
    const Vector_d*   label_priors,
    uint32_t          num_components
)
{
    uint32_t  l, num_labels;
    uint32_t  s, num_samples;
    uint32_t  m, num_markers;
    uint32_t  n, N;

    Matrix_d**   mu;
    Matblock_d** sigma;
    Vector_d**   pi;
    Vector_d*    priors;
    Matrix_d*    markers_d = NULL;

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

    num_markers = markers->num_cols;
    num_samples = markers->num_rows;
    num_labels  = get_num_haplo_groups();

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

    assert(mu    = malloc(num_labels*sizeof(Matrix_d*)));
    assert(sigma = malloc(num_labels*sizeof(Matblock_d*)));
    assert(pi    = malloc(num_labels*sizeof(Vector_d*)));
    priors = NULL;

    for (l = 0; l < num_labels; l++)
    {
        mu   [ l ] = NULL;
        sigma[ l ] = NULL;
        pi   [ l ] = NULL;

        N = 0;
        for (s = 0; s < num_samples; s++)
        {
            if (labels->elts[ s ] == l)
            {
                N++;
            }
        }

        // TODO we should probably handle the case of N=1 by creating a
        // distribution with the sample as the mean and a very small diagonal
        // covariance matrix.
        if (N < num_components*2)
            continue;

        create_matrix_d(&markers_d, N, num_markers);

        n = 0;
        for (s = 0; s < num_samples; s++)
        {
            if (labels->elts[ s ] == l)
            {
                for (m = 0; m < num_markers; m++)
                {
                    markers_d->elts[ n ][ m ] = markers->elts[ s ][ m ];
                }
                n++;
            }
        }

        train_gmm_d(&(mu[l]), &(sigma[l]), &(pi[l]), NULL, markers_d,
                num_components);
    }

    free_matrix_d(markers_d);

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

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

    assert(*model_out = malloc(sizeof(MV_gmm_model)));
    (*model_out)->num_markers    = num_markers;
    (*model_out)->num_components = num_components;
    (*model_out)->mu             = mu;
    (*model_out)->sigma          = sigma;
    (*model_out)->pi             = pi;
    (*model_out)->priors         = priors;
}


Error* predict_label_with_mv_gmm_model
(
    uint32_t*           label_out,
    double*             confidence_out,
    const Vector_i32*   markers,
    const MV_gmm_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   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;
    Matrix_d*  marker_vals      = NULL;

    const Matrix_d*   mu;
    const Matblock_d* sigma;
    const Vector_d*   pi;

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

    log_map = log(0);
    posterior_sum = 0;

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

    create_matrix_d(&marker_vals, 1, num_markers);

    for (label = 0; label < num_labels; label++)
    {
        if (!model->mu[ label ])
            continue;

        for (marker = 0; marker < num_markers; marker++)
        {
            if (markers->elts[ marker ] <= 0)
            {
                // Impute missing value with the mean.
                marker_vals->elts[0][ marker ] = 0;
                for (i = 0; i < model->num_components; i++)
                {
                    marker_vals->elts[0][ marker ] += 
                        model->mu[ label ]->elts[ i ][ marker ] * 
                        model->pi[ label ]->elts[ i ];
                }
            }
            else
            {
                marker_vals->elts[0][ marker ] = markers->elts[ marker ];
            }
        }

        mu = model->mu[ label ];
        sigma = model->sigma[ label ];
        pi = model->pi[ label ];

        log_ll = gmm_log_likelihood_d(mu, sigma, pi, marker_vals);

        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);
    free_matrix_d(marker_vals);

    return NULL;
}


Error* predict_labels_with_mv_gmm_model
(
    Vector_u32**        labels_out,
    Vector_d**          confidence_out,
    const Matrix_i32*   markers,
    const MV_gmm_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_mv_gmm_model(&(labels->elts[ sample ]), 
                &(confidence->elts[ sample ]), markers_v, model, 0)) != NULL)
        {
            return e;
        }
    }

    free_vector_i32(markers_v);

    return NULL;
}


Error* read_mv_gmm_model(MV_gmm_model** model_out, const char* fname)
{
    uint32_t      num_labels;
    uint32_t      i;
    FILE*         fp;
    Vector_u32*   v = NULL;
    MV_gmm_model* model;
    Error*        e;
    int           slen = 256;
    char          str[slen];

    num_labels = get_num_haplo_groups();

    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(MV_gmm_model))) != NULL);
        model = *model_out;
        model->priors = NULL;
    }

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

    assert(model->mu    = malloc(num_labels*sizeof(void*)));
    assert(model->sigma = malloc(num_labels*sizeof(void*)));
    assert(model->pi    = malloc(num_labels*sizeof(void*)));

    if ((e = read_vector_fp_u32(&v, fp)) || (v->num_elts != num_labels))
    {
        return JWSC_EARG("Improperly formatted model file");
    }

    for (i = 0; i < num_labels; i++)
    {
        model->mu[ i ]    = NULL;
        model->sigma[ i ] = NULL;
        model->pi[ i ]    = NULL;

        if (v->elts[ i ])
        {
            if ((e = read_matrix_with_header_fp_d(&(model->mu[i]), fp)) ||
                (e = read_matblock_with_header_fp_d(&(model->sigma[i]), fp)) ||
                (e = read_vector_with_header_fp_d(&(model->pi[i]), fp)))
            {
                return JWSC_EARG("Improperly formatted model file");
            }
        }
    }

    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);
    }

    return NULL;
}


Error* write_mv_gmm_model(MV_gmm_model* model, const char* fname)
{
    uint32_t num_labels;
    uint32_t i;
    FILE*    fp;
    int      slen = 256;
    char     str[slen];

    num_labels = get_num_haplo_groups();

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

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

    for (i = 0; i < num_labels; i++)
    {
        fprintf(fp, "%u ", (model->mu[ i ] != NULL));
    }
    fprintf(fp, "\n");

    for (i = 0; i < num_labels; i++)
    {
        if (model->mu[ i ])
        {
            write_matrix_with_header_fp_d(model->mu[ i ], fp);
            write_matblock_with_header_fp_d(model->sigma[ i ], fp);
            write_vector_with_header_fp_d(model->pi[ i ], fp);
        }
    }

    write_vector_with_header_fp_d(model->priors, fp);

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

    return NULL;
}


void free_mv_gmm_model(MV_gmm_model* model)
{
    uint32_t i;
    uint32_t num_labels;

    if (!model)
        return;

    num_labels = get_num_haplo_groups();

    for (i = 0; i < num_labels; i++)
    {
        free_matrix_d(model->mu[ i ]);
        free_matblock_d(model->sigma[ i ]);
        free_vector_d(model->pi[ i ]);
    }
    free(model->mu);
    free(model->sigma);
    free(model->pi);
    free_vector_d(model->priors);
    free(model);
}


static Error* read_mv_gmm_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_mv_gmm_model_tree_from_xml_node
(
     MV_gmm_model_tree**      tree_out, 
     const MV_gmm_model_tree* parent,
     uint32_t                 parent_label,
     xmlNode*                 xml_node
)
{
    MV_gmm_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_mv_gmm_model_tree(*tree_out);
    }

    assert(*tree_out = malloc(sizeof(MV_gmm_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, "mixture-components"))
        {
            if (sscanf((char*)it->children->content, "%d", 
                        &(tree->num_components)) != 1 || 
                    tree->num_components < 1)
            {
                return JWSC_EARG("Invalid number of mixture components");
            }
        }
        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_mv_gmm_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_mv_gmm_model_node
(
    MV_gmm_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_mv_gmm_model(&(node->model), node_labels, node_markers, 
            node->priors, node->num_components);

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

    free_vector_u32(node_labels);
    free_matrix_i32(node_markers);

    return NULL;
}


Error* train_mv_gmm_model_tree
(
    MV_gmm_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_mv_gmm_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_mv_gmm_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_mv_gmm_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 MV_gmm_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_mv_gmm_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_mv_gmm_model_tree
(
    Vector_u32**             labels_out,
    Vector_d**               confidence_out,
    const Matrix_i32*        markers,
    const MV_gmm_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_mv_gmm_model_node
(
     MV_gmm_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_mv_gmm_model(&(node->model), buf)))
    {
        return err;
    }

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

    return NULL;
}


Error* read_mv_gmm_model_tree
(
    MV_gmm_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_mv_gmm_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_mv_gmm_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_mv_gmm_model_node(*tree_out, model_dirname)))
    {
        snprintf(str, slen, "%s: %s", tree_xml_fname, err->msg);
        return JWSC_EARG(str);
    }

    return NULL;
}


Error* write_mv_gmm_model_tree
(
    const MV_gmm_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_mv_gmm_model(tree->model, buf)))
    {
        return err;
    }

    for (i = 0; i < tree->num_groups; i++)
    {
        if (tree->subtrees[ i ])
        {
            if ((err = write_mv_gmm_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_mv_gmm_model_tree(MV_gmm_model_tree* tree)
{
    uint32_t i;

    if (!tree)
        return;

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