text:
# -*- coding: utf-8 -*-
"""
L-functions from PARI

This is a wrapper around the general PARI L-functions functionality.

AUTHORS:

- Frédéric Chapoton (2018) interface

"""
# ****************************************************************************
#       Copyright (C) 2018 Frédéric Chapoton <chapoton@unistra.fr>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#                  https://www.gnu.org/licenses/
# ****************************************************************************

from operator import index as PyNumber_Index
from cypari2.gen import Gen
from sage.libs.pari import pari
from sage.structure.sage_object import SageObject
from sage.rings.all import (ZZ, RealField, ComplexField, PowerSeriesRing)

class lfun_generic(object):
    r"""
    Create a PARI `L`-function (:pari:`lfun` instance).

    The arguments are::

        lfun_generic(conductor, gammaV, weight, eps, poles, residues, init)

    where

    - ``conductor`` -- integer, the conductor

    - ``gammaV`` -- list of Gamma-factor parameters, e.g. [0] for
      Riemann zeta, [0,1] for ell.curves, (see examples).

    - ``weight`` -- positive real number, usually an integer e.g. 1 for
      Riemann zeta, 2 for `H^1` of curves/`\QQ`

    - ``eps`` -- complex number; sign in functional equation

    - ``poles`` -- (default: []) list of points where `L^*(s)` has
      (simple) poles; only poles with `Re(s)>weight/2` should be
      included

    - ``residues`` -- vector of residues of `L^*(s)` in those poles or
      set residues='automatic' (default value)

    - ``init`` -- list of coefficients

    RIEMANN ZETA FUNCTION:

    We compute with the Riemann Zeta function::

        sage: from sage.lfunctions.pari import lfun_generic, LFunction
        sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[1])
        sage: lf.init_coeffs([1]*2000)

    Now we can wrap this PARI L-function into one Sage L-function::

        sage: L = LFunction(lf); L
        L-series of conductor 1 and weight 1
        sage: L(1)
        Traceback (most recent call last):
        ...
        ArithmeticError: pole here
        sage: L(2)
        1.64493406684823
        sage: L.derivative(2)
        -0.937548254315844
        sage: h = RR('0.0000000000001')
        sage: (zeta(2+h) - zeta(2.))/h
        -0.937028232783632
        sage: L.taylor_series(2, k=5)
        1.64493406684823 - 0.937548254315844*z + 0.994640117149451*z^2 - 1.00002430047384*z^3 + 1.00006193307...*z^4 + O(z^5)
    """
    def __init__(self, conductor, gammaV, weight, eps, poles=[],
                 residues='automatic', prec=None, *args, **kwds):
        """
        Initialisation of a :pari:`lfun` from motivic data.

        This can happen either in one stage or in two stages: the coefficients
        can be given using the ``init`` keyword, or entered using
        the :meth:`init_coeffs`.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_generic, LFunction
            sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[1])
        """
        # before entering the coefficients, this attribute is None
        self._L = None

        self.conductor = conductor
        self.gammaV = gammaV
        self.weight = weight
        self.eps = eps
        self.poles = poles
        self.residues = residues
        self.prec = prec

        if args or kwds:
            self.init_coeffs(*args, **kwds)

    def init_coeffs(self, v, cutoff=None, w=1, *args, **kwds):
        """
        Set the coefficients `a_n` of the `L`-series.

        If `L(s)` is not equal to its dual, pass the coefficients of
        the dual as the second optional argument.

        INPUT:

        -  ``v`` -- list of complex numbers or unary function

        -  ``cutoff`` -- unused

        -  ``w`` --  list of complex numbers or unary function

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_generic, LFunction
            sage: lf = lfun_generic(conductor=1, gammaV=[0,1], weight=12, eps=1)
            sage: pari_coeffs = pari('k->vector(k,n,(5*sigma(n,3)+7*sigma(n,5))*n/12 - 35*sum(k=1,n-1,(6*k-4*(n-k))*sigma(k,3)*sigma(n-k,5)))')
            sage: lf.init_coeffs(pari_coeffs)

        Evaluate the resulting L-function at a point, and compare with
        the answer that one gets "by definition" (of L-function
        attached to a modular form)::

            sage: L = LFunction(lf)
            sage: L(14)
            0.998583063162746
            sage: a = delta_qexp(1000)
            sage: sum(a[n]/float(n)^14 for n in range(1,1000))
            0.9985830631627459

        Illustrate that one can give a list of complex numbers for v
        (see :trac:`10937`)::

            sage: l2 = lfun_generic(conductor=1, gammaV=[0,1], weight=12, eps=1)
            sage: l2.init_coeffs(list(delta_qexp(1000))[1:])
            sage: L2 = LFunction(l2)
            sage: L2(14)
            0.998583063162746

        TESTS:

        Verify that setting the `w` parameter does not raise an error
        (see :trac:`10937`)::

            sage: L2 = lfun_generic(conductor=1, gammaV=[0,1], weight=12, eps=1)
            sage: L2.init_coeffs(list(delta_qexp(1000))[1:], w=[1..1000])

        Additional arguments are ignored for compatibility with the old
        Dokchitser script::

            sage: L2.init_coeffs(list(delta_qexp(1000))[1:], foo="bar")
            doctest:...: DeprecationWarning: additional arguments for initializing an lfun_generic are ignored
            See https://trac.sagemath.org/26098 for details.
        """
        if args or kwds:
            from sage.misc.superseded import deprecation
            deprecation(26098, "additional arguments for initializing an lfun_generic are ignored")

        v = pari(v)
        if v.type() not in ('t_CLOSURE', 't_VEC'):
            raise TypeError("v (coefficients) must be a list or a function")

        # w = 0 means a*_n = a_n
        # w = 1 means a*_n = complex conjugate of a_n
        # otherwise w must be a list of coefficients
        w = pari(w)
        if w.type() not in ('t_INT', 't_CLOSURE', 't_VEC'):
            raise TypeError("w (dual coefficients) must be a list or a function or the special value 0 or 1")

        if not self.poles:
            self._L = pari.lfuncreate([v, w, self.gammaV, self.weight,
                                       self.conductor, self.eps])
        else:
            # TODO
            # poles = list(zip(poles, residues))
            self._L = pari.lfuncreate([v, w, self.gammaV, self.weight,
                                       self.conductor, self.eps,
                                       self.poles[0]])

    def __pari__(self):
        """
        Return the PARI L-function object.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_generic, LFunction
            sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[1])
            sage: lf.__pari__()
            Traceback (most recent call last):
            ...
            ValueError: call init_coeffs on the L-function first

            sage: lf.init_coeffs([1]*2000)
            sage: X = lf.__pari__()
            sage: X.type()
            't_VEC'
        """
        if self._L is None:
            raise ValueError("call init_coeffs on the L-function first")
        return self._L

def lfun_character(chi):
    """
    Create the L-function of a primitive Dirichlet character.

    If the given character is not primitive, it is replaced by its
    associated primitive character.

    OUTPUT:

    one :pari:`lfun` object

    EXAMPLES::

        sage: from sage.lfunctions.pari import lfun_character, LFunction
        sage: chi = DirichletGroup(6).gen().primitive_character()
        sage: L = LFunction(lfun_character(chi))
        sage: L(3)
        0.884023811750080

    TESTS:

    A non-primitive one::

        sage: L = LFunction(lfun_character(DirichletGroup(6).gen()))
        sage: L(4)
        0.940025680877124

    With complex arguments::

        sage: from sage.lfunctions.pari import lfun_character, LFunction
        sage: chi = DirichletGroup(6, CC).gen().primitive_character()
        sage: L = LFunction(lfun_character(chi))
        sage: L(3)
        0.884023811750080

    Check the values::

        sage: chi = DirichletGroup(24)([1,-1,-1]); chi
        Dirichlet character modulo 24 of conductor 24
        mapping 7 |--> 1, 13 |--> -1, 17 |--> -1
        sage: Lchi = lfun_character(chi)
        sage: v = [0] + Lchi.lfunan(30).sage()
        sage: all(v[i] == chi(i) for i in (7,13,17))
        True
    """
    if not chi.is_primitive():
        chi = chi.primitive_character()
    G, v = chi._pari_init_()
    return pari.lfuncreate([G, v])

def lfun_elliptic_curve(E):
    """
    Create the L-function of an elliptic curve.

    OUTPUT:

    one :pari:`lfun` object

    EXAMPLES::

        sage: from sage.lfunctions.pari import lfun_elliptic_curve, LFunction
        sage: E = EllipticCurve('11a1')
        sage: L = LFunction(lfun_elliptic_curve(E))
        sage: L(3)
        0.752723147389051
        sage: L(1)
        0.253841860855911

    Over number fields::

        sage: K.<a> = QuadraticField(2)
        sage: E = EllipticCurve([1,a])
        sage: L = LFunction(lfun_elliptic_curve(E))
        sage: L(3)
        1.00412346717019
    """
    return pari.lfuncreate(E)

def lfun_number_field(K):
    """
    Create the Dedekind zeta function of a number field.

    OUTPUT:

    one :pari:`lfun` object

    EXAMPLES::

        sage: from sage.lfunctions.pari import lfun_number_field, LFunction

        sage: L = LFunction(lfun_number_field(QQ))
        sage: L(3)
        1.20205690315959

        sage: K = NumberField(x**2-2, 'a')
        sage: L = LFunction(lfun_number_field(K))
        sage: L(3)
        1.15202784126080
        sage: L(0)
        0.000000000000000
    """
    return pari.lfuncreate(K)

def lfun_eta_quotient(scalings, exponents):
    """
    Return the L-function of an eta-quotient.

    INPUT:

    - scalings -- a list of integers, the scaling factors

    - exponents -- a list of integers, the exponents

    EXAMPLES::

        sage: from sage.lfunctions.pari import lfun_eta_quotient, LFunction
        sage: L = LFunction(lfun_eta_quotient([1], [24]))
        sage: L(1)
        0.0374412812685155

        sage: lfun_eta_quotient([6],[4])
        [[Vecsmall([7]), [Vecsmall([6]), Vecsmall([4])]], 0, [0, 1], 2, 36, 1]

        sage: lfun_eta_quotient([2,1,4], [5,-2,-2])
        Traceback (most recent call last):
        ...
        PariError: sorry, noncuspidal eta quotient is not yet implemented
    """
    from sage.matrix.constructor import matrix
    N = len(scalings)
    if N != len(exponents):
        raise ValueError('arguments should have the same length')
    m = matrix(ZZ, N, 2, [(x, y) for x, y in zip(scalings, exponents)])
    return pari.lfunetaquo(m)

def lfun_delta():
    """
    Return the L-function of Ramanujan's Delta modular form.

    EXAMPLES::

        sage: from sage.lfunctions.pari import lfun_delta, LFunction
        sage: L = LFunction(lfun_delta())
        sage: L(1)
        0.0374412812685155
    """
    return lfun_eta_quotient([1], [24])

def lfun_quadratic_form(qf):
    """
    Return the L-function of a positive definite quadratic form.

    EXAMPLES::

        sage: from sage.lfunctions.pari import lfun_quadratic_form, LFunction
        sage: Q = QuadraticForm(ZZ,2,[2,3,4])
        sage: L = LFunction(lfun_quadratic_form(Q))
        sage: L(3)
        0.377597233183583
    """
    if not qf.is_positive_definite():
        raise ValueError('quadratic form must be positive definite')
    return pari.lfunqf(qf.matrix())

class LFunction(SageObject):
    r"""
    Build the L-function from a PARI L-function.

    .. RUBRIC:: Rank 1 elliptic curve

    We compute with the `L`-series of a rank `1` curve. ::

        sage: E = EllipticCurve('37a')
        sage: L = E.lseries().dokchitser(algorithm="pari"); L
        PARI L-function associated to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
        sage: L(1)
        0.000000000000000
        sage: L.derivative(1)
        0.305999773834052
        sage: L.derivative(1,2)
        0.373095594536324
        sage: L.num_coeffs()
        50
        sage: L.taylor_series(1,4)
        0.000000000000000 + 0.305999773834052*z + 0.186547797268162*z^2 - 0.136791463097188*z^3 + O(z^4)
        sage: L.check_functional_equation()  # abs tol 4e-19
        1.08420217248550e-19

    .. RUBRIC:: Rank 2 elliptic curve

    We compute the leading coefficient and Taylor expansion of the
    `L`-series of a rank `2` elliptic curve::

        sage: E = EllipticCurve('389a')
        sage: L = E.lseries().dokchitser(algorithm="pari")
        sage: L.num_coeffs()
        163
        sage: L.derivative(1,E.rank())
        1.51863300057685
        sage: L.taylor_series(1,4)
        ...e-19 + (...e-19)*z + 0.759316500288427*z^2 - 0.430302337583362*z^3 + O(z^4)

    .. RUBRIC:: Number field

    We compute with the Dedekind zeta function of a number field::

        sage: x = var('x')
        sage: K = NumberField(x**4 - x**2 - 1,'a')
        sage: L = K.zeta_function(algorithm="pari")
        sage: L.conductor
        400
        sage: L.num_coeffs()
        348
        sage: L(2)
        1.10398438736918
        sage: L.taylor_series(2,3)
        1.10398438736918 - 0.215822638498759*z + 0.279836437522536*z^2 + O(z^3)

    .. RUBRIC:: Ramanujan `\Delta` L-function

    The coefficients are given by Ramanujan's tau function::

        sage: from sage.lfunctions.pari import lfun_generic, LFunction
        sage: lf = lfun_generic(conductor=1, gammaV=[0,1], weight=12, eps=1)
        sage: tau = pari('k->vector(k,n,(5*sigma(n,3)+7*sigma(n,5))*n/12 - 35*sum(k=1,n-1,(6*k-4*(n-k))*sigma(k,3)*sigma(n-k,5)))')
        sage: lf.init_coeffs(tau)
        sage: L = LFunction(lf)

    Now we are ready to evaluate, etc. ::

        sage: L(1)
        0.0374412812685155
        sage: L.taylor_series(1,3)
        0.0374412812685155 + 0.0709221123619322*z + 0.0380744761270520*z^2 + O(z^3)
    """
    def __init__(self, lfun, prec=None):
        """
        Initialization of the L-function from a PARI L-function.

        INPUT:

        - lfun -- a PARI :pari:`lfun` object or an instance of :class:`lfun_generic`
        - prec -- integer (default: 53) number of *bits* of precision

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_generic, LFunction
            sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], v=pari('k->vector(k,n,1)'))
            sage: L = LFunction(lf)
            sage: L.num_coeffs()
            4
        """
        if isinstance(lfun, lfun_generic):
            # preparation using motivic data
            self._L = lfun.__pari__()
            if prec is None:
                prec = lfun.prec
        elif isinstance(lfun, Gen):
            # already some PARI lfun
            self._L = lfun
        else:
            # create a PARI lfunction from other input data
            self._L = pari.lfuncreate(lfun)

        self._conductor = ZZ(self._L[4])  # needs check
        self._weight = ZZ(self._L[3])  # needs check
        if prec is None:
            self.prec = 53
        else:
            self.prec = PyNumber_Index(prec)
        self._RR = RealField(self.prec)
        self._CC = ComplexField(self.prec)
        # Complex field used for inputs, which ensures exact 1-to-1
        # conversion to/from PARI. Since PARI objects have a precision
        # in machine words (not bits), this is typically higher. For
        # example, the default of 53 bits of precision would become 64.
        self._CCin = ComplexField(pari.bitprecision(self._RR(1)))

    def _repr_(self):
        """
        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: L = LFunction(lfun_number_field(QQ)); L
            L-series of conductor 1 and weight 1
        """
        return "L-series of conductor %s and weight %s" % (self._conductor,
                                                           self._weight)

    @property
    def conductor(self):
        """
        Return the conductor.

        EXAMPLES::

            sage: from sage.lfunctions.pari import *
            sage: L = LFunction(lfun_number_field(QQ)); L.conductor
            1
            sage: E = EllipticCurve('11a')
            sage: L = LFunction(lfun_number_field(E)); L.conductor
            11
        """
        return self._conductor

    def num_coeffs(self, T=1):
        """
        Return number of coefficients `a_n` that are needed in
        order to perform most relevant `L`-function computations to
        the desired precision.

        EXAMPLES::

            sage: E = EllipticCurve('11a')
            sage: L = E.lseries().dokchitser(algorithm="pari")
            sage: L.num_coeffs()
            27
            sage: E = EllipticCurve('5077a')
            sage: L = E.lseries().dokchitser(algorithm="pari")
            sage: L.num_coeffs()
            591

            sage: from sage.lfunctions.pari import lfun_generic, LFunction
            sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], v=pari('k->vector(k,n,1)'))
            sage: L = LFunction(lf)
            sage: L.num_coeffs()
            4
        """
        # lfuncost takes a domain
        domain = [pari(1) / 2, 1, 1]
        return ZZ(pari.lfuncost(self._L, domain)[0])

    def Lambda(self, s):
        """
        Evaluate the completed L-function at s.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: L = LFunction(lfun_number_field(QQ))
            sage: L.Lambda(2)
            0.523598775598299
            sage: L.Lambda(1-2)
            0.523598775598299
        """
        s = self._CCin(s)
        R = self._CC
        return R(pari.lfunlambda(self._L, s, precision=self.prec))

    def hardy(self, t):
        """
        Evaluate the associated Hardy function at t.

        This only works for real t.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: L = LFunction(lfun_number_field(QQ))
            sage: L.hardy(2)
            -0.962008487244041

        TESTS::

            sage: L.hardy(.4+.3*I)
            Traceback (most recent call last):
            ...
            PariError: incorrect type in lfunhardy (t_COMPLEX)
        """
        t = self._CCin(t)
        R = self._RR
        return R(pari.lfunhardy(self._L, t, precision=self.prec))

    def derivative(self, s, D=1):
        """
        Return the derivative of the L-function at point s and order D.

        INPUT:

        -  ``s`` -- complex number

        - ``D`` -- optional integer (default 1)

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: L = LFunction(lfun_number_field(QQ))
            sage: L.derivative(2)
            -0.937548254315844
        """
        s = self._CCin(s)
        R = self._CC
        return R(pari.lfun(self._L, s, D, precision=self.prec))

    def taylor_series(self, s, k=6, var='z'):
        r"""
        Return the first `k` terms of the Taylor series expansion
        of the `L`-series about `s`.

        This is returned as a formal power series in ``var``.

        INPUT:

        -  ``s`` -- complex number; point about which to expand

        -  ``k`` -- optional integer (default: 6), series is
           `O(``var``^k)`

        -  ``var`` -- optional string (default: 'z'), variable of power
           series

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: lf = lfun_number_field(QQ)
            sage: L = LFunction(lf)
            sage: L.taylor_series(2, 3)
            1.64493406684823 - 0.937548254315844*z + 0.994640117149451*z^2 + O(z^3)
            sage: E = EllipticCurve('37a')
            sage: L = E.lseries().dokchitser(algorithm="pari")
            sage: L.taylor_series(1)
            0.000000000000000 + 0.305999773834052*z + 0.186547797268162*z^2 - 0.136791463097188*z^3 + 0.0161066468496401*z^4 + 0.0185955175398802*z^5 + O(z^6)

        We compute a Taylor series where each coefficient is to high
        precision::

            sage: E = EllipticCurve('389a')
            sage: L = E.lseries().dokchitser(200,algorithm="pari")
            sage: L.taylor_series(1,3)
            2...e-63 + (...e-63)*z + 0.75931650028842677023019260789472201907809751649492435158581*z^2 + O(z^3)

        Check that :trac:`25402` is fixed::

            sage: L = EllipticCurve("24a1").modular_form().lseries()
            sage: L.taylor_series(-1, 3)
            0.000000000000000 - 0.702565506265199*z + 0.638929001045535*z^2 + O(z^3)
        """
        pt = pari.Ser([s, 1], d=k)  # s + x + O(x^k)
        B = PowerSeriesRing(self._CC, var)
        return B(pari.lfun(self._L, pt, precision=self.prec))

    def zeros(self, maxi):
        """
        Return the zeros with imaginary part bounded by maxi.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: lf = lfun_number_field(QQ)
            sage: L = LFunction(lf)
            sage: L.zeros(20)
            [14.1347251417347]
        """
        R = self._CC
        return [R(z) for z in pari.lfunzeros(self._L, maxi)]

    def _clear_value_cache(self):
        """
        Clear the cache where values of the function are stored.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_number_field, LFunction
            sage: lf = lfun_number_field(QQ)
            sage: L = LFunction(lf)
            sage: L(4)
            1.08232323371114
            sage: L._clear_value_cache()
        """
        del self.__values

    def __call__(self, s):
        r"""
        Return the value of the L-function at point ``s``.

        INPUT:

        -  ``s`` -- complex number

        .. NOTE::

            Evaluation of the function takes a long time, so each
            evaluation is cached. Call :meth:`_clear_value_cache` to
            clear the evaluation cache.

        EXAMPLES::

            sage: E = EllipticCurve('5077a')
            sage: L = E.lseries().dokchitser(100, algorithm="pari")
            sage: L(1)
            0.00000000000000000000000000000
            sage: L(1+I)
            -1.3085436607849493358323930438 + 0.81298000036784359634835412129*I
        """
        s = self._CC(s)
        try:
            return self.__values[s]
        except AttributeError:
            self.__values = {}
        except KeyError:
            pass
        try:
            value = self._CC(pari.lfun(self._L, s, precision=self.prec))
        except NameError:
            raise ArithmeticError('pole here')
        else:
            return value

    def check_functional_equation(self):
        r"""
        Verify how well numerically the functional equation is satisfied.

        If what this function returns does not look like 0 at all,
        probably the functional equation is wrong, i.e. some of the
        parameters gammaV, conductor, etc., or the coefficients are wrong.

        EXAMPLES::

            sage: from sage.lfunctions.pari import lfun_generic, LFunction
            sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[1], v=pari('k->vector(k,n,1)'))
            sage: L = LFunction(lf)
            sage: L.check_functional_equation()
            4.33680868994202e-19

        If we choose the sign in functional equation for the
        `\zeta` function incorrectly, the functional equation
        does not check out::

            sage: lf = lfun_generic(conductor=1, gammaV=[0], weight=1, eps=-1, poles=[1], residues=[1])
            sage: lf.init_coeffs([1]*2000)
            sage: L = LFunction(lf)
            sage: L.check_functional_equation()
            16.0000000000000
        """
        return self._RR(2) ** pari.lfuncheckfeq(self._L)

max_stars_repo_name: yzpopulation/sage

max_issues_repo_name: yzpopulation/sage

max_forks_repo_name: yzpopulation/sage

reason: from sage


text:
# Copyright 2018 by Frank Trampe.
# Preprocess to monochrome and get vectors from img2vec.
# Inputs for training include the positive directory, the negative directory, and the training data destination.

import numpy as np
import math
import sys
import argparse
import os
import random
import pickle
import sklearn.preprocessing
import sklearn.linear_model
import sklearn.svm
from img_to_vec import Img2Vec
from PIL import Image

# Initialize img2vec. This triggers downloading the neural net model if not present.
img2vec = Img2Vec(cuda=False, model='resnet-18')

def image_make_rgb(img):
	# https://stackoverflow.com/questions/9166400/convert-rgba-png-to-rgb-with-pil
	# img must be loaded first.
	flatimg = Image.new("RGB", img.size, (255, 255, 255))
	splitimg = img.split()
	if len(splitimg) > 3:
		flatimg.paste(img, mask=img.split()[3])
	else:
		flatimg.paste(img)
	return flatimg

def get_vec_from_file(fpath):
	# Read in an image from the path and get a vector from it.
	img = Image.open(fpath)
	img.load()
	return img2vec.get_vec(image_make_rgb(img))

def train_or_test_on_images(model, mode, negative_directory, positive_directory):
	# List the directories and flag the y values.
	flist_neg = map(lambda f : [negative_directory + f, 0], [f for f in os.listdir(negative_directory)])
	flist_pos = map(lambda f : [positive_directory + f, 1], [f for f in os.listdir(positive_directory)])
	# Vectorize and shuffle.
	training_data_vectorized = list(map(lambda fv : [get_vec_from_file(fv[0]), fv[1]], list(flist_neg) + list(flist_pos)))
	random.shuffle(training_data_vectorized)
	# Split x and y into separate arrays.
	training_data_separated = [list(map(lambda fv : fv[0], training_data_vectorized)), list(map(lambda fv : fv[1], training_data_vectorized))]
	# Normalize.
	training_data_normalized = [sklearn.preprocessing.normalize(training_data_separated[0]), training_data_separated[1]]
	# Train if specified.
	skmodel = model
	if (not skmodel):
		# skmodel = sklearn.linear_model.LogisticRegression()
		skmodel = sklearn.svm.SVC(probability = True)
		skmodel.fit(training_data_normalized[0], training_data_normalized[1])
	# Score.
	skscore = skmodel.score(training_data_normalized[0], training_data_normalized[1])
	# Report.
	print("Accuracy: %.3f." % skscore)
	if mode == 0:
		return skmodel
	elif mode == 1:
		return skscore
	return None

def predict_on_image_file_with_logic(model, fpath, slice_width, slice_height, opts):
	rv = None
	img = Image.open(fpath)
	img.load()
	imgw = img.size[0]
	imgh = img.size[1]
	slices = []
	if slice_width is None and slice_height is None:
		# No slices, so we just work on the whole image.
		slices += [img]
	else:
		# Set slice sizes if unset.
		slice_width_t = slice_width
		if slice_width_t is None:
			slice_width_t = imgw
		if type(slice_width_t) == str:
			slice_width_t = int(slice_width_t)
		slice_height_t = slice_height
		if slice_height_t is None:
			slice_height_t = imgh
		if type(slice_height_t) == str:
			slice_height_t = int(slice_height_t)
		slices = []
		# Cut the image into slices.
		slice_pos_y = 0
		while slice_pos_y * slice_height_t < imgh:
			slice_pos_x = 0
			while slice_pos_x * slice_width_t < imgw:
				cropbox = [slice_pos_x * slice_width_t, slice_pos_y * slice_height_t, \
				min((slice_pos_x + 1) * slice_width_t, imgw), min((slice_pos_y + 1) * slice_height_t, imgh)]
				slices += [img.crop(cropbox)]
				slice_pos_x += 1
			slice_pos_y += 1
	# Run predictions.
	predictor = model.predict_proba
	if opts["classify"]:
		predictor = model.predict
	slice_results = predictor(list(map(lambda tslice: img2vec.get_vec(image_make_rgb(tslice)), slices)))
	# print(slice_results)
	if slice_results is None or len(slice_results) < len(slices):
		return None
	# Apply thresholds.
	if opts["slice_threshold"] is not None and not opts["classify"]:
		slice_results = list(map(lambda tv: tv[1] > float(opts["slice_threshold"]), slice_results))
	# print(slice_results)
	# Invert individual results.
	if opts["slice_invert"] and (type(opts["slice_invert"]) == int or int(opts["slice_invert"])):
		slice_results = list(map(lambda tv: not tv, slice_results))
	# Run an or filter on all results to get a composite result.
	rv = False
	positive_count = len(list(filter(lambda tv: tv, slice_results)))
	if positive_count:
		rv = True
	# Check for an absolute threshold. If present, check that the count meets it.
	if type(opts["slice_count_threshold_proportional"]) is not None and \
	float(opts["slice_count_threshold_proportional"]) * len(slices) > positive_count:
		rv = False
	# Check for a proportional threshold. If present, check that the count meets it.
	if type(opts["slice_count_threshold"]) is not None and \
	float(opts["slice_count_threshold"]) > positive_count:
		rv = False
	# Invert if requested.
	if opts["group_invert"] and (type(opts["group_invert"]) == int or int(opts["group_invert"])):
		rv = not rv
	return int(rv)

def main():
	ap = argparse.ArgumentParser()
	ap.add_argument("--train-negative", required=False, help="directory of training images with desired negative output")
	ap.add_argument("--train-positive", required=False, help="directory of training images with desired positive output")
	ap.add_argument("--test-negative", required=False, help="directory of test images with desired negative output")
	ap.add_argument("--test-positive", required=False, help="directory of test images with desired positive output")
	ap.add_argument("--predict-file", required=False, help="an image on which to predict using the model")
	ap.add_argument("--classify", required=False, help="whether to classify rather than give a probability")
	ap.add_argument("--slice-width", required=False, help="width of slices on which to perform the prediction")
	ap.add_argument("--slice-height", required=False, help="height of slices on which to perform the prediction")
	ap.add_argument("--slice-threshold", required=False, default=0.5, help="probability cut-off per slice")
	ap.add_argument("--slice-count-threshold", required=False, default=0.5, help="positive slice count threshold")
	ap.add_argument("--slice-count-threshold-proportional", required=False, default=0.5, help="positive slice count threshold relative to total")
	ap.add_argument("--slice-result-invert", required=False, help="whether to invert the result of each slice before composing")
	ap.add_argument("--final-result-invert", required=False, help="whether to invert the final result")
	ap.add_argument("--model-input", required=False, help="a model to load")
	ap.add_argument("--model-output", required=False, help="a model to save")
	args = vars(ap.parse_args())
	model = None
	score = None
	if args["train_negative"] is not None and args["train_positive"] is not None:
		model = train_or_test_on_images(None, 0, args["train_negative"], args["train_positive"])
	elif args["model_input"] is not None:
		modelfile = open(args["model_input"], "rb")
		model = pickle.load(modelfile)
		modelfile.close()
	if args["test_negative"] is not None and args["test_positive"] is not None:
		score = train_or_test_on_images(model, 1, args["test_negative"], args["test_positive"])
	if args["predict_file"] is not None:
		predict_opts = {"slice_invert": args["slice_result_invert"], "group_invert": args["final_result_invert"], \
		"classify": args["classify"], "slice_threshold": args["slice_threshold"], \
		"slice_count_threshold": args["slice_count_threshold"], \
		"slice_count_threshold_proportional": args["slice_count_threshold_proportional"]}
		print(predict_on_image_file_with_logic(model, args["predict_file"], args["slice_width"], args["slice_height"], predict_opts))
	if args["model_output"] is not None:
		modelfile = open(args["model_output"], "wb")
		pickle.dump(model, modelfile)
		modelfile.close()

main()

ext: py

max_stars_repo_name: frank-trampe/binclass

max_issues_repo_name: frank-trampe/binclass

max_forks_repo_name: frank-trampe/binclass

reason: import numpy


text:
# -*- coding: utf-8 -*-

# Copyright 2021 Reo Yoneyama (Nagoya University)
#  MIT License (https://opensource.org/licenses/MIT)

"""Spectral envelope estimation module based on the idea of CheapTrick.
   Please see https://www.sciencedirect.com/science/article/pii/S0167639314000697 for details."""

import math

import torch
import torch.nn as nn
import torch.fft

class AdaptiveWindowing(nn.Module):

    def __init__(self,
                 sampling_rate,
                 hop_size,
                 fft_size,
                 f0_floor,
                 f0_ceil,
        ):
        """Initilize AdaptiveWindowing module.
        Args:
            sampling_rate (int): Sampling rate.
            hop_size (int): Hop size.
            fft_size (int): FFT size.
            f0_floor (int): Minimum value of F0.
            f0_ceil (int): Maximum value of F0.
        """
        super(AdaptiveWindowing, self).__init__()

        self.sampling_rate = sampling_rate
        self.hop_size = hop_size
        self.fft_size = fft_size
        self.window = torch.zeros((f0_ceil+1, fft_size)).cuda()
        self.zero_padding = nn.ConstantPad2d((fft_size // 2, fft_size // 2, 0, 0), 0)

        # Pre-calculation of the window functions
        for f0 in range(f0_floor, f0_ceil + 1):
            half_win_len = round(1.5 * self.sampling_rate / f0)
            base_index = torch.arange(-half_win_len, half_win_len + 1, dtype=torch.int64)
            position = base_index / 1.5 / self.sampling_rate
            left = fft_size // 2 - half_win_len
            right = fft_size // 2 + half_win_len + 1
            window = torch.zeros(fft_size)
            window[left: right] = 0.5 * torch.cos(math.pi * position * f0) + 0.5
            average = torch.sum(window * window).pow(0.5)
            self.window[f0] = (window / average)

    def forward(self, x, f):
        """Calculate forward propagation.
        Args:
            x (Tensor): Waveform (B, fft_size // 2 + 1, T).
            f (Tensor): F0 sequence (B, T').
        Returns:
            Tensor: Power spectrogram (B, bin_size, T').
        """
        # Get the matrix of window functions corresponding to F0
        x = self.zero_padding(x).unfold(1, self.fft_size, self.hop_size)
        windows = self.window[f]
        # Adaptive windowing and calculate power spectrogram.
        # In test, change x[:, : -1, :] to x.
        x = torch.abs(torch.fft.rfft(x[:, : -1, :] * windows)).pow(2)

        return x

class AdaptiveLiftering(nn.Module):

    def __init__(self,
                 sampling_rate,
                 fft_size,
                 f0_floor,
                 f0_ceil,
                 q1=-0.15,
        ):
        """Initilize AdaptiveLiftering module.
        Args:
            sampling_rate (int): Sampling rate.
            fft_size (int): FFT size.
            f0_floor (int): Minimum value of F0.
            f0_ceil (int): Maximum value of F0.
            q1 (float): Parameter to remove effect of adjacent harmonics.
        """
        super(AdaptiveLiftering, self).__init__()

        self.sampling_rate = sampling_rate
        self.bin_size = fft_size // 2 + 1
        self.q1 = q1
        self.q0 = 1.0 - 2.0 * q1
        self.smoothing_lifter = torch.zeros((f0_ceil+1, self.bin_size)).cuda()
        self.compensation_lifter = torch.zeros((f0_ceil+1, self.bin_size)).cuda()

        # Pre-calculation of the smoothing lifters and compensation lifters
        for f0 in range(f0_floor, f0_ceil + 1):
            smoothing_lifter = torch.zeros(self.bin_size)
            compensation_lifter = torch.zeros(self.bin_size)
            quefrency = torch.arange(1, self.bin_size) / sampling_rate
            smoothing_lifter[0] = 1.0
            smoothing_lifter[1:] = torch.sin(math.pi * f0 * quefrency) / (math.pi * f0 * quefrency)
            compensation_lifter[0] = self.q0 + 2.0 * self.q1
            compensation_lifter[1:] = self.q0 + 2.0 * self.q1 * torch.cos(2.0 * math.pi * f0 * quefrency)
            self.smoothing_lifter[f0] = smoothing_lifter
            self.compensation_lifter[f0] = compensation_lifter

    def forward(self, x, f):
        """Calculate forward propagation.
        Args:
            x (Tensor): Power spectrogram (B, bin_size, T').
            f (Tensor): F0 sequence (B, T').
        Returns:
            Tensor: Estimated spectral envelope (B, bin_size, T').
        """
        # Setting the smoothing lifter and compensation lifter
        smoothing_lifter = self.smoothing_lifter[f]
        compensation_lifter = self.compensation_lifter[f]
        # Calculating cepstrum
        tmp = torch.cat((x, torch.flip(x[:, :, 1:-1], [2])), dim=2)
        cepstrum = torch.fft.rfft(
            torch.log(torch.clamp(tmp, min=1e-7))
            ).real
        # Liftering cepstrum with the lifters
        liftered_cepstrum = cepstrum * smoothing_lifter * compensation_lifter
        # Return the result to the spectral domain
        x = torch.fft.irfft(liftered_cepstrum)[:, :, : self.bin_size]

        return x

class CheapTrick(nn.Module):

    def __init__(self,
                 sampling_rate,
                 hop_size,
                 fft_size,
                 f0_floor,
                 f0_ceil,
                 uv_threshold=0,
                 q1=-0.15,
        ):
        """Initilize AdaptiveLiftering module.
        Args:
            sampling_rate (int): Sampling rate.
            hop_size (int): Hop size.
            fft_size (int): FFT size.
            f0_floor (int): Minimum value of F0.
            f0_ceil (int): Maximum value of F0.
            uv_threshold (float): V/UV determining threshold.
            q1 (float): Parameter to remove effect of adjacent harmonics.
        """
        super(CheapTrick, self).__init__()

        # fft_size must be larger than 3.0 * sampling_rate / f0_floor
        assert fft_size > 3.0 * sampling_rate / f0_floor
        self.f0_floor = f0_floor
        self.f0_ceil = f0_ceil
        self.uv_threshold = uv_threshold

        self.ada_wind = AdaptiveWindowing(
                            sampling_rate,
                            hop_size,
                            fft_size,
                            f0_floor,
                            f0_ceil,
        )
        self.ada_lift = AdaptiveLiftering(
                            sampling_rate,
                            fft_size,
                            f0_floor,
                            f0_ceil,
                            q1,
        )

    def forward(self, x, f):
        """Calculate forward propagation.
        Args:
            x (Tensor): Power spectrogram (B, T).
            f (Tensor): F0 sequence (B, T').
        Returns:
            Tensor: Estimated spectral envelope (B, bin_size, T').
        """
        # Step0: Round F0 values to integers.
        voiced = (f > self.uv_threshold) * torch.ones_like(f)
        f = voiced * f + (1 - voiced) * self.f0_ceil
        f = torch.round(
                torch.clamp(f, min=self.f0_floor, max=self.f0_ceil)
            ).to(torch.int64)
        # Step1: Adaptive windowing and calculate power spectrogram.
        x = self.ada_wind(x, f)
        # Step3: Smoothing (log axis) and spectral recovery on the cepstrum domain.
        x = self.ada_lift(x, f)

        return x

if __name__ == "__main__":
    """Test of spectral envelope extraction."""
    import numpy as np
    import pyworld as pw
    import soundfile as sf
    import librosa.display
    import matplotlib.pyplot as plt

    config = {
        'sampling_rate': 16000,
        'hop_size': 80,
        'fft_size': 1024,
        'f0_floor': 50,
        'f0_ceil': 500
    }

    device = 'cuda' if torch.cuda.is_available() else 'cpu'
    cheaptrick = CheapTrick(**config)
    cheaptrick.to(device)

    file_name = "../../egs/arctic/data/wav/arctic_evaluation/bdl/bdl_arctic_b0474.wav"
    x, sr = sf.read(file_name)
    x = x[:config['sampling_rate']]
    _f0, t = pw.dio(x, config['sampling_rate'], frame_period=config['hop_size'] * 1000 / config['sampling_rate'])
    f0 = pw.stonemask(x, _f0, t, config['sampling_rate'])
    ap = pw.d4c(x, f0, t, config['sampling_rate'])

    x = torch.from_numpy(np.array(x[np.newaxis, :])).clone().to(device)
    f0 = torch.from_numpy(np.array(f0[np.newaxis, :])).clone().to(device)
    sp = torch.exp(cheaptrick.forward(x, f0))
    sp = sp.to('cpu').numpy().copy()[0]
    f0 = f0.to('cpu').numpy().copy()[0]

    # confirm whether the signal is resynthesized properly
    y = pw.synthesize(f0, sp, ap, config['sampling_rate'], config['hop_size'] * 1000 / config['sampling_rate'])
    save_name = 'resynthesized.wav'
    sf.write(save_name, y, config['sampling_rate'])

    # confirm whether reasonable spectral envelopes are extracted
    sp_db = librosa.power_to_db(sp)
    librosa.display.specshow(data=sp_db.T, sr=config['sampling_rate'],
                             hop_length=config['hop_size'], y_axis='linear', x_axis='time')
    plt.colorbar(format="%+2.f dB")
    save_name = 'spectrogram.png'
    plt.savefig(save_name)

max_stars_repo_name: chomeyama/UnifiedSourceFilterGAN

max_issues_repo_name: chomeyama/UnifiedSourceFilterGAN

max_forks_repo_name: chomeyama/UnifiedSourceFilterGAN


text:
import tensorflow as tf
import numpy as np
from PIL import Image

with open('mnist_200.pb', 'rb') as f:
    a = tf.parse_tensor(f.read(), tf.float32)

with tf.Session() as sess:
    b = a.eval(session=sess)

b = np.squeeze(b)*256
b = b.astype(np.uint8)

img = Image.fromarray(b, 'L')
img.show()


text:
#! /usr/bin/env python
# Copyright 2014 Uri Laserson
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import sys
import optparse

# HACK
# importing vdj.clustering will import scipy.cluster, which will import pylab,
# if it's available.  By default, pylab will try to open some kind of
# windowing device.  If I dispatch jobs to the batch system, this causes an
# error.  But if I preempt this import by calling matplotlib and connecting it
# to the non-windowing Agg backend, it should avoid the problem.
import matplotlib as mpl
mpl.use('agg')

import vdj
import vdj.clustering

parser = optparse.OptionParser()
parser.add_option('-c','--cutoff',type='float',default=4.5)
parser.add_option('-t','--tag',default='')
parser.add_option('-l','--linkage',type='choice',choices=['single','complete'],default='single')
(options, args) = parser.parse_args()

if len(args) == 2:
    inhandle = open(args[0],'r')
    outhandle = open(args[1],'w')
elif len(args) == 1:
    inhandle = open(args[0],'r')
    outhandle = sys.stdout
elif len(args) == 0:
    inhandle = sys.stdin
    outhandle = sys.stdout

# NOTE: this script requires there to be a well-defined junction
#       sequence.  It raises an exception if not.  Therefore, seqs
#       must be pre-filtered for having legit junctions
# NOTE: this script must hold all chains in memory in order to
#       perform the clustering and then assign cluster names

# load data
chains = []
junctions = []
for chain in vdj.parse_imgt(inhandle):
    # check for presence of V, J, and non-trivial junction
    if not hasattr(chain,'v') or not hasattr(chain,'j') or not hasattr(chain,'junction'):
        raise ValueError, "Chain %s has no junction of V-J aln." % chain.descr
    chains.append(chain)
    junctions.append(chain.junction)

# perform the sequence clustering
(T,seq_idxs) = vdj.clustering.cluster_seqs(junctions,options.cutoff,options.linkage)

# tag chains with unique cluster IDs
if options.tag == '':
    tag = ''
else:
    tag = options.tag+'|'

for (i,chain) in enumerate(chains):
    cloneID = '%s%s' % (tag,T[seq_idxs[chain.junction]])
    chain.annotations['clone'] = cloneID
    print >>outhandle, chain

max_stars_repo_name: laserson/vdj

max_issues_repo_name: laserson/vdj

max_forks_repo_name: churchlab/vdj

reason: import scipy


text:
import nrrd
import tensorflow as tf
import os
import numpy as np

flags = tf.app.flags
FLAGS = flags.FLAGS
flags.DEFINE_string('directory', '../data/gan/V24/', 'Data directory.')
flags.DEFINE_string('saving_directory', '../data/gan/tfrecords/', 'Data directory.')
flags.DEFINE_string('filename', 'train.tfrecords', 'filename')

filename = FLAGS.filename
directory = FLAGS.directory
saving_directory = FLAGS.saving_directory

def _bytes_feature(value):
    return tf.train.Feature(bytes_list=tf.train.BytesList(value=[value]))

def _int64_feature(value):
    return tf.train.Feature(int64_list=tf.train.Int64List(value=[value]))

def get_im(path):
    images = sorted(os.listdir(path))
    nb_images = int(len(images))
    data = np.ndarray(shape=(int(nb_images/2), 195, 233, 2), dtype=np.float32)

    num = 0
    ni = 0

    for ima in images:
        ima = os.path.join(path, ima)
        # print(ima)
        # print('ni : ', ni)
        # print('num image : ', num)

        im_nrrd, _ = nrrd.read(ima)
        im_nrrd = im_nrrd.astype(np.float32)
        im_nrrd = im_nrrd.reshape((195,233))

        data[num, :, :, ni] = im_nrrd

        if ni == 0: ni = 1
        else:
            ni = 0
            num += 1
    return data, int(nb_images/2)

def get_la(path):
    labs = sorted(os.listdir(path))
    nb_labs = int(len(labs))
    data = np.ndarray(shape=(int(nb_labs), 195, 233, 1), dtype=np.float32)

    num = 0

    for ima in labs:
        ima = os.path.join(path, ima)
        # print(ima)
        # print('num label : ', num)

        im_nrrd, _ = nrrd.read(ima)
        im_nrrd = im_nrrd.astype(np.float32)
        im_nrrd = im_nrrd.reshape((195, 233))

        data[num, :, :, 0] = im_nrrd
        num += 1

    return data

im, nb_im = get_im(os.path.join(directory, 'V06'))
la = get_la(os.path.join(directory, 'V24_t1'))

writer = tf.python_io.TFRecordWriter(os.path.join(saving_directory, filename))

for i in xrange(nb_im):
    img = im[i]
    labl = la[i]

    # The reason to store image sizes was demonstrated
    # in the previous example -- we have to know sizes
    # of images to later read raw serialized string,
    # convert to 1d array and convert to respective
    # shape that image used to have.
    height = img.shape[0]
    width = img.shape[1]
    depth = img.shape[2]

    img_raw = img.tostring()
    labl_raw = labl.tostring()

    example = tf.train.Example(features=tf.train.Features(feature={
        'height': _int64_feature(height),
        'width': _int64_feature(width),
        'depth': _int64_feature(depth),
        'V24_raw': _bytes_feature(labl_raw),
        'V06_raw': _bytes_feature(img_raw)}))

    writer.write(example.SerializeToString())

writer.close()

def read_and_decode(filename_queue):
    reader = tf.TFRecordReader()

    _, serialized_example = reader.read(filename_queue)

    features = tf.parse_single_example(
        serialized_example,
        # Defaults are not specified since both keys are required.
        features={
            'height': tf.FixedLenFeature([], tf.int64),
            'width': tf.FixedLenFeature([], tf.int64),
            'depth': tf.FixedLenFeature([], tf.int64),
            'V24_raw': tf.FixedLenFeature([], tf.string),
            'V06_raw': tf.FixedLenFeature([], tf.string)
        })

    # Convert from a scalar string tensor (whose single string has
    # length mnist.IMAGE_PIXELS) to a uint8 tensor with shape
    # [mnist.IMAGE_PIXELS].
    image = tf.decode_raw(features['V06_raw'], tf.float32)
    label = tf.decode_raw(features['V24_raw'], tf.float32)

    # Set image and label shapes
    height = tf.cast(features['height'], tf.int32)
    width = tf.cast(features['width'], tf.int32)
    depth = tf.cast(features['depth'], tf.int32)

    image_shape = tf.pack([195, 233, 2])
    label_shape = tf.pack([195, 233,1])

    image = tf.reshape(image, image_shape)
    label = tf.reshape(label, label_shape)

    # Random transformations can be put here: right before you crop images
    # to predefined size. To get more information look at the stackoverflow
    # question linked above.
    images, labels = tf.train.shuffle_batch([image, label],
                                            batch_size=3,
                                            capacity=300,
                                            num_threads=2,
                                            min_after_dequeue=5)

    return images, labels

filename_queue = tf.train.string_input_producer(
    [os.path.join(directory,filename)], num_epochs=1)

# Even when reading in multiple threads, share the filename
# queue.
image, label = read_and_decode(filename_queue)

# The op for initializing the variables.
init_op = tf.group(tf.initialize_all_variables(),
                   tf.initialize_local_variables())

with tf.Session() as sess:
    sess.run(init_op)

    coord = tf.train.Coordinator()
    threads = tf.train.start_queue_runners(coord=coord)

    # Let's read off 3 batches just for example
    for i in xrange(3):
        img, lab = sess.run([image, label])
        print(img[0, :, :, :].shape)
        print(lab.shape)
        print(np.amin(lab))
        print('current batch')

        # We selected the batch size of two
        # So we should get two image pairs in each batch
        # Let's make sure it is random

        # nrrd.write('prediction_a_' + i.__str__() + '.nrrd', img[0, :, :, 0])
        # nrrd.write('l_a_' + i.__str__() + '.nrrd', lab[0, :, :, 0])

    coord.request_stop()
    coord.join(threads)

max_stars_repo_name: juanprietob/gan-brain

max_issues_repo_name: juanprietob/gan-brain

max_forks_repo_name: juanprietob/gan-brain


text:
#!/usr/bin/env python

import os, sys
# Add parent of parent folder to path
sys.path.append(os.path.dirname(os.path.dirname(os.path.realpath(__file__))))

import argparse
import pickle
import inspect
import pathlib
from datetime import datetime

import numpy as np
import matplotlib.pyplot as plt

import torch
import torch.optim as optim
import torch.nn.functional as F

from ignite.engine import Events, create_supervised_trainer, create_supervised_evaluator
from ignite.metrics import Loss, RunningAverage
from ignite.handlers import ModelCheckpoint
from ignite.contrib.handlers.tqdm_logger import ProgressBar

from ignite.contrib.handlers.tensorboard_logger import (
    TensorboardLogger, OutputHandler, OptimizerParamsHandler, WeightsHistHandler)

from termcolor import colored
from tqdm import tqdm

from utils.utilities import (create_folder, load_hdf5, scale, np_mean_absolute_error,
    pad_with_border, log_sp, mat_2d_to_3d, inverse_scale, get_stft_window_func,
    write_audio, read_audio, calculate_spectrogram, load_features, Range)
import models
import utils.config as config
from utils.stft import real_to_complex, istft, get_cola_constant, overlap_add

from audio_datasets import NoisySpeechFeaturesDataset
import wandb

def get_models():
    ''' Get dict 'name':value of all models by inspecting models.py '''
    # is_module uses short circuiting to check if it is a class _before_ checking
    # if it is a subclass, as issubclass throws exception if passed something that
    # is _not_ a class
    is_module = lambda x : inspect.isclass(x) and issubclass(x,torch.nn.Module)
    # Returns [(key,value), (key,value)], so turn into dict
    models_list =  inspect.getmembers(models,predicate=is_module)
    return dict(models_list)

def get_model_kwargs(model_class: torch.nn.Module):
    ''' Get keyword arguments of model, handles float or int based on default
        Returns dict(keyword, default)
    '''
    kwargs = dict()
    sig = inspect.signature(model_class)
    for name,p in sig.parameters.items():
        if p.kind == inspect.Parameter.KEYWORD_ONLY:
            assert(p.default is not inspect.Parameter.empty)
            kwargs[name] = p.default
    return kwargs

def model_kwargs_to_argparser(kwargs):
    ''' Create argument parser from output of get_model_kwargs '''
    parser = argparse.ArgumentParser()
    for name,default in kwargs.items():
        parser.add_argument(f'--{name}',type=type(default), default=default)
    return parser

def plot_spectrogram(x,y,n_concat):
    ''' Given Tensors x,y, returns a figure of the first 1000 items '''
    fig, axs = plt.subplots(2,1, sharex=True)
    axs[0].matshow(x[0 : 1000, n_concat // 2, :].numpy().T, origin='lower', aspect='auto', cmap='jet')
    axs[1].matshow(y[0 : 1000, :].numpy().T, origin='lower', aspect='auto', cmap='jet')
    axs[0].set_title('Mixture')
    axs[1].set_title('Clean')
    plt.tight_layout()
    return fig

def split_train_val_sequential(dataset: torch.utils.data.Dataset, ratio):
    train_len = len(dataset)
    validation_len = int(ratio * train_len)
    train_len = train_len - validation_len
    val_idx = list(range(train_len,len(dataset)))
    train_dataset = torch.utils.data.Subset(dataset, list(range(train_len)))
    val_dataset = torch.utils.data.Subset(dataset, val_idx)
    print("Split dataset, {} train {} validation".format(train_len, validation_len))
    return train_dataset, val_dataset

def train(args, model_vars):

    # Arugments & parameters
    workspace = args.workspace
    tr_snr = args.tr_snr
    te_snr = args.te_snr
    batch_size = args.batch_size
    lr = args.lr
    device = 'cpu'
    validate_every = 100
    # We already "shuffled" when choosing audio files.
    # Shuffling again would break any kind of RNN from working, as this would
    # randomly rearrange each spectrogram.
    shuffle_data = False
    # Disable parallel GPU usage as breaks sequential_sampler
    use_parallel = False

    # Initialise W&B, ensuring we remove the tricky "model" attribute from args
    safe_config = vars(args).copy() # Convert to dictionary and get a copy
    safe_config.pop('model') # Remove model attribute
    safe_config.pop('mode') # Remove needless mode attribute, will always be train
    safe_config.pop('loss') # Remove loss attribute
    safe_config['shuffle_data'] = shuffle_data
    safe_config.update(model_vars) # Add model hyperparams
    wandb.init(project="mss_speech_sep_lstm",config=safe_config,sync_tensorboard=True)

    if torch.cuda.is_available():
        device = 'cuda'
    else:
        print(colored('No GPU Found!','red'))

    # Paths
    train_hdf5_path = os.path.join(workspace, 'packed_features', 'spectrogram',
        'train', '{}db'.format(int(tr_snr)), 'data.h5')

    test_hdf5_path = os.path.join(workspace, 'packed_features', 'spectrogram',
        'test', '%ddb' % int(te_snr), 'data.h5')

    scaler_path = os.path.join(workspace, 'packed_features', 'spectrogram',
        'train', '{}db'.format(int(tr_snr)), 'scaler.p')

    models_dir = os.path.join(workspace, 'models', '{}db'.format(int(tr_snr)))

    # Save model to cloud when done
    wandb.save(str(models_dir) + '*.pth')

    # Datasets & dataloader
    print("Loading Dataset, Scaling data...")
    train_dataset = NoisySpeechFeaturesDataset(train_hdf5_path, scaler_path)
    test_dataset = NoisySpeechFeaturesDataset(train_hdf5_path, scaler_path)

    # Get properties embedded in dataset, e.g. n_concat
    data_prop = train_dataset.get_properties()

    # Split train_dataset into train and validation set
    train_dataset, validation_dataset = split_train_val_sequential(
                                            train_dataset,
                                            args.validation_split)

    # Load Model and pass in properties of dataset
    print(f"Using Model {args.model_name}")
    model = args.model(**data_prop, **model_vars)

    if use_parallel==True and torch.cuda.device_count() > 1:
        print(colored('Using {} GPUs'.format(torch.cuda.device_count()),'green'))
        model = torch.nn.DataParallel(model).to(device)

    # Data loader
    # TODO: setting num_workers > 1 breaks CUDA? Needs restart to fix!
    # See https://github.com/pytorch/pytorch/issues/2517
    print("Dataloader...")

    if shuffle_data:
        sampler = None
    else:
        sampler = torch.utils.data.SequentialSampler

    train_loader = torch.utils.data.DataLoader(train_dataset,
                                            batch_size=batch_size,
                                            shuffle=shuffle_data,
                                            sampler=sampler(train_dataset),
                                            pin_memory=True)
    validation_loader = torch.utils.data.DataLoader(validation_dataset,
                                            batch_size=batch_size,
                                            shuffle=shuffle_data,
                                            sampler=sampler(validation_dataset),
                                            pin_memory=True)
    test_loader = torch.utils.data.DataLoader(test_dataset,
                                            batch_size=batch_size,
                                            shuffle=shuffle_data,
                                            sampler=sampler(test_dataset),
                                            pin_memory=True)

    # Optimizer
    optimizer = optim.Adam(model.parameters(), lr=lr)

    # Loss function, specified by args
    loss = args.loss

    # Setup trainer and evaluator
    trainer = create_supervised_trainer(model, optimizer, loss, device=device)
    evaluator = create_supervised_evaluator(model,
                                            device=device,
                                            metrics={
                                                # Accuracy() # TODO find replacement for accuracy
                                                'loss': Loss(loss)
                                                })

    # TODO make difference between training and test loss clear

    # Progress bar attach
    RunningAverage(output_transform=lambda x: x).attach(trainer, 'loss')
    pbar = ProgressBar(persist=True)
    pbar.attach(trainer, ['loss'])

    # Tensorboard attach, training loss and optimizer params
    date_str = datetime.now().strftime("%H:%M:%S %d:%m:%y")
    tb_logger = TensorboardLogger(os.path.join(workspace, f'tensorboard/{args.model_name}/{date_str}'))

    # Add graph to tensorboard
    # FIXME Doesnt work
    # freq_bins = config.window_size // 2 + 1
    # tb_logger.writer.add_graph(model, torch.zeros(5,5,freq_bins), True)

    # Plot spectrogram before/after cleaning
    # TODO plot after!
    # n_concat = data_prop['n_concat']
    # fig = plot_spectrogram(train_dataset.x,train_dataset.y,n_concat)
    # tb_logger.writer.add_figure('Spectrogram',fig,global_step=0,close=True)

    # Attach tensorboard loggers
    tb_logger.attach(trainer,
                 log_handler=OutputHandler(tag="training", output_transform=lambda loss: {'loss': loss}),
                 event_name=Events.ITERATION_COMPLETED)

    tb_logger.attach(trainer,
                 log_handler=OptimizerParamsHandler(optimizer),
                 event_name=Events.ITERATION_STARTED)

    tb_logger.attach(trainer,
                log_handler=WeightsHistHandler(model),
                event_name=Events.EPOCH_COMPLETED)

    @trainer.on(Events.ITERATION_COMPLETED)
    def validate(trainer):
        ''' Show validation_loss results every validate_every '''
        if trainer.state.iteration % validate_every == 0:
            evaluator.run(validation_loader)
            metrics = evaluator.state.metrics
            # Log to tensorboard
            tb_logger.writer.add_scalar('training/val_loss', metrics['loss'],
                                        global_step=trainer.state.iteration)
            # Log to w&b
            wandb.log({"val_loss": metrics['loss']}, step=trainer.state.iteration)

    @trainer.on(Events.EPOCH_COMPLETED)
    def log_test_results(trainer):
        ''' Show test results every epoch '''
        state = evaluator.run(test_loader)
        #pylint: disable=no-member
        metrics = state.metrics
        loss = metrics['loss']
        # Show in progress bar
        pbar.log_message("Test Set Results - Epoch: {}  Avg loss: {:.2f}"
            .format(trainer.state.epoch, loss))
        # Show in train set
        wandb.log({"Test Loss": loss}, step=trainer.state.iteration)

    # Handler to save checkpoints
    chkpoint_handler = ModelCheckpoint(models_dir,
                                        filename_prefix=f"{args.model_name}_chkpoint_",
                                        save_interval=2, # Save every 2nd epoch
                                        require_empty=False # Overwrite
                                        )
    trainer.add_event_handler(Events.EPOCH_COMPLETED, chkpoint_handler,
                                {
                                    'ig_model': model,
                                    'ig_optimizer': optimizer
                                })

    # Watch model with W&B, logging both gradients and parameters
    wandb.watch(model, log='all')

    print(f"Starting training {args.model_name}")
    try:
        trainer.run(train_loader, max_epochs=args.max_epochs)
    except KeyboardInterrupt:
        print(colored("Interrupted, exiting..",'red'))
    else:
        # Print done
        print(colored('FINISHED','green'))
    finally:
        # We need to close the logger with we are done
        tb_logger.close()

    # TODO run inference and upload PESQ, STOI results

def infer_test_audio(device, model, features_dir, audio_name, n_concat, scaler):
    ''' Run inference on an audio file using model '''

    window_size = config.window_size
    freq_bins = window_size // 2 + 1
    window = get_stft_window_func(config.window_type)(window_size)
    overlap = config.overlap
    hop_size = window_size - overlap

    # Load features
    feature_path = os.path.join(features_dir, audio_name)
    data = load_features(feature_path)
    mixed_complx_x, speech_x, noise_x, alpha, mrcg = data

    # Convert complex input x to real
    mixed_x = np.abs(mixed_complx_x)

    # Pad and log data
    n_pad = (n_concat - 1) // 2
    mixed_x = pad_with_border(mixed_x, n_pad)
    mixed_x = log_sp(mixed_x)
    speech_x = log_sp(speech_x)

    # Scale data
    mixed_x = scale(mixed_x, scaler)

    # Cut input spectrogram to 3D segments with n_concat
    mixed_x_3d = mat_2d_to_3d(mixed_x, agg_num=n_concat, hop=1)

    # Covert data to torch.Tensor
    mixed_x_3d = torch.Tensor(mixed_x_3d).to(device)

    # Predict
    prediction = model(mixed_x_3d)

    # Move back to CPU
    prediction = prediction.data.cpu().numpy()

    # Inverse scale
    mixed_x = inverse_scale(mixed_x, scaler)
    prediction = inverse_scale(prediction, scaler)

    # Recover enhanced wav
    prediction_sp = np.exp(prediction)
    complex_sp = real_to_complex(prediction_sp, mixed_complx_x)
    frames = istft(complex_sp)

    # Overlap add
    cola_constant = get_cola_constant(hop_size, window)
    enh_audio = overlap_add(frames, hop_size, cola_constant)

    return enh_audio

def inference(args, model_vars):
    """Inference all test data, write out recovered wavs to disk.

    Args:
      workspace: str, path of workspace.
      tr_snr: float, training SNR.
      te_snr: float, testing SNR.
      n_concat: int, number of frames to concatenta, should equal to n_concat
          in the training stage.
      iter: int, iteration of model to load.
    """
    # This function is from original sednn codebase, lightly modified

    # Arguments & parameters
    workspace = args.workspace
    tr_snr = args.tr_snr
    te_snr = args.te_snr
    n_concat = args.n_concat
    data_type = 'test'
    device = 'cuda' if torch.cuda.is_available() else 'cpu'
    sample_rate = config.sample_rate
    window_size = config.window_size
    freq_bins = window_size // 2 + 1

    # Paths
    mixed_audios_dir = os.path.join(workspace, 'mixed_audios', data_type,
        '{}db'.format(int(te_snr)))

    features_dir = os.path.join(workspace, 'features', 'spectrogram', data_type,
        '{}db'.format(int(te_snr)))

    scaler_path = os.path.join(workspace, 'packed_features', 'spectrogram',
        'train', '{}db'.format(int(tr_snr)), 'scaler.p')

    # TODO model path changes when number of epochs change
    model_path = os.path.join(workspace, 'models', '{}db'.format(int(tr_snr)),
        f'{args.model_name}_chkpoint__ig_model_10.pth')

    enh_audios_dir = args.enhanced_dir

    # Load model
    model = args.model(n_concat, freq_bins, **model_vars)
    state_dict = torch.load(model_path)

    # Hack to remove results of DataParallel
    # See https://discuss.pytorch.org/t/solved-keyerror-unexpected-key-module-encoder-embedding-weight-in-state-dict/1686/4
    # TODO detect if required, e.g. might be trained on only 1 GPU
    if "ght" in state_dict:
        from collections import OrderedDict
        new_state_dict = OrderedDict()
        for k, v in state_dict.items():
            name = k[7:] # remove `module.`
            new_state_dict[name] = v
        model.load_state_dict(new_state_dict)
    else:
        model.load_state_dict(state_dict)

    # Move model to device if needed
    model.to(device)

    # Load scaler
    scaler = pickle.load(open(scaler_path, 'rb'))

    # Get list of audio names
    audio_names = os.listdir(features_dir)

    # Create output folder
    create_folder(enh_audios_dir)

    for name in tqdm(audio_names):
        # Infer the audio
        enh_audio = infer_test_audio(device, model, features_dir, name, n_concat, scaler)

        # Write out enhanced wav
        bare_name = os.path.splitext(name)[0]
        out_path = os.path.join(enh_audios_dir, '{}.enh.wav'.format(bare_name))
        write_audio(out_path, enh_audio, sample_rate)

if __name__ == '__main__':

    models = get_models() # Get dict of models in models.py

    loss_functions = {
        'L1': F.l1_loss,
        'MSE': F.mse_loss
    }

    parser = argparse.ArgumentParser(description='Run training/inference',
                                    # Can't use abbrev for parse_known_args()
                                     allow_abbrev=False)

    parser.add_argument("model_name", help="Model from models.py to use",
                        choices=models.keys())

    subparsers = parser.add_subparsers(dest='mode')

    parser_train = subparsers.add_parser('train')
    parser_train.add_argument('--workspace', type=str, required=True)
    parser_train.add_argument('--tr_snr', type=float, required=True)
    parser_train.add_argument('--te_snr', type=float, required=True)
    parser_train.add_argument('--lr', type=float, default=1e-5)
    parser_train.add_argument('--batch_size', type=int, default=1000)
    parser_train.add_argument('--max_epochs', type=int, default=10)
    parser_train.add_argument('--loss_func', dest="loss_name",
                              choices=loss_functions, default='MSE')
    parser_train.add_argument('--validation_split', type=float, default=0.1,
                              choices=[Range(0.0, 0.2)])

    parser_inference = subparsers.add_parser('inference')
    parser_inference.add_argument('--workspace', type=str, required=True)
    parser_inference.add_argument('--enhanced_dir', type=str, required=True)
    parser_inference.add_argument('--tr_snr', type=float, required=True)
    parser_inference.add_argument('--te_snr', type=float, required=True)
    parser_inference.add_argument('--n_concat', type=int, required=True)

    args, remaining_args = parser.parse_known_args()

    # Add model to arguments
    args.model = models[args.model_name]

    # Add extra model arguments if given
    extra_kwargs = get_model_kwargs(args.model)
    extra_parser = model_kwargs_to_argparser(extra_kwargs)
    model_vars = vars(extra_parser.parse_args(remaining_args))

    if args.mode == 'train':
        args.loss = loss_functions[args.loss_name]
        train(args, model_vars)
    elif args.mode == 'inference':
        inference(args, model_vars)
    else:
        raise Exception('Error argument!')


text:
import numpy as np
import matplotlib.pyplot as pl

def load_data():
    """
    loads training and testing data
    """

    cls = ['3', '4', '6'] # 3: Happy, 4: Sad, 6: Neutral
    cls_names = ['Happy', 'Sad', 'Neutral']
    # training data
    loc = '../data/train/'
    train_K = 500
    test_K = 500
    x_training = np.zeros((train_K * len(cls), 24 * 24), dtype=np.float)
    y_training = np.zeros((train_K * len(cls), len(cls)), dtype=np.float)
    k = 0
    for n in np.arange(train_K):
        for c in np.arange(len(cls)):
            # load image
            I = pl.imread(loc + cls[c] + '_' + str(n) + '.jpg')[:,:,0]
            x_training[k,:] = np.reshape(I, (I.shape[0] * I.shape[1]))
            y_training[k,c] = 1.0
            k = k + 1
    # add all one column to the first axis
    x_training = np.hstack((np.ones((x_training.shape[0],1)), x_training))

    # testing data
    loc = '../data/test/'
    x_testing = np.zeros((test_K * len(cls), 24 * 24), dtype=np.float)
    y_testing  = np.zeros((test_K * len(cls), len(cls)), dtype=np.float)
    k = 0
    for n in np.arange(test_K):
        for c in np.arange(len(cls)):
            # load image
            I = pl.imread(loc + cls[c] + '_' + str(n) + '.jpg')[:,:,0]
            x_testing[k,:] = np.reshape(I, (I.shape[0] * I.shape[1]))
            y_testing[k,c] = 1.0
            k = k + 1
    # add all one column to the first axis
    x_testing = np.hstack((np.ones((x_testing.shape[0],1)), x_testing))
    return x_training, y_training, x_testing, y_testing

def normalization_parameters(x):
    """
    computes and returns normalization parameters (mean and std vectors of each
    column) from data x where it is assumed that the very first column is all
    1s
    """
    mean_vector = np.mean(x, axis=0)
    mean_vector[0] = 0
    std_vector = np.std(x, axis=0)
    std_vector[0] = 1.0
    return mean_vector, std_vector

def normalize_data(x, mean_vector, std_vector):
    """
    normalizes the input data x by parameters mean_vector and std_vector
    """
    x = x - mean_vector
    x = x / std_vector
    return x

def act(s):
    """
    activation function
    """
    u = 1.0 / (1.0 + np.exp(-s))
    return u

def is_converged(theta, theta_p, epsilon=1e-8):
    """
    checks convergence
    theta: theta at iteration t
    thetap: theta at iteration t-1
    epsilon: a small number used for thresholdin
    """
    L = len(theta)
    d = np.max(np.abs(theta[0]-theta_p[0]))
    for l in np.arange(L):
        d_temp = np.max(np.abs(theta[l]-theta_p[l]))
        if d_temp > d:
            d = d_temp
    if  d <= epsilon:
        return True
    else:
        return False

def ffnn_learn(x, y, alpha, theta):
    """
    learns and returns the parameters of ffnn using SGD
    """
    N = x.shape[0]
    L = len(theta)
    converged = False
    epoch = 0
    while not converged:
        mse = 0.0
        for n in np.arange(N): # for each observation
            xn = x[n, :][np.newaxis, :]
            yn = y[n, :][np.newaxis, :]
            s = dict()
            u = dict()
            s[0] = xn.copy()
            u[0] = s[0].copy()
            for l in np.arange(1, L+1): # apply forward propagation
                s[l] = np.dot(u[l-1], theta[l-1])
                u[l] = act(s[l])
                if l != L:
                    u[l][0,0] = 1.0
            en = u[l] - yn # error of the network
            mse += np.sum(en**2)
            theta_p = theta.copy()
            for l in np.arange(L-1, -1, -1): # back propagation
                delta = (u[l+1] * (1.0 - u[l+1]) * en)
                theta[l] = theta[l] - alpha * np.dot(u[l].T, delta)
                if l != L-1:
                    theta[l][0,0] = 0.0
                # accumulate the error at the previous layer
                en = np.dot(en, theta[l].T)
            converged = is_converged(theta, theta_p)
            if converged:
                break;
        epoch += 1
        mse = mse / n
        msg = 'epoch: {:7d}, MSE: {:.6f}'.format(epoch, mse)
        print(msg)
    return theta

def ffnn_classify(x, theta):
    """
    classifies x according ffnn architecture
    """
    L = len(theta)
    s = dict()
    u = dict()
    s[0] = x.copy()
    u[0] = s[0].copy()
    for l in np.arange(1, L+1): # apply forward propagation
        s[l] = np.dot(u[l-1], theta[l-1])
        u[l] = act(s[l])
        if l != L:
            u[l][0,0] = 1.0
    return u[l]


text:
import numpy as np
from sklearn.neighbors import KNeighborsClassifier as KNN
import matplotlib.pyplot as plt
from sklearn import svm
from sklearn.metrics import accuracy_score
from sklearn.metrics import roc_auc_score
from sklearn import preprocessing

class Solution():
    # structure of the solution
    def __init__(self):
        self.num_features = None
        self.num_agents = None
        self.max_iter = None
        self.obj_function = None
        self.execution_time = None
        self.convergence_curve = {}
        self.best_agent = None
        self.best_fitness = None
        self.best_accuracy = None
        self.final_population = None
        self.final_fitness = None
        self.final_accuracy = None

class Data():
    # structure of the training data
    def __init__(self):
        self.train_X = None
        self.train_Y = None
        self.val_X = None
        self.val_Y = None

def initialize(num_agents, num_features):
    # define min and max number of features  #Jake modify 1109 原先最大0.6最小0.3
    max_features = int(0.8 * num_features)
    min_features = int(0.3 * num_features)  # max_features - 1  #  avoid 沒有顯著差異
    # initialize the agents with zeros
    agents = np.zeros((num_agents, num_features))

    # select random features for each agent
    for agent_no in range(num_agents):

        # find random indices
        cur_count = np.random.randint(min_features, max_features)
        temp_vec = np.random.rand(1, num_features)
        temp_idx = np.argsort(temp_vec)[0][0:cur_count]

        # select the features with the ranom indices
        agents[agent_no][temp_idx] = 1

    return agents

def sort_agents(agents, obj, data, fitness=None):
    # sort the agents according to fitness
    train_X, val_X, train_Y, val_Y = data.train_X, data.val_X, data.train_Y, data.val_Y
    (obj_function, weight_acc) = obj

    if fitness is None:
        # if there is only one agent
        if len(agents.shape) == 1:
            num_agents = 1
            fitness = obj_function(
                agents, train_X, val_X, train_Y, val_Y, weight_acc)
            return agents, fitness

        # for multiple agents
        else:
            num_agents = agents.shape[0]
            fitness = np.zeros(num_agents)
            for id, agent in enumerate(agents):
                fitness[id] = obj_function(
                    agent, train_X, val_X, train_Y, val_Y, weight_acc)

    idx = np.argsort(-fitness)
    sorted_agents = agents[idx].copy()
    sorted_fitness = fitness[idx].copy()

    return sorted_agents, sorted_fitness

def display(dataframe, agents, fitness, agent_name='Agent'):
    # display the population
    print('\nNumber of agents: {}'.format(agents.shape[0]))
    print('\n------------- Best Agent ---------------')
    print('Fitness: {}'.format(fitness[0]))
    print(agents[0])
    # Jake add 1109
    print([column[0]
          for column in zip(dataframe.columns, agents[0]) if column[1]])
    col = [column[0] for column in zip(
        dataframe.columns, agents[0]) if column[1]]  # Jake add 1109
    print('Number of Features: {}'.format(int(np.sum(agents[0]))))
    print('----------------------------------------\n')

    for id, agent in enumerate(agents):
        print('{} {} - Fitness: {}, Number of Features: {}'.format(agent_name,
              id+1, fitness[id], int(np.sum(agent))))

    print('================================================================================\n')
    return(col)

def compute_accuracy(agent, train_X, test_X, train_Y, test_Y):
    # compute classification accuracy of the given agents
    cols = np.flatnonzero(agent)
    if(cols.shape[0] == 0):
        return 0
    clf = svm.SVC(kernel='linear', decision_function_shape='ovr')

    train_data = train_X[:, cols]
    train_label = train_Y
    test_data = test_X[:, cols]
    test_label = test_Y

    clf.fit(train_data, train_label)
    acc = clf.score(test_data, test_label)

    test_predict = clf.predict(test_data)
    acc2 = accuracy_score(test_label, test_predict)

    # 多分類
    roc = multiclass_roc_auc_score(test_label, test_predict, average='macro')

    print("acc=", acc)
    print("accuracy score=", acc2)
    print("roc=", roc)
    return acc

def compute_fitness(agent, train_X, test_X, train_Y, test_Y, weight_acc=0.9):
    # compute a basic fitness measure
    if(weight_acc == None):
        weight_acc = 0.9

    weight_feat = 1 - weight_acc
    num_features = agent.shape[0]

    acc = compute_accuracy(agent, train_X, test_X, train_Y, test_Y)
    feat = (num_features - np.sum(agent))/num_features

    fitness = weight_acc * acc + weight_feat * feat
    return fitness

def Conv_plot(convergence_curve):
    # plot convergence curves
    num_iter = len(convergence_curve['fitness'])
    iters = np.arange(num_iter) + 1
    fig, axes = plt.subplots(1)
    fig.tight_layout(pad=5)
    fig.suptitle('Convergence Curves')

    axes.set_title('Convergence of Fitness over Iterations')
    axes.set_xlabel('Iteration')
    axes.set_ylabel('Avg. Fitness')
    axes.plot(iters, convergence_curve['fitness'])

    return fig, axes

def multiclass_roc_auc_score(y_test, y_pred, average="macro"):
    lb = preprocessing.LabelBinarizer().fit(y_test)
    y_test = lb.transform(y_test)
    y_pred = lb.transform(y_pred)
    return roc_auc_score(y_test, y_pred, average=average)


text:
import numpy as np

def mse(input_array: np.ndarray, reference_array: np.ndarray) -> float:
    r"""
    The mse function compute the mean square error of predication array.
        .. math::
            mse=\frac{1}{N}\sum_i^N (p_i-r_i)^2

    where mse is the mean square error measurement, p is the predication array, r is the reference array.
    Note:reference array shape must be equal to input array shape

    :param input_array: A numpy array of any shape
    :param reference_array: A numpy array of any shape
    :return: a floating point number that represent the mean square error measurement
    """
    return float(np.mean(np.power(input_array - reference_array, 2)))

def nmse(input_array: np.ndarray, reference_array: np.ndarray, epsilon: float = 0.00001) -> float:
    r"""
    The nmse function compute the normalized mean square error of predication array.
        .. math::
            nmse=\frac{1}{N}\sum_i^N \frac{(p_i-r_i)^2}{r_i^2+\epsilon}

    where nmse is the normalized mean square error measurement, p is the predication array, r is the reference array
    and epsilon is a floating point number fo numeric stability.
    Note:reference array shape must be equal to input array shape

    :param input_array: A numpy array of any shape
    :param reference_array: A numpy array of any shape
    :param epsilon: a floating point number fo numric stabiliy
    :return: a floating point number that represent the normalized mean square error measurement
    """
    return float(np.mean(np.power(input_array - reference_array, 2) / (epsilon + np.power(reference_array, 2))))

def rmse(input_array: np.ndarray, reference_array: np.ndarray) -> float:
    r"""
    The rmse function compute the mean square error of predication array.
        .. math::
            mse=\sqrt{\frac{1}{N}\sum_i^N (p_i-r_i)^2}

    where mse is the mean square error measurement, p is the predication array, r is the reference array.
    Note:reference array shape must be equal to input array shape

    :param input_array: A numpy array of any shape
    :param reference_array: A numpy array of any shape
    :return: a floating point number that represent the mean square error measurement
    """
    return float(np.sqrt(np.mean(np.power(input_array - reference_array, 2))))

max_stars_repo_name: haihabi/PyNNcml

max_issues_repo_name: haihabi/PyNNcml

max_forks_repo_name: haihabi/PyNNcml


text:
# -*- coding: utf-8 -*-
"""
Purpose: To load ARPES MAP data from the ZIP file produced by Scienta SES program.
Version: 20211003
@author: Pranab Das (GitHub: @pranabdas)
[data, energy, theta, phi] = load_map_data("map_data.zip")
"""
def load_ses_map(filename):
    import numpy as np
    import zipfile
    from io import BytesIO
    from urllib.request import urlopen

    # Check the path contains the zip file
    if (filename[-4:].lower() != '.zip'):
        print('Please provide a .zip file!')

    if (filename[:7]=='http://') or (filename[:8]=='https://'):
        web=True
    else:
        web=False

    if (web):
        try:
            filename = urlopen(filename)
            map_data = zipfile.ZipFile(BytesIO(filename.read()))
        except:
            print('Could not read url.')
    else:
        map_data = zipfile.ZipFile(filename, 'r')

    ini_file = map_data.open('viewer.ini')
    contents = ini_file.readlines()
    ini_file.close()

    for line in range(len(contents)):
        contents[line] = contents[line].decode().strip()

    for line in range(len(contents)):
        energy_length = contents[line].split("=")
        if energy_length[0].strip()=="width":
            energy_length = int(energy_length[-1])
            break

    for line in range(len(contents)):
        theta_length = contents[line].split("=")
        if theta_length[0].strip()=="height":
            theta_length = int(theta_length[-1])
            break

    for line in range(len(contents)):
        phi_length = contents[line].split("=")
        if phi_length[0].strip()=="depth":
            phi_length = int(phi_length[-1])
            break

    for line in range(len(contents)):
        map_data_path = contents[line].split("=")
        if map_data_path[0].strip()=="path":
            map_data_path = str(map_data_path[-1].strip())
            break

    fid = map_data.open(map_data_path, 'r')
    data = np.frombuffer(fid.read(), dtype=np.dtype(np.float32))
    data = np.float64(data)
    data = data.reshape((energy_length, theta_length, phi_length), order='F')
    fid.close()

    for line in range(len(contents)):
        energy_offset = contents[line].split("=")
        if energy_offset[0].strip()=="width_offset":
            energy_offset = float(energy_offset[-1])
            break

    for line in range(len(contents)):
        energy_delta = contents[line].split("=")
        if energy_delta[0].strip()=="width_delta":
            energy_delta = float(energy_delta[-1])
            break

    energy_end = energy_offset + (energy_length - 1)*energy_delta
    energy = np.linspace(energy_offset, energy_end, num = energy_length)

    for line in range(len(contents)):
        theta_offset = contents[line].split("=")
        if theta_offset[0].strip()=="height_offset":
            theta_offset = float(theta_offset[-1])
            break

    for line in range(len(contents)):
        theta_delta = contents[line].split("=")
        if theta_delta[0].strip()=="height_delta":
            theta_delta = float(theta_delta[-1])
            break

    theta_end = theta_offset + (theta_length - 1)*theta_delta
    theta = np.linspace(theta_offset, theta_end, num = theta_length)

    for line in range(len(contents)):
        phi_offset = contents[line].split("=")
        if phi_offset[0].strip()=="depth_offset":
            phi_offset = float(phi_offset[-1])
            break

    for line in range(len(contents)):
        phi_delta = contents[line].split("=")
        if phi_delta[0].strip()=="depth_delta":
            phi_delta = float(phi_delta[-1])
            break

    phi_end = phi_offset + (phi_length - 1)*phi_delta
    phi = np.linspace(phi_offset, phi_end, num = phi_length)

    return data, energy, theta, phi

max_stars_repo_name: thePranab/arpespythontools

max_issues_repo_name: thePranab/arpespythontools

max_forks_repo_name: thePranab/arpespythontools


text:
import os
import pdb
import pickle
import re
from collections import defaultdict
from itertools import chain

import chakin
import numpy as np
import scipy.sparse as sp
from gensim.models import KeyedVectors
from gensim.models.wrappers.fasttext import FastText
from sklearn.metrics import classification_report, f1_score
from sklearn.preprocessing import MultiLabelBinarizer
from tqdm import tqdm

# sequence operation
# =========================================================
def clean_str(string):
    string = re.sub(r"[^A-Za-z0-9(),!?\'\`]", " ", string)
    string = re.sub(r"\'s", " \'s", string)
    string = re.sub(r"\'ve", " \'ve", string)
    string = re.sub(r"n\'t", " n\'t", string)
    string = re.sub(r"\'re", " \'re", string)
    string = re.sub(r"\'d", " \'d", string)
    string = re.sub(r"\'ll", " \'ll", string)
    string = re.sub(r",", " , ", string)
    string = re.sub(r"!", " ! ", string)
    string = re.sub(r"\(", " \( ", string)
    string = re.sub(r"\)", " \) ", string)
    string = re.sub(r"\?", " \? ", string)
    string = re.sub(r"\s{2,}", " ", string)
    return string.strip().lower()

# read data from text file
# =========================================================
def  make_data_list(data, kind_of_data, tree_info, max_sen_len, vocab, catgy, article_id, useWords):
    data_list = []
    for line in tqdm(data,desc="Loading " + kind_of_data + " data"):
        tmp_dict = dict()
        line = line[:-1]
        tmp_dict['text'] = ' '.join(clean_str(' '.join(line.split("\t")[1].split(" "))).split(" ")[:useWords])
        [vocab[word] for word in tmp_dict['text'].split(" ")]
        tmp_dict['num_words'] = len(tmp_dict['text'].split(" "))
        max_sen_len = max(max_sen_len, tmp_dict['num_words'])
        tmp_dict['split'] = kind_of_data
        tmp_dict['hie_info'] = list(set([tree_info[cat] for cat in line.split("\t")[0].split(",")]))
        tmp_dict['catgy'] = [cat for cat in line.split("\t")[0].split(",")]
        [catgy[cat] for cat in line.split("\t")[0].split(",")]
        tmp_dict['id'] = str(article_id)
        article_id += 1
        data_list.append(tmp_dict)
        del tmp_dict
    return data_list, max_sen_len, vocab, catgy, article_id

# read data
# =========================================================
def data_load(train, valid, test, tree_info, use_words):
    vocab = defaultdict( lambda: len(vocab) )
    catgy = defaultdict( lambda: len(catgy) )
    article_id = 0
    max_sen_len = 0

    train_list, max_sen_len, vocab, catgy, article_id = make_data_list(train, 'train', tree_info, max_sen_len, vocab, catgy, article_id, use_words)
    valid_list, max_sen_len, vocab, catgy, article_id = make_data_list(valid, 'valid', tree_info, max_sen_len, vocab, catgy, article_id, use_words)
    test_list, max_sen_len, vocab, catgy, article_id = make_data_list(test, 'test', tree_info, max_sen_len, vocab, catgy, article_id, use_words)
    class_dim = len(catgy)

    data = {}
    data['train'] = train_list
    data['test'] = test_list
    data['valid'] = valid_list
    data['vocab'] = vocab
    data['catgy'] = catgy
    data['max_sen_len'] = max_sen_len
    data['class_dim'] = class_dim
    return data

# read word embedding
# =========================================================
def embedding_weights_load(words_map,embedding_weights_path):
    pre_trained_embedding = None
    try:
        model = FastText.load_fasttext_format(embedding_weights_path)
        pre_trained_embedding = "bin"
    except:
        print ("fastText binary file (.bin) is not found!")
        if os.path.exists("./Word_embedding/wiki.en.vec"):
            print ("Using wikipedia(en) pre-trained word vectors.")
        else:
            print ("Downloading wikipedia(en) pre-trained word vectors.")
            chakin.download(number=2, save_dir="./Word_embedding")
        print ("Loading vectors...")
        if os.path.exists("./Word_embedding_model.pkl"):
            with open("./Word_embedding_model.pkl", mode="rb") as f:
                model = pickle.load(f)
        else:
            model =  KeyedVectors.load_word2vec_format('./Word_embedding/wiki.en.vec')
            with open("Word_embedding_model.pkl", mode="wb") as f:
                pickle.dump(model, f)
        pre_trained_embedding = "txt"

    vocab_size = len(words_map)
    word_dimension = model['a'].shape[0]
    w = np.zeros((vocab_size,word_dimension),dtype=np.float32)

    for k,v in words_map.items():
        word = k
        word_number = v

        try:
                w[word_number][:] = model[word]
        except KeyError as e:
                if pre_trained_embedding == "bin":
                    w[word_number][:] = model.seeded_vector(word)
                else:
                    np.random.seed(word_number)
                    w[word_number][:] = np.random.uniform(-0.25, 0.25, word_dimension)
    return w

# Conversion from network output to label
# =========================================================
def get_catgy_mapping(network_output_order_list, test_labels, prediction,current_depth):

    predict_result = []
    grand_labels = []

    for i in range(len(test_labels)):
        predict_result.append([])
        grand_labels.append([])

    class_dim = prediction.shape[1]

    row_idx, col_idx, val_idx = [], [], []
    for i in range(len(test_labels)):
        l_list = list(set(test_labels[i]))
        for y in l_list:
            row_idx.append(i)
            col_idx.append(y)
            val_idx.append(1)
    m = max(row_idx) + 1
    n = max(col_idx) + 1
    n = max(class_dim, n)
    test_labels = sp.csr_matrix((val_idx, (row_idx, col_idx)), shape=(m, n), dtype=np.int8).todense()

    np_orderList = np.array(network_output_order_list)

    for i,j in tqdm(enumerate(prediction), desc="Generating predict labels..."):
        one_hots = np.where(j == 1)[0]
        if len(one_hots) >= 1:
            predict_result[i] = np_orderList[one_hots].tolist()

    output_grand_truth_file_name = "CNN/RESULT/grand_truth_" + current_depth + ".csv"
    with open(output_grand_truth_file_name, 'w') as f:
        f.write(','.join(network_output_order_list)+"\n")

    with open(output_grand_truth_file_name, 'a') as f:
        for i,j in tqdm(enumerate(test_labels), desc="Generating grand truth labels..."):
            one_hots = np.where(j == 1)[1]
            if len(one_hots) >= 1:
                grand_labels[i] = np_orderList[one_hots].tolist()
                f.write(",".join(grand_labels[i])+"\n")
            else:
                f.write("\n")

    return grand_labels,predict_result

# Write results to a file
# =========================================================
def write_out_prediction(GrandLabels, PredResult, input_data_dic):

    # Writing out prediction
    # ===================================================
    print ("-"*50)
    print ("Writing out prediction...")
    test_data = input_data_dic['test']
    result_file = open("./CNN/RESULT/Prediction.txt", mode="w")
    result_file.write("Grand-truth-label\tPrediction-labels\tInput-text\n")
    for g,p,t in zip(GrandLabels, PredResult, test_data):
        result_file.write("{}\t{}\t{}\n".format(','.join(sorted(g)), ','.join(sorted(p)), t['text']))
    result_file.close()

# conversion of data
#========================================================

# conversion from text data to ndarray
# =========================================================
def build_input_sentence_data(sentences):
    x = np.array(sentences)
    return x

# conversion from sequence label to the number
# =========================================================
def build_input_label_data(labels, class_order):
    from sklearn.preprocessing import MultiLabelBinarizer
    from itertools import chain

    bml = MultiLabelBinarizer(classes=class_order, sparse_output=True)
    indexes = sp.find(bml.fit_transform(labels))
    y = []

    for i in range(len(labels)):
        y.append([])
    for i,j in zip(indexes[0], indexes[1]):
        y[i].append(j)
    return y

# padding operation
# =========================================================
def pad_sentences(sentences, padding_word=-1, max_length=50):
    sequence_length = max(max(len(x) for x in sentences), max_length)
    padded_sentences = []
    for i in range(len(sentences)):
        sentence = sentences[i]
        if len(sentence) < max_length:
            num_padding = sequence_length - len(sentence)
            new_sentence = sentence + [padding_word] * num_padding
        else:
            new_sentence = sentence[:max_length]
        padded_sentences.append(new_sentence)
    return padded_sentences

# conversion from documents and labels to the numbers
# =========================================================
def build_problem(learning_categories, depth, input_data_dic):

    train_data = input_data_dic['train']
    validation_data = input_data_dic['valid']
    test_data = input_data_dic['test']
    vocab = input_data_dic['vocab']
    max_sen_len = input_data_dic['max_sen_len']

    if depth == "flat":
        trn_text = [[vocab[word] for word in doc['text'].split()] for doc in train_data]
        trn_labels = [doc['catgy'] for doc in train_data]
        val_text = [[vocab[word] for word in doc['text'].split()] for doc in validation_data]
        val_labels = [doc['catgy'] for doc in validation_data]
        tst_text = [[vocab[word] for word in doc['text'].split()] for doc in test_data]
        tst_labels = [doc['catgy'] for doc in test_data]

    else:
        layer = int(depth[:-2])
        trn_text = [[vocab[word] for word in doc['text'].split()] for doc in train_data if (layer in doc['hie_info']) or ((layer-1) in doc['hie_info'])]
        trn_labels = [list( set(doc['catgy']) & set(learning_categories)) for doc in train_data if (layer in doc['hie_info']) or ((layer-1) in doc['hie_info'])]
        val_text = [[vocab[word] for word in doc['text'].split()] for doc in validation_data if (layer in doc['hie_info']) or ((layer-1) in doc['hie_info'])]
        val_labels = [list( set(doc['catgy']) & set(learning_categories)) for doc in validation_data if (layer in doc['hie_info']) or ((layer-1) in doc['hie_info'])]
        tst_text = [[vocab[word] for word in doc['text'].split()] for doc in test_data]
        tst_labels = [list( set(doc['catgy']) & set(learning_categories)) if layer in doc['hie_info'] else [] for doc in test_data]

    trn_padded = pad_sentences(trn_text, max_length=max_sen_len)
    val_padded = pad_sentences(val_text, max_length=max_sen_len)
    tst_padded = pad_sentences(tst_text, max_length=max_sen_len)
    x_trn = build_input_sentence_data(trn_padded)
    x_val = build_input_sentence_data(val_padded)
    x_tst = build_input_sentence_data(tst_padded)
    y_trn = build_input_label_data(trn_labels,learning_categories)
    y_val = build_input_label_data(val_labels, learning_categories)
    y_tst = build_input_label_data(tst_labels, learning_categories)

    return x_trn, y_trn, x_val, y_val, x_tst, y_tst

# conversion from the number to an ordinal number
# =========================================================
def order_n(i): return {1:"1st", 2:"2nd", 3:"3rd"}.get(i) or "%dth"%i

max_issues_repo_name: CharlotteSean/HFT-CNN

max_forks_repo_name: CharlotteSean/HFT-CNN

reason: import numpy,import scipy