"""
Data providers for genome visualizations.
"""

import abc
import itertools
import logging
import math
import os
import random
import re
import sys
from collections.abc import Iterator
from contextlib import contextmanager
from json import loads
from typing import (
    Any,
    IO,
    Optional,
    Union,
)

import pysam
from bx.bbi.bigbed_file import BigBedFile
from bx.bbi.bigwig_file import BigWigFile
from bx.interval_index_file import Indexes

from galaxy.datatypes.interval import (
    Bed,
    Gff,
    Gtf,
)
from galaxy.datatypes.util.gff_util import (
    convert_gff_coords_to_bed,
    GFFFeature,
    GFFInterval,
    GFFReaderWrapper,
    parse_gff_attributes,
)
from galaxy.exceptions import MessageException
from galaxy.model import DatasetInstance
from galaxy.visualization.data_providers.basic import BaseDataProvider
from galaxy.visualization.data_providers.cigar import get_ref_based_read_seq_and_cigar

log = logging.getLogger(__name__)

IntWebParam = Union[str, int]

#
# Utility functions.
#

# pysam 0.16.0.1 emits logs containing the word 'Error', this can confuse the stdout/stderr checkers.
# Can be removed once https://github.com/pysam-developers/pysam/issues/939 is resolved.
pysam.set_verbosity(0)

PAYLOAD_LIST_TYPE = list[Optional[Union[str, int, float, list[tuple[int, int]]]]]


def float_nan(n):
    """
    Return None instead of NaN to pass jQuery 1.4's strict JSON
    """
    if n != n:  # NaN != NaN
        return None
    else:
        return float(n)


def get_bounds(reads, start_pos_index, end_pos_index):
    """
    Returns the minimum and maximum position for a set of reads.
    """
    max_low = sys.maxsize
    max_high = -sys.maxsize
    for read in reads:
        if read[start_pos_index] < max_low:
            max_low = read[start_pos_index]
        if read[end_pos_index] > max_high:
            max_high = read[end_pos_index]
    return max_low, max_high


def _convert_between_ucsc_and_ensemble_naming(chrom):
    """
    Convert between UCSC chromosome ('chr1') naming conventions and Ensembl
    naming conventions ('1')
    """
    if chrom.startswith("chr"):
        # Convert from UCSC to Ensembl
        return chrom[3:]
    else:
        # Convert from Ensembl to UCSC
        return "chr" + chrom


def _chrom_naming_matches(chrom1, chrom2):
    return (chrom1.startswith("chr") and chrom2.startswith("chr")) or (
        not chrom1.startswith("chr") and not chrom2.startswith("chr")
    )


# Was previously a BaseDataProvider but it doesn't have the same interface or use any of
# the parent methods. It is also only used explicitly, constructed by name, and access with
# get_data directly using a specified query.
class FeatureLocationIndexDataProvider:
    """
    Reads/writes/queries feature location index (FLI) datasets.
    """

    def __init__(self, converted_dataset):
        self.converted_dataset = converted_dataset

    def get_data(self, query):
        if self.converted_dataset is None or not self.converted_dataset.is_ok:
            raise MessageException("The dataset is not available or is in an error state.")
        # Init.
        result = []
        with open(self.converted_dataset.get_file_name()) as textloc_file:
            line = textloc_file.readline()
            if not line:
                raise MessageException("The dataset is empty.")
            try:
                line_len = int(line)
            except ValueError:
                raise MessageException(f"Expected an integer at first line, but found: '{line}'")
            if line_len < 1:
                raise MessageException(f"The first line must be a positive integer, but found: {line_len}")

            file_len = os.path.getsize(self.converted_dataset.get_file_name())
            query = query.lower()

            # Find query in file using binary search.
            low = 0
            high = int(file_len / line_len)
            while low < high:
                mid: int = (low + high) // 2
                position = mid * line_len
                textloc_file.seek(position)

                # Compare line with query and update low, high.
                line = textloc_file.readline()
                if line < query:
                    low = mid + 1
                else:
                    high = mid

            # Need to move back one line because last line read may be included in
            # results.
            position = low * line_len
            textloc_file.seek(position)

            # At right point in file, generate hits.
            while True:
                line = textloc_file.readline()
                if not line.startswith(query):
                    break
                if line[-1:] == "\n":
                    line = line[:-1]
                result.append(line.split()[1:])

        return result


class GenomeDataProvider(BaseDataProvider):
    """
    Base class for genome data providers. All genome providers use BED coordinate
    format (0-based, half-open coordinates) for both queries and returned data.
    """

    dataset_type: str

    # Mapping from column name to payload data; this mapping is used to create
    # filters. Key is column name, value is a dict with mandatory key 'index'
    # and optional key 'name'. E.g. this defines column 4:
    # col_name_data_attr_mapping = {4 : { index: 5, name: 'Score' } }
    col_name_data_attr_mapping: dict[Union[str, int], dict] = {}

    def __init__(
        self,
        converted_dataset=None,
        original_dataset=None,
        dependencies=None,
        error_max_vals="Only the first %i %s in this region are displayed.",
    ):
        super().__init__(
            converted_dataset=converted_dataset,
            original_dataset=original_dataset,
            dependencies=dependencies,
            error_max_vals=error_max_vals,
        )

    def valid_chroms(self):
        """
        Returns chroms/contigs that the dataset contains
        """
        return None  # by default

    def has_data(self, chrom):
        """
        Returns true if dataset has data in the specified genome window, false
        otherwise.
        """
        raise Exception("Unimplemented Function")

    @contextmanager
    def open_data_file(self):
        """
        Open data file for reading data.
        """
        raise Exception("Unimplemented Function")

    def get_iterator(self, data_file, chrom, start, end, **kwargs) -> Iterator[str]:
        """
        Returns an iterator that provides data in the region chrom:start-end
        """
        raise Exception("Unimplemented Function")

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Process data from an iterator to a format that can be provided to client.
        """
        raise Exception("Unimplemented Function")

    def get_data(self, chrom: str, start: IntWebParam, end: IntWebParam, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Returns data in region defined by chrom, start, and end. start_val and
        max_vals are used to denote the data to return: start_val is the first element to
        return and max_vals indicates the number of values to return.

        Return value must be a dictionary with the following attributes:
            dataset_type, data
        """
        start, end = int(start), int(end)
        with self.open_data_file() as data_file:
            iterator = self.get_iterator(data_file, chrom, start, end, **kwargs)
            try:
                data = self.process_data(iterator, start_val, max_vals, start=start, end=end, **kwargs)
            except ValueError as e:
                err_msg = f"Could not return data, error was '{e}'"
                log.warning(err_msg, exc_info=True)
                raise MessageException(err_msg)
        return data

    def get_genome_data(self, chroms_info, **kwargs):
        """
        Returns data for complete genome.
        """
        genome_data = []
        for chrom_info in chroms_info["chrom_info"]:
            chrom = chrom_info["chrom"]
            chrom_len = chrom_info["len"]
            chrom_data = self.get_data(chrom, 0, chrom_len, **kwargs)
            # FIXME: data providers probably should never return None.
            # Some data providers return None when there's no data, so
            # create a dummy dict if necessary.
            if not chrom_data:
                chrom_data = {"data": None}
            chrom_data["region"] = f"{chrom}:{0}-{chrom_len}"
            genome_data.append(chrom_data)

        return {"data": genome_data, "dataset_type": self.dataset_type}

    def get_filters(self):
        """
        Returns filters for provider's data. Return value is a list of
        filters; each filter is a dictionary with the keys 'name', 'index', 'type'.
        NOTE: This method uses the original dataset's datatype and metadata to
        create the filters.
        """
        # Get column names.
        try:
            column_names = self.original_dataset.datatype.column_names  # type: ignore[attr-defined]
        except AttributeError:
            try:
                column_names = list(range(self.original_dataset.metadata.columns))
            except Exception:  # Give up
                return []

        # Dataset must have column types; if not, cannot create filters.
        try:
            column_types = self.original_dataset.metadata.column_types
        except AttributeError:
            return []

        # Create and return filters.
        filters = []
        if self.original_dataset.metadata.viz_filter_cols:
            for viz_col_index in self.original_dataset.metadata.viz_filter_cols:
                # Some columns are optional, so can't assume that a filter
                # column is in dataset.
                if viz_col_index >= len(column_names):
                    continue
                col_name = column_names[viz_col_index]
                # Make sure that column has a mapped index. If not, do not add filter.
                try:
                    attrs = self.col_name_data_attr_mapping[col_name]
                except KeyError:
                    continue
                filters.append({"name": attrs["name"], "type": column_types[viz_col_index], "index": attrs["index"]})
        return filters

    def get_default_max_vals(self):
        return 5000


#
# -- Base mixins and providers --
#


class FilterableMixin:
    original_dataset: DatasetInstance

    def get_filters(self):
        """Returns a dataset's filters."""
        # Get filters.
        # TODOs:
        # (a) might be useful to move this into each datatype's set_meta method;
        # (b) could look at first N lines to ensure GTF attribute types are consistent.
        #
        filters = []
        # HACK: first 8 fields are for drawing, so start filter column index at 9.
        filter_col = 8
        if isinstance(self.original_dataset.datatype, Gff):
            # Can filter by score and GTF attributes.
            filters = [
                {"name": "Score", "type": "number", "index": filter_col, "tool_id": "Filter1", "tool_exp_name": "c6"}
            ]
            filter_col += 1
            if isinstance(self.original_dataset.datatype, Gtf):
                # Create filters based on dataset metadata.
                for name, a_type in self.original_dataset.metadata.attribute_types.items():
                    if a_type in ["int", "float"]:
                        filters.append(
                            {
                                "name": name,
                                "type": "number",
                                "index": filter_col,
                                "tool_id": "gff_filter_by_attribute",
                                "tool_exp_name": name,
                            }
                        )
                        filter_col += 1
        elif isinstance(self.original_dataset.datatype, Bed):
            # Can filter by score column only.
            filters = [
                {"name": "Score", "type": "number", "index": filter_col, "tool_id": "Filter1", "tool_exp_name": "c5"}
            ]

        return filters


class TabixDataProvider(GenomeDataProvider, FilterableMixin):
    """
    Tabix index data provider for the Galaxy track browser.
    """

    dataset_type = "tabix"

    col_name_data_attr_mapping: dict[Union[str, int], dict] = {4: {"index": 4, "name": "Score"}}

    @contextmanager
    def open_data_file(self):
        # We create a symlink to the index file. This is
        # required until https://github.com/pysam-developers/pysam/pull/586 is merged.
        index_path = self.converted_dataset.get_file_name()
        with pysam.TabixFile(self.dependencies["bgzip"].get_file_name(), index=index_path) as f:
            yield f

    def get_iterator(self, data_file, chrom, start, end, **kwargs) -> Iterator[str]:
        # chrom must be a string, start/end integers.
        # in previous versions of pysam, unicode was accepted for chrom, but not in 8.4
        chrom = str(chrom)
        start = int(start)
        end = int(end)
        if end >= (2 << 29):
            end = 2 << 29 - 1  # Tabix-enforced maximum
        # Get iterator using either naming scheme.
        iterator: Iterator[str] = iter([])
        if chrom in data_file.contigs:
            iterator = data_file.fetch(reference=chrom, start=start, end=end)
        else:
            # Try alternative naming scheme.
            chrom = _convert_between_ucsc_and_ensemble_naming(chrom)
            if chrom in data_file.contigs:
                iterator = data_file.fetch(reference=chrom, start=start, end=end)

        return iterator


#
# -- Interval data providers --
#


class IntervalDataProvider(GenomeDataProvider):
    """
    Processes interval data from native format to payload format.

    Payload format: [ uid (offset), start, end, name, strand, thick_start, thick_end, blocks ]
    """

    dataset_type = "interval_index"

    def get_iterator(self, data_file, chrom, start, end, **kwargs):
        raise Exception("Unimplemented Function")

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Provides
        """
        # Build data to return. Payload format is:
        # [ <guid/offset>, <start>, <end>, <name>, <strand> ]
        #
        # First three entries are mandatory, others are optional.
        #
        filter_cols = loads(kwargs.get("filter_cols", "[]"))
        no_detail = "no_detail" in kwargs
        rval = []
        message = None

        # Subtract one b/c columns are 1-based but indices are 0-based.
        def col_fn(col):
            return None if col is None else col - 1

        start_col = self.original_dataset.metadata.startCol - 1
        end_col = self.original_dataset.metadata.endCol - 1
        strand_col = col_fn(self.original_dataset.metadata.strandCol)
        name_col = col_fn(self.original_dataset.metadata.nameCol)
        for count, line in enumerate(iterator):
            if count < start_val:
                continue
            if max_vals and count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "features")
                break

            feature = line.split()
            length = len(feature)
            # Unique id is just a hash of the line
            payload: PAYLOAD_LIST_TYPE = [hash(line), int(feature[start_col]), int(feature[end_col])]

            if no_detail:
                rval.append(payload)
                continue

            # Name, strand.
            if name_col:
                payload.append(feature[name_col])
            if strand_col:
                # Put empty name as placeholder.
                if not name_col:
                    payload.append("")
                payload.append(feature[strand_col])

            # Score (filter data)
            if length >= 5 and filter_cols and filter_cols[0] == "Score":
                try:
                    payload.append(float(feature[4]))
                except Exception:
                    payload.append(feature[4])

            rval.append(payload)

        return {"data": rval, "message": message}


class IntervalTabixDataProvider(TabixDataProvider, IntervalDataProvider):
    """
    Provides data from a BED file indexed via tabix.
    """


#
# -- BED data providers --
#


class BedDataProvider(GenomeDataProvider):
    """
    Processes BED data from native format to payload format.

    Payload format: [ uid (offset), start, end, name, strand, thick_start, thick_end, blocks ]
    """

    dataset_type = "interval_index"

    def get_iterator(self, data_file, chrom, start, end, **kwargs):
        raise Exception("Unimplemented Method")

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Provides
        """
        # Build data to return. Payload format is:
        # [ <guid/offset>, <start>, <end>, <name>, <strand>, <thick_start>,
        #   <thick_end>, <blocks> ]
        #
        # First three entries are mandatory, others are optional.
        #
        filter_cols = loads(kwargs.get("filter_cols", "[]"))
        no_detail = "no_detail" in kwargs
        rval = []
        message = None
        for count, line in enumerate(iterator):
            if count < start_val:
                continue
            if max_vals and count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "features")
                break
            # TODO: can we use column metadata to fill out payload?
            # TODO: use function to set payload data

            feature = line.split()
            length = len(feature)
            # Unique id is just a hash of the line
            payload: PAYLOAD_LIST_TYPE = [hash(line), int(feature[1]), int(feature[2])]

            if no_detail:
                rval.append(payload)
                continue

            # Name, strand, thick start, thick end.
            if length >= 4:
                payload.append(feature[3])
            if length >= 6:
                payload.append(feature[5])
            if length >= 8:
                payload.append(int(feature[6]))
                payload.append(int(feature[7]))

            # Blocks.
            if length >= 12:
                block_sizes = [int(n) for n in feature[10].split(",") if n != ""]
                block_starts = [int(n) for n in feature[11].split(",") if n != ""]
                blocks = list(zip(block_sizes, block_starts))
                payload.append(
                    [(int(feature[1]) + block[1], int(feature[1]) + block[1] + block[0]) for block in blocks]
                )

            # Score (filter data)
            if length >= 5 and filter_cols and filter_cols[0] == "Score":
                # If dataset doesn't have name/strand/thick start/thick end/blocks,
                # add placeholders. There should be 8 entries if all attributes
                # are present.
                payload.extend([None for i in range(8 - len(payload))])

                try:
                    payload.append(float(feature[4]))
                except Exception:
                    payload.append(feature[4])

            rval.append(payload)

        return {"data": rval, "dataset_type": self.dataset_type, "message": message}


class BedTabixDataProvider(TabixDataProvider, BedDataProvider):
    """
    Provides data from a BED file indexed via tabix.
    """


class RawBedDataProvider(BedDataProvider):
    """
    Provide data from BED file.

    NOTE: this data provider does not use indices, and hence will be very slow
    for large datasets.
    """

    def get_iterator(self, data_file, chrom, start, end, **kwargs):
        # Read first line in order to match chrom naming format.
        line = data_file.readline()
        dataset_chrom = line.split()[0]
        if not _chrom_naming_matches(chrom, dataset_chrom):
            chrom = _convert_between_ucsc_and_ensemble_naming(chrom)
        # Undo read.
        data_file.seek(0)

        def line_filter_iter():
            with open(self.original_dataset.get_file_name()) as data_file:
                for line in data_file:
                    if line.startswith("track") or line.startswith("browser"):
                        continue
                    feature = line.split()
                    feature_chrom = feature[0]
                    feature_start = int(feature[1])
                    feature_end = int(feature[2])
                    if (
                        (chrom is not None and feature_chrom != chrom)
                        or (start is not None and feature_start > end)
                        or (end is not None and feature_end < start)
                    ):
                        continue
                    yield line

        return line_filter_iter()


#
# -- VCF data providers --
#


class VcfDataProvider(GenomeDataProvider):
    """
    Abstract class that processes VCF data from native format to payload format.

    Payload format: An array of entries for each locus in the file. Each array
    has the following entries:

        1. GUID (unused)
        2. location (0-based)
        3. reference base(s)
        4. alternative base(s)
        5. quality score
        6. whether variant passed filter
        7. sample genotypes -- a single string with samples separated by commas; empty string
           denotes the reference genotype
        8. allele counts for each alternative

    """

    col_name_data_attr_mapping: dict[Union[str, int], dict] = {"Qual": {"index": 6, "name": "Qual"}}

    dataset_type = "variant"

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Returns a dict with the following attributes::

            data - a list of variants with the format

            .. raw:: text

                [<guid>, <start>, <end>, <name>, cigar, seq]

            message - error/informative message

        """
        data = []
        message = None

        def get_mapping(ref, alt):
            """
            Returns ( offset, new_seq, cigar ) tuple that defines mapping of
            alt to ref. Cigar format is an array of [ op_index, length ] pairs
            where op_index is the 0-based index into the string "MIDNSHP=X"
            """

            cig_ops = "MIDNSHP=X"

            ref_len = len(ref)
            alt_len = len(alt)

            # Substitutions?
            if ref_len == alt_len:
                return 0, alt, [[cig_ops.find("M"), ref_len]]

            # Deletions?
            alt_in_ref_index = ref.find(alt)
            if alt_in_ref_index != -1:
                return alt_in_ref_index, ref[alt_in_ref_index + 1 :], [[cig_ops.find("D"), ref_len - alt_len]]

            # Insertions?
            ref_in_alt_index = alt.find(ref)
            if ref_in_alt_index != -1:
                return ref_in_alt_index, alt[ref_in_alt_index + 1 :], [[cig_ops.find("I"), alt_len - ref_len]]

        # Pack data.
        genotype_re = re.compile(r"/|\|")
        for count, line in enumerate(iterator):
            if count < start_val:
                continue
            if max_vals and count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "features")
                break

            # Split line and aggregate data.
            feature = line.split()
            pos, c_id, ref, alt, qual, c_filter, info = feature[1:8]

            # Format and samples data are optional.
            samples_data = []
            if len(feature) > 8:
                samples_data = feature[9:]

            # VCF is 1-based but provided position is 0-based.
            pos = int(pos) - 1

            # FIXME: OK to skip?
            if alt == ".":
                count -= 1
                continue

            # Set up array to track allele counts.
            allele_counts = [0 for i in range(alt.count(",") + 1)]
            sample_gts = []

            if samples_data:
                # Process and pack samples' genotype and count alleles across samples.
                alleles_seen: dict[int, bool] = {}
                has_alleles = False

                for sample in samples_data:
                    # Parse and count alleles.
                    genotype = sample.split(":")[0]
                    has_alleles = False
                    alleles_seen.clear()
                    for allele_str in genotype_re.split(genotype):
                        try:
                            # This may throw a ValueError if allele is missing.
                            allele = int(allele_str)

                            # Only count allele if it hasn't been seen yet.
                            if allele != 0 and allele not in alleles_seen:
                                allele_counts[allele - 1] += 1
                                alleles_seen[allele] = True
                                has_alleles = True
                        except ValueError:
                            pass

                    # If no alleles, use empty string as proxy.
                    if not has_alleles:
                        genotype = ""

                    sample_gts.append(genotype)
            else:
                # No samples, so set allele count and sample genotype manually.
                allele_counts = [1]
                sample_gts = ["1/1"]

            # Add locus data.
            locus_data = [-1, pos, c_id, ref, alt, qual, c_filter, ",".join(sample_gts)]
            locus_data.extend(allele_counts)
            data.append(locus_data)

        return {"data": data, "message": message}


class VcfTabixDataProvider(TabixDataProvider, VcfDataProvider):
    """
    Provides data from a VCF file indexed via tabix.
    """

    dataset_type = "variant"


class RawVcfDataProvider(VcfDataProvider):
    """
    Provide data from VCF file.

    NOTE: this data provider does not use indices, and hence will be very slow
    for large datasets.
    """

    @contextmanager
    def open_data_file(self):
        with open(self.original_dataset.get_file_name()) as f:
            yield f

    def get_iterator(self, data_file, chrom, start, end, **kwargs) -> Iterator[str]:
        # Skip comments.
        line: Optional[str] = None
        for line in data_file:
            if not line.startswith("#"):
                break

        # If last line is a comment, there are no data lines.
        if not line or line.startswith("#"):
            return iter([])

        # Match chrom naming format.
        dataset_chrom = line.split()[0]
        if not _chrom_naming_matches(chrom, dataset_chrom):
            chrom = _convert_between_ucsc_and_ensemble_naming(chrom)

        def line_in_region(vcf_line, chrom, start, end):
            """Returns true if line is in region."""
            variant_chrom, variant_start = vcf_line.split()[0:2]
            # VCF format is 1-based.
            variant_start = int(variant_start) - 1
            return variant_chrom == chrom and variant_start >= start and variant_start <= end

        def line_filter_iter():
            """Yields lines in data that are in region chrom:start-end"""
            # Yield data line read above.
            if line_in_region(line, chrom, start, end):
                yield line

            # Search for and yield other data lines.
            for data_line in data_file:
                if line_in_region(data_line, chrom, start, end):
                    yield data_line

        return line_filter_iter()


class BamDataProvider(GenomeDataProvider, FilterableMixin):
    """
    Provides access to intervals from a sorted indexed BAM file. Coordinate
    data is reported in BED format: 0-based, half-open.
    """

    dataset_type = "bai"

    def get_filters(self):
        """
        Returns filters for dataset.
        """
        # HACK: first 7 fields are for drawing, so start filter column index at 7.
        filter_col = 7
        filters = []
        filters.append({"name": "Mapping Quality", "type": "number", "index": filter_col})
        return filters

    @contextmanager
    def open_data_file(self):
        # Attempt to open the BAM file with index
        with pysam.AlignmentFile(
            self.original_dataset.get_file_name(), mode="rb", index_filename=self.converted_dataset.get_file_name()
        ) as f:
            yield f

    def get_iterator(self, data_file, chrom, start, end, **kwargs) -> Iterator[str]:
        """
        Returns an iterator that provides data in the region chrom:start-end
        """

        # Fetch and return data.
        chrom = str(chrom)
        start = int(start)
        end = int(end)
        try:
            data = data_file.fetch(start=start, end=end, reference=chrom)
        except ValueError:
            # Try alternative chrom naming.
            chrom = _convert_between_ucsc_and_ensemble_naming(chrom)
            try:
                data = data_file.fetch(start=start, end=end, reference=chrom)
            except ValueError:
                return iter([])
        return data

    def process_data(
        self,
        iterator,
        start_val=0,
        max_vals=sys.maxsize,
        ref_seq=None,
        iterator_type="nth",
        mean_depth=None,
        start=0,
        end=0,
        **kwargs,
    ):
        """
        Returns a dict with the following attributes::

            data - a list of reads with the format
                [<guid>, <start>, <end>, <name>, <read_1>, <read_2>, [empty], <mapq_scores>]

                where <read_1> has the format
                    [<start>, <end>, <cigar>, <strand>, <read_seq>]

                and <read_2> has the format
                    [<start>, <end>, <cigar>, <strand>, <read_seq>]

                Field 7 is empty so that mapq scores' location matches that in single-end reads.
                For single-end reads, read has format:
                    [<guid>, <start>, <end>, <name>, <cigar>, <strand>, <seq>, <mapq_score>]

                NOTE: read end and sequence data are not valid for reads outside of
                requested region and should not be used.

            max_low - lowest coordinate for the returned reads
            max_high - highest coordinate for the returned reads
            message - error/informative message
        """
        # No iterator indicates no reads.
        if iterator is None:
            return {"data": [], "message": None}

        #
        # Helper functions.
        #

        def decode_strand(read_flag, mask):
            """Decode strand from read flag."""

            strand_flag = read_flag & mask == 0
            if strand_flag:
                return "+"
            else:
                return "-"

        def _random_read_iterator(read_iterator, threshold):
            """
            An iterator that returns a random stream of reads from the read_iterator
            as well as corresponding pairs for returned reads.
            threshold is a value in [0,1] that denotes the percentage of reads
            to return.
            """
            for e in read_iterator:
                if e.qname in paired_pending or random.uniform(0, 1) <= threshold:
                    yield e

        def _nth_read_iterator(read_iterator, threshold):
            """
            An iterator that returns every nth read.
            """

            # Convert threshold to N for stepping through iterator.
            n = int(1 / threshold)
            return itertools.islice(read_iterator, None, None, n)

        # -- Choose iterator. --

        # Calculate threshold for non-sequential iterators based on mean_depth and read length.
        try:
            first_read = next(iterator)
        except StopIteration:
            # no reads.
            return {"data": [], "message": None, "max_low": start, "max_high": start}

        read_len = len(first_read.seq)
        num_reads = max((end - start) * mean_depth / float(read_len), 1)
        threshold = float(max_vals) / num_reads
        iterator = itertools.chain(iter([first_read]), iterator)

        # Use specified iterator type, save for when threshold is >= 1.
        # A threshold of >= 1 indicates all reads are to be returned, so no
        # sampling needed and seqential iterator will be used.
        if iterator_type == "sequential" or threshold >= 1:
            read_iterator = iterator
        elif iterator_type == "random":
            read_iterator = _random_read_iterator(iterator, threshold)
        elif iterator_type == "nth":
            read_iterator = _nth_read_iterator(iterator, threshold)

        #
        # Encode reads as list of lists.
        #
        results = []
        paired_pending: dict[str, dict[str, Any]] = {}
        unmapped = 0
        message = None
        count = 0
        for read in read_iterator:
            if count < start_val:
                continue
            if (count - start_val - unmapped) >= max_vals:
                message = self.error_max_vals % (max_vals, "reads")
                break

            # If not mapped, skip read.
            is_mapped = read.flag & 0x0004 == 0
            if not is_mapped:
                unmapped += 1
                continue

            qname = read.qname
            seq = read.seq
            strand = decode_strand(read.flag, 0x0010)
            if read.cigar is not None:
                read_len = sum(cig[1] for cig in read.cigar)  # Use cigar to determine length
            else:
                read_len = len(seq)  # If no cigar, just use sequence length

            if read.is_proper_pair:
                if qname in paired_pending:
                    # Found pair.
                    pair = paired_pending[qname]
                    results.append(
                        [
                            hash("{}_{}".format(pair["start"], qname)),
                            pair["start"],
                            read.pos + read_len,
                            qname,
                            [pair["start"], pair["end"], pair["cigar"], pair["strand"], pair["seq"]],
                            [read.pos, read.pos + read_len, read.cigar, strand, seq],
                            None,
                            [pair["mapq"], read.mapq],
                        ]
                    )
                    del paired_pending[qname]
                else:
                    # Insert first of pair.
                    paired_pending[qname] = {
                        "start": read.pos,
                        "end": read.pos + read_len,
                        "seq": seq,
                        "mate_start": read.mpos,
                        "rlen": read_len,
                        "strand": strand,
                        "cigar": read.cigar,
                        "mapq": read.mapq,
                    }
                    count += 1
            else:
                results.append(
                    [
                        hash(f"{read.pos}_{qname}"),
                        read.pos,
                        read.pos + read_len,
                        qname,
                        read.cigar,
                        strand,
                        read.seq,
                        read.mapq,
                    ]
                )
                count += 1

        # Take care of reads whose mates are out of range.
        for qname, read in paired_pending.items():
            if read["mate_start"] < read["start"]:
                # Mate is before read.
                read_start = read["mate_start"]
                read_end = read["end"]
                # Make read_1 start=end so that length is 0 b/c we don't know
                # read length.
                r1 = [read["mate_start"], read["mate_start"]]
                r2 = [read["start"], read["end"], read["cigar"], read["strand"], read["seq"]]
            else:
                # Mate is after read.
                read_start = read["start"]
                # Make read_2 start=end so that length is 0 b/c we don't know
                # read length. Hence, end of read is start of read_2.
                read_end = read["mate_start"]
                r1 = [read["start"], read["end"], read["cigar"], read["strand"], read["seq"]]
                r2 = [read["mate_start"], read["mate_start"]]

            results.append([hash(f"{read_start}_{qname}"), read_start, read_end, qname, r1, r2, [read["mapq"], 125]])

        # Clean up. TODO: is this needed? If so, we'll need a cleanup function after processing the data.
        # bamfile.close()

        def compress_seq_and_cigar(read, start_field, cigar_field, seq_field):
            """
            Use reference-based compression to compress read sequence and cigar.
            """
            read_seq, read_cigar = get_ref_based_read_seq_and_cigar(
                read[seq_field].upper(), read[start_field], ref_seq.sequence, ref_seq.start, read[cigar_field]
            )
            read[seq_field] = read_seq
            read[cigar_field] = read_cigar

        def convert_cigar(read, start_field, cigar_field, seq_field):
            """
            Convert read cigar from pysam format to string format.
            """
            cigar_ops = "MIDNSHP=X"
            read_cigar = ""
            for op_tuple in read[cigar_field]:
                read_cigar += f"{op_tuple[1]}{cigar_ops[op_tuple[0]]}"
            read[cigar_field] = read_cigar

        # Choose method for processing reads. Use reference-based compression
        # if possible. Otherwise, convert cigar.
        if ref_seq:
            # Uppercase for easy comparison.
            ref_seq.sequence = ref_seq.sequence.upper()
            process_read = compress_seq_and_cigar
        else:
            process_read = convert_cigar

        # Process reads.
        for read in results:
            if isinstance(read[5], list):
                # Paired-end read.
                if len(read[4]) > 2:
                    process_read(read[4], 0, 2, 4)
                if len(read[5]) > 2:
                    process_read(read[5], 0, 2, 4)
            else:
                # Single-end read.
                process_read(read, 1, 4, 6)

        max_low, max_high = get_bounds(results, 1, 2)

        return {"data": results, "message": message, "max_low": max_low, "max_high": max_high}


class SamDataProvider(BamDataProvider):
    dataset_type = "bai"

    def __init__(self, converted_dataset=None, original_dataset=None, dependencies=None):
        """Create SamDataProvider."""
        super().__init__(
            converted_dataset=converted_dataset, original_dataset=original_dataset, dependencies=dependencies
        )

        # To use BamDataProvider, original dataset must be BAM and
        # converted dataset must be BAI. Use BAI from BAM metadata.
        if converted_dataset:
            self.original_dataset = converted_dataset
            self.converted_dataset = converted_dataset.metadata.bam_index


class BBIDataProvider(GenomeDataProvider):
    """
    BBI data provider for the Galaxy track browser.
    """

    dataset_type = "bigwig"

    @abc.abstractmethod
    def _get_dataset(self) -> tuple[IO[bytes], Union[BigBedFile, BigWigFile]]: ...

    def valid_chroms(self):
        # No way to return this info as of now
        return None

    def has_data(self, chrom):
        f, bbi = self._get_dataset()
        all_dat = bbi.query(chrom, 0, 2147483647, 1) or bbi.query(
            _convert_between_ucsc_and_ensemble_naming(chrom), 0, 2147483647, 1
        )
        f.close()
        return all_dat is not None

    def get_data(self, chrom: str, start, end, start_val=0, max_vals=sys.maxsize, **kwargs):
        start = int(start)
        end = int(end)

        # Helper function for getting summary data regardless of chromosome
        # naming convention.
        def _summarize_bbi(bbi, chrom, start, end, num_points):
            return bbi.summarize(chrom, start, end, num_points) or bbi.summarize(
                _convert_between_ucsc_and_ensemble_naming(chrom), start, end, num_points
            )

        # Bigwig can be a standalone bigwig file, in which case we use
        # original_dataset, or coming from wig->bigwig conversion in
        # which we use converted_dataset
        f, bbi = self._get_dataset()

        # If stats requested, compute overall summary data for the range
        # start:endbut no reduced data. This is currently used by client
        # to determine the default range.
        if "stats" in kwargs:
            summary = _summarize_bbi(bbi, chrom, start, end, 1)
            f.close()

            min_val = 0
            max_val = 0
            mean = 0
            sd = 0.0
            if summary is not None:
                # Does the summary contain any defined values?
                valid_count = summary.valid_count[0]
                if summary.valid_count > 0:
                    # Compute $\mu \pm 2\sigma$ to provide an estimate for upper and lower
                    # bounds that contain ~95% of the data.
                    mean = summary.sum_data[0] / valid_count
                    var = max(summary.sum_squares[0] - mean, 0)  # Prevent variance underflow.
                    if valid_count > 1:
                        var /= valid_count - 1
                    sd = math.sqrt(var)
                    min_val = summary.min_val[0]
                    max_val = summary.max_val[0]

            return dict(data=dict(min=min_val, max=max_val, mean=mean, sd=sd))

        def summarize_region(bbi, chrom, start, end, num_points):
            """
            Returns results from summarizing a region using num_points.
            NOTE: num_points cannot be greater than end - start or BBI
            will return None for all positions.
            """
            result = []

            # Get summary; this samples at intervals of length
            # (end - start)/num_points -- i.e. drops any fractional component
            # of interval length.
            if summary := _summarize_bbi(bbi, chrom, start, end, num_points):
                # mean = summary.sum_data / summary.valid_count

                # Standard deviation by bin, not yet used
                # var = summary.sum_squares - mean
                # var /= minimum( valid_count - 1, 1 )
                # sd = sqrt( var )

                pos = start
                step_size = (end - start) / num_points

                for i in range(num_points):
                    result.append((pos, float_nan(summary.sum_data[i] / summary.valid_count[i])))
                    pos += step_size

            return result

        # Approach is different depending on region size.
        num_samples = int(kwargs.get("num_samples", 100))
        if end - start < num_samples:
            # Get values for individual bases in region, including start and end.
            # To do this, need to increase end to next base and request number of points.
            num_points = end - start + 1
            end += 1
        else:
            #
            # The goal is to sample the region between start and end uniformly
            # using ~N (num_samples) data points. The challenge is that the size of
            # sampled intervals rarely is full bases, so sampling using N points
            # will leave the end of the region unsampled due to remainders for
            # each interval. To recitify this, a new N is calculated based on the
            # step size that covers as much of the region as possible.
            #
            # However, this still leaves some of the region unsampled. This
            # could be addressed by repeatedly sampling remainder using a
            # smaller and smaller step_size, but that would require iteratively
            # going to BBI, which could be time consuming.
            #

            # Start with N samples.
            num_points = num_samples
            step_size = (end - start) // num_points
            # Add additional points to sample in the remainder not covered by
            # the initial N samples.
            remainder_start = start + step_size * num_points
            additional_points = (end - remainder_start) // step_size
            num_points += additional_points

        result = summarize_region(bbi, chrom, start, end, num_points)

        # Cleanup and return.
        f.close()
        return {"data": result, "dataset_type": self.dataset_type}


class BigBedDataProvider(BBIDataProvider):
    def _get_dataset(self):
        # Nothing converts to bigBed so we don't consider converted dataset
        f = open(self.original_dataset.get_file_name(), "rb")
        return f, BigBedFile(file=f)


class BigWigDataProvider(BBIDataProvider):
    """
    Provides data from BigWig files; position data is reported in 1-based
    coordinate system, i.e. wiggle format.
    """

    def _get_dataset(self):
        if self.converted_dataset is not None:
            f = open(self.converted_dataset.get_file_name(), "rb")
        else:
            f = open(self.original_dataset.get_file_name(), "rb")
        return f, BigWigFile(file=f)


class IntervalIndexDataProvider(GenomeDataProvider, FilterableMixin):
    """
    Interval index files used for GFF, Pileup files.
    """

    col_name_data_attr_mapping = {4: {"index": 4, "name": "Score"}}

    dataset_type = "interval_index"

    @contextmanager
    def open_data_file(self):
        i = Indexes(self.converted_dataset.get_file_name())
        yield i

    def get_iterator(self, data_file, chrom, start, end, **kwargs) -> Iterator[str]:
        """
        Returns an iterator for data in data_file in chrom:start-end
        """
        if chrom not in data_file.indexes:
            # Try alternative naming.
            chrom = _convert_between_ucsc_and_ensemble_naming(chrom)

        return data_file.find(chrom, start, end)

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        results = []
        message = None
        with open(self.original_dataset.get_file_name()) as source:
            # Build data to return. Payload format is:
            # [ <guid/offset>, <start>, <end>, <name>, <score>, <strand>, <thick_start>,
            #   <thick_end>, <blocks> ]
            #
            # First three entries are mandatory, others are optional.
            filter_cols = loads(kwargs.get("filter_cols", "[]"))
            no_detail = "no_detail" in kwargs
            for count, val in enumerate(iterator):
                offset = val[2]
                if count < start_val:
                    continue
                if count - start_val >= max_vals:
                    message = self.error_max_vals % (max_vals, "features")
                    break
                source.seek(offset)
                # TODO: can we use column metadata to fill out payload?

                # GFF dataset.
                reader = GFFReaderWrapper(source, fix_strand=True)
                feature = next(reader)
                payload = package_gff_feature(feature, no_detail, filter_cols)
                payload.insert(0, offset)

                results.append(payload)

        return {"data": results, "message": message}


class RawGFFDataProvider(GenomeDataProvider):
    """
    Provide data from GFF file that has not been indexed.

    NOTE: this data provider does not use indices, and hence will be very slow
    for large datasets.
    """

    dataset_type = "interval_index"

    def get_iterator(self, data_file, chrom, start, end, **kwargs):
        """
        Returns an iterator that provides data in the region chrom:start-end as well as
        a file offset.
        """
        source = open(self.original_dataset.get_file_name())

        # Read first line in order to match chrom naming format.
        line = source.readline()

        # If line empty, assume file is empty and return empty iterator.
        if len(line) == 0:
            return iter([])

        # Determine chromosome naming format.
        dataset_chrom = line.split()[0]
        if not _chrom_naming_matches(chrom, dataset_chrom):
            chrom = _convert_between_ucsc_and_ensemble_naming(chrom)
        # Undo read.
        source.seek(0)

        def features_in_region_iter():
            offset = 0
            for feature in GFFReaderWrapper(source, fix_strand=True):
                # Only provide features that are in region.
                feature_start, feature_end = convert_gff_coords_to_bed([feature.start, feature.end])
                if feature.chrom == chrom and feature_end > start and feature_start < end:
                    yield feature, offset
                offset += feature.raw_size

        return features_in_region_iter()

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Process data from an iterator to a format that can be provided to client.
        """
        filter_cols = loads(kwargs.get("filter_cols", "[]"))
        no_detail = "no_detail" in kwargs
        results = []
        message = None

        for count, (feature, offset) in enumerate(iterator):
            if count < start_val:
                continue
            if count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "reads")
                break

            payload = package_gff_feature(feature, no_detail=no_detail, filter_cols=filter_cols)
            payload.insert(0, offset)
            results.append(payload)

        return {"data": results, "dataset_type": self.dataset_type, "message": message}


class GtfTabixDataProvider(TabixDataProvider):
    """
    Returns data from GTF datasets that are indexed via tabix.
    """

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        # Loop through lines and group by transcript_id; each group is a feature.

        # TODO: extend this code or use code in gff_util to process GFF/3 as well
        # and then create a generic GFFDataProvider that can be used with both
        # raw and tabix datasets.
        features: dict[str, list[GFFInterval]] = {}

        for line in iterator:
            line_attrs = parse_gff_attributes(line.split("\t")[8])
            transcript_id = line_attrs["transcript_id"]
            feature_list: list[GFFInterval]
            if transcript_id in features:
                feature_list = features[transcript_id]
            else:
                feature_list = []
                features[transcript_id] = feature_list
            feature_list.append(GFFInterval(None, line.split("\t")))

        # Process data.
        filter_cols = loads(kwargs.get("filter_cols", "[]"))
        no_detail = "no_detail" in kwargs
        results = []
        message = None

        for count, intervals in enumerate(features.values()):
            if count < start_val:
                continue
            if count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "reads")
                break

            feature = GFFFeature(None, intervals=intervals)
            payload = package_gff_feature(feature, no_detail=no_detail, filter_cols=filter_cols)
            payload.insert(0, feature.intervals[0].attributes["transcript_id"])
            results.append(payload)

        return {"data": results, "message": message}


#
# -- ENCODE Peak data providers.
#


class ENCODEPeakDataProvider(GenomeDataProvider):
    """
    Abstract class that processes ENCODEPeak data from native format to payload format.

    Payload format: [ uid (offset), start, end, name, strand, thick_start, thick_end, blocks ]
    """

    def get_iterator(self, data_file, chrom, start, end, **kwargs):
        raise Exception("Unimplemented Method")

    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Provides
        """

        # FIXMEs:
        # (1) should be able to unify some of this code with BedDataProvider.process_data
        # (2) are optional number of parameters supported?

        # Build data to return. Payload format is:
        # [ <guid/offset>, <start>, <end>, <name>, <strand>, <thick_start>,
        #   <thick_end>, <blocks> ]
        #
        # First three entries are mandatory, others are optional.
        #
        no_detail = "no_detail" in kwargs
        rval = []
        message = None
        for count, line in enumerate(iterator):
            if count < start_val:
                continue
            if max_vals and count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "features")
                break

            feature = line.split()

            # Feature initialization.
            payload: PAYLOAD_LIST_TYPE = [
                # GUID is just a hash of the line
                hash(line),
                # Add start, end.
                int(feature[1]),
                int(feature[2]),
            ]

            if no_detail:
                rval.append(payload)
                continue

            # Extend with additional data.
            payload.extend(
                [
                    # Add name, strand.
                    feature[3],
                    feature[5],
                    # Thick start, end are feature start, end for now.
                    int(feature[1]),
                    int(feature[2]),
                    # No blocks.
                    None,
                    # Filtering data: Score, signalValue, pValue, qValue.
                    float(feature[4]),
                    float(feature[6]),
                    float(feature[7]),
                    float(feature[8]),
                ]
            )

            rval.append(payload)

        return {"data": rval, "message": message}


class ENCODEPeakTabixDataProvider(TabixDataProvider, ENCODEPeakDataProvider):
    """
    Provides data from an ENCODEPeak dataset indexed via tabix.
    """

    def get_filters(self):
        """
        Returns filters for dataset.
        """
        # HACK: first 8 fields are for drawing, so start filter column index at 9.
        filter_col = 8
        filters = []
        filters.append(
            {"name": "Score", "type": "number", "index": filter_col, "tool_id": "Filter1", "tool_exp_name": "c6"}
        )
        filter_col += 1
        filters.append(
            {"name": "Signal Value", "type": "number", "index": filter_col, "tool_id": "Filter1", "tool_exp_name": "c7"}
        )
        filter_col += 1
        filters.append(
            {"name": "pValue", "type": "number", "index": filter_col, "tool_id": "Filter1", "tool_exp_name": "c8"}
        )
        filter_col += 1
        filters.append(
            {"name": "qValue", "type": "number", "index": filter_col, "tool_id": "Filter1", "tool_exp_name": "c9"}
        )
        return filters


#
# -- ChromatinInteraction data providers --
#


class ChromatinInteractionsDataProvider(GenomeDataProvider):
    def process_data(self, iterator, start_val=0, max_vals=sys.maxsize, **kwargs):
        """
        Provides
        """

        rval = []
        message = None
        for count, line in enumerate(iterator):
            if count < start_val:
                continue
            if max_vals and count - start_val >= max_vals:
                message = self.error_max_vals % (max_vals, "interactions")
                break

            feature = line.split()

            s1 = int(feature[1])
            e1 = int(feature[2])
            c = feature[3]
            s2 = int(feature[4])
            e2 = int(feature[5])
            v = float(feature[6])

            # Feature initialization.
            payload = [
                # GUID is just a hash of the line
                hash(line),
                # Add start1, end1, chr2, start2, end2, value.
                s1,
                e1,
                c,
                s2,
                e2,
                v,
            ]

            rval.append(payload)

        return {"data": rval, "message": message}

    def get_default_max_vals(self):
        return 100000


class ChromatinInteractionsTabixDataProvider(TabixDataProvider, ChromatinInteractionsDataProvider):
    def get_iterator(
        self, data_file, chrom, start=0, end=sys.maxsize, interchromosomal=False, **kwargs
    ) -> Iterator[str]:
        """ """
        # Modify start as needed to get earlier interactions with start region.
        span = int(end) - int(start)
        filter_start = max(0, int(start) - span - span / 2)

        def filter(iter):
            for line in iter:
                feature = line.split()
                s1 = int(feature[1])
                e1 = int(feature[2])
                c = feature[3]
                s2 = int(feature[4])
                e2 = int(feature[5])
                # Check for intrachromosal interactions.
                if ((s1 + s2) / 2 <= end) and ((e1 + e2) / 2 >= start) and (c == chrom):
                    yield line
                # Check for interchromosal interactions.
                if interchromosomal and c != chrom:
                    yield line

        return filter(TabixDataProvider.get_iterator(self, data_file, chrom, filter_start, end))


#
# -- Helper methods. --
#


def package_gff_feature(feature, no_detail=False, filter_cols=None) -> PAYLOAD_LIST_TYPE:
    """Package a GFF feature in an array for data providers."""
    filter_cols = filter_cols or []
    feature = convert_gff_coords_to_bed(feature)

    # No detail means only start, end.
    if no_detail:
        return [feature.start, feature.end]

    # Return full feature.
    payload = [
        feature.start,
        feature.end,
        feature.name(),
        feature.strand,
        # No notion of thick start, end in GFF, so make everything
        # thick.
        feature.start,
        feature.end,
    ]

    # Add blocks.
    block_sizes = [(interval.end - interval.start) for interval in feature.intervals]
    block_starts = [(interval.start - feature.start) for interval in feature.intervals]
    blocks = list(zip(block_sizes, block_starts))
    payload.append([(feature.start + block[1], feature.start + block[1] + block[0]) for block in blocks])

    # Add filter data to payload.
    for col in filter_cols:
        if col == "Score":
            if feature.score == "nan":
                payload.append(feature.score)
            else:
                try:
                    f = float(feature.score)
                    payload.append(f)
                except Exception:
                    payload.append(feature.score)
        elif col in feature.attributes:
            if feature.attributes[col] == "nan":
                payload.append(feature.attributes[col])
            else:
                try:
                    f = float(feature.attributes[col])
                    payload.append(f)
                except Exception:
                    payload.append(feature.attributes[col])
        else:
            # Dummy value.
            payload.append(0)
    return payload
