input.c

Go to the documentation of this file.

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


#include <config.h>

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <inttypes.h>
#include <assert.h>
#include <errno.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/vector/vector.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matblock/matblock.h>

#include <xml.h>
#include <haplo_groups.h>
#include <options.h>
#include <input.h>


#define  MAX_COL_SIZE  256


static uint32_t read_block_into_buffer(char* buf, uint32_t size, FILE* fp)
{
    uint32_t n;
    char     c;

    assert(size > 0);
    n = 0;

    while (n < size && (c = fgetc(fp)) != EOF)
    {
        *buf = c;
        buf++;
        n++;
    }

    return n;
}


static Error* read_file_into_buffer
(
    char**      buf_out, 
    uint32_t*   buf_len_out, 
    const char* fname
)
{
    FILE*    fp;
    char*    buf;
    char*    buf_p;
    uint32_t buf_len;
    uint32_t block_size;
    uint32_t n;
    int      slen = 256;
    char     str[slen];

    assert(*buf_out == NULL);

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

    block_size = 1024;
    buf_len = block_size;
    assert(buf = malloc(buf_len*sizeof(char)));
    buf_p = buf;

    while ((n = read_block_into_buffer(buf_p, block_size, fp)) == block_size)
    {
        buf_len += block_size;
        assert(buf = realloc(buf, buf_len*sizeof(char)));
        buf_p = buf + (buf_len - block_size);
    }
    buf_len -= block_size - n;

    if (buf_len == 0)
    {
        snprintf(str, slen, "%s: %s", fname, "No data to read");
        return JWSC_EIO(str);
    }

    assert(buf = realloc(buf, buf_len*sizeof(char)));

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

    *buf_out = buf;
    *buf_len_out = buf_len;

    return NULL;
}


static uint8_t skip_to_csv_column
(
    const char** buf_in_out, 
    uint32_t     buf_len, 
    uint32_t     col
)
{
    const char* buf_in;
    const char* buf_p;
    uint32_t cur_col;

    assert(buf_len > 0);

    buf_in = *buf_in_out;
    buf_p  = *buf_in_out;
    cur_col = 1;

    while ((cur_col != col) && ((buf_p - buf_in) < buf_len-1) && 
            (*buf_p != '\n'))
    {
        if (*buf_p == ',')
        {
            cur_col++;
        }
        buf_p++;
    }

    *buf_in_out = buf_p;

    if (cur_col != col)
    {
        return 0;
    }

    return 1;
}


static const char* read_csv_column
(
    const char* src, 
    uint32_t    src_len, 
    char*       dst, 
    uint32_t    dst_len
)
{
    const char* src_p;
    char* dst_p;

    assert(src_len > 0 && dst_len > 0);

    src_p = src;
    dst_p = dst;

    while (((src_p - src) < src_len-1) && (*src_p != ','))
    {
        if ((dst_p - dst) == dst_len-1)
        {
            break;
        }

        *dst_p = *src_p;
        dst_p++;
        src_p++;
    }

    *dst_p = '\0';

    if (*src_p == ',')
    {
        src_p++;
    }

    return src_p;
}


static uint8_t skip_to_txt_column
(
    const char** buf_in_out, 
    uint32_t     buf_len, 
    uint32_t     col
)
{
    const char* buf_in;
    const char* buf_p;
    uint32_t cur_col;

    assert(buf_len > 0);

    buf_in = *buf_in_out;
    buf_p  = *buf_in_out;
    cur_col = 1;

    while ((cur_col != col) && ((buf_p - buf_in) < buf_len-1) && 
            (*buf_p != '\n'))
    {
        if (*buf_p == ' ')
        {
            cur_col++;
            while (*buf_p == ' ') 
                buf_p++;
        }
        else
        {
            buf_p++;
        }
    }

    *buf_in_out = buf_p;

    if (cur_col != col)
    {
        return 0;
    }

    return 1;
}


static const char* read_txt_column
(
    const char* src, 
    uint32_t    src_len, 
    char*       dst, 
    uint32_t    dst_len
)
{
    const char* src_p;
    char* dst_p;

    assert(src_len > 0 && dst_len > 0);

    src_p = src;
    dst_p = dst;

    while (((src_p - src) < src_len-1) && (*src_p != ' '))
    {
        if ((dst_p - dst) == dst_len-1)
        {
            break;
        }

        *dst_p = *src_p;
        dst_p++;
        src_p++;
    }

    *dst_p = '\0';

    while (*src_p == ' ')
    {
        src_p++;
    }

    return src_p;
}


static uint8_t skip_line
(
    const char** buf_in_out, 
    const char*  buf, 
    uint32_t     buf_len
)
{
    const char* buf_p;

    buf_p = *buf_in_out;

    while (((buf_p - buf) < buf_len) && (*buf_p != '\n')) 
    {
        buf_p++;
    }
    buf_p++;

    *buf_in_out = buf_p;

    return ((buf_p - buf) < buf_len);
}


static Error* get_num_samples_from_txt_csv
(
    uint32_t*   num_samples_out, 
    const char* file_buf, 
    uint32_t    file_len
)
{
    const char* file_buf_p;
    uint32_t num_samples;

    file_buf_p = file_buf;

    if (opts.header_in) 
    {
        if (!skip_line(&file_buf_p, file_buf, file_len))
        {
            return JWSC_EIO("No samples to read");
        }
    }

    num_samples = 1;
    while (skip_line(&file_buf_p, file_buf, file_len))
    {
        num_samples++;
    }

    *num_samples_out = num_samples;

    return NULL;
}


static Error* read_ids_from_txt_csv
(
    Matblock_u8** ids_out,
    const char*   file_buf,
    uint32_t      file_len,
    uint32_t      num_samples,
    uint8_t       aux
)
{
    const char* file_buf_p;
    const char* file_buf_p_saved;
    uint32_t sample;
    uint32_t num_ids;
    uint32_t id;
    uint32_t col;
    uint32_t len;
    int      slen = 256;
    char     s[slen];

    if (!aux && !opts.id_cols)
        return NULL;
    else if (aux && !opts.aux_id_cols)
        return NULL;

    file_buf_p = file_buf;

    if (opts.header_in) 
    {
        skip_line(&file_buf_p, file_buf, file_len);
    }

    num_ids = (aux) ? opts.aux_id_cols->num_elts : opts.id_cols->num_elts;
    create_zero_matblock_u8(ids_out, num_samples, num_ids, MAX_COL_SIZE);

    for (sample = 0; sample < num_samples; sample++)
    {
        file_buf_p_saved = file_buf_p;

        for (id = 0; id < num_ids; id++)
        {
            file_buf_p = file_buf_p_saved;

            col = (aux) ? opts.aux_id_cols->elts[ id ] : 
                opts.id_cols->elts[ id ];
            len = file_len - (file_buf_p - file_buf);

            switch (opts.input_format)
            {
                case HAPLO_INPUT_TXT:
                    if (!skip_to_txt_column(&file_buf_p, len, col))
                    {
                        snprintf(s, slen, "ID %d missing in sample %d", 
                                id+1, sample+1);
                        return JWSC_EIO(s);
                    }
                    read_txt_column(file_buf_p, len, (char*)
                            ((*ids_out)->elts[ sample ][ id ]), MAX_COL_SIZE);
                    break;
                case HAPLO_INPUT_CSV:
                    if (!skip_to_csv_column(&file_buf_p, len, col))
                    {
                        snprintf(s, slen, "ID %d missing in sample %d", 
                                id+1, sample+1);
                        return JWSC_EIO(s);
                    }
                    read_csv_column(file_buf_p, len, (char*)
                            ((*ids_out)->elts[ sample ][ id ]), MAX_COL_SIZE);
                    break;
                case HAPLO_INPUT_XML:
                    abort();
            }
        }

        skip_line(&file_buf_p, file_buf, file_len);
    }

    return NULL;
}


static Error* read_labels_from_txt_csv
(
    Vector_u32** labels_out,
    const char*  file_buf,
    uint32_t     file_len,
    uint32_t     num_samples,
    uint8_t      aux
)
{
    const char* file_buf_p;
    uint32_t label_col;
    uint32_t sample;
    uint32_t len;
    char     col_buf[ MAX_COL_SIZE ] = {0};
    int      slen = 256;
    char     s[slen];

    if (!aux && opts.label_col == 0)
        return NULL;
    else if (aux && opts.aux_label_col == 0)
        return NULL;

    label_col = (aux) ? opts.aux_label_col : opts.label_col;

    file_buf_p = file_buf;

    if (opts.header_in) 
    {
        skip_line(&file_buf_p, file_buf, file_len);
    }

    create_zero_vector_u32(labels_out, num_samples);

    for (sample = 0; sample < num_samples; sample++)
    {
        len = file_len - (file_buf_p - file_buf);

        switch (opts.input_format)
        {
            case HAPLO_INPUT_TXT:
                if (!skip_to_txt_column(&file_buf_p, len, label_col))
                {
                    snprintf(s, slen, "No label in sample %d", sample+1);
                    return JWSC_EIO(s);
                }
                read_txt_column(file_buf_p, len, col_buf, MAX_COL_SIZE);
                break;
            case HAPLO_INPUT_CSV:
                if (!skip_to_csv_column(&file_buf_p, len, label_col))
                {
                    snprintf(s, slen, "No label in sample %d", sample+1);
                    return JWSC_EIO(s);
                }
                read_csv_column(file_buf_p, len, col_buf, MAX_COL_SIZE);
                break;
            case HAPLO_INPUT_XML:
                abort();
        }

        if (lookup_haplo_group_index_from_label(
                    &((*labels_out)->elts[ sample ]),
                    (const char*) col_buf) != NULL)
        {
            snprintf(s, slen, "Label '%s' in sample %d not recognized", 
                    col_buf, sample+1);
            return JWSC_EIO(s);
        }

        skip_line(&file_buf_p, file_buf, file_len);
    }

    return NULL;
}


static Error* read_markers_from_txt_csv
(
    Matrix_i32** markers_out,
    const char*  file_buf,
    uint32_t     file_len,
    uint32_t     num_samples,
    uint8_t      aux
)
{
    const char* file_buf_p;
    const char* file_buf_p_saved;
    uint32_t num_markers;
    uint32_t first_marker_col;
    uint32_t sample;
    uint32_t marker;
    uint32_t col;
    uint32_t len;
    int32_t  marker_val;
    char     col_buf[ MAX_COL_SIZE ] = {0};
    const Vector_u32* marker_cols;
    int       slen = 256;
    char      s[slen];

    if (aux)
        assert(opts.aux_first_marker_col > 0);
    else
        assert(opts.first_marker_col > 0);

    file_buf_p = file_buf;

    if (opts.header_in) 
    {
        skip_line(&file_buf_p, file_buf, file_len);
    }

    num_markers = (aux) ? opts.aux_num_markers : opts.num_markers;
    marker_cols = (aux) ? opts.aux_marker_cols : opts.marker_cols;
    first_marker_col = (aux) ? opts.aux_first_marker_col :
        opts.first_marker_col;

    create_matrix_i32(markers_out, num_samples, num_markers);

    for (sample = 0; sample < num_samples; sample++)
    {
        file_buf_p_saved = file_buf_p;

        for (marker = 0; marker < num_markers; marker++)
        {
            file_buf_p = file_buf_p_saved;

            col = (marker_cols) ? marker_cols->elts[ marker ] 
                                : first_marker_col + marker;

            len = file_len - (file_buf_p - file_buf);

            switch (opts.input_format)
            {
                case HAPLO_INPUT_TXT:
                    if (!skip_to_txt_column(&file_buf_p, len, col))
                    {
                        snprintf(s, slen, "Marker %d missing in sample %d",
                                marker+1, sample+1);
                        return JWSC_EIO(s);
                    }
                    read_txt_column(file_buf_p, len, col_buf, MAX_COL_SIZE);
                    break;
                case HAPLO_INPUT_CSV:
                    if (!skip_to_csv_column(&file_buf_p, len, col))
                    {
                        snprintf(s, slen, "Marker %d missing in sample %d",
                                marker+1, sample+1);
                        return JWSC_EIO(s);
                    }
                    read_csv_column(file_buf_p, len, col_buf, MAX_COL_SIZE);
                    break;
                case HAPLO_INPUT_XML:
                    abort();
            }

            if (sscanf(col_buf, "%d", &marker_val) != 1)
            {
                snprintf(s, slen, "Marker %d in sample %d not valid",
                        marker+1, sample+1);
                return JWSC_EIO(s);
            }
            (*markers_out)->elts[ sample ][ marker ] = marker_val;
        }

        skip_line(&file_buf_p, file_buf, file_len);
    }

    return NULL;
}


static Error* read_input_from_txt_csv
(
    Matblock_u8** ids_out,
    Vector_u32**  labels_out,
    Matrix_i32**  markers_out,
    const char*   file_buf,
    uint32_t      file_len,
    uint8_t       aux
)
{
    uint32_t num_samples;
    Error* err;

    if ((err = get_num_samples_from_txt_csv(&num_samples, file_buf, file_len)))
    {
        return err;
    }

    if ((err = read_ids_from_txt_csv(ids_out, file_buf, file_len, num_samples,
                    aux)))
    {
        return err;
    }
    
    if ((err = read_labels_from_txt_csv(labels_out, file_buf, file_len, 
                    num_samples, aux)))
    {
        return err;
    }

    if ((err = read_markers_from_txt_csv(markers_out, file_buf, file_len, 
                    num_samples, aux)))
    {
        return err;
    }

    return NULL;
}


static uint32_t get_num_samples_from_xml(xmlDoc* doc)
{
    const xmlNode* root;
    const xmlNode* it;
    uint32_t N = 0;

    root = xmlDocGetRootElement(doc);

    for (it = root->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "sample"))
        {
            N++;
        }
    }

    return N;
}


static Error* read_ids_from_xml
(
    Matblock_u8** ids_out, 
    xmlDoc*       doc, 
    uint32_t      num_samples,
    uint8_t       aux
)
{
    uint32_t i, j;
    uint32_t num_ids;
    const xmlNode* root;
    const xmlNode* it;
    const xmlNode* itt;
    int  slen = 256;
    char str[slen];

    if (!aux && !opts.id_cols)
        return NULL;
    else if (aux && !opts.aux_id_cols)
        return NULL;

    num_ids = (aux) ? opts.aux_id_cols->num_elts : opts.id_cols->num_elts;
    create_zero_matblock_u8(ids_out, num_samples, num_ids, MAX_COL_SIZE);

    root = xmlDocGetRootElement(doc);

    i = 0;
    for (it = root->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "sample"))
        {
            j = 0;
            for (itt = it->children; itt; itt = itt->next)
            {
                if (itt->type == XML_ELEMENT_NODE && 
                        XMLStrEqual(itt->name, "id"))
                {
                    if (j == num_ids)
                    {
                        snprintf(str, slen, 
                                "Sample %d has too many IDs", i+1);
                        return JWSC_EARG(str);
                    }
                    strncpy((char*)(*ids_out)->elts[ i ][ j ], 
                            (const char*)itt->children->content, MAX_COL_SIZE);
                    j++;
                }
            }
            if (j < num_ids)
            {
                snprintf(str, slen, "Sample %d missing IDs", i+1);
                return JWSC_EARG(str);
            }
            i++;
        }
    }
    assert(i == num_samples);

    return NULL;
}


static Error* read_labels_from_xml
(
    Vector_u32** labels_out, 
    xmlDoc*      doc, 
    uint32_t     num_samples,
    uint8_t      aux
)
{
    uint32_t i;
    const char*    label;
    const xmlNode* root;
    const xmlNode* it;
    const xmlNode* itt;
    Error* err;
    int    slen = 256;
    char   str[slen];

    if (!aux && opts.label_col == 0)
        return NULL;
    else if (aux && opts.aux_label_col == 0)
        return NULL;

    create_zero_vector_u32(labels_out, num_samples);

    root = xmlDocGetRootElement(doc);

    i = 0;
    for (it = root->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "sample"))
        {
            label = 0;
            for (itt = it->children; itt; itt = itt->next)
            {
                if (itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "label"))
                {
                    label = (const char*) itt->children->content;

                    if ((err = lookup_haplo_group_index_from_label(
                                &((*labels_out)->elts[ i ]), label)))
                    {
                        return err;
                    }
                    break;
                }
            }
            if (!label)
            {
                snprintf(str, slen, "Sample %d missing label", i+1);
                return JWSC_EARG(str);
            }
            i++;
        }
    }
    assert(i == num_samples);

    return NULL;
}


static Error* read_markers_from_xml
(
    Matrix_i32** markers_out, 
    xmlDoc*      doc, 
    uint32_t     num_samples,
    uint8_t      aux
)
{
    uint32_t i, j;
    uint32_t num_markers;
    int32_t  marker_val;
    const xmlNode* root;
    const xmlNode* it;
    const xmlNode* itt;
    int  slen = 256;
    char str[slen];

    if (aux)
        assert(opts.aux_first_marker_col > 0);
    else
        assert(opts.first_marker_col > 0);

    num_markers = (aux) ? opts.aux_num_markers : opts.num_markers;

    create_matrix_i32(markers_out, num_samples, num_markers);

    root = xmlDocGetRootElement(doc);

    i = 0;
    for (it = root->children; it; it = it->next)
    {
        if (it->type == XML_ELEMENT_NODE && XMLStrEqual(it->name, "sample"))
        {
            j = 0;
            for (itt = it->children; itt; itt = itt->next)
            {
                if (itt->type == XML_ELEMENT_NODE &&
                        XMLStrEqual(itt->name, "marker"))
                {
                    if (j == num_markers)
                    {
                        snprintf(str, slen, 
                                "Sample %d has too many markers", i+1);
                        return JWSC_EARG(str);
                    }
                    if (sscanf((const char*)itt->children->content, "%d",  
                            &marker_val) != 1)
                    {
                        snprintf(str, slen, 
                                "Marker %d in sample %d not valid", j+1, i+1);
                        return JWSC_EARG(str);
                    }
                    (*markers_out)->elts[ i ][ j ] = marker_val;
                    j++;
                }
            }
            if (j < num_markers)
            {
                snprintf(str, slen, "Sample %d missing markers", i+1);
                return JWSC_EARG(str);
            }
            i++;
        }
    }
    assert(i == num_samples);

    return NULL;
}


static Error* read_input_from_xml
(
    Matblock_u8** ids_out,
    Vector_u32**  labels_out,
    Matrix_i32**  markers_out,
    const char*   file_buf,
    uint32_t      file_len,
    const char*   fname,
    uint8_t       aux
)
{
    uint32_t       num_samples;
    xmlParserCtxt* xml_parse_ctxt;
    xmlValidCtxt*  xml_valid_ctxt;
    xmlDoc*        xml_doc;
    xmlDtd*        xml_dtd;
    Error*         err;
    int            slen = 256;
    char           str[slen];

    assert(xml_parse_ctxt = xmlNewParserCtxt());

    if (!(xml_doc = xmlCtxtReadMemory(xml_parse_ctxt, file_buf, file_len, fname,
                    NULL, 0)))
    {
        return JWSC_EARG("Could not parse file");
    } 

    xmlFreeParserCtxt(xml_parse_ctxt);

    if (opts.input_dtd_fname)
    {
        assert(xml_valid_ctxt = xmlNewValidCtxt());

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

        if (!xmlValidateDtd(xml_valid_ctxt, xml_doc, xml_dtd))
        {
            return JWSC_EARG("XML file not valid");
        }

        xmlFreeValidCtxt(xml_valid_ctxt);
        xmlFreeDtd(xml_dtd);
    }

    num_samples = get_num_samples_from_xml(xml_doc);

    if ((err = read_ids_from_xml(ids_out, xml_doc, num_samples, aux)))
    {
        return err;
    }

    if ((err = read_labels_from_xml(labels_out, xml_doc, num_samples, aux)))
    {
        return err;
    }

    if ((err = read_markers_from_xml(markers_out, xml_doc, num_samples, aux)))
    {
        return err;
    }

    return NULL;
}


Error* read_input
(
    Matblock_u8** ids_out,
    Vector_u32**  labels_out,
    Matrix_i32**  markers_out,
    const char*   fname
)
{
    uint32_t file_len;
    char*    file_buf = NULL;
    Error*   err;
    int  slen = 256;
    char str[slen];

    if ((err = read_file_into_buffer(&file_buf, &file_len, fname)))
    {
        return err;
    }

    switch (opts.input_format)
    {
        case HAPLO_INPUT_TXT:
        case HAPLO_INPUT_CSV:
            if ((err = read_input_from_txt_csv(ids_out, labels_out, markers_out,
                            file_buf, file_len, 0)))
            {
                snprintf(str, slen, "%s: %s", fname, err->msg);
                return JWSC_EIO(str);
            }
            break;
        case HAPLO_INPUT_XML:
            if ((err = read_input_from_xml(ids_out, labels_out, markers_out,
                            file_buf, file_len, fname, 0)))
            {
                snprintf(str, slen, "%s: %s", fname, err->msg);
                return JWSC_EIO(str);
            }
            break;
    }

    free(file_buf);

    return err;
}


Error* read_aux_input
(
    Matblock_u8** ids_out,
    Vector_u32**  labels_out,
    Matrix_i32**  markers_out,
    const char*   fname
)
{
    uint32_t file_len;
    char*    file_buf = NULL;
    Error*   err;
    int  slen = 256;
    char str[slen];

    if ((err = read_file_into_buffer(&file_buf, &file_len, fname)))
    {
        return err;
    }

    switch (opts.input_format)
    {
        case HAPLO_INPUT_TXT:
        case HAPLO_INPUT_CSV:
            if ((err = read_input_from_txt_csv(ids_out, labels_out,
                            markers_out, file_buf, file_len, 1)))
            {
                snprintf(str, slen, "%s: %s", fname, err->msg);
                return JWSC_EIO(str);
            }
            break;
        case HAPLO_INPUT_XML:
            if ((err = read_input_from_xml(ids_out, labels_out, markers_out,
                            file_buf, file_len, fname, 1)))
            {
                snprintf(str, slen, "%s: %s", fname, err->msg);
                return JWSC_EIO(str);
            }
            break;
    }

    free(file_buf);

    return err;
}


Error* impute_marker_from_parent_of_haplogroup_index(
        uint32_t haplo_group_index,
        uint32_t marker_no,
        uint32_t sample_no,
        Matrix_i32* marker_sums,
        Matrix_i32* imp_markers
        ){
                
                uint32_t parent_i;  // where the index of the parent will be stored
                Error* e= lookup_haplo_group_ancestor_index_from_haplo_group_index(&parent_i,haplo_group_index);
                if(e){
                        return e;
                }
                
        
                if (marker_sums->elts[ parent_i ][ marker_no ] == 0)
            {
                        return impute_marker_from_parent_of_haplogroup_index(parent_i,marker_no,sample_no,marker_sums,imp_markers);
            }
            else
            {
                        
                imp_markers->elts[ sample_no ][ marker_no ] = marker_sums->elts[ parent_i ][ marker_no];
            }

                return NULL;
                
        }


Error* impute_missing_markers_from_avg
(
    const Vector_u32* imp_labels, 
    Matrix_i32*       imp_markers,
    const Vector_u32* src_labels, 
    const Matrix_i32* src_markers
)
{
    uint32_t s, num_imp_samples;
    uint32_t num_src_samples;
    uint32_t m, num_markers;
    uint32_t l, ll, num_labels;

    Matrix_i32* marker_sums = NULL;
    Matrix_u32* marker_cnts = NULL;
    int  slen = 256;
    char str[slen];

    assert(imp_markers->num_cols == src_markers->num_cols);

    num_imp_samples = imp_markers->num_rows;
    num_src_samples = src_markers->num_rows;
    num_markers = imp_markers->num_cols;
    num_labels  = get_num_haplo_groups();

    create_zero_matrix_i32(&marker_sums, num_labels, num_markers);
    create_zero_matrix_u32(&marker_cnts, num_labels, num_markers);

    for (s = 0; s < num_src_samples; s++)
    {
        l = src_labels->elts[ s ];
        for (ll = 0; ll < num_labels; ll++)
        {
            if (is_ancestor(l, ll))
            {
                for (m = 0; m < num_markers; m++)
                {
                    if (src_markers->elts[ s ][ m ])
                    {
                        marker_sums->elts[ ll ][ m ] += 
                            src_markers->elts[ s ][ m ];
                        marker_cnts->elts[ ll ][ m ]++;
                    }
                }
            }
        }
    }

    for (l = 0; l < num_labels; l++)
    {
        for (m = 0; m < num_markers; m++)
        {
            if (marker_cnts->elts[ l ][ m ])
            {
                marker_sums->elts[ l ][ m ] /= marker_cnts->elts[ l ][ m ];
            }
        }       
    }

    for (s = 0; s < num_imp_samples; s++)
    {
        l = imp_labels->elts[ s ];
        for (m = 0; m < num_markers; m++)
        {
            if (!imp_markers->elts[ s ][ m ] && 
                    marker_cnts->elts[ l ][ m ] == 0)
            {
                if(impute_marker_from_parent_of_haplogroup_index(l,m,s,marker_sums,imp_markers)){
                        snprintf(str, slen, "Sample %d: %s", s+1, 
                                "No data to impute missing markers");
                        return JWSC_EARG(str);
                  }
            }
            else if (!imp_markers->elts[ s ][ m ])
            {
                imp_markers->elts[ s ][ m ] = marker_sums->elts[ l ][ m ];
            }
        }
    }

    free_matrix_u32(marker_cnts);
    free_matrix_i32(marker_sums);

    return NULL;
}


Error* impute_missing_markers_from_nn
(
    Matrix_i32*       imp_markers,
    const Matrix_i32* src_markers
)
{
    // TODO
    abort();

    return NULL;
}