text:
# -*- coding: utf-8 -*-
"""This module provides alternative implementations of array and table
classes defined in the :mod:`allel.model.ndarray` module, using
chunked arrays for data storage. Chunked arrays can be compressed and
optionally stored on disk, providing a means for working with data too
large to fit uncompressed in main memory.

Either HDF5 (via `h5py <http://www.h5py.org/>`_) or `bcolz
<http://bcolz.blosc.org>`_ can be used as the underlying storage
layer. Choice of storage layer can be made via the `storage` keyword
argument which all class methods accept.  This argument can either be
a string identifying one of the predefined storage layer
configurations, or an object implementing the chunked storage API. For more
information about controlling storage see the :mod:`allel.chunked` module.

"""
from __future__ import absolute_import, print_function, division
import itertools

import numpy as np

from allel.compat import copy_method_doc, string_types
from allel.model import ndarray as _ndarray
from allel import chunked as _chunked
from allel.io import write_vcf_header, write_vcf_data, iter_gff3

__all__ = ['GenotypeChunkedArray', 'HaplotypeChunkedArray',
           'AlleleCountsChunkedArray', 'VariantChunkedTable',
           'FeatureChunkedTable', 'AlleleCountsChunkedTable']

class GenotypeChunkedArray(_chunked.ChunkedArray):
    """Alternative implementation of the
    :class:`allel.model.ndarray.GenotypeArray` class, wrapping a
    chunked array as the backing store.

    Parameters
    ----------
    data : array_like, int, shape (n_variants, n_samples, ploidy)
        Genotype data to be wrapped. May be a bcolz carray, h5py dataset, or
        anything providing a similar interface.

    Examples
    --------

    Wrap an HDF5 dataset::

        >>> import h5py
        >>> with h5py.File('callset.h5', mode='w') as h5f:
        ...     h5g = h5f.create_group('/3L/calldata')
        ...     h5g.create_dataset('genotype',
        ...                        data=[[[0, 0], [0, 1]],
        ...                              [[0, 1], [1, 1]],
        ...                              [[0, 2], [-1, -1]]],
        ...                        dtype='i1',
        ...                        chunks=(2, 2, 2))
        ...
        <HDF5 dataset "genotype": shape (3, 2, 2), type "|i1">
        >>> import allel
        >>> callset = h5py.File('callset.h5', mode='r')
        >>> g = allel.GenotypeChunkedArray(callset['/3L/calldata/genotype'])
        >>> g
        GenotypeChunkedArray((3, 2, 2), int8, nbytes=12, cbytes=16, cratio=0.8, shuffle=False, chunks=(2, 2, 2), data=h5py._hl.dataset.Dataset)
        >>> g.data
        <HDF5 dataset "genotype": shape (3, 2, 2), type "|i1">

    Obtain a numpy array by slicing, e.g.::

        >>> g[:]
        GenotypeArray((3, 2, 2), dtype=int8)
        [[[ 0  0]
          [ 0  1]]
         [[ 0  1]
          [ 1  1]]
         [[ 0  2]
          [-1 -1]]]

    Note that most methods will return a chunked array, using whatever
    chunked storage is set as default (bcolz carray) or specified
    directly via the `storage` keyword argument. E.g.::

        >>> g.copy()
        GenotypeChunkedArray((3, 2, 2), int8, nbytes=12, cbytes=16.0K, cratio=0.0, cname=blosclz, clevel=5, shuffle=True, chunks=(4096, 2, 2), data=bcolz.carray_ext.carray)
        >>> g.copy(storage='hdf5mem_zlib1')
        GenotypeChunkedArray((3, 2, 2), int8, nbytes=12, cbytes=4.5K, cratio=0.0, cname=gzip, clevel=1, shuffle=False, chunks=(262144, 2, 2), data=h5py._hl.dataset.Dataset)

    """  # flake8: noqa

    def __init__(self, data):
        super(GenotypeChunkedArray, self).__init__(data)
        self._check_input_data(self.data)
        self._mask = None

    @staticmethod
    def _check_input_data(data):
        if len(data.shape) != 3:
            raise ValueError('expected 3 dimensions')
        if data.dtype.kind not in 'ui':
            raise TypeError('expected integer dtype')

    def __getitem__(self, *args):
        out = super(GenotypeChunkedArray, self).__getitem__(*args)
        if hasattr(out, 'shape') \
                and len(self.shape) == len(out.shape) \
                and self.shape[2] == out.shape[2]:
            # dimensionality and ploidy preserved
            out = _ndarray.GenotypeArray(out)
            if self.mask is not None:
                # attempt to slice mask too
                m = self.mask.__getitem__(*args)
                out.mask = m
        return out

    def _repr_html_(self):
        return self[:6].to_html_str(caption=repr(self))

    @property
    def n_variants(self):
        return self.shape[0]

    @property
    def n_samples(self):
        return self.shape[1]

    @property
    def ploidy(self):
        return self.shape[2]

    @property
    def n_calls(self):
        return self.shape[0] * self.shape[1]

    @property
    def n_allele_calls(self):
        return self.shape[0] * self.shape[1] * self.shape[2]

    @property
    def mask(self):
        return self._mask

    @mask.setter
    def mask(self, mask):

        # check input
        if not hasattr(mask, 'shape'):
            mask = np.asarray(mask)
        if mask.shape != self.shape[:2]:
            raise ValueError('mask has incorrect shape')

        # store
        self._mask = _chunked.ChunkedArray(mask)

    def compress(self, condition, axis=0, **storage_kwargs):
        out = super(GenotypeChunkedArray, self).compress(condition, axis=axis,
                                                         **storage_kwargs)
        if self.mask is not None:
            out.mask = self.mask.compress(condition, axis=axis,
                                          **storage_kwargs)
        return out

    def take(self, indices, axis=0, **storage_kwargs):
        out = super(GenotypeChunkedArray, self).take(indices, axis=axis,
                                                     **storage_kwargs)
        if self.mask is not None:
            out.mask = self.mask.take(indices, axis=axis, **storage_kwargs)
        return out

    def subset(self, sel0, sel1, **storage_kwargs):
        out = super(GenotypeChunkedArray, self).subset(sel0, sel1,
                                                       **storage_kwargs)
        if self.mask is not None:
            out.mask = self.mask.subset(sel0, sel1, **storage_kwargs)
        return out

    def fill_masked(self, value=-1, **storage_kwargs):
        out = self.apply_method('fill_masked', kwargs=dict(value=value),
                                **storage_kwargs)
        return GenotypeChunkedArray(out)

    def is_called(self, **storage_kwargs):
        return self.apply_method('is_called', **storage_kwargs)

    def is_missing(self, **storage_kwargs):
        return self.apply_method('is_missing', **storage_kwargs)

    def is_hom(self, allele=None, **storage_kwargs):
        return self.apply_method('is_hom', kwargs=dict(allele=allele),
                                 **storage_kwargs)

    def is_hom_ref(self, **storage_kwargs):
        return self.apply_method('is_hom_ref', **storage_kwargs)

    def is_hom_alt(self, **storage_kwargs):
        return self.apply_method('is_hom_alt', **storage_kwargs)

    def is_het(self, allele=None, **storage_kwargs):
        return self.apply_method('is_het', kwargs=dict(allele=allele),
                                 **storage_kwargs)

    def is_call(self, call, **storage_kwargs):
        return self.apply_method('is_call', kwargs=dict(call=call),
                                 **storage_kwargs)

    def _count(self, method_name, axis, kwargs=None, **storage_kwargs):
        if kwargs is None:
            kwargs = dict()

        def mapper(block):
            method = getattr(block, method_name)
            return method(**kwargs)
        out = self.sum(axis=axis, mapper=mapper, **storage_kwargs)
        return out

    def count_called(self, axis=None, **storage_kwargs):
        return self._count('is_called', axis, **storage_kwargs)

    def count_missing(self, axis=None, **storage_kwargs):
        return self._count('is_missing', axis, **storage_kwargs)

    def count_hom(self, allele=None, axis=None, **storage_kwargs):
        return self._count('is_hom', axis, kwargs=dict(allele=allele),
                           **storage_kwargs)

    def count_hom_ref(self, axis=None, **storage_kwargs):
        return self._count('is_hom_ref', axis, **storage_kwargs)

    def count_hom_alt(self, axis=None, **storage_kwargs):
        return self._count('is_hom_alt', axis, **storage_kwargs)

    def count_het(self, allele=None, axis=None, **storage_kwargs):
        return self._count('is_het', axis, kwargs=dict(allele=allele),
                           **storage_kwargs)

    def count_call(self, call, axis=None, **storage_kwargs):
        return self._count('is_call', axis, kwargs=dict(call=call),
                           **storage_kwargs)

    def to_haplotypes(self, **storage_kwargs):
        out = self.apply_method('to_haplotypes', **storage_kwargs)
        return HaplotypeChunkedArray(out)

    def to_n_ref(self, fill=0, dtype='i1', **storage_kwargs):
        out = self.apply_method('to_n_ref', kwargs=dict(fill=fill,
                                                        dtype=dtype),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def to_n_alt(self, fill=0, dtype='i1', **storage_kwargs):
        out = self.apply_method('to_n_alt', kwargs=dict(fill=fill,
                                                        dtype=dtype),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def to_allele_counts(self, alleles=None, **storage_kwargs):
        # determine alleles to count
        if alleles is None:
            m = self.max()
            alleles = list(range(m+1))

        out = self.apply_method('to_allele_counts',
                                kwargs=dict(alleles=alleles),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def to_packed(self, boundscheck=True, **storage_kwargs):
        out = self.apply_method('to_packed',
                                kwargs=dict(boundscheck=boundscheck),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    @staticmethod
    def from_packed(packed, **storage_kwargs):
        def f(block):
            return _ndarray.GenotypeArray.from_packed(block)
        out = _chunked.apply(packed, f, **storage_kwargs)
        return GenotypeChunkedArray(out)

    def count_alleles(self, max_allele=None, subpop=None, **storage_kwargs):
        # if max_allele not specified, count all alleles
        if max_allele is None:
            max_allele = self.max()

        out = self.apply_method('count_alleles',
                                kwargs=dict(max_allele=max_allele,
                                            subpop=subpop),
                                **storage_kwargs)
        return AlleleCountsChunkedArray(out)

    def count_alleles_subpops(self, subpops, max_allele=None,
                              **storage_kwargs):
        if max_allele is None:
            max_allele = self.max()

        def f(block):
            return block.count_alleles_subpops(subpops, max_allele=max_allele)
        out = _chunked.apply(self, f, create='table', **storage_kwargs)
        return AlleleCountsChunkedTable(out)

    def to_gt(self, phased=False, max_allele=None, **storage_kwargs):
        out = self.apply_method('to_gt',
                                kwargs=dict(phased=phased,
                                            max_allele=max_allele),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def map_alleles(self, mapping, **storage_kwargs):
        def f(block, bmapping):
            return block.map_alleles(bmapping, copy=False)
        domain = (self, mapping)
        out = _chunked.apply(domain, f, **storage_kwargs)
        return GenotypeChunkedArray(out)

# copy docstrings
copy_method_doc(GenotypeChunkedArray.fill_masked,
                _ndarray.GenotypeArray.fill_masked)
copy_method_doc(GenotypeChunkedArray.subset,
                _ndarray.GenotypeArray.subset)
copy_method_doc(GenotypeChunkedArray.is_called,
                _ndarray.GenotypeArray.is_called)
copy_method_doc(GenotypeChunkedArray.is_missing,
                _ndarray.GenotypeArray.is_missing)
copy_method_doc(GenotypeChunkedArray.is_hom,
                _ndarray.GenotypeArray.is_hom)
copy_method_doc(GenotypeChunkedArray.is_hom_ref,
                _ndarray.GenotypeArray.is_hom_ref)
copy_method_doc(GenotypeChunkedArray.is_hom_alt,
                _ndarray.GenotypeArray.is_hom_alt)
copy_method_doc(GenotypeChunkedArray.is_het,
                _ndarray.GenotypeArray.is_het)
copy_method_doc(GenotypeChunkedArray.is_call,
                _ndarray.GenotypeArray.is_call)
copy_method_doc(GenotypeChunkedArray.to_haplotypes,
                _ndarray.GenotypeArray.to_haplotypes)
copy_method_doc(GenotypeChunkedArray.to_n_ref,
                _ndarray.GenotypeArray.to_n_ref)
copy_method_doc(GenotypeChunkedArray.to_n_alt,
                _ndarray.GenotypeArray.to_n_alt)
copy_method_doc(GenotypeChunkedArray.to_allele_counts,
                _ndarray.GenotypeArray.to_allele_counts)
copy_method_doc(GenotypeChunkedArray.to_packed,
                _ndarray.GenotypeArray.to_packed)
GenotypeChunkedArray.from_packed.__doc__ = \
    _ndarray.GenotypeArray.from_packed.__doc__
copy_method_doc(GenotypeChunkedArray.count_alleles,
                _ndarray.GenotypeArray.count_alleles)
copy_method_doc(GenotypeChunkedArray.count_alleles_subpops,
                _ndarray.GenotypeArray.count_alleles_subpops)
copy_method_doc(GenotypeChunkedArray.to_gt, _ndarray.GenotypeArray.to_gt)
copy_method_doc(GenotypeChunkedArray.map_alleles,
                _ndarray.GenotypeArray.map_alleles)
copy_method_doc(GenotypeChunkedArray.hstack, _ndarray.GenotypeArray.hstack)
copy_method_doc(GenotypeChunkedArray.vstack, _ndarray.GenotypeArray.vstack)

class HaplotypeChunkedArray(_chunked.ChunkedArray):
    """Alternative implementation of the
    :class:`allel.model.ndarray.HaplotypeArray` class, using a chunked array as
    the backing store.

    Parameters
    ----------
    data : array_like, int, shape (n_variants, n_haplotypes)
        Haplotype data to be wrapped. May be a bcolz carray, h5py dataset, or
        anything providing a similar interface.

    """

    def __init__(self, data):
        super(HaplotypeChunkedArray, self).__init__(data)
        self._check_input_data(self.data)

    @staticmethod
    def _check_input_data(data):
        if len(data.shape) != 2:
            raise ValueError('expected 2 dimensions')
        if data.dtype.kind not in 'ui':
            raise TypeError('expected integer dtype')

    def __getitem__(self, *args):
        out = super(HaplotypeChunkedArray, self).__getitem__(*args)
        if hasattr(out, 'shape') and len(self.shape) == len(out.shape):
            # dimensionality preserved
            out = _ndarray.HaplotypeArray(out)
        return out

    def _repr_html_(self):
        return self[:6].to_html_str(caption=repr(self))

    @property
    def n_variants(self):
        return self.shape[0]

    @property
    def n_haplotypes(self):
        return self.shape[1]

    def to_genotypes(self, ploidy=2, **storage_kwargs):

        # check ploidy is compatible
        if (self.n_haplotypes % ploidy) > 0:
            raise ValueError('incompatible ploidy')

        # build output
        def f(block):
            return block.to_genotypes(ploidy)

        out = self.apply(f, **storage_kwargs)
        return GenotypeChunkedArray(out)

    def is_called(self, **storage_kwargs):
        return self.__ge__(0, **storage_kwargs)

    def is_missing(self, **storage_kwargs):
        return self.__lt__(0, **storage_kwargs)

    def is_ref(self, **storage_kwargs):
        return self.__eq__(0, **storage_kwargs)

    def is_alt(self, **storage_kwargs):
        return self.__gt__(0, **storage_kwargs)

    def is_call(self, allele, **storage_kwargs):
        return self.__eq__(allele, **storage_kwargs)

    def _count(self, method_name, axis, kwargs=None, **storage_kwargs):
        if kwargs is None:
            kwargs = dict()

        def mapper(block):
            method = getattr(block, method_name)
            return method(**kwargs)
        out = self.sum(axis=axis, mapper=mapper, **storage_kwargs)
        return out

    def count_called(self, axis=None, **storage_kwargs):
        return self._count('is_called', axis=axis, **storage_kwargs)

    def count_missing(self, axis=None, **storage_kwargs):
        return self._count('is_missing', axis=axis, **storage_kwargs)

    def count_ref(self, axis=None, **storage_kwargs):
        return self._count('is_ref', axis=axis, **storage_kwargs)

    def count_alt(self, axis=None, **storage_kwargs):
        return self._count('is_alt', axis=axis, **storage_kwargs)

    def count_call(self, allele, axis=None, **storage_kwargs):
        return self._count('is_call', axis=axis,
                           kwargs=dict(allele=allele),
                           **storage_kwargs)

    def count_alleles(self, max_allele=None, subpop=None, **storage_kwargs):
        # if max_allele not specified, count all alleles
        if max_allele is None:
            max_allele = self.max()

        def f(block):
            return block.count_alleles(max_allele=max_allele, subpop=subpop)
        out = self.apply(f, **storage_kwargs)
        return AlleleCountsChunkedArray(out)

    def count_alleles_subpops(self, subpops, max_allele=None,
                              **storage_kwargs):
        if max_allele is None:
            max_allele = self.max()

        def f(block):
            return block.count_alleles_subpops(subpops, max_allele=max_allele)
        out = _chunked.apply(self, f, create='table', **storage_kwargs)
        return AlleleCountsChunkedTable(out)

    def map_alleles(self, mapping, **storage_kwargs):
        def f(block, bmapping):
            return block.map_alleles(bmapping, copy=False)
        domain = (self, mapping)
        out = _chunked.apply(domain, f, **storage_kwargs)
        return HaplotypeChunkedArray(out)

# copy docstrings
copy_method_doc(HaplotypeChunkedArray.to_genotypes,
                _ndarray.HaplotypeArray.to_genotypes)
copy_method_doc(HaplotypeChunkedArray.count_alleles,
                _ndarray.HaplotypeArray.count_alleles)
copy_method_doc(HaplotypeChunkedArray.count_alleles_subpops,
                _ndarray.HaplotypeArray.count_alleles_subpops)
copy_method_doc(HaplotypeChunkedArray.map_alleles,
                _ndarray.HaplotypeArray.map_alleles)

class AlleleCountsChunkedArray(_chunked.ChunkedArray):
    """Alternative implementation of the
    :class:`allel.model.ndarray.AlleleCountsArray` class, using a chunked
    array as the backing store.

    Parameters
    ----------
    data : array_like, int, shape (n_variants, n_alleles)
        Allele counts data to be wrapped. May be a bcolz carray,
        h5py dataset, or anything providing a similar interface.

    """

    def __init__(self, data):
        super(AlleleCountsChunkedArray, self).__init__(data)
        self._check_input_data(self.data)

    @staticmethod
    def _check_input_data(data):
        if len(data.shape) != 2:
            raise ValueError('expected 2 dimensions')
        if data.dtype.kind not in 'ui':
            raise TypeError('expected integer dtype')

    def __getitem__(self, *args):
        out = super(AlleleCountsChunkedArray, self).__getitem__(*args)
        if hasattr(out, 'shape') and len(self.shape) == len(out.shape) and \
                out.shape[1] == self.shape[1]:
            # dimensionality and allele indices preserved
            out = _ndarray.AlleleCountsArray(out)
        return out

    def _repr_html_(self):
        return self[:6].to_html_str(caption=repr(self))

    @property
    def n_variants(self):
        return self.shape[0]

    @property
    def n_alleles(self):
        return self.shape[1]

    def to_frequencies(self, fill=np.nan, **storage_kwargs):
        out = self.apply_method('to_frequencies',
                                kwargs=dict(fill=fill),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def allelism(self, **storage_kwargs):
        out = self.apply_method('allelism', **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def max_allele(self, **storage_kwargs):
        out = self.apply_method('max_allele', **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def is_variant(self, **storage_kwargs):
        out = self.apply_method('is_variant', **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def is_non_variant(self, **storage_kwargs):
        out = self.apply_method('is_non_variant', **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def is_segregating(self, **storage_kwargs):
        out = self.apply_method('is_segregating', **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def is_non_segregating(self, allele=None, **storage_kwargs):
        out = self.apply_method('is_non_segregating',
                                kwargs=dict(allele=allele),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def is_singleton(self, allele=1, **storage_kwargs):
        out = self.apply_method('is_singleton',
                                kwargs=dict(allele=allele),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def is_doubleton(self, allele=1, **storage_kwargs):
        out = self.apply_method('is_doubleton',
                                kwargs=dict(allele=allele),
                                **storage_kwargs)
        return _chunked.ChunkedArray(out)

    def _count(self, method_name, kwargs=None, **storage_kwargs):
        if kwargs is None:
            kwargs = dict()

        def mapper(block):
            method = getattr(block, method_name)
            return method(**kwargs)
        out = self.sum(mapper=mapper, **storage_kwargs)
        return out

    def count_variant(self, **storage_kwargs):
        return self._count('is_variant', **storage_kwargs)

    def count_non_variant(self, **storage_kwargs):
        return self._count('is_non_variant', **storage_kwargs)

    def count_segregating(self, **storage_kwargs):
        return self._count('is_segregating', **storage_kwargs)

    def count_non_segregating(self, allele=None, **storage_kwargs):
        return self._count('is_non_segregating', kwargs=dict(allele=allele),
                           **storage_kwargs)

    def count_singleton(self, allele=1, **storage_kwargs):
        return self._count('is_singleton', kwargs=dict(allele=allele),
                           **storage_kwargs)

    def count_doubleton(self, allele=1, **storage_kwargs):
        return self._count('is_doubleton', kwargs=dict(allele=allele),
                           **storage_kwargs)

    def map_alleles(self, mapping, **storage_kwargs):
        def f(block, bmapping):
            return block.map_alleles(bmapping)
        domain = (self, mapping)
        out = _chunked.apply(domain, f, **storage_kwargs)
        return AlleleCountsChunkedArray(out)

copy_method_doc(AlleleCountsChunkedArray.allelism,
                _ndarray.AlleleCountsArray.allelism)
copy_method_doc(AlleleCountsChunkedArray.max_allele,
                _ndarray.AlleleCountsArray.max_allele)
copy_method_doc(AlleleCountsChunkedArray.map_alleles,
                _ndarray.AlleleCountsArray.map_alleles)

class VariantChunkedTable(_chunked.ChunkedTable):
    """Alternative implementation of the
    :class:`allel.model.ndarray.VariantTable` class, using a chunked table as
    the backing store.

    Parameters
    ----------
    data: table_like
        Data to be wrapped. May be a tuple or list of columns (array-like),
        a dict mapping names to columns, a bcolz ctable, h5py group,
        numpy recarray, or anything providing a similar interface.
    names : sequence of strings
        Column names.

    Examples
    --------

    Wrap columns stored as datasets within an HDF5 group::

        >>> import h5py
        >>> chrom = [b'chr1', b'chr1', b'chr2', b'chr2', b'chr3']
        >>> pos = [2, 7, 3, 9, 6]
        >>> dp = [35, 12, 78, 22, 99]
        >>> qd = [4.5, 6.7, 1.2, 4.4, 2.8]
        >>> ac = [(1, 2), (3, 4), (5, 6), (7, 8), (9, 10)]
        >>> with h5py.File('callset.h5', mode='w') as h5f:
        ...     h5g = h5f.create_group('/3L/variants')
        ...     h5g.create_dataset('CHROM', data=chrom, chunks=True)
        ...     h5g.create_dataset('POS', data=pos, chunks=True)
        ...     h5g.create_dataset('DP', data=dp, chunks=True)
        ...     h5g.create_dataset('QD', data=qd, chunks=True)
        ...     h5g.create_dataset('AC', data=ac, chunks=True)
        ...
        <HDF5 dataset "CHROM": shape (5,), type "|S4">
        <HDF5 dataset "POS": shape (5,), type "<i8">
        <HDF5 dataset "DP": shape (5,), type "<i8">
        <HDF5 dataset "QD": shape (5,), type "<f8">
        <HDF5 dataset "AC": shape (5, 2), type "<i8">
        >>> import allel
        >>> callset = h5py.File('callset.h5', mode='r')
        >>> vt = allel.VariantChunkedTable(callset['/3L/variants'],
        ...                                names=['CHROM', 'POS', 'AC', 'QD', 'DP'])
        >>> vt
        VariantChunkedTable(5, nbytes=220, cbytes=220, cratio=1.0, data=h5py._hl.group.Group)

    Obtain a single row::

        >>> vt[0]
        row(CHROM=b'chr1', POS=2, AC=array([1, 2]), QD=4.5, DP=35)

    Obtain a numpy array by slicing::

        >>> vt[:] # doctest: +ELLIPSIS
        VariantTable((5,), dtype=[('CHROM', 'S4'), ('POS', '<i8'), ('AC', ...
        [(b'chr1', 2, [1, 2], 4.5, 35) (b'chr1', 7, [3, 4], 6.7, 12)
         (b'chr2', 3, [5, 6], 1.2, 78) (b'chr2', 9, [7, 8], 4.4, 22)
         (b'chr3', 6, [9, 10], 2.8, 99)]

    Access a subset of columns::

        >>> vt[['CHROM', 'POS']]
        VariantChunkedTable(5, nbytes=60, cbytes=60, cratio=1.0, data=builtins.list)

    Note that most methods will return a chunked table, using whatever
    chunked storage is set as default (bcolz ctable) or specified
    directly via the `storage` keyword argument. E.g.::

        >>> vt.copy()
        VariantChunkedTable(5, nbytes=220, cbytes=80.0K, cratio=0.0, data=bcolz.ctable.ctable)
        >>> vt.copy(storage='hdf5mem_zlib1')
        VariantChunkedTable(5, nbytes=220, cbytes=22.5K, cratio=0.0, data=h5py._hl.files.File)

    """  # flake8: noqa

    view_cls = _ndarray.VariantTable

    def __init__(self, data, names=None, index=None):
        super(VariantChunkedTable, self).__init__(data, names=names)
        self.index = None
        if index is not None:
            self.set_index(index)

    @property
    def n_variants(self):
        return len(self)

    def set_index(self, spec):
        if isinstance(spec, string_types):
            self.index = _ndarray.SortedIndex(self[spec][:], copy=False)
        elif isinstance(spec, (tuple, list)) and len(spec) == 2:
            self.index = _ndarray.SortedMultiIndex(self[spec[0]][:],
                                                   self[spec[1]][:],
                                                   copy=False)
        else:
            raise ValueError('invalid index argument, expected string or '
                             'pair of strings, found %s' % repr(spec))

    def to_vcf(self, path, rename=None, number=None, description=None,
               fill=None, blen=None, write_header=True):
        with open(path, 'w') as vcf_file:
            if write_header:
                write_vcf_header(vcf_file, self, rename=rename, number=number,
                                 description=description)
            blen = _chunked.get_blen_table(self)
            for i in range(0, len(self), blen):
                j = min(i+blen, len(self))
                block = self[i:j]
                write_vcf_data(vcf_file, block, rename=rename, fill=fill)

class FeatureChunkedTable(_chunked.ChunkedTable):
    """Alternative implementation of the
    :class:`allel.model.ndarray.FeatureTable` class, using a chunked table as
    the backing store.

    Parameters
    ----------
    data: table_like
        Data to be wrapped. May be a tuple or list of columns (array-like),
        a dict mapping names to columns, a bcolz ctable, h5py group,
        numpy recarray, or anything providing a similar interface.
    names : sequence of strings
        Column names.

    """

    view_cls = _ndarray.FeatureTable

    @property
    def n_features(self):
        return len(self)

    def to_mask(self, size, start_name='start', stop_name='end'):
        m = np.zeros(size, dtype=bool)
        start = self[start_name]
        stop = self[stop_name]
        for i, j in zip(start, stop):
            # assume 1-based inclusive coords
            m[i-1:j] = True
        return m

    @staticmethod
    def from_gff3(path, attributes=None, region=None, score_fill=-1,
                  phase_fill=-1, attributes_fill=b'.', dtype=None,
                  blen=None, storage=None, create='table', expectedlen=200000,
                  **storage_kwargs):

        # setup iterator
        recs = iter_gff3(path, attributes=attributes, region=region,
                         score_fill=score_fill, phase_fill=phase_fill,
                         attributes_fill=attributes_fill)

        # read a sample to determine dtype, blen
        recs_sample = list(itertools.islice(recs, 1000))
        names = 'seqid', 'source', 'type', 'start', 'end', 'score', 'strand', \
                'phase'
        if attributes:
            names += tuple(attributes)
        ra = np.rec.array(recs_sample, names=names, dtype=dtype)
        dtype = ra.dtype

        # setup output
        storage = _chunked.get_storage(storage)
        out = getattr(storage, create)(ra, expectedlen=expectedlen,
                                       **storage_kwargs)
        blen = _chunked.get_blen_table(out, blen=blen)

        # read block-wise
        block = list(itertools.islice(recs, 0, blen))
        while block:
            a = np.asarray(block, dtype=dtype)
            out.append(a)
            block = list(itertools.islice(recs, 0, blen))

        out = FeatureChunkedTable(out)
        return out

class AlleleCountsChunkedTable(_chunked.ChunkedTable):

    def __getitem__(self, item):
        out = super(AlleleCountsChunkedTable, self).__getitem__(item)
        if isinstance(item, string_types):
            # rewrap
            out = AlleleCountsChunkedArray(out.data)
        return out

max_stars_repo_name: podpearson/scikit-allel

max_issues_repo_name: podpearson/scikit-allel

max_forks_repo_name: podpearson/scikit-allel


text:
from cmath import cos, sin
import numpy as np
import scipy.integrate
import math
from scipy.spatial.transform import Rotation
from numpy.linalg import inv

class Quadrotor_linear():
    def __init__(self):
        """
        Define model parameters and state space for linear quadrotor dynamics
        """

        self.mass            =  0.90  # kg
        self.Ixx             = 30*1.43e-5  # kg*m^2
        self.Iyy             = 30*1.43e-5  # kg*m^2
        self.Izz             = 30*2.89e-5  # kg*m^2
        self.arm_length      = 0.046  # meters
        g = 9.81 # m/s^2
        L = self.arm_length
        self.to_TM = np.array([[1,  1,  1,  1],
                               [ 0,  L,  0, -L],
                               [-L,  0,  L,  0]])
        self.t_step = 0.1

        """
        Continuous state space without considering yawing(assuming yaw angle is 0 for all the time)
        state = [x y z dx dy dz phi theta phi_dot theta_dot]
        dot_state = [dx dy dz ddx ddy ddz phi_dot theta_dot phi_ddot theta_ddot]
        u = [F1 F2 F3 F4]
        """
        #                     x  y  z  dx dy dz phi theta phi_dot theta_dot
        self.A_c = np.array([[0, 0, 0, 1, 0, 0, 0,  0,    0,      0        ],    #dx
                             [0, 0, 0, 0, 1, 0, 0,  0,    0,      0        ],    #dy
                             [0, 0, 0, 0, 0, 1, 0,  0,    0,      0        ],    #dz
                             [0, 0, 0, 0, 0, 0, 0,  g,    0,      0        ],    #ddx
                             [0, 0, 0, 0, 0, 0, -g, 0,    0,      0        ],    #ddy
                             [0, 0, 0, 0, 0, 0, 0,  0,    0,      0        ],    #ddz
                             [0, 0, 0, 0, 0, 0, 0,  0,    1,      0        ],    #phi_dot
                             [0, 0, 0, 0, 0, 0, 0,  0,    0,      1        ],    #theta_dot
                             [0, 0, 0, 0, 0, 0, 0,  0,    0,      0        ],    #phi_ddot
                             [0, 0, 0, 0, 0, 0, 0,  0,    0,      0        ]])   #theta_ddot
        self.B_c = np.array([[0, 0, 0], #dx
                             [0, 0, 0], #dy
                             [0, 0, 0], #dz
                             [0, 0, 0], #ddx
                             [0, 0, 0], #ddy
                             [1/self.mass, 0, 0], #ddz
                             [0, 0, 0],           #phi_dot
                             [0, 0, 0],           #theta_dot
                             [0, 1/self.Ixx, 0],  #phi_ddot
                             [0, 0, 1/self.Iyy]  #theta_ddot
                             ]) @ self.to_TM
        self.C_c = np.identity(10)
        self.D_c = np.zeros((1,4))

        # Discretization state space
        self.A = np.eye(10) + self.A_c * self.t_step
        self.B = self.B_c * self.t_step
        self.C = self.C_c
        self.D = self.D_c

    def next_x(self, x, u):
        return self.A.dot(x).reshape(-1,1) + self.B.dot(u)

    def disturbed_next_x(self,x,u,real_d,Bd):
        return self.A @ x.reshape(-1,1) + self.B @ u.reshape(-1,1) + Bd @ real_d.reshape(-1,1)

    def disturbed_output(self,x,real_d, Cd, sigma_noise):
        return self.C @ x.reshape(-1,1) + Cd @ real_d.reshape(-1,1) + np.random.normal(loc=np.zeros((1,10)),scale=sigma_noise).reshape(10,1)

    def from_nonlinear(self, quadrotor):
        linearized_state = np.zeros(10).reshape(-1,1)
        linearized_state[:8] = quadrotor.state[:8].reshape(-1,1)
        linearized_state[8:] = quadrotor.state[9:11].reshape(-1,1)
        return linearized_state

### Non linear quadrotor model. Wasnt used in the end ###
class Quadrotor():
    """
    Quadrotor forward dynamics model.
    """
    def __init__(self):
        self.mass            =  0.030  # kg
        self.Ixx             = 1.43e-5  # kg*m^2
        self.Iyy             = 1.43e-5  # kg*m^2
        self.Izz             = 2.89e-5  # kg*m^2
        self.arm_length      = 0.046  # meters

        #TODO: change force limit
        self.cmd_rotor_forces_min = -0.7  # N
        self.cmd_rotor_forces_max = 0.7  # N
        self.k_thrust        = 2.3e-08  # N/(rad/s)**2
        self.k_drag          = 7.8e-11   # Nm/(rad/s)**2
        self.inertia = np.diag(np.array([self.Ixx, self.Iyy, self.Izz])) # kg*m^2
        self.g = 9.81 # m/s^2
        # Precomputes
        k = self.k_drag/self.k_thrust
        L = self.arm_length
        self.to_TM = np.array([[1,  1,  1,  1],
                               [ 0,  L,  0, -L],
                               [-L,  0,  L,  0],
                               [ k, -k,  k, -k]])
        self.inv_inertia = np.linalg.inv(self.inertia)
        #self.weight = np.array([0, 0, -self.mass*self.g])
        self.weight = np.array([0, 0, 0])
        self.t_step = 0.05
        # self.inv_inertia = inv(self.inertia)
        # self.weight = np.array([0, 0, 0])
        # self.t_step = 0.05
    def reset(self, position=[0, 0, 0], yaw =0, pitch=0, roll=0):
        '''
        state is a 12 dimensional vector
            postion*3 velocity*3 attitude(rpy)*3 angular velocity*3
        state = [x y z dx dy dz psi theta phi p q r]
        dot_state = [dx dy dz ddx ddy ddz dqw dqx dqy dqz dr dp dq]
        '''
        s = np.zeros(12)
        s[0] = position[0]
        s[1] = position[1]
        s[2] = position[2]
        s[6] = pitch
        s[7] = roll
        s[8] = yaw
        self.state = s
        return self.state

    def step(self, cmd_rotor_forces):
        '''
        Considering the max and min of rotor forces
        action is a 4 dimensional vector: conmmand rotor forces
        action = [F1, F2, F3, F4]
        '''
        rotor_thrusts = np.clip(cmd_rotor_forces, self.cmd_rotor_forces_min, self.cmd_rotor_forces_max)
        TM = self.to_TM @ cmd_rotor_forces
        T = TM[0]
        M = TM[1:4].flatten()
        # Form autonomous ODE for constant inputs and integrate one time step.
        def s_dot_fn(t, s):
            return self._s_dot_fn(t, s, T, M)
        '''
        The next state can be obtained through integration （Runge-Kutta）
        '''
        s = self.state
        sol = scipy.integrate.solve_ivp(s_dot_fn, (0, self.t_step), s, first_step=self.t_step)
        self.state = sol['y'][:,-1]
        return self.state

    def _s_dot_fn(self, t, s, u1, u2):
        """
        Compute derivative of state for quadrotor given fixed control inputs as
        an autonomous ODE.
        """
        state = self.state
        x_dot = state[3:6]
        v_dot = np.array([u1[0]/self.mass * (cos(state[8]) * sin(state[7]) + cos(state[7]) * sin(state[6]) * sin(state[8])),
                          u1[0]/self.mass * (sin(state[8]) * sin(state[7]) - cos(state[8]) * cos(state[7]) * sin(state[6])),
                          u1[0]/self.mass * (cos(state[6]) * cos(state[7]))])
        q_dot = state[9:]
        q_dot_hat = Quadrotor.hat_map(q_dot)
        w_dot = self.inv_inertia @ (u2 - q_dot_hat @ (self.inertia @ q_dot))

        # Pack into vector of derivatives.
        s_dot = np.zeros((12,))
        s_dot[0:3] = x_dot
        s_dot[3:6] = v_dot
        s_dot[6:9] = q_dot
        s_dot[9:] = w_dot
        return s_dot
    @classmethod
    def hat_map(cls, s):
        """
        Given vector s in R^3, return associate skew symmetric matrix S in R^3x3
        """
        return np.array([[0, -s[2], s[1]],
                         [s[2], 0, -s[0]],
                         [-s[1], s[0], 0]])

ext: py

max_stars_repo_name: smoggy-P/MPC-Collision-Avoidance

max_issues_repo_name: smoggy-P/MPC-Collision-Avoidance

max_forks_repo_name: smoggy-P/MPC-Collision-Avoidance

reason: import numpy,from numpy,import scipy,from scipy


text:
import torch
import torch.nn.functional as F
import torch.nn as nn
import numpy as np
import random
import math
from matplotlib import pyplot as plt
from sklearn import metrics
from torch.utils.data import DataLoader
from dataset.torch_dataset import TorchDataset
from dataset.tumor_dataset import TumorDataset
from metrics.evaluate_cls import evaluate_multi_cls
from models.gcn_conv import GCNConv
import os.path as osp
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.nn import Parameter
from tqdm import tqdm
import pandas as pd
from sklearn.metrics import confusion_matrix

# 2 layer mlp
# translayer -> adj
# 2 layer gcn
# linear pred

class CNN_Model(nn.Module):
    def __init__(self, output_dim):
        super(CNN_Model, self).__init__()
        # 32x32x3 image
        self.conv1 = nn.Conv2d(3, 64, 3, padding=1)
        # 16x16x16
        self.conv2 = nn.Conv2d(64, 128, 3, padding=1)
        # 8x8x32
        self.conv3 = nn.Conv2d(128, 256, 3, padding=1)
        ## max pool
        self.pool = nn.MaxPool2d(2, 2)
        self.bn1 = nn.BatchNorm2d(64)
        self.bn2 = nn.BatchNorm2d(128)
        self.bn3 = nn.BatchNorm2d(256)
        # linear layer (64x4x4->500)

        self.fc1 = nn.Linear(64 * 16 * 16, 500)
        #        self.bn = nn.BatchNorm1d(500)
        #        self.fc1 = nn.Linear(3*128*128,500)
        # linear layer (500->10)
        self.fc2 = nn.Linear(500, output_dim)
        # dropout
        self.dropout = nn.Dropout(0.5)

    def forward(self, x):
        # add sequence of convolutional and max pooling layers

        x = x.permute(0, 3, 1, 2)
        x = self.pool(F.relu(self.conv1(x)))
        #        x = self.conv1(x)
        #        x = self.pool(x)
        ##        x = self.pool(x)
        #        x = self.bn1(x)
        #
        #        x = F.leaky_relu(x)
        #        x = self.pool(F.relu(self.conv2(x)))
        #        x = self.conv2(x)
        #        x = self.pool(x)
        #        x = self.bn2(x)
        #        x = F.leaky_relu(x)
        #        x = self.pool(F.relu(self.conv3(x)))
        #        x = self.conv3(x)
        #        x = self.pool(x)
        #        x = self.bn3(x)
        #        x = F.leaky_relu(x)

        x = x.reshape(-1, 64 * 16 * 16)

        x = self.dropout(x)
        # add 1st hidden layer, with relu activation function
        x = F.leaky_relu(self.fc1(x))
        #        x = self.bn(x)

        x = self.fc2(x)
        return x

class Model(nn.Module):
    def __init__(self, **param_dict):
        super(Model, self).__init__()
        self.param_dict = param_dict
        self.input_dim = self.param_dict['ft_dim']
        self.out_dim = self.param_dict['label_num']
        self.h_dim = self.param_dict['h_dim']
        self.dropout_num = self.param_dict['dropout_num']
        self.add_res = self.param_dict['add_res']
        self.img_input_dim = self.param_dict['img_input_dim']
        self.img_output_dim = self.param_dict['img_output_dim']

        self.linear1 = nn.Linear(self.input_dim, self.h_dim)
        self.linear2 = nn.Linear(self.h_dim, self.h_dim)

        self.adj_trans_linear = nn.Linear(self.h_dim, self.h_dim)
        self.gcn_layer1 = GCNConv(self.h_dim, self.h_dim)
        self.gcn_layer2 = GCNConv(self.h_dim, self.h_dim)

        self.linear_pred = nn.Linear(self.h_dim + self.img_output_dim, self.out_dim)
        self.activation = nn.ELU()
        self.dropout_layer = nn.Dropout(p=self.dropout_num)
        self.CNN = CNN_Model(self.img_output_dim)

    def forward(self, node_ft, img):

        res_mat = torch.zeros(node_ft.size()[0], self.h_dim).to(node_ft.device)

        node_ft = self.activation(self.linear1(node_ft))
        node_ft = self.dropout_layer(node_ft)

        node_ft = self.activation(self.linear2(node_ft))
        node_ft = self.dropout_layer(node_ft)

        res_mat += node_ft

        # adj
        trans_adj_ft = self.adj_trans_linear(node_ft)
        trans_adj_ft = torch.tanh(trans_adj_ft)
        w = torch.norm(trans_adj_ft, p=2, dim=-1).view(-1, 1)
        w_mat = w * w.t()
        adj = torch.mm(trans_adj_ft, trans_adj_ft.t()) / w_mat

        node_ft = self.activation(self.gcn_layer1(node_ft, adj))
        node_ft = self.dropout_layer(node_ft)
        res_mat += node_ft

        node_ft = self.activation(self.gcn_layer2(node_ft, adj))
        node_ft = self.dropout_layer(node_ft)
        res_mat += node_ft
        image = self.CNN(img)
        if self.add_res:
            res_mat = torch.cat((res_mat, image), 1)
            node_embedding = res_mat
            pred = self.linear_pred(res_mat)
        else:
            node_embedding = node_ft
            pred = self.linear_pred(node_ft)
        pred = F.log_softmax(pred, dim=-1)
        return pred, adj, node_embedding

model_save_dir = 'save_model_param_new'
current_path = osp.dirname(osp.realpath(__file__))

class Trainer(object):
    def __init__(self, **param_dict):
        self.param_dict = param_dict
        self.setup_seed(self.param_dict['seed'])
        self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
        self.dataset = TumorDataset(
            train_split_param=self.param_dict['train_split_param'],
            ft_stand=self.param_dict['ft_stand']
        )
        self.dataset.generate_dataset_info()
        self.param_dict.update(self.dataset.dataset_info)
        # self.dataset.to_tensor(self.device)
        self.file_name = __file__.split('/')[-1].replace('.py', '')
        self.trainer_info = '{}_seed={}_batch={}'.format(self.file_name, self.param_dict['seed'],
                                                         self.param_dict['batch_size'])
        # self.save_model_path = osp.join(current_path, model_save_dir, self.trainer_info)
        self.loss_op = torch.nn.NLLLoss()
        self.build_model()

    def build_model(self):
        self.model = Model(**self.param_dict).to(self.device)
        print(self.model)
        self.optimizer = torch.optim.Adam(self.model.parameters(), lr=self.param_dict['lr'])
        self.best_res = None
        self.min_dif = -1e10

    def setup_seed(self, seed):
        torch.manual_seed(seed)
        torch.cuda.manual_seed_all(seed)
        np.random.seed(seed)
        random.seed(seed)
        torch.backends.cudnn.deterministic = True

    def iteration(self, epoch, dataloader, is_training=False):
        if is_training:
            self.model.train()
        else:
            self.model.eval()

        all_pred = []
        all_label = []
        all_loss = []
        node_embedding = None
        for ft_mat, img, label_mat in dataloader:
            ft_mat = ft_mat.cuda().float()
            label_mat = label_mat.cuda().long()
            img = img.cuda().float()
            pred, adj, node_embedding0 = self.model(ft_mat, img)
            if is_training:
                # print(pred.size(), label_mat.size())
                c_loss = self.loss_op(pred, label_mat)
                param_l2_loss = 0
                param_l1_loss = 0
                for name, param in self.model.named_parameters():
                    if 'bias' not in name:
                        param_l2_loss += torch.norm(param, p=2)
                        param_l1_loss += torch.norm(param, p=1)

                param_l2_loss = self.param_dict['param_l2_coef'] * param_l2_loss
                adj_l1_loss = self.param_dict['adj_loss_coef'] * torch.norm(adj)
                loss = c_loss + param_l2_loss + adj_l1_loss
                # print('c_loss = ', c_loss.detach().to('cpu').item(),
                #       ' adj_l1_loss = ', adj_l1_loss.detach().to('cpu').item(),
                #       ' param_l2_loss = ', param_l2_loss.detach().to('cpu').item()
                #       )
                loss.backward()
                self.optimizer.step()
                self.optimizer.zero_grad()
                all_loss.append(loss.detach().to('cpu').item())

            max_value, max_pos = torch.max(pred, dim=1)
            pred = max_pos.detach().to('cpu').numpy()
            label_mat = label_mat.detach().to('cpu').numpy()
            all_pred = np.hstack([all_pred, pred])
            all_label = np.hstack([all_label, label_mat])
            node_embedding = node_embedding0

        return all_pred, all_label, all_loss, node_embedding

    def print_res(self, res_list, epoch):
        train_acc, valid_acc, test_primary_acc, test_transfer_acc, \
        train_macro_f1, valid_macro_f1, test_primary_macro_f1, test_transfer_macro_f1, \
        train_micro_p, valid_micro_p, test_primary_micro_p, test_transfer_micro_p, \
        train_micro_r, valid_micro_r, test_primary_micro_r, test_transfer_micro_r = res_list
        """
        msg_log = 'Epoch: {:03d}, Acc Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f}, Test transfer: {:.4f} ' \
                  'Macro F1 Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f}, Test transfer: {:.4f} ' \
                  'Micro P Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f}, Test transfer: {:.4f} ' \
                  'Micro R Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f}, Test transfer: {:.4f} '.format(
            epoch, train_acc, valid_acc, test_primary_acc, test_transfer_acc, \
            train_macro_f1, valid_macro_f1, test_primary_macro_f1, test_transfer_macro_f1, \
            train_micro_p, valid_micro_p, test_primary_micro_p, test_transfer_micro_p, \
            train_micro_r, valid_micro_r, test_primary_micro_r, test_transfer_micro_r)
        """
        msg_log = 'Epoch: {:03d}, Acc Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f} ' \
                  'Macro F1 Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f} ' \
                  'Micro P Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f} ' \
                  'Micro R Train: {:.4f}, Val: {:.4f}, Test primary: {:.4f} '.format(
            epoch, train_acc, valid_acc, test_primary_acc, \
            train_macro_f1, valid_macro_f1, test_primary_macro_f1, \
            train_micro_p, valid_micro_p, test_primary_micro_p, \
            train_micro_r, valid_micro_r, test_primary_micro_r)
        print(msg_log)

    def start(self, display=True):
        train_dataset = TorchDataset(dataset=self.dataset, split_type='train')
        train_dataloader = DataLoader(train_dataset, batch_size=self.param_dict['batch_size'], shuffle=True,
                                      num_workers=2)

        valid_dataset = TorchDataset(self.dataset, split_type='valid')
        valid_dataloader = DataLoader(valid_dataset, batch_size=self.param_dict['batch_size'], shuffle=True,
                                      num_workers=2)

        test_primary_dataset = TorchDataset(self.dataset, split_type='test_is_primary')
        test_primary_dataloader = DataLoader(test_primary_dataset, batch_size=self.param_dict['batch_size'],
                                             shuffle=False, num_workers=2)

        test_transfer_dataset = TorchDataset(self.dataset, split_type='test_is_transfer')
        test_transfer_dataloader = DataLoader(test_transfer_dataset, batch_size=self.param_dict['batch_size'],
                                              shuffle=False, num_workers=2)

        for epoch in range(1, self.param_dict['epoch_num'] + 1):
            train_pred, train_label, train_loss, node_embedding = self.iteration(epoch=epoch,
                                                                                dataloader=train_dataloader,
                                                                                is_training=True)
            train_acc, train_micro_f1, train_macro_f1, train_micro_p, train_micro_r = evaluate_multi_cls(train_pred,
                                                                                train_label)

            valid_pred, valid_label, valid_loss, node_embedding = self.iteration(epoch=epoch,
                                                                                dataloader=valid_dataloader,
                                                                                is_training=False)
            valid_acc, valid_micro_f1, valid_macro_f1, valid_micro_p, valid_micro_r = evaluate_multi_cls(valid_pred,
                                                                                valid_label)

            test_primary_pred, test_primary_label, test_primary_loss, node_embedding = self.iteration(epoch=epoch,
                                                                                dataloader=test_primary_dataloader,
                                                                                is_training=False)
            test_primary_acc, test_primary_micro_f1, test_primary_macro_f1, test_primary_micro_p, test_primary_micro_r \
                = evaluate_multi_cls(test_primary_pred, test_primary_label)

            test_transfer_pred, test_transfer_label, test_transfer_loss, node_embedding = \
                self.iteration(epoch=epoch, dataloader=test_transfer_dataloader, is_training=False)
            test_transfer_acc, test_transfer_micro_f1, test_transfer_macro_f1, test_transfer_micro_p, test_transfer_micro_r \
                = evaluate_multi_cls(test_transfer_pred, test_transfer_label)

            res_list = [
                train_acc, valid_acc, test_primary_acc, test_transfer_acc,
                train_macro_f1, valid_macro_f1, test_primary_macro_f1, test_transfer_macro_f1,
                train_micro_p, valid_micro_p, test_primary_micro_p, test_transfer_micro_p,
                train_micro_r, valid_micro_r, test_primary_micro_r, test_transfer_micro_r
            ]

            if valid_acc > self.min_dif:
                self.min_dif = valid_acc
                self.best_res = res_list
                self.best_epoch = epoch
                # save model
                save_complete_model_path = osp.join(current_path, model_save_dir, self.trainer_info + '_complete.pth')
                torch.save(self.model, save_complete_model_path)
                same_model_param_path = osp.join(current_path, model_save_dir, self.trainer_info + '_param.pkl')
                torch.save(self.model.state_dict(), same_model_param_path)

            if display:
                self.print_res(res_list, epoch)

            if epoch % 50 == 0 and epoch > 0:
                print('Best res')
                self.print_res(self.best_res, self.best_epoch)

    def whole_graph_evaluate(self):
        # load_params
        model_param_path = osp.join(current_path, model_save_dir, self.trainer_info + '_param.pkl')
        print('load params ', model_param_path)
        self.model.load_state_dict(torch.load(model_param_path))
        train_dataset = TorchDataset(dataset=self.dataset, split_type='train')
        train_dataloader = DataLoader(train_dataset, num_workers=2, batch_size=self.param_dict['batch_size'],
                                      shuffle=False)

        valid_dataset = TorchDataset(self.dataset, split_type='valid')
        valid_dataloader = DataLoader(valid_dataset, num_workers=2, batch_size=self.param_dict['batch_size'],
                                      shuffle=False)

        test_primary_dataset = TorchDataset(self.dataset, split_type='test_is_primary')
        test_primary_dataloader = DataLoader(test_primary_dataset, num_workers=2,
                                             batch_size=self.param_dict['batch_size'], shuffle=False)

        test_transfer_dataset = TorchDataset(self.dataset, split_type='test_is_transfer')
        test_transfer_dataloader = DataLoader(test_transfer_dataset, num_workers=2,
                                              batch_size=self.param_dict['batch_size'], shuffle=False)

        self.model.eval()

        train_pred, train_label, train_loss, train_embedding = self.iteration(epoch=1, dataloader=train_dataloader,
                                                                              is_training=True)
        train_acc, train_micro_f1, train_macro_f1, train_micro_p, train_micro_r = evaluate_multi_cls(train_pred,
                                                                                train_label)

        valid_pred, valid_label, valid_loss, valid_embedding = self.iteration(epoch=1, dataloader=valid_dataloader,
                                                                              is_training=False)
        valid_acc, valid_micro_f1, valid_macro_f1, valid_micro_p, valid_micro_r = evaluate_multi_cls(valid_pred,
                                                                                valid_label)

        test_primary_pred, test_primary_label, test_primary_loss, test_primary_embedding = self.iteration(epoch=1,
                                                                                dataloader=test_primary_dataloader,
                                                                                is_training=False)
        test_primary_acc, test_primary_micro_f1, test_primary_macro_f1, test_primary_micro_p, test_primary_micro_r \
            = evaluate_multi_cls(test_primary_pred, test_primary_label)

        test_transfer_pred, test_transfer_label, test_transfer_loss, test_transfer_embedding = \
            self.iteration(epoch=1, dataloader=test_transfer_dataloader, is_training=False)
        test_transfer_acc, test_transfer_micro_f1, test_transfer_macro_f1, test_transfer_micro_p, test_transfer_micro_r \
            = evaluate_multi_cls(test_transfer_pred, test_transfer_label)

        node_embedding = test_primary_embedding.detach().to('cpu').numpy()

        np.save(
            osp.join(current_path, 'save_embedding', self.trainer_info + '_embedding.npy'),
            node_embedding
        )

        all_pred = train_pred
        all_label = train_label
        print('all_pred = ', all_pred)
        print('all_label = ', all_label)

        cm = confusion_matrix(all_label, all_pred)
        print("confusion_matrix: ", cm)
        np.savetxt('cm_new.csv', cm, delimiter = ',')

        #        all_pred_prob, adj, node_embedding = self.model(self.dataset.ft_mat)
        #        max_value, all_pred = torch.max(all_pred_prob, dim=1)

        #        node_embedding = node_embedding.detach().to('cpu').numpy()
        #
        #        all_pred = all_pred.detach().to('cpu').numpy()
        #        all_label = self.dataset.label_mat.detach().to('cpu').numpy()
        #        print('all_pred = ', all_pred)
        #        print('all_label = ', all_label)

        # train_acc, train_micro_f1, train_macro_f1, train_micro_p, train_micro_r = evaluate_multi_cls(train_pred,
        #                                                                                train_label)
        #
        # valid_acc, valid_micro_f1, valid_macro_f1, valid_micro_p, valid_micro_r = evaluate_multi_cls(valid_pred,
        #                                                                                valid_label)
        #
        # test_primary_acc, test_primary_micro_f1, test_primary_macro_f1, test_primary_micro_p, test_primary_micro_r \
        #     = evaluate_multi_cls(test_primary_pred, test_primary_label)
        #
        # test_transfer_acc, test_transfer_micro_f1, test_transfer_macro_f1, test_transfer_micro_p, test_transfer_micro_r \
        #     = evaluate_multi_cls(test_transfer_pred, test_transfer_label)

        res_list = [
            train_acc, valid_acc, test_primary_acc, test_transfer_acc,
            train_macro_f1, valid_macro_f1, test_primary_macro_f1, test_transfer_macro_f1,
            train_micro_p, valid_micro_p, test_primary_micro_p, test_transfer_micro_p,
            train_micro_r, valid_micro_r, test_primary_micro_r, test_transfer_micro_r
        ]
        print('whole_graph_evaluate')
        self.print_res(res_list, 0)

if __name__ == '__main__':
    for seed in [2, 3, 5]:
        param_dict = {
            'seed': seed,
            'train_split_param': [0.95, 0.05],
            'ft_stand': False,
            'dropout_num': 0.3,
            'layer_num': 4,
            'epoch_num': 200,
            'lr': 1e-4,
            'param_l2_coef': 1e-2,
            'batch_size': 1024,
            'h_dim': 512,
            'adj_loss_coef': 1e-2,
            'add_res': True,
            'img_input_dim': 1024,
            'img_output_dim': 1,
        }
        trainer = Trainer(**param_dict)
        trainer.start()
        # trainer.whole_graph_evaluate()

reason: import numpy


text:
import numpy as np
from sklearn.datasets import load_digits
from sklearn.svm import SVC
#loaing digit images
digit = load_digits()
#only feature data
training_data = digit.data
#only target data
training_target = digit.target
#training data extract from original data
td_original = np.delete(training_data,-1,axis=0)
#training target extract from original data
tt_original = np.delete(training_target,-1)
#now calling support vector machine
clf = SVC()
#Now time for training algo
trained = clf.fit(td_original,tt_original)
#Now time forv predicting data
output = trained.predict(digit.data[-1].reshape(1,64))
print(output)
#plotting that testing image
import matplotlib.pyplot as plt
plt.imshow(digit.images[-1])
plt.show()


text:
import pandas as pa
import numpy as np
import matplotlib.pyplot as plt

datset=pa.read_csv("Ads_CTR_Optimisation.csv")
#datset contions a matrix of 0s and 1s, which are nothing but the rewards corresponding to each option at each round

#%%implementing the ucb reinforcement method

import math

N=10000         #no of rounds
d=10            #no of options available
                #n is the ierator of the number of round
                #i is the iterator of the number of option
number_of_selections= [0]*d         #vector of lenght =options containing the number of times
                                    #each of the option is selected
sum_of_rewards=[0]*d                #vector containing the sum of rewards for each of the options available
ads_selected=[]                     #contains the option nmber selected at each round,
                                    #i.e at the end of the algo, it will be a list of 10000 elements
total_reward=0                      #contains the total reward provided at the end of the algo(i.e. 10000 rounds)

                                    #as initially each of the option is selected 0 times and also the sum of the rewards for each of the option is 0
                                    #both the vectors are initialised to 0
for n in range(0,N):            #for the first o=10 rounds, at each round o, oth option is selected,providing the algo some initial innfo to start wwith
    upper_bound_list=[]         #contains the ucb for each option at each round
    for i in range(0,d):
        if number_of_selections[i]>0:           #since we have no info in the beginning, each option has to be selected atleast ones to provide the algo with some info
            average_reward=sum_of_rewards[i]/number_of_selections[i]
            delta_i=math.sqrt((3/2)*(math.log(n+1)/number_of_selections[i]))
            upper_bound=average_reward+delta_i
            upper_bound_list.append(upper_bound)
        else:                                       #if th round is no tselected even ones
            upper_bound=1e400                #upper bound to a very large value so thaat this option gets selected(value = 10^400........ math.pow(10,400) could not work due to range limit)
            upper_bound_list.append(upper_bound)
    x=upper_bound_list.index(max(upper_bound_list))     #find the index of option with the max ucb
    ads_selected.append(x)                              #append the selected option to the list
    number_of_selections[x]=number_of_selections[x]+1   #update the selection count for the selectedd option
    reward=datset.values[n,x]                                  #select the reward for the selected opt(index x) from the dataset
    sum_of_rewards[x]=sum_of_rewards[x]+reward          #update the reward vector for the sum of reward of each option
    total_reward=total_reward+reward

#%%
plt.hist(ads_selected)
plt.show()

max_stars_repo_path: other ML algos/ucb.py

max_issues_repo_path: other ML algos/ucb.py

max_forks_repo_path: other ML algos/ucb.py


text:
import numpy as np
from .base import BaseLayer
import logging

logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)

class Dense(BaseLayer):
	name = "Dense"

	def forward(self, A_m1):
		# * Using 'np.dot()' or '@' ?
		#	https://stackoverflow.com/questions/34142485/difference-between-numpy-dot-and-python-3-5-matrix-multiplication
		logger.debug(f"{self.__class__.__name__}: forward()")
		logger.debug(f"\t{self.params['W'].shape = }")
		logger.debug(f"\t{A_m1.shape = }")

		self.cache['Z'] = np.dot(self.params['W'], A_m1) + self.params['b']

		logger.debug(f"\t{self.cache['Z'].shape = }")

		self.cache['A'] = self.g.forward(self.cache['Z'])

		logger.debug(f"\t{self.cache['A'].shape = }")

		# logging
		# logger.debug(f"\t{self.cache['A'] = }")
		return self.cache['A']

	def backward(self, params_p1, grads_p1, A_m1):
		logger.debug(f"{self.__class__.__name__}: Backward()")
		logger.debug(f"\t{self.n_units =}")

		logger.debug(f"\t{grads_p1['dZ'].shape = }")
		logger.debug(f"\t{params_p1['W'].T.shape = }")
		logger.debug(f"\t{self.cache['Z'].shape = }")
		logger.debug(f"\t{self.g.backward(self.cache['Z']).shape = }")

		self.grads['dZ'] = np.dot(params_p1['W'].T, grads_p1['dZ']) * self.g.backward(self.cache['Z'])

		m = self.n_units
		self.grads['dW'] = (1 / m) * np.dot(self.grads['dZ'], A_m1.T)
		self.grads['db'] = (1 / m) * np.sum(self.grads['dZ'], axis=1, keepdims=True)

		logger.debug(f"\tdZ: {self.grads['dZ'].shape}")

		return self.grads['dZ']

if __name__ == "__main__":
	np.random.seed(seed=42)

	prev = 2
	curr = 4
	nb_examples = 10
	X = np.random.randn(prev, nb_examples)

	print(f"Prev layers size: {prev}")
	print(f"Curr layers size: {curr}")
	print(f"X values:\n{X}")
	print()

	layer = Dense(curr, prev, debug=True)
	print()
	print(layer)
	print()
	Y = layer.forward(X)
	print(f"Y= {Y}")
	print()
	dA = layer.backward(X, Y)
	print(f"dA= {dA}")

max_stars_repo_name: ezalos/NotYetSelfAware

max_issues_repo_name: ezalos/NotYetSelfAware

max_forks_repo_name: ezalos/NotYetSelfAware


text:
import networkx as nx
import random

from tools import safe_sample, log
from nature import Regulon

MIN_P_INSERT, MAX_P_INSERT = 0.1618, 0.42

def insert_motif(G, id):
    M, new_ids = Regulon(parent=id)
    G = nx.compose(G, M)
    predecessors, successors = G.predecessors(id), new_ids[0]
    for predecessor in predecessors:
        for successor in successors:
            G.add_edge(predecessor, successor)
    predecessors, successors = new_ids[-1], list(G.successors(id))
    for predecessor in predecessors:
        for successor in successors:
            G.add_edge(predecessor, successor)
    G.remove_node(id)
    return G

def insert_motifs(G):
    p_insert = safe_sample(MIN_P_INSERT, MAX_P_INSERT)
    nodes = G.nodes(data=True)
    for node in nodes:
        try:
            if node[1]["shape"] is "square" and random.random() < p_insert:
                insert_motif(G, id=node[0])
        except Exception as e:
            log('insert_motifs Exception', e)
    return G

max_stars_repo_name: bionicles/neuromax

max_issues_repo_name: bionicles/neuromax

max_forks_repo_name: bionicles/neuromax

reason: import networkx


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Oct  8 12:55:31 2020

@author: xxx
"""

import pickle
import os
import numpy as np
from collections import Counter

import shutil
import re

EGCN_strategies = {
'xyz':'Single Modal Position',
'ang':'Single Modal Orientation',
'xyzang':'MultiModal SLE',
'cat':'MultiModal FLE-1',
'lit':'MultiModal FLE-2',
'add':'MultiModal DLE-1',
'add_1':'MultiModal DLE-2',
'mul':'MultiModal MLE'
}

ui_prmd_exercises = {
1:'E1, ',
2:'E2, ',
3:'E3, ',
4:'E4, ',
5:'E5, ',
6:'E6, ',
7:'E7, ',
8:'E8, ',
9:'E9, ',
10:'E10,'
}

eval_criterias = ['cv_cs','cv_rd']
for eval_criteria in eval_criterias:

    label_path = './data/UI_PRMD/'+ eval_criteria +'/xyz'
    for index, method in enumerate(['xyz','ang','xyzang','cat','lit','add','add_1','mul']):# 'xyz','ang','xyzang','cat','lit','add','add_1','mul'
        print('\nThe cross-validation ('+eval_criteria+') results of ' + EGCN_strategies[method] + ':')
        result_path = './work_dir/ui_prmd/'+ eval_criteria +'/' + method

        if method != 'add_1':
            folds = ['1', '2', '3', '4', '5']
            actions = [1,2,3,4,5,6,7,8,9,10]
            for action in actions:
                labels = []
                results = []
                fold_acc = []
                for fold in folds:
                    label = open(os.path.join(label_path,str(action),fold, 'eval_label.pkl'), 'rb')

                    label = np.array(pickle.load(label))
                    label = label.T.tolist()
                    if not os.path.exists(os.path.join(result_path,str(action),fold, 'result.pkl')):
                        fold_acc.append(0)
                        #print('no result: ', os.path.join(result_path,str(action),fold))
                        continue
                    r1 = open(os.path.join(result_path,str(action),fold, 'result.pkl'), 'rb')

                    r1 = list(pickle.load(r1).items())
                    right_num_11 = 0
                    for i in range(len(label)):
                        l = label[i][1]
                        _, r11 = r1[i]

                        r11 = np.argmax(r11)
                        right_num_11 += int(r11 == int(l))

                    acc = right_num_11/len(label)
                    fold_acc.append(acc)

                    labels.extend(label)
                    results.extend(r1)

                right_num_action = 0
                for i in range(len(labels)):
                    l = labels[i][1]
                    _, r_a = results[i]

                    r_a = np.argmax(r_a)
                    right_num_action += int(r_a == int(l))

                acc_action = right_num_action/len(labels)
                print("Exercise "+ui_prmd_exercises[action]+" f1: {:0.4f}, f2: {:0.4f}, f3: {:0.4f}, f4: {:0.4f}, f5: {:0.4f}, all:{:0.4f}".format(fold_acc[0],fold_acc[1],fold_acc[2],fold_acc[3],fold_acc[4],acc_action))

        else:
            result_path_xyz = './work_dir/ui_prmd/'+ eval_criteria +'/xyz'
            result_path_ang = './work_dir/ui_prmd/'+ eval_criteria +'/ang'

            folds = ['1', '2', '3', '4', '5']
            actions = [1,2,3,4,5,6,7,8,9,10]
            for action in actions:
                labels = []
                results_1 = []
                results_2 = []
                fold_acc = []
                for fold in folds:
                    label = open(os.path.join(label_path,str(action),fold, 'eval_label.pkl'), 'rb')
                    label = np.array(pickle.load(label))
                    label = label.T.tolist()
                    if not os.path.exists(os.path.join(result_path_xyz,str(action),fold, 'result.pkl')):
                        fold_acc.append(0)
                        #print('no result: ', os.path.join(result_path,str(action),fold))
                        continue
                    r1 = open(os.path.join(result_path_xyz,str(action),fold, 'result.pkl'), 'rb')
                    r2 = open(os.path.join(result_path_ang,str(action),fold, 'result.pkl'), 'rb')
                    r1 = list(pickle.load(r1).items())
                    r2 = list(pickle.load(r2).items())
                    right_num_11 = 0
                    for i in range(len(label)):
                        l = label[i][1]
                        _, r11 = r1[i]
                        _, r22 = r2[i]

                        r11 = r11 + r22
                        r11 = np.argmax(r11)
                        right_num_11 += int(r11 == int(l))

                    acc = right_num_11/len(label)
                    fold_acc.append(acc)

                    labels.extend(label)
                    results_1.extend(r1)
                    results_2.extend(r2)

                right_num_action = 0
                for i in range(len(labels)):
                    l = labels[i][1]
                    _, r_a_1 = results_1[i]
                    _, r_a_2 = results_2[i]
                    r_a_1 = r_a_1 + r_a_2
                    r_a_1 = np.argmax(r_a_1)
                    right_num_action += int(r_a_1 == int(l))

                acc_action = right_num_action/len(labels)
                print("Exercise "+ui_prmd_exercises[action]+" f1: {:0.4f}, f2: {:0.4f}, f3: {:0.4f}, f4: {:0.4f}, f5: {:0.4f}, all:{:0.4f}".format(fold_acc[0],fold_acc[1],fold_acc[2],fold_acc[3],fold_acc[4],acc_action))

max_stars_repo_name: bruceyo/EGCN

max_issues_repo_name: bruceyo/EGCN

max_forks_repo_name: bruceyo/EGCN


text:
import os
import sys
import numpy as np
import statistics as stats
from astropy.io import fits
from sklearn.neighbors import NearestNeighbors
from matplotlib.collections import PatchCollection
import matplotlib.pyplot as plt
import matplotlib as mpl

def plot_Bins(Bins,x_min,x_max,y_min,y_max,StN_Target,file_dir,filename):
    fig = plt.figure()
    #fig.set_size_inches(7, 6.5)
    ax = plt.axes(xlim=(x_min,x_max), ylim=(y_min,y_max))
    N = len(Bins)
    StN_list = []
    SNR_list = []
    bin_nums = []
    max_StN = max([bin.StN[0] for bin in Bins])
    StN_list = [bin.StN[0] for bin in Bins]
    median_StN = np.median(StN_list)
    stand_dev = stats.stdev(StN_list)
    mini_pallete = ['mediumspringgreen','salmon','cyan','orchid','yellow','blue','red','magenta','black','white']
    binNumber  = 0
    for bin in Bins:
        bin_nums.append(bin.bin_number)
        SNR = bin.StN[0]/median_StN
        SNR_list.append(SNR)
        for pixel in bin.pixels:
            x_coord = pixel.pix_x
            y_coord = pixel.pix_y
            #patches.append(Rectangle((x_coord,y_coord),1,1))
            if binNumber%10 == 0:
                color = mini_pallete[0]
            if binNumber%10 == 1:
                color = mini_pallete[1]
            if binNumber%10 == 2:
                color = mini_pallete[2]
            if binNumber%10 == 3:
                color = mini_pallete[3]
            if binNumber%10 == 4:
                color = mini_pallete[4]
            if binNumber%10 == 5:
                color = mini_pallete[5]
            if binNumber%10 == 6:
                color = mini_pallete[6]
            if binNumber%10 == 7:
                color = mini_pallete[7]
            if binNumber%10 == 8:
                color = mini_pallete[8]
            if binNumber%10 == 9:
                color = mini_pallete[9]
            #Shift because x_coord,y_coord are the center points
            rectangle = plt.Rectangle((x_coord,y_coord),1,1, fc=color)
            ax.add_patch(rectangle)
        binNumber += 1
    centroids_x = [Bins[i].centroidx[0] for i in range(len(Bins))]
    centroids_y = [Bins[i].centroidy[0] for i in range(len(Bins))]
    ax.scatter(centroids_x,centroids_y,marker='+',c="black")
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.title("Bin Mosaic")
    plt.savefig(file_dir+'/'+filename+".png")
    plt.clf()
    plt.hist(StN_list)
    plt.xlim((median_StN-3*stand_dev,median_StN+3*stand_dev))
    plt.ylabel("Number of Bins")
    plt.xlabel("Signal-to-Noise")
    plt.title("Signal-to-Noise per Bin")
    plt.savefig(file_dir+'/histograms/'+filename+".png")
    plt.clf()
    SNR_std = stats.stdev(SNR_list)
    SNR_med = np.median(SNR_list)
    SNR_nel = len(SNR_list)
    n_el = SNR_nel #Just change this
    plt.scatter(np.arange(n_el),SNR_list, marker = '+', color='salmon',label="Data Points")
    plt.plot(np.arange(n_el),[SNR_med for i in range(n_el)], linestyle='solid', color= 'forestgreen', label="Median")
    plt.plot(np.arange(n_el),[SNR_med+SNR_std for i in range(n_el)], linestyle='--', color= 'black')
    plt.plot(np.arange(n_el),[SNR_med-SNR_std for i in range(n_el)], linestyle='--', color= 'black', label="sigma")
    plt.title("Signal to Noise Ratio")
    plt.ylim((min(SNR_list),max(SNR_list)))
    plt.xlabel("Bin Number")
    plt.ylabel("Signal-to-Noise Normalized by Median Value")
    plt.legend(loc='center left',bbox_to_anchor=(1.0, 0.5),
          ncol=1, fancybox=True, shadow=True)
    plt.savefig(file_dir+'/'+filename+"_scatter.png", bbox_inches="tight")
    plt.clf()
    return None

class Bin:
    """
    Bin Class Information
    Everything in here should be self-explanatory... if not let me know and I
    will most happily comment it! :)
    """
    def __init__(self, number):
        self.bin_number = number
        self.pixels = []
        self.pixel_neighbors = []
        self.centroidx = [0] #So I can pass by value
        self.centroidy = [0]
        self.centroidx_prev = [0]
        self.centroidy_prev = [0]
        self.StN = [0]
        self.StN_prev = [0]
        #self.Signal = [0]
        #self.Noise = [0]
        self.Area = [0]
        self.Area_prev = [0]
        self.scale_length = [0]
        self.scale_length_prev = [0]
        self.successful = False
        self.WVT_successful = False
        self.avail_reassign = True #Can a pixel be reassigned to you?

    def add_pixel(self, Pixel):
        #self.StN[0] = 0
        self.pixels.append(Pixel)
        self.StN[0] += Pixel.StN
        #print(self.StN[0])
        #self.Signal[0] += Pixel.Signal
        #self.Noise[0] += Pixel.Noise
        #if self.Noise[0] != 0:
        #    self.StN[0] = self.Signal[0]/(np.sqrt(self.Noise[0]))
        for neigh in Pixel.neighbors:
            if (neigh not in self.pixel_neighbors):
                self.pixel_neighbors.append(neigh)

    def clear_pixels(self):
        self.pixels = []
        self.StN[0] = 0
        #self.Signal[0] = 0
        #self.Noise[0] = 0
        self.successful = False
        self.WVT_successful = False
        self.avail_reassign = True

    def remove_pixel(self, Pixel):
        self.pixels.remove(Pixel)
        self.StN[0] -= Pixel.StN
        #self.Signal[0] -= Pixel.Signal
        #self.Noise[0] -= Pixel.Noise
        #if self.Noise[0] != 0:
        #    self.StN[0] = self.Signal[0]/(np.sqrt(self.Noise[0]))

    def success(self):
        self.successful = True

    def WVT_success(self):
        self.WVT_successful = True

    def availabilty(self):
        self.avail_reassign = False

    def update_StN_prev(self):
        self.StN_prev[0] = self.StN[0]

    def CalcCentroid(self):
        self.centroidx_prev[0] = self.centroidx[0]
        self.centroidy_prev[0] = self.centroidy[0]
        self.centroidx[0] = 0
        self.centroidy[0] = 0
        n_cent = 0
        for pixel in self.pixels:
            self.centroidx[0] += pixel.pix_x
            self.centroidy[0] += pixel.pix_y
            n_cent += 1
        self.centroidx[0] *= (1/n_cent)
        self.centroidy[0] *= (1/n_cent)

    def CalcArea(self,pixel_length):
        self.Area_prev[0] = self.Area[0]
        self.Area[0] = 0
        self.Area[0] += pixel_length**2*len(self.pixels)

    def CalcScaleLength(self,StN_Target):
        self.scale_length_prev[0] = self.scale_length[0]
        self.scale_length[0] = 0
        self.scale_length[0] += np.sqrt((self.Area[0]/np.pi)*(StN_Target/self.StN[0]))

class Pixel:
    def __init__(self, number, pix_x, pix_y, SNR):
        self.pix_number = number
        self.pix_x = pix_x
        self.pix_y = pix_y
        #self.Signal = signal
        #self.Noise = noise
        self.StN = SNR
        #if self.Noise != 0:
        #    self.StN = self.Signal/np.sqrt(self.Noise)
        self.neighbors = []
        self.neighbors_x = []
        self.neighbors_y = []
        self.assigned_to_bin = False
        self.assigned_bin = None
    def add_neighbor(self, pixel,x_pixel,y_pixel):
        self.neighbors.append(pixel)
        self.neighbors_x.append(x_pixel)
        self.neighbors_y.append(y_pixel)
    def add_to_bin(self,bin):
        self.assigned_to_bin = True
        self.assigned_bin = bin
    def clear_bin(self):
        self.assigned_to_bin = False
        self.assigned_bin = None

def read_in(SNR_map):
    #Collect Pixel Data
    print(os.getcwd())
    hdu_list = fits.open(SNR_map, memmap=True)
    counts = hdu_list[0].data
    y_len = counts.shape[0]
    x_len = counts.shape[1]
    hdu_list.close()
    x_min = 0; y_min = 0;
    x_max = x_len; y_max = y_len;
    Pixels = []
    pixel_count = 0
    for col in range(int(x_len)):
        for row in range(int(y_len)):
            SNR = counts[row][col]
            Pixels.append(Pixel(pixel_count,x_min+col,y_min+row,SNR)) #Bottom Left Corner!
            pixel_count += 1
    print("We have "+str(pixel_count)+" Pixels! :)")
    return Pixels, x_min, x_max, y_min, y_max

def Nearest_Neighbors(pixel_list):
    print("Running Nearest Neighbor Algorithm")
    xvals = []
    yvals = []
    num_neigh = 9
    for pixel in pixel_list:
        xvals.append(pixel.pix_x)
        yvals.append(pixel.pix_y)
    X = np.column_stack((xvals,yvals))
    nbrs = NearestNeighbors(n_neighbors=num_neigh, algorithm='ball_tree').fit(X)
    distances, indices = nbrs.kneighbors(X)
    pix_num = 0
    for pixel in pixel_list:
        for j in range(num_neigh-1):
            if distances[pix_num][j+1] == 1:
                index = indices[pix_num][j+1]
                pixel.add_neighbor(pixel_list[index],xvals[index],yvals[index])
            else:
                pass #not adjacent
        pix_num += 1
    print("Finished Nearest Neighbor Algorithm")
    return None

def dist(p1x,p1y,p2x,p2y):
    return np.sqrt((p1x-p2x)**2+(p1y-p2y)**2)
def closest_node(Bin_current,unassigned_pixels):
    closest_val = 1e16 #just some big number
    p1x = Bin_current.centroidx[0]
    p1y = Bin_current.centroidy[0]
    for pix_neigh in Bin_current.pixel_neighbors:
        if pix_neigh.assigned_to_bin == False: #Dont bother with already assigned pixels!
            p2x = pix_neigh.pix_x
            p2y = pix_neigh.pix_y
            new_dist = dist(p1x,p1y,p2x,p2y)
            if new_dist < closest_val:
                closest_val = new_dist
                closest_pixel = pix_neigh
    if closest_val == 1e16:
        #None of the neighbors work so just pick a random unassigned pixel
        dist_list = [dist(p1x,p1y,pixel.pix_x,pixel.pix_y) for pixel in unassigned_pixels]
        closest_pixel = unassigned_pixels[dist_list.index(min(dist_list))]
    return closest_pixel

def adjacency(current_bin,closest_node):
    if closest_node in current_bin.pixel_neighbors:
        return True
    else:
        return False

def Roundness(current_bin,closest_pixel,pixel_length):
    pixel_list = current_bin.pixels
    n = len(pixel_list)+1 #Including new
    rad_equiv = np.sqrt(n/np.pi)*(pixel_length)
    cen_x_new = 0 ; cen_y_new = 0
    xvals = [pixel.pix_x for pixel in current_bin.pixels]
    yvals = [pixel.pix_y for pixel in current_bin.pixels]
    cen_x_new = (sum(xvals)+closest_pixel.pix_x)/n
    cen_y_new = (sum(yvals)+closest_pixel.pix_y)/n
    dists = []
    for i in range(n-1):
        dists.append(np.sqrt((xvals[i] - cen_x_new)**2+(yvals[i] - cen_y_new)**2))
    dists.append(np.sqrt((closest_pixel.pix_x-cen_x_new)**2+(closest_pixel.pix_y-cen_y_new)**2))
    rad_max = max(dists) # maximum distance between the centroid of the bin and any of the bin pixels
    roundness = rad_max/rad_equiv - 1.
    return roundness

def Potential_SN(Current_bin,closest_pix):
    Current_bin.add_pixel(closest_pix)
    new_StN = Current_bin.StN[0]
    Current_bin.remove_pixel(closest_pix)
    return new_StN

def Bin_data(Bins,missing_pixels,min_x,min_y, output_directory, filename):
    Bins.sort(key=lambda bin: bin.bin_number)
    file = open(output_directory+'/'+filename+'.txt',"w+")
    file.write("This text file contains information necessary for chandra to bin the pixels appropriately for image.fits \n")
    file.write("pixel_x pixel_y bin \n")
    file2 = open(output_directory+'/'+filename+'_bins.txt','w+')
    file2.write("Bin data for paraview script to plot Weighted Voronoi Diagram \n")
    file2.write("centroidx centroidy weight \n")
    binCount = 0
    for bin in Bins:
        for pixel in bin.pixels:
            file.write(str(pixel.pix_x-min_x)+" "+str(pixel.pix_y-min_y)+" "+str(binCount)+' \n')
        file2.write(str(bin.centroidx[0])+" "+str(bin.centroidy[0])+" "+str(bin.scale_length[0])+" \n")
        binCount += 1
    file.close()
    file2.close()
    return None

def bin_num_pixel(binList,currentPixel):
    bin_number_interest = None
    for binNumber in range(len(binList)):
        if currentPixel in binList[binNumber].pixels:
            bin_number_interest = binNumber
            break
    return bin_number_interest

def assigned_missing_pixels(pixels):
    #get all pixels in domain but not already assigned
    pos_x = [pix.pix_x for pix in pixels]
    pos_y = [pix.pix_y for pix in pixels]
    already_pixel = {}
    for i in range(len(pos_x)):
        already_pixel[(pos_x[i],pos_y[i])] = True
    x_range = [i for i in range(min(pos_x), max(pos_x)+1)]
    y_range = [i for i in range(min(pos_y), max(pos_y)+1)]
    Pixels_unbinned = []
    unbinned_num = max([pix.pix_number for pix in pixels])+1
    for x in x_range:
        for y in y_range:
            if (x,y) not in already_pixel.keys():
                new_pixel = Pixel(unbinned_num,x,y,0)
                Pixels_unbinned.append(new_pixel)
                unbinned_num += 1
            else:
                pass #Pixel already binned
    print("We have "+str(len(Pixels_unbinned))+" unbinned pixels")
    return Pixels_unbinned

def reassign_pixels(bin,bins_successful,sucessful_centroids):
    for pixel in bin.pixels:
        pixel.clear_bin()
        potential_bins = bins_successful[:]
        potential_centroids = sucessful_centroids[:]
        while pixel.assigned_to_bin == False:
            try:
                dists = [dist(pixel.pix_x,pixel.pix_y,centx,centy) for (centx,centy) in potential_centroids]
                closest_bin_index = dists.index(min(dists))
                closest_bin = potential_bins[closest_bin_index]
                if closest_bin.availabilty == False and len(potential_centroids) > 1:
                    del potential_centroids[closest_bin_index]
                    potential_bins.remove(closest_bin)
                else:
                    pixel.add_to_bin(closest_bin)
                    closest_bin.add_pixel(pixel)
            except:
                pass#print(potential_centroids)
    return None

def Bin_Acc(Pixels,pixel_length,StN_Target,roundness_crit):
    #step 1:setup list of bin objects
    print("Starting Bin Accretion Algorithm")
    unassigned_pixels = Pixels[:]
    binCount = 0
    Bin_list = []
    bins_successful = []
    criteria_a = False; criteria_b = False; criteria_c = False
    StN_list_pixels = [pixel.StN for pixel in Pixels]
    max_StN_ind = StN_list_pixels.index(max(StN_list_pixels))
    max_StN_pix = Pixels[max_StN_ind]
    Current_bin = Bin(binCount)
    Bin_list.append(Current_bin)
    Current_bin.add_pixel(max_StN_pix)
    Current_bin.CalcCentroid()
    max_StN_pix.add_to_bin(binCount)
    unassigned_pixels.remove(max_StN_pix)
    closest_pix = None
    closest_not_in_pix = None
    while len(unassigned_pixels) != 0:
        criteria_a = True; criteria_b = True; criteria_c = True;
        while (criteria_a == True and criteria_b == True and criteria_c == True) and len(unassigned_pixels) != 0:
            closest_pix = closest_node(Current_bin,unassigned_pixels)
            criteria_a = adjacency(Current_bin,closest_pix)
            criteria_b = True if (Roundness(Current_bin,closest_pix,pixel_length) < roundness_crit) else False
            criteria_c = True if (Potential_SN(Current_bin,closest_pix) < 0.75*StN_Target) else False
            if (criteria_a == True and criteria_b == True and criteria_c == True):
                Current_bin.add_pixel(closest_pix)
                Current_bin.CalcCentroid()
                closest_pix.add_to_bin(binCount)
                unassigned_pixels.remove(closest_pix)
            else:
                closest_not_in_pix = closest_pix
        if (Current_bin.StN[0] > 0.5*StN_Target):
            Current_bin.success()
            bins_successful.append(Current_bin)
        #else:
            #print(Current_bin.StN[0])
        if len(unassigned_pixels) == 0:
            break #All pixels assigned so dont create a new bin. that would be silly
        else:
            binCount += 1
            Current_bin = Bin(binCount)
            Bin_list.append(Current_bin)
            Current_bin.add_pixel(closest_not_in_pix)
            Current_bin.CalcCentroid()
            closest_not_in_pix.add_to_bin(binCount)
            unassigned_pixels.remove(closest_not_in_pix)
    print('We have '+str(len(Bin_list))+' bins.')
    print('We have '+str(len(unassigned_pixels))+' unassigned pixels.')
    print("Reassigning unsuccessful bins")
    for bin in bins_successful:
        bin.CalcCentroid()
        bin.CalcArea(pixel_length)
        bin.CalcScaleLength(StN_Target)
    sucessful_centroids = [(bin.centroidx[0],bin.centroidy[0]) for bin in bins_successful]
    for bin in Bin_list:
        if bin.successful == False:
            reassign_pixels(bin,bins_successful,sucessful_centroids)
            bin.clear_pixels()
        else:
            bin.CalcCentroid()
            bin.CalcArea(pixel_length)
            bin.CalcScaleLength(StN_Target)
    for bin in bins_successful: bin.update_StN_prev()
    print("Completed Bin Accretion Algorithm")
    print("There are a total of "+str(len(bins_successful)+1)+" bins!")
    return bins_successful

def converged_met(Bins,ToL):
    True_count = 0
    for bin in Bins:
        StN_old = bin.StN_prev[0]
        StN_new = bin.StN[0]
        if abs(StN_new-StN_old)/StN_old < ToL:
            True_count += 1
        else:
            pass
    if True_count/len(Bins) > 0.9:
        return True
    else:
        return False

def Rebin_Pixels(binList,pixel_list,pixel_length,StN_Target):
    for bin in binList:
        bin.clear_pixels()
    WVT_successful_bins = []
    for pixel in pixel_list:
        pixel.clear_bin()
        distances_bins = [dist(pixel.pix_x,pixel.pix_y,bin.centroidx_prev[0],bin.centroidy_prev[0])/bin.scale_length_prev[0] for bin in binList]
        closest_bin = binList[distances_bins.index(min(distances_bins))]
        pixel.add_to_bin(closest_bin)
        closest_bin.add_pixel(pixel)
        if closest_bin.StN[0] > 0 and closest_bin.WVT_successful == False:
            WVT_successful_bins.append(closest_bin)
            closest_bin.WVT_success()
    sucessful_centroids = [(bin.centroidx[0],bin.centroidy[0]) for bin in WVT_successful_bins]
    for bin in binList:
        if bin not in WVT_successful_bins:
            reassign_pixels(bin,WVT_successful_bins,sucessful_centroids)
    for bin in WVT_successful_bins:
        bin.CalcCentroid()
        bin.CalcArea(pixel_length)
        bin.CalcScaleLength(StN_Target)
    return WVT_successful_bins

def WVT(Bin_list_init,Pixel_Full,StN_Target,ToL,pixel_length,image_dir):
    print("Beginning WVT")
    Bin_list_prev = Bin_list_init[:]
    converged = False
    its_to_conv = 0
    while converged == False and its_to_conv<5:
        print("We are on step "+str(its_to_conv+1))
        bins_with_SN = Rebin_Pixels(Bin_list_prev,Pixel_Full,pixel_length,StN_Target)[:]
        converged = converged_met(bins_with_SN,ToL)
        Bin_list_prev = bins_with_SN[:]
        for bin in bins_with_SN: bin.update_StN_prev()
        bin_SN_List = [bin.StN[0] for bin in bins_with_SN]
        plt.hist(bin_SN_List)
        plt.ylabel("Number of Bins")
        plt.xlabel("Signal-to-Noise")
        its_to_conv += 1
        plt.savefig(image_dir+'/histograms/iteration_'+str(its_to_conv)+".png")
        plt.clf()
    print("Completed WVT in "+str(its_to_conv)+" step(s)!!")
    print("There are a total of "+str(len(bins_with_SN)+1)+" bins!")
    return bins_with_SN

max_stars_repo_name: benjaminvigneron/LUCI

max_issues_repo_name: benjaminvigneron/LUCI

max_forks_repo_name: benjaminvigneron/LUCI

reason: import numpy,from astropy


text:
import numpy as np
import torchvision.transforms.functional as F
import pyjet.data as data
import pyjet.augmenters as aug
import pyjet.backend as J
import logging
from scipy.misc import imresize
from skimage.util import random_noise

"""This file contains the various data objects and manipulators
I'll use to handle the data. The RandomCropper implements a pyjet
augmenter and crops the given images to
"""

class DataTransformer(object):
    """This class is designed to transform a stream
    of data into some new output data. More specifically
    it takes the output of a DatasetGenerator and yields
    an output of its own
    """

    def __init__(self, labels=False):
        self.labels = labels

    def transform(self, batch):
        raise NotImplementedError

    def __call__(self, generator):
        return TransformerGenerator(self, generator)

class TransformerGenerator(data.BatchGenerator):
    def __init__(self, transformer, generator):
        # Copy the steps per epoch and batch size if it has one
        if hasattr(generator, "steps_per_epoch") and hasattr(generator, "batch_size"):
            super(TransformerGenerator, self).__init__(
                steps_per_epoch=generator.steps_per_epoch,
                batch_size=generator.batch_size,
            )
        else:
            logging.warning(
                "Input generator does not have a "
                "steps_per_epoch or batch_size "
                "attribute. Continuing without them."
            )

        self.transformer = transformer
        self.generator = generator

    def __next__(self):
        return self.transformer.transform(next(self.generator))

def batch_apply(func, x):
    x_out = np.stack([func(xi) for xi in x], axis=0)
    return x_out

def crappfier(factor, noise=True):
    """A simple function to get a bilinear resizer with given factor"""

    # Do nothing if we're not crappifying
    if factor == 1.0:
        return lambda x: x

    def crappify_img(npimg):
        npimg = imresize(
            imresize(npimg, factor, interp="bilinear"), 1 / factor, interp="bilinear"
        )
        if noise:
            npimg = random_noise(npimg)
        else:
            # convert to float
            npimg = npimg / 255.0

        return npimg

    return lambda x: batch_apply(crappify_img, x)

def cropper(crop_size):
    def crop_img(npimg):
        h, w, _ = npimg.shape
        ch, cw = crop_size
        # If the image is too small, we'll make it bigger
        # This copies the entire image which might not be necessary
        # If speed becomes and issue look here.
        if h < ch or w < cw:
            newimg_shape = (min(ch, h), min(cw, w))
            new_npimg = np.zeros(newimg_shape + npimg.shape[2:])
            new_npimg[[slice(i) for i in npimg.shape]] = npimg
            npimg = new_npimg
            h, w, _ = npimg.shape

        cropy = np.random.randint(h - ch + 1)
        cropx = np.random.randint(w - cw + 1)
        return npimg[cropy : cropy + ch, cropx : cropx + cw]

    return lambda x: batch_apply(crop_img, x)

def channels_first(x):
    return x.transpose(0, 3, 1, 2).astype(np.float32)

def transform(x, crop_size, factor, noise=True):
    crap = crappfier(factor, noise=noise)
    crop = cropper(crop_size)
    # Crop the image and crappify
    x = crop(x)
    y = x / 255.0  # get the label and convert to float
    x = crap(x)
    # Convert to channels first
    return channels_first(x), channels_first(y)

class PrepareData(DataTransformer):
    def __init__(self, crop_size, factor, noise):
        super().__init__(labels=False)
        self.crop_size = crop_size
        self.factor = factor
        self.noise = noise

    def transform(self, x):
        return transform(x, self.crop_size, self.factor, noise=self.noise)

reason: import numpy,from scipy