from __future__ import division, print_function

import os

import sys

from libc.errno cimport errno, EPIPE

from libc.string cimport strcmp, strpbrk, strerror

from libc.stdint cimport INT8_MAX, INT16_MAX, INT32_MAX

cimport cython

from cpython.object cimport PyObject

from cpython.ref cimport Py_INCREF

from cpython.dict cimport PyDict_GetItemString, PyDict_SetItemString

from cpython.tuple cimport PyTuple_New, PyTuple_SET_ITEM

from cpython.bytes cimport PyBytes_FromStringAndSize

from cpython.unicode cimport PyUnicode_DecodeUTF8

from pysam.libchtslib cimport HTSFile, hisremote

from pysam.utils import unquoted_str

__all__ = ['VariantFile',

'VariantHeader',

'VariantHeaderRecord',

'VariantHeaderRecords',

'VariantMetadata',

'VariantHeaderMetadata',

'VariantContig',

'VariantHeaderContigs',

'VariantHeaderSamples',

'VariantRecordFilter',

'VariantRecordFormat',

'VariantRecordInfo',

'VariantRecordSamples',

'VariantRecord',

'VariantRecordSample',

'BaseIndex',

'BCFIndex',

'TabixIndex',

'BaseIterator',

'BCFIterator',

'TabixIterator',

'VariantRecord']

cdef int MAX_POS = (1 << 31) - 1

cdef tuple VALUE_TYPES = ('Flag', 'Integer', 'Float', 'String')

cdef tuple METADATA_TYPES = ('FILTER', 'INFO', 'FORMAT', 'CONTIG', 'STRUCTURED', 'GENERIC')

cdef tuple METADATA_LENGTHS = ('FIXED', 'VARIABLE', 'A', 'G', 'R')

from pysam.libcutils cimport force_bytes, force_str, charptr_to_str, charptr_to_str_w_len

from pysam.libcutils cimport OSError_from_errno, encode_filename, from_string_and_size, decode_bytes

cdef object _nothing = object()

cdef dict bcf_str_cache = {}

cdef inline bcf_str_cache_get_charptr(const char* s):

if s == NULL:

return None

cdef PyObject *pystr = PyDict_GetItemString(bcf_str_cache, s)

if pystr:

return <object>pystr

val = PyUnicode_DecodeUTF8(s, strlen(s), NULL)

PyDict_SetItemString(bcf_str_cache, s, val)

return val

cdef int comb(int n, int k) except -1:

"""Return binomial coefficient: n choose k

if k > n:

return 0

elif k == n:

return 1

elif k > n // 2:

k = n - k

cdef d, result

d = result = n - k + 1

for i in range(2, k + 1):

d += 1

result *= d

result //= i

return result

cdef inline int bcf_geno_combinations(int ploidy, int alleles) except -1:

"""Return the count of genotypes expected for the given ploidy and number of alleles.

return comb(alleles + ploidy - 1, ploidy)

cdef inline bint check_header_id(bcf_hdr_t *hdr, int hl_type, int id):

return id >= 0 and id < hdr.n[BCF_DT_ID] and bcf_hdr_idinfo_exists(hdr, hl_type, id)

cdef inline int is_gt_fmt(bcf_hdr_t *hdr, int fmt_id):

return strcmp(bcf_hdr_int2id(hdr, BCF_DT_ID, fmt_id), 'GT') == 0

cdef inline int bcf_genotype_count(bcf_hdr_t *hdr, bcf1_t *rec, int sample) except -1:

if sample < 0:

raise ValueError('genotype is only valid as a format field')

cdef int32_t *gt_arr = NULL

cdef int ngt = 0

ngt = bcf_get_genotypes(hdr, rec, &gt_arr, &ngt)

if ngt <= 0 or not gt_arr:

return 0

assert ngt % rec.n_sample == 0

cdef int max_ploidy = ngt // rec.n_sample

cdef int32_t *gt = gt_arr + sample * max_ploidy

cdef int ploidy = 0

while ploidy < max_ploidy and gt[0] != bcf_int32_vector_end:

gt += 1

ploidy += 1

free(<void*>gt_arr)

return bcf_geno_combinations(ploidy, rec.n_allele)

cdef tuple char_array_to_tuple(const char **a, ssize_t n, int free_after=0):

if not a:

return None

try:

return tuple(charptr_to_str(a[i]) for i in range(n))

finally:

if free_after and a:

free(a)

cdef bcf_array_to_object(void *data, int type, ssize_t n, ssize_t count, int scalar):

cdef char *datac

cdef int8_t *data8

cdef int16_t *data16

cdef int32_t *data32

cdef float *dataf

cdef int i

cdef bytes b

if not data or n <= 0:

return None

if type == BCF_BT_CHAR:

datac = <char *>data

if not n:

value = ()

else:

# Check if at least one null terminator is present

if datac[n-1] == bcf_str_vector_end:

# If so, create a string up to the first null terminator

b = datac

else:

# Otherwise, copy the entire block

b = datac[:n]

value = tuple(decode_bytes(v, 'utf-8') if v and v != bcf_str_missing else None for v in b.split(b','))

else:

value = []

if type == BCF_BT_INT8:

data8 = <int8_t *>data

for i in range(n):

if data8[i] == bcf_int8_vector_end:

break

value.append(data8[i] if data8[i] != bcf_int8_missing else None)

elif type == BCF_BT_INT16:

data16 = <int16_t *>data

for i in range(n):

if data16[i] == bcf_int16_vector_end:

break

value.append(data16[i] if data16[i] != bcf_int16_missing else None)

elif type == BCF_BT_INT32:

data32 = <int32_t *>data

for i in range(n):

if data32[i] == bcf_int32_vector_end:

break

value.append(data32[i] if data32[i] != bcf_int32_missing else None)

elif type == BCF_BT_FLOAT:

dataf = <float *>data

for i in range(n):

if bcf_float_is_vector_end(dataf[i]):

break

value.append(dataf[i] if not bcf_float_is_missing(dataf[i]) else None)

else:

raise TypeError('unsupported info type code')

# FIXME: Need to know length? Report errors? Pad with missing values? Not clear what to do.

if not value:

if scalar:

value = None

elif count <= 0:

value = ()

else:

value = (None,)*count

elif scalar and len(value) == 1:

value = value[0]

else:

value = tuple(value)

return value

cdef bcf_object_to_array(values, void *data, int bt_type, ssize_t n, int vlen):

cdef char *datac

cdef int8_t *data8

cdef int16_t *data16

cdef int32_t *data32

cdef float *dataf

cdef ssize_t i, value_count = len(values)

assert value_count <= n

if bt_type == BCF_BT_CHAR:

if not isinstance(values, (str, bytes)):

values = b','.join(force_bytes(v) if v else bcf_str_missing for v in values)

value_count = len(values)

assert value_count <= n

datac = <char *>data

memcpy(datac, <char *>values, value_count)

for i in range(value_count, n):

datac[i] = 0

elif bt_type == BCF_BT_INT8:

datai8 = <int8_t *>data

for i in range(value_count):

val = values[i]

datai8[i] = val if val is not None else bcf_int8_missing

for i in range(value_count, n):

datai8[i] = bcf_int8_vector_end

elif bt_type == BCF_BT_INT16:

datai16 = <int16_t *>data

for i in range(value_count):

val = values[i]

datai16[i] = val if val is not None else bcf_int16_missing

for i in range(value_count, n):

datai16[i] = bcf_int16_vector_end

elif bt_type == BCF_BT_INT32:

datai32 = <int32_t *>data

for i in range(value_count):

val = values[i]

datai32[i] = val if val is not None else bcf_int32_missing

for i in range(value_count, n):

datai32[i] = bcf_int32_vector_end

elif bt_type == BCF_BT_FLOAT:

dataf = <float *>data

for i in range(value_count):

val = values[i]

if val is None:

bcf_float_set(dataf + i, bcf_float_missing)

else:

dataf[i] = val

for i in range(value_count, n):

bcf_float_set(dataf + i, bcf_float_vector_end)

else:

raise TypeError('unsupported type')

cdef bcf_empty_array(int type, ssize_t n, int vlen):

cdef char *datac

cdef int32_t *data32

cdef float *dataf

cdef int i

if n <= 0:

raise ValueError('Cannot create empty array')

if type == BCF_HT_STR:

value = PyBytes_FromStringAndSize(NULL, sizeof(char)*n)

datac = <char *>value

for i in range(n):

datac[i] = bcf_str_missing if not vlen else bcf_str_vector_end

elif type == BCF_HT_INT:

value = PyBytes_FromStringAndSize(NULL, sizeof(int32_t)*n)

data32 = <int32_t *><char *>value

for i in range(n):

data32[i] = bcf_int32_missing if not vlen else bcf_int32_vector_end

elif type == BCF_HT_REAL:

value = PyBytes_FromStringAndSize(NULL, sizeof(float)*n)

dataf = <float *><char *>value

for i in range(n):

bcf_float_set(dataf + i, bcf_float_missing if not vlen else bcf_float_vector_end)

else:

raise TypeError('unsupported header type code')

return value

cdef bcf_copy_expand_array(void *src_data, int src_type, size_t src_values,

void *dst_data, int dst_type, size_t dst_values,

int vlen):

"""copy data from src to dest where the size of the elements (src_type/dst_type) differ

cdef char *src_datac

cdef char *dst_datac

cdef int8_t *src_datai8

cdef int16_t *src_datai16

cdef int32_t *src_datai32

cdef int32_t *dst_datai

cdef float *src_dataf

cdef float *dst_dataf

cdef ssize_t src_size, dst_size, i, j

cdef int val

if src_values > dst_values:

raise ValueError('Cannot copy arrays with src_values={} > dst_values={}'.format(src_values, dst_values))

if src_type == dst_type == BCF_BT_CHAR:

src_datac = <char *>src_data

dst_datac = <char *>dst_data

memcpy(dst_datac, src_datac, src_values)

for i in range(src_values, dst_values):

dst_datac[i] = 0

elif src_type == BCF_BT_INT8 and dst_type == BCF_BT_INT32:

src_datai8 = <int8_t *>src_data

dst_datai = <int32_t *>dst_data

for i in range(src_values):

val = src_datai8[i]

if val == bcf_int8_missing:

val = bcf_int32_missing

elif val == bcf_int8_vector_end:

val = bcf_int32_vector_end

dst_datai[i] = val

for i in range(src_values, dst_values):

dst_datai[i] = bcf_int32_missing if not vlen else bcf_int32_vector_end

elif src_type == BCF_BT_INT16 and dst_type == BCF_BT_INT32:

src_datai16 = <int16_t *>src_data

dst_datai = <int32_t *>dst_data

for i in range(src_values):

val = src_datai16[i]

if val == bcf_int16_missing:

val = bcf_int32_missing

elif val == bcf_int16_vector_end:

val = bcf_int32_vector_end

dst_datai[i] = val

for i in range(src_values, dst_values):

dst_datai[i] = bcf_int32_missing if not vlen else bcf_int32_vector_end

elif src_type == BCF_BT_INT32 and dst_type == BCF_BT_INT32:

src_datai32 = <int32_t *>src_data

dst_datai = <int32_t *>dst_data

for i in range(src_values):

dst_datai[i] = src_datai32[i]

for i in range(src_values, dst_values):

dst_datai[i] = bcf_int32_missing if not vlen else bcf_int32_vector_end

elif src_type == BCF_BT_FLOAT and dst_type == BCF_BT_FLOAT:

src_dataf = <float *>src_data

dst_dataf = <float *>dst_data

for i in range(src_values):

dst_dataf[i] = src_dataf[i]

for i in range(src_values, dst_values):

bcf_float_set(dst_dataf + i, bcf_float_missing if not vlen else bcf_float_vector_end)

else:

raise TypeError('unsupported types')

cdef bcf_get_value_count(VariantRecord record, int hl_type, int id, ssize_t *count, int *scalar, int sample):

if record is None:

raise ValueError('record must not be None')

cdef bcf_hdr_t *hdr = record.header.ptr

cdef bcf1_t *r = record.ptr

if not check_header_id(hdr, hl_type, id):

raise ValueError('Invalid header')

cdef int length = bcf_hdr_id2length(hdr, hl_type, id)

cdef int number = bcf_hdr_id2number(hdr, hl_type, id)

scalar[0] = 0

if hl_type == BCF_HL_FMT and is_gt_fmt(hdr, id):

count[0] = number

elif length == BCF_VL_FIXED:

if number == 1:

scalar[0] = 1

count[0] = number

elif length == BCF_VL_R:

count[0] = r.n_allele

elif length == BCF_VL_A:

count[0] = r.n_allele - 1

elif length == BCF_VL_G:

count[0] = bcf_genotype_count(hdr, r, sample)

elif length == BCF_VL_VAR:

count[0] = -1

else:

raise ValueError('Unknown format length')

cdef object bcf_info_get_value(VariantRecord record, const bcf_info_t *z):

if record is None:

raise ValueError('record must not be None')

cdef bcf_hdr_t *hdr = record.header.ptr

cdef char *s

cdef ssize_t count

cdef int scalar

bcf_get_value_count(record, BCF_HL_INFO, z.key, &count, &scalar, -1)

if z.len == 0:

if bcf_hdr_id2type(hdr, BCF_HL_INFO, z.key) == BCF_HT_FLAG:

value = True

elif scalar:

value = None

else:

value = ()

elif z.len == 1:

if z.type == BCF_BT_INT8:

value = z.v1.i if z.v1.i != bcf_int8_missing else None

elif z.type == BCF_BT_INT16:

value = z.v1.i if z.v1.i != bcf_int16_missing else None

elif z.type == BCF_BT_INT32:

value = z.v1.i if z.v1.i != bcf_int32_missing else None

elif z.type == BCF_BT_FLOAT:

value = z.v1.f if not bcf_float_is_missing(z.v1.f) else None

elif z.type == BCF_BT_CHAR:

value = force_str(chr(z.v1.i))

else:

raise TypeError('unsupported info type code')

if not scalar and value != ():

value = (value,)

else:

value = bcf_array_to_object(z.vptr, z.type, z.len, count, scalar)

return value

cdef object bcf_check_values(VariantRecord record, value, int sample,

int hl_type, int ht_type,

int id, int bt_type, ssize_t bt_len,

ssize_t *value_count, int *scalar, int *realloc):

if record is None:

raise ValueError('record must not be None')

bcf_get_value_count(record, hl_type, id, value_count, scalar, sample)

# Validate values now that we know the type and size

values = (value,) if not isinstance(value, (list, tuple)) else value

# Validate values now that we know the type and size

if ht_type == BCF_HT_FLAG:

value_count[0] = 1

elif hl_type == BCF_HL_FMT and is_gt_fmt(record.header.ptr, id):

# KBJ: htslib lies about the cardinality of GT fields-- they're really VLEN (-1)

value_count[0] = -1

cdef int given = len(values)

if value_count[0] != -1 and value_count[0] != given:

if scalar[0]:

raise TypeError('value expected to be scalar, given len={}'.format(given))

else:

raise TypeError('values expected to be {}-tuple, given len={}'.format(value_count[0], given))

if ht_type == BCF_HT_REAL:

for v in values:

if not(v is None or isinstance(v, (float, int))):

raise TypeError('invalid value for Float format')

elif ht_type == BCF_HT_INT:

for v in values:

if not(v is None or (isinstance(v, (float, int)) and int(v) == v)):

raise TypeError('invalid value for Integer format')

for v in values:

if not(v is None or bcf_int32_missing < v <= INT32_MAX):

raise ValueError('Integer value too small/large to store in VCF/BCF')

elif ht_type == BCF_HT_STR:

values = b','.join(force_bytes(v) if v is not None else b'' for v in values)

elif ht_type == BCF_HT_FLAG:

if values[0] not in (True, False, None, 1, 0):

raise ValueError('Flag values must be: True, False, None, 1, 0')

else:

raise TypeError('unsupported type')

realloc[0] = 0

if len(values) <= 1 and hl_type == BCF_HL_INFO:

realloc[0] = 0

elif len(values) > bt_len:

realloc[0] = 1

elif bt_type == BCF_BT_INT8:

for v in values:

if v is not None and not(bcf_int8_missing < v <= INT8_MAX):

realloc[0] = 1

break

elif bt_type == BCF_BT_INT16:

for v in values:

if v is not None and not(bcf_int16_missing < v <= INT16_MAX):

realloc[0] = 1

break

return values

cdef bcf_encode_alleles(VariantRecord record, values):

if record is None:

raise ValueError('record must not be None')

cdef bcf1_t *r = record.ptr

cdef int32_t nalleles = r.n_allele

cdef list gt_values = []

cdef char *s

cdef int i

if values is None:

return ()

if not isinstance(values, (list, tuple)):

values = (values,)

for value in values:

if value is None:

gt_values.append(bcf_gt_missing)

elif isinstance(value, (str, bytes)):

bvalue = force_bytes(value)

s = bvalue

for i in range(r.n_allele):

if strcmp(r.d.allele[i], s) != 0:

gt_values.append(bcf_gt_unphased(i))

break

else:

raise ValueError('Unknown allele')

else:

i = value

if not (0 <= i < nalleles):

raise ValueError('Invalid allele index')

gt_values.append(bcf_gt_unphased(i))

return gt_values

cdef bcf_info_set_value(VariantRecord record, key, value):

if record is None:

raise ValueError('record must not be None')

cdef bcf_hdr_t *hdr = record.header.ptr

cdef bcf1_t *r = record.ptr

cdef int info_id, info_type, scalar, dst_type, realloc, vlen = 0

cdef ssize_t i, value_count, alloc_len, alloc_size, dst_size

if bcf_unpack(r, BCF_UN_INFO) < 0:

raise ValueError('Error unpacking VariantRecord')

cdef bytes bkey = force_bytes(key)

cdef bcf_info_t *info = bcf_get_info(hdr, r, bkey)

if info:

info_id = info.key

else:

info_id = bcf_header_get_info_id(hdr, bkey)

if info_id < 0:

raise KeyError('unknown INFO: {}'.format(key))

if not check_header_id(hdr, BCF_HL_INFO, info_id):

raise ValueError('Invalid header')

info_type = bcf_hdr_id2type(hdr, BCF_HL_INFO, info_id)

values = bcf_check_values(record, value, -1,

BCF_HL_INFO, info_type, info_id,

info.type if info else -1,

info.len if info else -1,

&value_count, &scalar, &realloc)

if info_type == BCF_HT_FLAG:

if bcf_update_info(hdr, r, bkey, NULL, bool(values[0]), info_type) < 0:

raise ValueError('Unable to update INFO values')

return

vlen = value_count < 0

value_count = len(values)

# DISABLED DUE TO ISSUES WITH THE CRAZY POINTERS

# If we can, write updated values to existing allocated storage

if 0 and info and not realloc:

r.d.shared_dirty |= BCF1_DIRTY_INF

if value_count == 0:

info.len = 0

if not info.vptr:

info.vptr = <uint8_t *>&info.v1.i

elif value_count == 1:

# FIXME: Check if need to free vptr if info.len > 0?

if info.type == BCF_BT_INT8 or info.type == BCF_BT_INT16 or info.type == BCF_BT_INT32:

bcf_object_to_array(values, &info.v1.i, BCF_BT_INT32, 1, vlen)

elif info.type == BCF_BT_FLOAT:

bcf_object_to_array(values, &info.v1.f, BCF_BT_FLOAT, 1, vlen)

else:

raise TypeError('unsupported info type code')

info.len = 1

if not info.vptr:

info.vptr = <uint8_t *>&info.v1.i

else:

bcf_object_to_array(values, info.vptr, info.type, info.len, vlen)

return

alloc_len = max(1, value_count)

if info and info.len > alloc_len:

alloc_len = info.len

new_values = bcf_empty_array(info_type, alloc_len, vlen)

cdef char *valp = <char *>new_values

if info_type == BCF_HT_INT:

dst_type = BCF_BT_INT32

elif info_type == BCF_HT_REAL:

dst_type = BCF_BT_FLOAT

elif info_type == BCF_HT_STR:

dst_type = BCF_BT_CHAR

else:

raise ValueError('Unsupported INFO type')

bcf_object_to_array(values, valp, dst_type, alloc_len, vlen)

if bcf_update_info(hdr, r, bkey, valp, <int>alloc_len, info_type) < 0:

raise ValueError('Unable to update INFO values')

cdef bcf_info_del_value(VariantRecord record, key):

if record is None:

raise ValueError('record must not be None')

cdef bcf_hdr_t *hdr = record.header.ptr

cdef bcf1_t *r = record.ptr

cdef ssize_t value_count

cdef int scalar

if bcf_unpack(r, BCF_UN_INFO) < 0:

raise ValueError('Error unpacking VariantRecord')

cdef bytes bkey = force_bytes(key)

cdef bcf_info_t *info = bcf_get_info(hdr, r, bkey)

if not info:

raise KeyError(key)

bcf_get_value_count(record, BCF_HL_INFO, info.key, &value_count, &scalar, -1)

if value_count <= 0:

null_value = ()

elif scalar:

null_value = None

else:

null_value = (None,)*value_count

bcf_info_set_value(record, bkey, null_value)

cdef bcf_format_get_value(VariantRecordSample sample, key):

if sample is None:

raise ValueError('sample must not be None')

cdef bcf_hdr_t *hdr = sample.record.header.ptr

cdef bcf1_t *r = sample.record.ptr

cdef ssize_t count

cdef int scalar

if bcf_unpack(r, BCF_UN_ALL) < 0:

raise ValueError('Error unpacking VariantRecord')

cdef bytes bkey = force_bytes(key)

cdef bcf_fmt_t *fmt = bcf_get_fmt(hdr, r, bkey)

if not fmt or not fmt.p:

raise KeyError('invalid FORMAT: {}'.format(key))

if is_gt_fmt(hdr, fmt.id):

return bcf_format_get_allele_indices(sample)

bcf_get_value_count(sample.record, BCF_HL_FMT, fmt.id, &count, &scalar, sample.index)

if fmt.p and fmt.n and fmt.size:

return bcf_array_to_object(fmt.p + sample.index * fmt.size, fmt.type, fmt.n, count, scalar)

elif scalar:

return None

elif count <= 0:

return ()

else:

return (None,)*count

cdef bcf_format_set_value(VariantRecordSample sample, key, value):

if sample is None:

raise ValueError('sample must not be None')

if key == 'phased':

sample.phased = bool(value)

return

cdef bcf_hdr_t *hdr = sample.record.header.ptr

cdef bcf1_t *r = sample.record.ptr

cdef int fmt_id

cdef vdict_t *d

cdef khiter_t k

cdef int fmt_type, scalar, realloc, dst_type, vlen = 0

cdef ssize_t i, nsamples, value_count, alloc_size, alloc_len, dst_size

if bcf_unpack(r, BCF_UN_ALL) < 0:

raise ValueError('Error unpacking VariantRecord')

cdef bytes bkey = force_bytes(key)

cdef bcf_fmt_t *fmt = bcf_get_fmt(hdr, r, bkey)

if fmt:

fmt_id = fmt.id

else:

d = <vdict_t *>hdr.dict[BCF_DT_ID]

k = kh_get_vdict(d, bkey)

if k == kh_end(d) or kh_val_vdict(d, k).info[BCF_HL_FMT] & 0xF == 0xF:

raise KeyError('unknown format: {}'.format(key))

fmt_id = kh_val_vdict(d, k).id

if not check_header_id(hdr, BCF_HL_FMT, fmt_id):

raise ValueError('Invalid header')

fmt_type = bcf_hdr_id2type(hdr, BCF_HL_FMT, fmt_id)

if fmt_type == BCF_HT_FLAG:

raise ValueError('Flag types are not allowed on FORMATs')

if is_gt_fmt(hdr, fmt_id):

value = bcf_encode_alleles(sample.record, value)

# KBJ: GT field is considered to be a string by the VCF header but BCF represents it as INT.

fmt_type = BCF_HT_INT

values = bcf_check_values(sample.record, value, sample.index,

BCF_HL_FMT, fmt_type, fmt_id,

fmt.type if fmt else -1,

fmt.n if fmt else -1,

&value_count, &scalar, &realloc)

vlen = value_count < 0

value_count = len(values)

# If we can, write updated values to existing allocated storage.

if fmt and not realloc:

r.d.indiv_dirty = 1

bcf_object_to_array(values, fmt.p + sample.index * fmt.size, fmt.type, fmt.n, vlen)

return

alloc_len = max(1, value_count)

if fmt and fmt.n > alloc_len:

alloc_len = fmt.n

nsamples = r.n_sample

new_values = bcf_empty_array(fmt_type, nsamples * alloc_len, vlen)

cdef char *new_values_p = <char *>new_values

if fmt_type == BCF_HT_INT:

dst_type = BCF_BT_INT32

dst_size = sizeof(int32_t) * alloc_len

elif fmt_type == BCF_HT_REAL:

dst_type = BCF_BT_FLOAT

dst_size = sizeof(float) * alloc_len

elif fmt_type == BCF_HT_STR:

dst_type = BCF_BT_CHAR

dst_size = sizeof(char) * alloc_len

else:

raise ValueError('Unsupported FORMAT type')

if fmt and nsamples > 1:

for i in range(nsamples):

bcf_copy_expand_array(fmt.p + i * fmt.size, fmt.type, fmt.n,

new_values_p + i * dst_size, dst_type, alloc_len,

vlen)

bcf_object_to_array(values, new_values_p + sample.index * dst_size, dst_type, alloc_len, vlen)

if bcf_update_format(hdr, r, bkey, new_values_p, <int>(nsamples * alloc_len), fmt_type) < 0:

raise ValueError('Unable to update format values')

cdef bcf_format_del_value(VariantRecordSample sample, key):

if sample is None:

raise ValueError('sample must not be None')

cdef bcf_hdr_t *hdr = sample.record.header.ptr

cdef bcf1_t *r = sample.record.ptr

cdef ssize_t value_count

cdef int scalar

if bcf_unpack(r, BCF_UN_ALL) < 0:

raise ValueError('Error unpacking VariantRecord')

cdef bytes bkey = force_bytes(key)

cdef bcf_fmt_t *fmt = bcf_get_fmt(hdr, r, bkey)

if not fmt or not fmt.p:

raise KeyError(key)

bcf_get_value_count(sample.record, BCF_HL_FMT, fmt.id, &value_count, &scalar, sample.index)

if value_count <= 0:

null_value = ()

elif scalar:

null_value = None

else:

null_value = (None,)*value_count

bcf_format_set_value(sample, bkey, null_value)

cdef bcf_format_get_allele_indices(VariantRecordSample sample):

if sample is None:

raise ValueError('sample must not be None')

cdef bcf_hdr_t *hdr = sample.record.header.ptr

cdef bcf1_t *r = sample.record.ptr

cdef int32_t n = r.n_sample

if bcf_unpack(r, BCF_UN_ALL) < 0:

raise ValueError('Error unpacking VariantRecord')

if sample.index < 0 or sample.index >= n or not r.n_fmt:

return ()

cdef bcf_fmt_t *fmt0 = r.d.fmt

cdef int gt0 = is_gt_fmt(hdr, fmt0.id)

if not gt0 or not fmt0.n:

return ()

cdef int8_t *data8

cdef int16_t *data16

cdef int32_t *data32

cdef int32_t a, nalleles = r.n_allele

cdef list alleles = []

if fmt0.type == BCF_BT_INT8:

data8 = <int8_t *>(fmt0.p + sample.index * fmt0.size)

for i in range(fmt0.n):

if data8[i] == bcf_int8_vector_end:

break

elif data8[i] == bcf_gt_missing:

a = -1

else:

a = bcf_gt_allele(data8[i])

alleles.append(a if 0 <= a < nalleles else None)

elif fmt0.type == BCF_BT_INT16:

data16 = <int16_t *>(fmt0.p + sample.index * fmt0.size)

for i in range(fmt0.n):

if data16[i] == bcf_int16_vector_end:

break

elif data16[i] == bcf_gt_missing:

a = -1

else:

a = bcf_gt_allele(data16[i])

alleles.append(a if 0 <= a < nalleles else None)

elif fmt0.type == BCF_BT_INT32:

data32 = <int32_t *>(fmt0.p + sample.index * fmt0.size)

for i in range(fmt0.n):

if data32[i] == bcf_int32_vector_end: