text:
"""Precision for ranking."""
import numpy as np

from matchzoo.engine.base_metric import (
    BaseMetric, sort_and_couple, RankingMetric
)

class Precision(RankingMetric):
    """Precision metric."""

    ALIAS = 'precision'

    def __init__(self, k: int = 1, threshold: float = 0.):
        """
        :class:`PrecisionMetric` constructor.

        :param k: Number of results to consider.
        :param threshold: the label threshold of relevance degree.
        """
        self._k = k
        self._threshold = threshold

    def __repr__(self) -> str:
        """:return: Formated string representation of the metric."""
        return f"{self.ALIAS}@{self._k}({self._threshold})"

    def __call__(self, y_true: np.array, y_pred: np.array) -> float:
        """
        Calculate precision@k.

        Example:
            >>> y_true = [0, 0, 0, 1]
            >>> y_pred = [0.2, 0.4, 0.3, 0.1]
            >>> Precision(k=1)(y_true, y_pred)
            0.0
            >>> Precision(k=2)(y_true, y_pred)
            0.0
            >>> Precision(k=4)(y_true, y_pred)
            0.25
            >>> Precision(k=5)(y_true, y_pred)
            0.2

        :param y_true: The ground true label of each document.
        :param y_pred: The predicted scores of each document.
        :return: Precision @ k
        :raises: ValueError: len(r) must be >= k.
        """
        if self._k <= 0:
            raise ValueError(f"k must be greater than 0."
                             f"{self._k} received.")
        coupled_pair = sort_and_couple(y_true, y_pred)
        precision = 0.0
        for idx, (label, score) in enumerate(coupled_pair):
            if idx >= self._k:
                break
            if label > self._threshold:
                precision += 1.
        return precision / self._k

max_stars_repo_name: ChrisRBXiong/MatchZoo-py

max_issues_repo_name: ChrisRBXiong/MatchZoo-py

max_forks_repo_name: ChrisRBXiong/MatchZoo-py


text:
r"""
The potential fields of a homogeneous prolate ellipsoid.
"""
from __future__ import division, absolute_import

import numpy as np

from fatiando.constants import CM, T2NT
from fatiando import utils

def tf(xp, yp, zp, ellipsoids, F, inc, dec, demag=True, pmag=None):
    r"""
    The total-field anomaly produced by prolate ellipsoids.

    .. math::

        \Delta T = |\mathbf{T}| - |\mathbf{F}|,

    where :math:`\mathbf{T}` is the measured field and :math:`\mathbf{F}` is
    the local-geomagnetic field.

    The anomaly of a homogeneous ellipsoid can be calculated as:

    .. math::

        \Delta T \approx \hat{\mathbf{F}}\cdot\mathbf{B}.

    where :math:`\mathbf{B}` is the magnetic induction produced by the
    ellipsoid.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> North, y -> East and
    z -> Down.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated.
    * ellipsoids : list of :class:`mesher.ProlateEllipsoid`
        The ellipsoids. Ellipsoids must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        Ellipsoids that are ``None`` will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, will include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of calculating the magnetization of the ellipsoid. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * tf : array
        The total-field anomaly

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """
    fx, fy, fz = utils.dircos(inc, dec)
    Bx = bx(xp, yp, zp, ellipsoids, F, inc, dec, demag, pmag)
    By = by(xp, yp, zp, ellipsoids, F, inc, dec, demag, pmag)
    Bz = bz(xp, yp, zp, ellipsoids, F, inc, dec, demag, pmag)

    return fx*Bx + fy*By + fz*Bz

def bx(xp, yp, zp, ellipsoids, F, inc, dec, demag=True, pmag=None):
    r"""
    The x component of the magnetic induction produced by prolate
    ellipsoids.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> North, y -> East and
    z -> Down.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * ellipsoids : list of :class:`mesher.ProlateEllipsoid`
        The ellipsoids. Ellipsoids must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        Ellipsoids that are ``None`` will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, will include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of the resultant magnetization of the ellipsoids. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * bx: array
        The x component of the magnetic induction

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """

    res = 0
    for ellipsoid in ellipsoids:
        if ellipsoid is None:
            continue
        b1 = _bx(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        b2 = _by(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        b3 = _bz(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        res += ellipsoid.transf_matrix[0, 0]*b1 \
            + ellipsoid.transf_matrix[0, 1]*b2 \
            + ellipsoid.transf_matrix[0, 2]*b3

    return res

def by(xp, yp, zp, ellipsoids, F, inc, dec, demag=True, pmag=None):
    r"""
    The y component of the magnetic induction produced by prolate
    ellipsoids.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> North, y -> East and
    z -> Down.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * ellipsoids : list of :class:`mesher.ProlateEllipsoid`
        The ellipsoids. Ellipsoids must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        Ellipsoids that are ``None`` will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, will include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of the resultant magnetization of the ellipsoids. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * by: array
        The y component of the magnetic induction

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """

    res = 0
    for ellipsoid in ellipsoids:
        if ellipsoid is None:
            continue
        b1 = _bx(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        b2 = _by(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        b3 = _bz(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        res += ellipsoid.transf_matrix[1, 0]*b1 \
            + ellipsoid.transf_matrix[1, 1]*b2 \
            + ellipsoid.transf_matrix[1, 2]*b3

    return res

def bz(xp, yp, zp, ellipsoids, F, inc, dec, demag=True, pmag=None):
    r"""
    The z component of the magnetic induction produced by prolate
    ellipsoids.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> North, y -> East and
    z -> Down.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * ellipsoids : list of :class:`mesher.ProlateEllipsoid`
        The ellipsoids. Ellipsoids must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        Ellipsoids that are ``None`` will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, will include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of the resultant magnetization of the ellipsoids. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * bz: array
        The z component of the magnetic induction

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """

    res = 0
    for ellipsoid in ellipsoids:
        if ellipsoid is None:
            continue
        b1 = _bx(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        b2 = _by(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        b3 = _bz(xp, yp, zp, ellipsoid, F, inc, dec, demag, pmag)
        res += ellipsoid.transf_matrix[2, 0]*b1 \
            + ellipsoid.transf_matrix[2, 1]*b2 \
            + ellipsoid.transf_matrix[2, 2]*b3

    return res

def _bx(xp, yp, zp, ellipsoid, F, inc, dec, demag=True, pmag=None):
    r"""
    The x component of the magnetic induction produced by prolate
    ellipsoids in the ellipsoid system.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> semi-axis a,
    y -> semi-axis b and z -> semi-axis c.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.
        The ellipsoid. It must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        If the ellipsoid is ``None``, it will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of the resultant magnetization of the ellipsoid. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * bx : array
        The x component of the magnetic induction in the ellipsoid system.

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """

    if pmag is None:
        mx, my, mz = magnetization(ellipsoid, F, inc, dec, demag)
    else:
        assert demag is not True, 'the use of a forced magnetization \
impedes the computation of self-demagnetization'
        mx, my, mz = pmag

    # Transform the magnetization to the local coordinate system
    V = ellipsoid.transf_matrix
    mx_local = V[0, 0]*mx + V[1, 0]*my + V[2, 0]*mz
    my_local = V[0, 1]*mx + V[1, 1]*my + V[2, 1]*mz
    mz_local = V[0, 2]*mx + V[1, 2]*my + V[2, 2]*mz

    x1, x2, x3 = x1x2x3(xp, yp, zp, ellipsoid)
    lamb = _lamb(x1, x2, x3, ellipsoid)
    denominator = _dlamb_aux(x1, x2, x3, ellipsoid, lamb)
    dlamb = _dlamb(x1, x2, x3, ellipsoid, lamb, denominator, deriv='x')
    h1 = _hv(ellipsoid, lamb, v='x')
    h2 = _hv(ellipsoid, lamb, v='y')
    h3 = _hv(ellipsoid, lamb, v='z')
    g = _gv(ellipsoid, lamb, v='x')

    res = dlamb*(h1*x1*mx_local + h2*x2*my_local + h3*x3*mz_local)
    res += g*mx_local

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis
    volume = 4*np.pi*a*b*b/3

    res *= -1.5*volume*CM*T2NT

    return res

def _by(xp, yp, zp, ellipsoid, F, inc, dec, demag=True, pmag=None):
    r"""
    The y component of the magnetic induction produced by prolate
    ellipsoids in the ellipsoid system.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> semi-axis a,
    y -> semi-axis b and z -> semi-axis c.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`
        The ellipsoid. It must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        If the ellipsoid is ``None``, it will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of the resultant magnetization of the ellipsoid. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * by : array
        The y component of the magnetic induction in the ellipsoid system.

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """

    if pmag is None:
        mx, my, mz = magnetization(ellipsoid, F, inc, dec, demag)
    else:
        assert demag is not True, 'the use of a forced magnetization \
impedes the computation of self-demagnetization'
        mx, my, mz = pmag

    # Transform the magnetization to the local coordinate system
    V = ellipsoid.transf_matrix
    mx_local = V[0, 0]*mx + V[1, 0]*my + V[2, 0]*mz
    my_local = V[0, 1]*mx + V[1, 1]*my + V[2, 1]*mz
    mz_local = V[0, 2]*mx + V[1, 2]*my + V[2, 2]*mz

    x1, x2, x3 = x1x2x3(xp, yp, zp, ellipsoid)
    lamb = _lamb(x1, x2, x3, ellipsoid)
    denominator = _dlamb_aux(x1, x2, x3, ellipsoid, lamb)
    dlamb = _dlamb(x1, x2, x3, ellipsoid, lamb, denominator, deriv='y')
    h1 = _hv(ellipsoid, lamb, v='x')
    h2 = _hv(ellipsoid, lamb, v='y')
    h3 = _hv(ellipsoid, lamb, v='z')
    g = _gv(ellipsoid, lamb, v='y')

    res = dlamb*(h1*x1*mx_local + h2*x2*my_local + h3*x3*mz_local)
    res += g*my_local

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis
    volume = 4*np.pi*a*b*b/3

    res *= -1.5*volume*CM*T2NT

    return res

def _bz(xp, yp, zp, ellipsoid, F, inc, dec, demag=True, pmag=None):
    r"""
    The z component of the magnetic induction produced by prolate
    ellipsoids in the ellipsoid system.

    This code follows the approach presented by Emerson et al. (1985).

    The coordinate system of the input parameters is x -> semi-axis a,
    y -> semi-axis b and z -> semi-axis c.

    Input units should be SI. Output is in nT.

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`
        The ellipsoid. It must have the physical properties
        ``'principal susceptibilities'`` and ``'susceptibility angles'``
        as prerequisite to calculate the self-demagnetization.
        If the ellipsoid is ``None``, it will be ignored.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, will include the self-demagnetization.
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead
        of the resultant magnetization of the ellipsoid. Use this, e.g.,
        for sensitivity matrix building.

    Returns:

    * bz : array
        The z component of the magnetic induction in the ellipsoid system.

    References:

    Emerson, D. W., Clark, D., and Saul, S.: Magnetic exploration models
    incorporating remanence, demagnetization and anisotropy: HP 41C handheld
    computer algorithms, Exploration Geophysics, 16, 1-122, 1985.

    """

    if pmag is None:
        mx, my, mz = magnetization(ellipsoid, F, inc, dec, demag)
    else:
        assert demag is not True, 'the use of a forced magnetization \
impedes the computation of self-demagnetization'
        mx, my, mz = pmag

    # Transform the magnetization to the local coordinate system
    V = ellipsoid.transf_matrix
    mx_local = V[0, 0]*mx + V[1, 0]*my + V[2, 0]*mz
    my_local = V[0, 1]*mx + V[1, 1]*my + V[2, 1]*mz
    mz_local = V[0, 2]*mx + V[1, 2]*my + V[2, 2]*mz

    x1, x2, x3 = x1x2x3(xp, yp, zp, ellipsoid)
    lamb = _lamb(x1, x2, x3, ellipsoid)
    denominator = _dlamb_aux(x1, x2, x3, ellipsoid, lamb)
    dlamb = _dlamb(x1, x2, x3, ellipsoid, lamb, denominator, deriv='z')
    h1 = _hv(ellipsoid, lamb, v='x')
    h2 = _hv(ellipsoid, lamb, v='y')
    h3 = _hv(ellipsoid, lamb, v='z')
    g = _gv(ellipsoid, lamb, v='z')

    res = dlamb*(h1*x1*mx_local + h2*x2*my_local + h3*x3*mz_local)
    res += g*mz_local

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis
    volume = 4*np.pi*a*b*b/3

    res *= -1.5*volume*CM*T2NT

    return res

def x1x2x3(xp, yp, zp, ellipsoid):
    '''
    Calculates the x, y and z coordinates referred to the
    ellipsoid coordinate system.

    Parameters:

    * xp, yp, zp: numpy arrays 1D
        x, y and z coordinates of points referred to the main
        system (in meters).
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.

    Returns:

    * x1, x2, x3: numpy arrays 1D
        x, y and z coordinates of points referred to the ellipsoid
        system (in meters).
    '''

    assert xp.size == yp.size == zp.size, \
        'xp, yp and zp must have the same size'

    dx = xp - ellipsoid.x
    dy = yp - ellipsoid.y
    dz = zp - ellipsoid.z

    x1 = ellipsoid.transf_matrix[0, 0]*dx + \
        ellipsoid.transf_matrix[1, 0]*dy + \
        ellipsoid.transf_matrix[2, 0]*dz
    x2 = ellipsoid.transf_matrix[0, 1]*dx + \
        ellipsoid.transf_matrix[1, 1]*dy + \
        ellipsoid.transf_matrix[2, 1]*dz
    x3 = ellipsoid.transf_matrix[0, 2]*dx + \
        ellipsoid.transf_matrix[1, 2]*dy + \
        ellipsoid.transf_matrix[2, 2]*dz

    return x1, x2, x3

def _lamb(x1, x2, x3, ellipsoid):
    '''
    Calculates the parameter lambda for a prolate ellipsoid.

    The parameter lambda is defined as the largest root of
    the quadratic equation defining the surface of the prolate
    ellipsoid.

    Parameters:

    * x1, x2, x3: numpy arrays 1D
        x, y and z coordinates of points referred to the ellipsoid
        system (in meters).
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.

    Returns:

    * lamb: numpy array 1D
        Parameter lambda for each point in the ellipsoid system.
    '''

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis

    # auxiliary variables (http://mathworld.wolfram.com/QuadraticFormula.html)
    p1 = a*a + b*b - x1*x1 - x2*x2 - x3*x3
    p0 = a*a*b*b - b*b*x1*x1 - a*a*(x2*x2 + x3*x3)

    delta = np.sqrt(p1*p1 - 4*p0)

    lamb = (-p1 + delta)/2.

    return lamb

def _dlamb(x1, x2, x3, ellipsoid, lamb, denominator, deriv='x'):
    '''
    Calculates the spatial derivative of the parameter lambda
    with respect to the coordinates x, y or z in the ellipsoid system.

    Parameters:

    * x1, x2, x3: numpy arrays 1D
        x, y and z coordinates of points referred to the ellipsoid
        system (in meters).
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.
    * lamb: numpy array 1D
        Parameter lambda for each point in the ellipsoid system.
    * denominator: numpy array 1D
        Denominator of the expression used to calculate the spatial
        derivative of the parameter lambda.
    * deriv: string
        Defines the coordinate with respect to which the
        derivative will be calculated. It must be 'x', 'y' or 'z'.

    Returns:

    * dlamb_dv: numpy array 1D
        Derivative of lambda with respect to the coordinate
        v = x, y, z in the ellipsoid system.
    '''

    assert deriv in ['x', 'y', 'z'], 'deriv must represent a coordinate \
        x, y or z'

    assert denominator.size == lamb.size == x1.size == x2.size == x3.size, \
        'x1, x2, x3, lamb and denominator must have the same size'

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis

    if deriv is 'x':
        dlamb_dv = (2*x1/(a*a + lamb))/denominator

    if deriv is 'y':
        dlamb_dv = (2*x2/(b*b + lamb))/denominator

    if deriv is 'z':
        dlamb_dv = (2*x3/(b*b + lamb))/denominator

    return dlamb_dv

def _dlamb_aux(x1, x2, x3, ellipsoid, lamb):
    '''
    Calculates an auxiliary variable used to calculate the spatial
    derivatives of the parameter lambda with respect to the
    coordinates x, y and z in the ellipsoid system.

    Parameters:

    * x1, x2, x3: numpy arrays 1D
        x, y and z coordinates of points referred to the ellipsoid
        system (in meters).
    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.
    * lamb: numpy array 1D
        Parameter lambda for each point in the ellipsoid system.

    Returns:

    * aux: numpy array 1D
        Denominator of the expression used to calculate the spatial
        derivative of the parameter lambda.
    '''

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis

    aux1 = x1/(a*a + lamb)
    aux2 = x2/(b*b + lamb)
    aux3 = x3/(b*b + lamb)
    aux = aux1*aux1 + aux2*aux2 + aux3*aux3

    return aux

def demag_factors(ellipsoid):
    '''
    Calculates the demagnetizing factors n11 and n22.

    Parameters:

    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.

    Returns:

    * n11, n22: floats
        Demagnetizing factors (in SI) along, respectively, the
        large and small axes of the prolate ellipsoid.
    '''

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis

    m = a/b

    aux1 = m*m - 1
    aux2 = np.sqrt(aux1)

    n11 = (1/aux1)*((m/aux2)*np.log(m + aux2) - 1)
    n22 = 0.5*(1 - n11)

    return n11, n22

def magnetization(ellipsoid, F, inc, dec, demag):
    '''
    Calculates the resultant magnetization corrected from
    demagnetizing in the main system.

    Parameters:

    * ellipsoid: element of :class:`mesher.ProlateEllipsoid`.
    * F, inc, dec : floats
       The intensity (in nT), inclination and declination (in degrees) of
       the local-geomagnetic field.
    * demag : boolean
        If True, will include the self-demagnetization.

    Returns:

    * resultant_mag: numpy array 1D
        Resultant magnetization (in A/m) in the main system.
    '''

    # Remanent magnetization
    if 'remanent magnetization' in ellipsoid.props:
        intensity = ellipsoid.props['remanent magnetization'][0]
        inclination = ellipsoid.props['remanent magnetization'][1]
        declination = ellipsoid.props['remanent magnetization'][2]
        remanent_mag = utils.ang2vec(intensity, inclination, declination)
    else:
        remanent_mag = np.zeros(3)

    suscep = ellipsoid.susceptibility_tensor

    # Induced magnetization
    if suscep is not None:
        geomag_field = utils.ang2vec(F/(4*np.pi*100), inc, dec)
        induced_mag = np.dot(suscep, geomag_field)
    else:
        induced_mag = np.zeros(3)

    # Self-demagnetization
    if demag is True:

        assert suscep is not None, 'self-demagnetization requires a \
susceptibility tensor'

        n11, n22 = demag_factors(ellipsoid)
        coord_transf_matrix = ellipsoid.transf_matrix
        suscep_tilde = np.dot(np.dot(coord_transf_matrix.T, suscep),
                              coord_transf_matrix)
        aux = np.linalg.inv(np.identity(3) + np.dot(suscep_tilde,
                                                    np.diag([n11, n22, n22])))
        Lambda = np.dot(np.dot(coord_transf_matrix, aux),
                        coord_transf_matrix.T)

        # resultant magnetization in the main system
        resultant_mag = np.dot(Lambda, induced_mag + remanent_mag)

    else:

        assert (suscep is not None) or ('remanent magnetization'
                                        in ellipsoid.props), 'neglecting \
            self-demagnetization requires a susceptibility tensor or a rem\
            anent magnetization'

        # resultant magnetization in the main system
        resultant_mag = induced_mag + remanent_mag

    return resultant_mag

def _hv(ellipsoid, lamb, v='x'):
    '''
    Calculates an auxiliary variable used to calculate the
    depolarization tensor outside the ellipsoidal body.

    Parameters:

    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.
    * lamb: numpy array 1D
        Parameter lambda for each point in the ellipsoid system.
    * v: string
        Defines the coordinate with respect to which the
        variable hv will be calculated. It must be 'x', 'y' or 'z'.

    Returns:

    * hv: numpy array 1D
        Auxiliary variable.
    '''

    assert v in ['x', 'y', 'z'], "v must be 'x', 'y' or 'z'"

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis

    aux1 = a*a + lamb
    aux2 = b*b + lamb

    if v is 'x':
        hv = -1./(np.sqrt(aux1*aux1*aux1)*aux2)

    if v is 'y' or 'z':
        hv = -1./(np.sqrt(aux1)*aux2*aux2)

    return hv

def _gv(ellipsoid, lamb, v='x'):
    '''
    Diagonal term of the depolarization tensor defined outside the
    ellipsoidal body.

    Parameters:

    * ellipsoid : element of :class:`mesher.ProlateEllipsoid`.
    * lamb: numpy array 1D
        Parameter lambda for each point in the ellipsoid system.
    * v: string
        Defines the coordinate with respect to which the
        variable gv will be calculated. It must be 'x', 'y' or 'z'.

    Returns:

    * gv: numpy array 1D
        Diagonal term of the depolarization tensor calculated for
        each lambda.
    '''

    assert v in ['x', 'y', 'z'], "v must be 'x', 'y' or 'z'"

    a = ellipsoid.large_axis
    b = ellipsoid.small_axis

    log = np.log((np.sqrt(a*a-b*b)+np.sqrt(a*a+lamb))/np.sqrt(b*b+lamb))
    aux1 = 1./np.sqrt((a*a - b*b)*(a*a - b*b)*(a*a - b*b))

    if v is 'x':
        aux2 = np.sqrt((a*a - b*b)/(a*a + lamb))
        gv = 2*aux1*(log - aux2)

    if v is 'y' or 'z':
        aux2 = np.sqrt((a*a - b*b)*(a*a + lamb))/(b*b + lamb)
        gv = aux1*(aux2 - log)

    return gv

ext: py

max_stars_repo_name: pinga-lab/magnetic-ellipsoid

max_issues_repo_name: pinga-lab/magnetic-ellipsoid

max_forks_repo_name: pinga-lab/magnetic-ellipsoid

reason: import numpy


text:
import math
import numpy as np
import tensorflow as tf
import scipy
import matplotlib.pyplot as plt
from tensorflow.python.framework import ops

np.random.seed(1)
tf.compat.v1.disable_eager_execution()

def random_mini_batches(X, Y, mini_batch_size = 64, seed = 0):
    """
    Creates a list of random minibatches from (X, Y)

    Arguments:
    X -- input data, of shape (input size, number of examples) (m, Hi, Wi, Ci)
    Y -- true "label" vector (containing 0 if cat, 1 if non-cat), of shape (1, number of examples) (m, n_y)
    mini_batch_size - size of the mini-batches, integer
    seed -- this is only for the purpose of grading, so that you're "random minibatches are the same as ours.

    Returns:
    mini_batches -- list of synchronous (mini_batch_X, mini_batch_Y)
    """

    m = X.shape[0]                  # number of training examples
    mini_batches = []
    np.random.seed(seed)

    # Step 1: Shuffle (X, Y)
    permutation = list(np.random.permutation(m))
    shuffled_X = X[permutation,:,:,:]
    shuffled_Y = Y[permutation,:]

    # Step 2: Partition (shuffled_X, shuffled_Y). Minus the end case.
    num_complete_minibatches = math.floor(m/mini_batch_size) # number of mini batches of size mini_batch_size in your partitioning
    for k in range(0, num_complete_minibatches):
        mini_batch_X = shuffled_X[k * mini_batch_size : k * mini_batch_size + mini_batch_size,:,:,:]
        mini_batch_Y = shuffled_Y[k * mini_batch_size : k * mini_batch_size + mini_batch_size,:]
        mini_batch = (mini_batch_X, mini_batch_Y)
        mini_batches.append(mini_batch)

    # end case (last mini-batch < mini_batch_size)
    if m % mini_batch_size != 0:
        mini_batch_X = shuffled_X[num_complete_minibatches * mini_batch_size : m,:,:,:]
        mini_batch_Y = shuffled_Y[num_complete_minibatches * mini_batch_size : m,:]
        mini_batch = (mini_batch_X, mini_batch_Y)
        mini_batches.append(mini_batch)

    return mini_batches

def create_placeholders(n_H0, n_W0, n_C0, n_y):
    X = tf.compat.v1.placeholder(tf.float32, shape = (None, n_H0, n_W0, n_C0), name='X')
    Y = tf.compat.v1.placeholder(tf.float32, shape = (None, n_y), name = 'Y')
    keep_prob = tf.compat.v1.placeholder(tf.float32, name="keep_prob")

    return X,Y,keep_prob

def initialize_parameters():

    # adjust this parameter
    W1 = tf.compat.v1.get_variable("W1", [4,4,1,8], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed = 0))
    W2 = tf.compat.v1.get_variable("W2", [2,2,8,16], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed = 0))

    parameters = {"W1": W1,
                "W2": W2}

    return parameters

def forward_propagation(X, keep_prob):
    tf.random.set_seed(1)

    # CONV 1
    W1 = tf.compat.v1.get_variable("W1", [4,4,1,64], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed = 0))
    b1 = tf.compat.v1.get_variable("b1", [64,1], initializer = tf.compat.v1.zeros_initializer())

    Z1 = tf.nn.conv2d(input = X, filters = W1, strides = [1,1,1,1], padding = "SAME")
    A1 = tf.nn.relu(Z1)
    print(A1)
    # MAX POOL 1
    P1 = tf.nn.max_pool2d(input=A1, ksize = [1,4,4,1], strides = [1,4,4,1], padding = "SAME")
    print(P1)

    # CONV 2
    W2 = tf.compat.v1.get_variable("W2", [2,2,64,128], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed = 0))
    b2 = tf.compat.v1.get_variable("b2", [128,1], initializer = tf.compat.v1.zeros_initializer())

    Z2 = tf.nn.conv2d(input = P1, filters = W2, strides = [1,1,1,1], padding = "SAME")
    A2 = tf.nn.relu(Z2)
    print(A2)
    # MAX POOL 2
    P2 = tf.nn.max_pool2d(input=A2, ksize = [1,4,4,1], strides = [1,4,4,1], padding = "SAME")
    print(P2)
    # FLATTEN
    P2 = tf.compat.v1.layers.flatten(P2)
    print(P2)
    # FULL CONNECT 1
    W3 = tf.compat.v1.get_variable('W3', [512, P2.shape[1:2].num_elements()], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed=0))
    b3 = tf.compat.v1.get_variable('b3', [512,1], initializer = tf.compat.v1.zeros_initializer())
    Z3 = tf.add(tf.matmul(W3,tf.linalg.matrix_transpose(P2)), b3)
    A3 = tf.nn.relu(Z3)
    print(A3)
    # FULL CONNECT 2
    W4 = tf.compat.v1.get_variable('W4', [256, 512], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed=0))
    b4 = tf.compat.v1.get_variable('b4', [256,1], initializer = tf.compat.v1.zeros_initializer())
    A4_drop = tf.nn.dropout(A3, 1 - (keep_prob))
    Z4 = tf.add(tf.matmul(W4,A4_drop), b4)
    A4 = tf.nn.relu(Z4)
    print(A4)
    # FULL CONNECT 3
    W5 = tf.compat.v1.get_variable('W5', [63,256], initializer = tf.compat.v1.keras.initializers.VarianceScaling(scale=1.0, mode="fan_avg", distribution="uniform", seed=0))
    b5 = tf.compat.v1.get_variable('b5', [63,1], initializer = tf.compat.v1.zeros_initializer())
    A5_drop = tf.nn.dropout(A4, 1 - (keep_prob))
    Z5 = tf.add(tf.matmul(W5,A5_drop), b5)
    print(Z5)

    Z5 = tf.linalg.matrix_transpose(Z5)

    return Z5

def compute_cost(Z3, Y):
    cost = tf.reduce_mean(input_tensor=tf.nn.softmax_cross_entropy_with_logits(logits = Z3, labels = Y))

    return cost

def model(X_train, Y_train, X_test, Y_test, learning_rate = 0.001, keep_prob_prob = 1,
        num_epochs = 20, minibatch_size = 64, print_cost = True):

	ops.reset_default_graph()
	tf.random.set_seed(1)
	seed = 3
	m, n_H0, n_W0, n_C0 = X_train.shape
	n_y = Y_train.shape[1]
	costs = []

	X, Y, keep_prob = create_placeholders(n_H0, n_W0, n_C0, n_y)

	Z3 = forward_propagation(X, keep_prob)

	cost = compute_cost(Z3, Y)

	optimizer = tf.compat.v1.train.AdamOptimizer(learning_rate = learning_rate).minimize(cost)

	init = tf.compat.v1.global_variables_initializer()

	# This is for accuracy testing
	y_pred = tf.nn.softmax(Z3)
	y_pred_class = tf.argmax(input=y_pred, axis = 1)
	y_true_class = tf.argmax(input=Y, axis = 1)

	correct_prediction = tf.equal(y_pred_class, y_true_class)
	accuracy = tf.reduce_mean(input_tensor=tf.cast(correct_prediction, tf.float32))

	with tf.compat.v1.Session() as sess:

		sess.run(init)

		for epoch in range(num_epochs):

			minibatch_cost = 0
			num_minibatches = int(m / minibatch_size)
			seed = seed + 1
			minibatches = random_mini_batches(X_train, Y_train, minibatch_size, seed)

			for minibatch in minibatches:

				minibatch_X, minibatch_Y = minibatch

				_, temp_cost = sess.run([optimizer, cost], {X:minibatch_X, Y:minibatch_Y, keep_prob :keep_prob_prob})
				minibatch_cost += temp_cost / num_minibatches

			if print_cost == True:
				print("Cost after epoch %i %f" %(epoch, minibatch_cost))
				costs.append(minibatch_cost)

			acc = sess.run(accuracy, feed_dict = {X: X_test, Y: Y_test, keep_prob : keep_prob_prob})
			print("Accuracy on test dataset after epoch %i: %f" %(epoch, acc))

		plt.plot(np.squeeze(costs))
		plt.ylabel('cost')
		plt.xlabel('iterations (per tens)')
		plt.title("Learning rate = " + str(learning_rate))
		plt.show()
		saver = tf.compat.v1.train.Saver()
		saver.save(sess, 'my-model.ckpt')

		predict_op = tf.argmax(input=Z3, axis=1)
		correct_prediction = tf.equal(predict_op, tf.argmax(input=Y,axis=1))

		accuracy = tf.reduce_mean(input_tensor=tf.cast(correct_prediction, "float"))
		print(accuracy)
		train_accuracy = accuracy.eval({X: X_train, Y: Y_train, keep_prob: keep_prob_prob})
		test_accuracy = accuracy.eval({X: X_test, Y: Y_test, keep_prob: keep_prob_prob})
		print("Train Accuracy:", train_accuracy)
		print("Test Accuracy:", test_accuracy)

		return None

if __name__ == "__main__":
	train_inputs, train_outputs , test_inputs, test_outputs = np.load('training_inputs_java.npy'), np.load('training_outputs_java.npy'), np.load('test_inputs_java.npy'), np.load('test_outputs_java.npy')

	X = np.array(train_inputs)
	Y = np.array(train_outputs)
	Y = Y.reshape(Y.shape[0],Y.shape[1])
	X_test = np.array(test_inputs)
	Y_test = np.array(test_outputs)
	Y_test = Y_test.reshape(Y_test.shape[0], Y_test.shape[1])
	print(X.shape)
	print(Y.shape)
	print(X_test.shape)
	print(Y_test.shape)

	parameters = model(X_train = X, Y_train = Y, X_test = X_test, Y_test = Y_test)
	print(parameters)

max_stars_repo_name: rvdrwn/javaOCR

max_issues_repo_name: rvdrwn/javaOCR

max_forks_repo_name: rvdrwn/javaOCR

reason: import numpy,import scipy


text:
#!/usr/bin/env python
# Copyright 2014-2021 The PySCF Developers. All Rights Reserved.
#
# 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.
#
# Author: Tianyu Zhu <zhutianyu1991@gmail.com>
#

"""
Spin-unrestricted random phase approximation (direct RPA/dRPA in chemistry)
with N^4 scaling

Method:
    Main routines are based on GW-AC method descirbed in:
    T. Zhu and G.K.-L. Chan, J. Chem. Theory. Comput. 17, 727-741 (2021)
    X. Ren et al., New J. Phys. 14, 053020 (2012)
"""

import numpy as np
from pyscf import lib
from pyscf.lib import logger
from pyscf.ao2mo import _ao2mo
from pyscf import df, scf
from pyscf.mp.ump2 import get_nocc, get_nmo, get_frozen_mask
from pyscf.gw.rpa import RPA, _get_scaled_legendre_roots

einsum = lib.einsum

# ****************************************************************************
# core routines, kernel, rpa_ecorr, rho_response
# ****************************************************************************

def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=None, verbose=logger.NOTE):
    """
    RPA correlation and total energy

    Args:
        Lpq : density fitting 3-center integral in MO basis.
        nw : number of frequency point on imaginary axis.
        vhf_df : using density fitting integral to compute HF exchange.

    Returns:
        e_tot : RPA total energy
        e_hf : EXX energy
        e_corr : RPA correlation energy
    """
    mf = rpa._scf
    # only support frozen core
    if rpa.frozen is not None:
        assert isinstance(rpa.frozen, int)
        assert (rpa.frozen < rpa.nocc[0] and rpa.frozen < rpa.nocc[1])

    if Lpq is None:
        Lpq = rpa.ao2mo(mo_coeff)

    # Grids for integration on imaginary axis
    freqs, wts = _get_scaled_legendre_roots(nw)

    # Compute HF exchange energy (EXX)
    dm = mf.make_rdm1()
    uhf = scf.UHF(rpa.mol)
    e_hf = uhf.energy_elec(dm=dm)[0]
    e_hf += mf.energy_nuc()

    # Compute RPA correlation energy
    e_corr = get_rpa_ecorr(rpa, Lpq, freqs, wts)

    # Compute totol energy
    e_tot = e_hf + e_corr

    logger.debug(rpa, '  RPA total energy = %s', e_tot)
    logger.debug(rpa, '  EXX energy = %s, RPA corr energy = %s', e_hf, e_corr)

    return e_tot, e_hf, e_corr

def get_rpa_ecorr(rpa, Lpq, freqs, wts):
    """
    Compute RPA correlation energy
    """
    mo_energy = _mo_energy_without_core(rpa, rpa._scf.mo_energy)
    nocca, noccb = rpa.nocc
    nw = len(freqs)
    naux = Lpq[0].shape[0]

    homo = max(mo_energy[0][nocca-1], mo_energy[1][noccb-1])
    lumo = min(mo_energy[0][nocca], mo_energy[1][noccb])
    if (lumo-homo) < 1e-3:
        logger.warn(rpa, 'Current RPA code not well-defined for degeneracy!')

    e_corr = 0.
    for w in range(nw):
        Pi = get_rho_response(freqs[w], mo_energy, Lpq[0,:,:nocca,nocca:], Lpq[1,:,:noccb,noccb:])
        ec_w = np.log(np.linalg.det(np.eye(naux) - Pi))
        ec_w += np.trace(Pi)
        e_corr += 1./(2.*np.pi) * ec_w * wts[w]

    return e_corr

def get_rho_response(omega, mo_energy, Lpqa, Lpqb):
    '''
    Compute density response function in auxiliary basis at freq iw
    '''
    naux, nocca, nvira = Lpqa.shape
    naux, noccb, nvirb = Lpqb.shape
    eia_a = mo_energy[0,:nocca,None] - mo_energy[0,None,nocca:]
    eia_b = mo_energy[1,:noccb,None] - mo_energy[1,None,noccb:]
    eia_a = eia_a / (omega**2 + eia_a*eia_a)
    eia_b = eia_b / (omega**2 + eia_b*eia_b)
    Pia_a = Lpqa * (eia_a * 2.0)
    Pia_b = Lpqb * (eia_b * 2.0)
    # Response from both spin-up and spin-down density
    Pi = einsum('Pia, Qia -> PQ', Pia_a, Lpqa) + einsum('Pia, Qia -> PQ', Pia_b, Lpqb)
    return Pi

def _mo_energy_without_core(rpa, mo_energy):
    moidx = get_frozen_mask(rpa)
    mo_energy = (mo_energy[0][moidx[0]], mo_energy[1][moidx[1]])
    return np.asarray(mo_energy)

def _mo_without_core(rpa, mo):
    moidx = get_frozen_mask(rpa)
    mo = (mo[0][:,moidx[0]], mo[1][:,moidx[1]])
    return np.asarray(mo)

class URPA(RPA):

    def dump_flags(self):
        log = logger.Logger(self.stdout, self.verbose)
        log.info('')
        log.info('******** %s ********', self.__class__)
        log.info('method = %s', self.__class__.__name__)
        nocca, noccb = self.nocc
        nmoa, nmob = self.nmo
        nvira = nmoa - nocca
        nvirb = nmob - noccb
        log.info('RPA (nocca, noccb) = (%d, %d), (nvira, nvirb) = (%d, %d)',
                 nocca, noccb, nvira, nvirb)
        if self.frozen is not None:
            log.info('frozen orbitals = %s', str(self.frozen))
        return self

    get_nocc = get_nocc
    get_nmo = get_nmo
    get_frozen_mask = get_frozen_mask

    def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40):
        """
        Args:
            mo_energy : 2D array (2, nmo), mean-field mo energy
            mo_coeff : 3D array (2, nmo, nmo), mean-field mo coefficient
            Lpq : 4D array (2, naux, nmo, nmo), 3-index ERI
            nw: interger, grid number

        Returns:
            self.e_tot : RPA total eenrgy
            self.e_hf : EXX energy
            self.e_corr : RPA correlation energy
        """
        if mo_coeff is None:
            mo_coeff = _mo_without_core(self, self._scf.mo_coeff)
        if mo_energy is None:
            mo_energy = _mo_energy_without_core(self, self._scf.mo_energy)

        cput0 = (logger.process_clock(), logger.perf_counter())
        self.dump_flags()
        self.e_tot, self.e_hf, self.e_corr = \
                        kernel(self, mo_energy, mo_coeff, Lpq=Lpq, nw=nw, verbose=self.verbose)

        logger.timer(self, 'RPA', *cput0)
        return self.e_corr

    def ao2mo(self, mo_coeff=None):
        nmoa, nmob = self.nmo
        nao = self.mo_coeff[0].shape[0]
        naux = self.with_df.get_naoaux()
        mem_incore = (nmoa**2*naux + nmob**2*naux + nao**2*naux) * 8/1e6
        mem_now = lib.current_memory()[0]

        moa = np.asarray(mo_coeff[0], order='F')
        mob = np.asarray(mo_coeff[1], order='F')
        ijslicea = (0, nmoa, 0, nmoa)
        ijsliceb = (0, nmob, 0, nmob)
        Lpqa = None
        Lpqb = None
        if (mem_incore + mem_now < 0.99*self.max_memory) or self.mol.incore_anyway:
            Lpqa = _ao2mo.nr_e2(self.with_df._cderi, moa, ijslicea, aosym='s2', out=Lpqa)
            Lpqb = _ao2mo.nr_e2(self.with_df._cderi, mob, ijsliceb, aosym='s2', out=Lpqb)
            return np.asarray((Lpqa.reshape(naux,nmoa,nmoa),Lpqb.reshape(naux,nmob,nmob)))
        else:
            logger.warn(self, 'Memory may not be enough!')
            raise NotImplementedError

if __name__ == '__main__':
    from pyscf import gto, dft
    mol = gto.Mole()
    mol.verbose = 4
    mol.atom = 'F 0 0 0'
    mol.basis = 'def2-svp'
    mol.spin = 1
    mol.build()

    mf = dft.UKS(mol)
    mf.xc = 'pbe0'
    mf.kernel()

    rpa = URPA(mf)
    rpa.kernel()
    print ('RPA e_tot, e_hf, e_corr = ', rpa.e_tot, rpa.e_hf, rpa.e_corr)
    assert(abs(rpa.e_corr- -0.20980646878974454) < 1e-6)
    assert(abs(rpa.e_tot- -99.49292565821425) < 1e-6)


text:
import numpy as np
from utils import load_data, zero_mean_unit_variance, vec_to_one_hot_matrix

class JuraEnv(object):
	def __init__(self):
		train_fn = 'data/jura_train_dataframe.pkl'
		test_fn = 'data/jura_val_dataframe.pkl'

		input_features = ['xloc', 'yloc']
		# input_features = ['xloc', 'yloc', 'landuse', 'rock']
		output_features = ['Cd', 'Ni', 'Zn']
		train_x, train_y, test_x, test_y = load_data(train_fn, test_fn, input_features, output_features)

		if len(input_features)==4:
			lu_train = vec_to_one_hot_matrix(train_x[:,2].astype(int)-1)
			ro_train = vec_to_one_hot_matrix(train_x[:,3].astype(int)-1)
			train_x = np.concatenate([train_x[:,:2], lu_train, ro_train], axis=1)
			lu_test = vec_to_one_hot_matrix(test_x[:,2].astype(int)-1)
			ro_test = vec_to_one_hot_matrix(test_x[:,3].astype(int)-1)
			test_x = np.concatenate([test_x[:,:2], lu_test, ro_test], axis=1)

		# arrange all data by category
		# predict Cd at test locations given Cd, Ni and Zn at training locations and Ni and Zn at test locations
		self.num_tasks = len(output_features)
		target_feature = 'Cd'
		self.target_task = output_features.index(target_feature)
		others = list(range(self.num_tasks))
		others.remove(self.target_task)

		# normalize data
		# take log of Cd and Zn measurements
		correct_skewness = ['Cd', 'Zn']
		for f in correct_skewness:
			idx = output_features.index(f)
			train_y[:, idx] = np.log(train_y[:, idx])
			test_y[:, idx] = np.log(test_y[:, idx])

		# zero mean unit variance normalization
		# in MOAL, this is done across train and test jointly
		all_y = zero_mean_unit_variance(np.concatenate([train_y, test_y], 0))
		train_y = all_y[:len(train_x)]
		test_y = all_y[len(train_x):]

		# test set consists of target_task at test locations
		# training set consists of rest all
		# train_ind = [0, 0, ..., 1, 1, ...]
		ind1 = np.arange(self.num_tasks).reshape(1,-1).repeat(len(train_x), 0).flatten()
		ind2 = np.array(others).reshape(1,-1).repeat(len(test_x), 0).flatten()
		x1 = np.repeat(train_x, self.num_tasks, 0)
		x2 = np.repeat(test_x, len(others), 0)
		y1 = train_y.flatten()
		y2 = test_y[:,others].flatten()

		# effective training set
		self.X = np.concatenate([x1, x2])
		self.Y = np.concatenate([y1, y2])
		self.ind = np.concatenate([ind1, ind2])

		# effective test set
		self.test_X = test_x
		self.test_Y = test_y[:,self.target_task]
		self.test_ind = np.full(len(test_x), self.target_task)

		# normalize X between -1 and 1
		# x = np.concatenate([self.X, self.test_X])
		# x = 2*((x-x.min(0))/x.max(0)) - 1
		# self.X = x[:len(self.Y), :]
		# self.test_X = x[len(self.Y):, :]

	@property
	def num_samples(self):
		return len(self.X)

max_stars_repo_name: sumitsk/MOGP-AL

max_issues_repo_name: sumitsk/MOGP-AL

max_forks_repo_name: sumitsk/MOGP-AL


text:
import os
import gym
from dqn_tf import DeepQNetwork, Agent
from utils import plotLearning
import numpy as np
from gym import wrappers
import matplotlib.pyplot as plt

def preprocess(observation):
    observation = observation / 255
    return np.mean(observation[30:,:], axis=2).reshape(180,160,1)

def stack_frames(stacked_frames, frame, buffer_size):
    if stacked_frames is None:
        stacked_frames = np.zeros((buffer_size, *frame.shape))
        for idx, _ in enumerate(stacked_frames):
            stacked_frames[idx,:] = frame
    else:
        stacked_frames[0:buffer_size-1,:] = stacked_frames[1:,:]
        stacked_frames[buffer_size-1, :] = frame

    stacked_frames = stacked_frames.reshape(1, *frame.shape[0:2], buffer_size)

    return stacked_frames

if __name__ == '__main__':
    #os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID"
    #os.environ["CUDA_VISIBLE_DEVICES"] = "0"
    #os.environ['TF_FORCE_GPU_ALLOW_GROWTH'] = 'true'

    env = gym.make('Breakout-v0')
    load_checkpoint = False
    agent = Agent(gamma=0.99, epsilon=1.0, alpha=0.000025, input_dims=(180,160,4),
                  n_actions=3, mem_size=25000, batch_size=64)
    if load_checkpoint:
        agent.load_models()
    filename = 'breakout-alpha0p000025-gamma0p9-only-one-fc-2.png'
    scores = []
    eps_history = []
    numGames = 50000
    stack_size = 4
    score = 0
    # uncomment the line below to record every episode.
    #env = wrappers.Monitor(env, "tmp/breakout-0",
    #                         video_callable=lambda episode_id: True, force=True)
    """
    print("Loading up the agent's memory with random gameplay")

    while agent.mem_cntr < 25000:
        done = False
        observation = env.reset()
        observation = preprocess(observation)
        stacked_frames = None
        observation = stack_frames(stacked_frames, observation, stack_size)
        while not done:
            action = np.random.choice([0, 1, 2])
            action += 1
            observation_, reward, done, info = env.step(action)
            observation_ = stack_frames(stacked_frames,
                                        preprocess(observation_), stack_size)
            action -= 1
            agent.store_transition(observation, action,
                                   reward, observation_, int(done))
            observation = observation_
    print("Done with random gameplay. Game on.")
    """
    n_steps = 0
    for i in range(numGames):
        done = False
        #if i % 100 == 0 and i > 0:
        #    x = [j+1 for j in range(i)]

        #    plotLearning(x, scores, eps_history, filename)
        observation = env.reset()
        observation = preprocess(observation)
        stacked_frames = None
        observation = stack_frames(stacked_frames, observation, stack_size)
        score = 0
        while not done:
            action = agent.choose_action(observation)
            action += 1
            observation_, reward, done, info = env.step(action)
            n_steps += 1
            observation_ = stack_frames(stacked_frames,
                                        preprocess(observation_), stack_size)
            score += reward
            action -= 1
            agent.store_transition(observation, action,
                                   reward, observation_, int(done))
            observation = observation_
            if n_steps % 4 == 0:
                agent.learn()
        if i % 12 == 0 and i > 0:
            avg_score = np.mean(scores[max(0, i-12):(i+1)])
            print('episode: ', i,'score: ', score,
                 ' average score %.3f' % avg_score,
                'epsilon %.3f' % agent.epsilon)
            agent.save_models()
        else:
            print('episode: ', i,'score: ', score)
        eps_history.append(agent.epsilon)
        scores.append(score)
    x = [i+1 for i in range(numGames)]
    plotLearning(x, scores, eps_history, filename)

max_stars_repo_name: SuperSaiyan-God/Reinforcement-Learning

max_issues_repo_name: SuperSaiyan-God/Reinforcement-Learning

max_forks_repo_name: SuperSaiyan-God/Reinforcement-Learning


text:
#!/usr/bin/env python3

"""
Evaluate the proximity of statements in pre-defined clusters.
"""

import json
import logging
import sys
import warnings
from pathlib import Path

import numpy as np
import tokenizers as tk
import torch

from theissues import clustering, training, utils

def distances_euclid(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    """
    Calculate the pair-wise distance between each vector in `x` and `y`, where
    the eucledian space is the last axis.
    """
    return np.linalg.norm(x - y, axis=-1)

def distances_cosine(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    """
    Calculate the pair-wise cosine distance between each vector in `x` and `y`,
    where the eucledian space is the last axis.
    """
    dot_product = np.sum(x * y, axis=-1)
    norm = np.linalg.norm(x, axis=-1) * np.linalg.norm(y, axis=-1)
    return 1.0 - dot_product / norm

def main(
    path_statements: Path,
    path_sequences: Path,
    path_embeddings: Path,
    path_tokenizer: Path,
    path_log: Path,
):
    logging_handlers = [logging.StreamHandler(sys.stdout)]
    if path_log is not None:
        # clear the log file
        with path_log.open("w"):
            pass
        logging_handlers.append(logging.FileHandler(path_log))
    logging.basicConfig(
        level=logging.INFO,
        format="%(message)s",
        handlers=logging_handlers,
    )

    tokenizer = tk.Tokenizer.from_file(str(path_tokenizer))
    special_tokens = utils.SpecialTokens(tokenizer)

    logging.info("Reading statements ...")
    with path_statements.open("r") as fio:
        statements = [utils.Statement(*json.loads(l)) for l in fio]

    logging.info("Reading sequences ...")
    with path_sequences.open("rb") as fio:
        sequences = np.load(fio)
    sequence_splits = training.build_sequences_splits(
        torch.from_numpy(sequences), special_tokens.bos_id
    ).numpy()

    with path_embeddings.open("rb") as fio:
        embeddings: np.ndarray = np.load(fio)

    num_statements = len(statements)
    if sequence_splits.shape[0] != num_statements:
        raise RuntimeError("sequences don't match statements")
    if embeddings.shape[0] != num_statements:
        raise RuntimeError("embeddings don't match statements")

    logging.info("Building clusters ...")
    topic_groups = clustering.build_topic_groups([s.text for s in statements])

    logging.info("Cluster sizes:")
    for topic, group in topic_groups.items():
        logging.info("%s: %d", topic, len(group))
        if not group:
            warnings.warn(f"empty topic group: {topic}")

    logging.info("Evaluating clusters ...")

    rng: np.random.Generator = np.random.default_rng()

    # NOTE: relying on dict insertion order
    pairs_intra = clustering.build_intra_topic_pairs(
        groups=topic_groups.values(),
        embeddings=embeddings,
        num_pairs=2 ** 15,
        rng=rng,
    )
    pairs_inter = clustering.build_inter_topic_pairs(
        groups=topic_groups.values(),
        embeddings=embeddings,
        num_pairs=2 ** 15,
        rng=rng,
    )

    logging.info("Euclidean distances:")

    distances_intra = distances_euclid(pairs_intra[:, 0], pairs_intra[:, 1])
    distances_inter = distances_euclid(pairs_inter[:, 0], pairs_inter[:, 1])

    intra_mu = np.mean(distances_intra)
    intra_sig = np.std(distances_intra)
    inter_mu = np.mean(distances_inter)
    inter_sig = np.std(distances_inter)
    separation = (inter_mu - intra_mu) / (inter_sig ** 2 + intra_sig ** 2) ** 0.5

    logging.info("Intra distance: %.2e ± %.2e", intra_mu, intra_sig)
    logging.info("Inter distance: %.2e ± %.2e", inter_mu, inter_sig)
    logging.info("Separation: %+.2e", separation)

    logging.info("Cosine distances:")

    distances_intra = distances_cosine(pairs_intra[:, 0], pairs_intra[:, 1])
    distances_inter = distances_cosine(pairs_inter[:, 0], pairs_inter[:, 1])

    intra_mu = np.mean(distances_intra)
    intra_sig = np.std(distances_intra)
    inter_mu = np.mean(distances_inter)
    inter_sig = np.std(distances_inter)
    separation = (inter_mu - intra_mu) / (inter_sig ** 2 + intra_sig ** 2) ** 0.5

    logging.info("Intra distance: %.2e ± %.2e", intra_mu, intra_sig)
    logging.info("Inter distance: %.2e ± %.2e", inter_mu, inter_sig)
    logging.info("Separation: %+.2e", separation)

    logging.info("Cosine neighbours:")
    # want to re-use the same examples everytime this is run
    rng: np.random.Generator = np.random.default_rng(123)

    for topic, group in topic_groups.items():
        if not group:
            continue
        ref_idx = rng.choice(group)
        ref_embedding = embeddings[ref_idx]
        distances = distances_cosine(ref_embedding[np.newaxis, :], embeddings)
        sort_idx = np.argsort(distances)
        logging.info("Topic %s", topic)
        for idx in sort_idx[:3]:
            logging.info("> %s", statements[idx].text)

if __name__ == "__main__":
    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument("path_statements", type=Path)
    parser.add_argument("path_sequences", type=Path)
    parser.add_argument("path_embeddings", type=Path)
    parser.add_argument("path_tokenizer", type=Path)
    parser.add_argument("--path_log", type=Path)

    main(**vars(parser.parse_args()))

max_stars_repo_name: gmcgoldr/theissues

max_issues_repo_name: gmcgoldr/theissues

max_forks_repo_name: gmcgoldr/theissues


text:
import numpy as np
import os
import keras
import cv2

model = keras.models.load_model("model.h5")
current_path = os.path.join("sim_data", "tl_00012.png")
image = np.array(cv2.imread(current_path))
image_array = np.asarray(image)
pred = model.predict(image_array[None,:,:,:], batch_size=1)
print pred
print np.argmax(pred)

max_stars_repo_name: mbmiroslaw/CarND-Capstone-master

max_issues_repo_name: mbmiroslaw/CarND-Capstone-master

max_forks_repo_name: mbmiroslaw/CarND-Capstone-master


text:
# Copyright 2019-2020 ETH Zurich and the DaCe authors. All rights reserved.
import dace
import numpy as np

sdfg = dace.SDFG('internal')
sdfg.add_array('inp', [2], dace.float32)
state = sdfg.add_state()
t = state.add_tasklet('p', {'i'}, set(), 'printf("hello world %f\\n", i)')
r = state.add_read('inp')
state.add_edge(r, None, t, 'i', dace.Memlet.simple('inp', '1'))

@dace.program
def caller(A: dace.float32[4]):
    sdfg(inp=A[1:3])

def test():
    A = np.random.rand(4).astype(np.float32)
    caller(A)
    print('Should print', A[2])

if __name__ == '__main__':
    test()

max_stars_repo_name: tobiasholenstein/dace

max_issues_repo_name: tobiasholenstein/dace

max_forks_repo_name: tobiasholenstein/dace


text:
import bson
import pymongo
import json
from bson import ObjectId
from pymongo import MongoClient
import string
import tangelo

import networkx as nx

def run(host,database,graphname,centername,degree=1):
    # Create an empty response object.
    response = {}
    collectionNames = []
    centerid = 0

    print 'looking for one-hop around:',centername
   # this method traverses the documents in the selected graph collection and builds a JSON object
   # that represents the graph to the application.  It might be faster to adopt to using a standard
   # networkX JSON description, but this is certainly simple and flexible for an initial prototype.

    client = MongoClient(host, 27017)
    db = client[database]
    # get a list of all collections (excluding system collections)
    collection = db[graphname]

    # create an empty graph to fill from the collection
    graph = nx.Graph()

    # loop through the records in the network and take the appropriate action for each type
    for record in collection.find({'type':'node'},{'_id':0}):
        # processing for nodes.  Add node to graph, then add attributes
            #print 'node:',record
        if 'name' in record['data']:
            recordname = record['data']['name']
        elif 'username' in record['data']:
            recordname = record['data']['username']
        elif 'name' in record:
            recordname = record['name']
        else:
            recordname = ''

        if recordname == centername:
                centerid = record['data']['id']
                centername = recordname
                graph.add_node(centerid)
                #print 'one hop: found node id:',centerid, 'node record:',record
                for attrib in record['data']:
                    # don't need the test here because of _id suppression at the query above
                    #if (attrib != '_id'):
                    graph.node[centerid][attrib] = record['data'][attrib]
                    #print 'added',attrib,'to ',centerid,' record now:',graph.node[centerid]
    #print nx.info(graph,n=centerid)

    # go through the links in the network and output any which are connected to te central node. Also output
    # the other node involved so a complete one-hop graph is generated

    query = {'$and': [{'type':'link'},{'$or': [{'data.source':centerid},{'data.target':centerid}]}]}
    #for record in collection.find({'type':'link'},{'_id':0}):
    for record in collection.find(query,{'_id':0}):
        source = record['data']['source']
        target = record['data']['target']
        print 'source,target, centerid:',source,target,centerid
        edgeToAdd = False
        # examine this edge to see if it connects to the center node, if so, add it to the output graph
        if (source == centerid):
            graph.add_edge(source,target)
            nodetoadd = target
            edgeToAdd = True
        if  (target == centerid):
            graph.add_edge(source,target)
            nodetoadd = source
            edgeToAdd = True

        # this is called when a node is entererd into the one-hop graph by edge reference. We need
        # to do a mongo query to get the full attribute information for the node.  There will only be
        # one node record, so we can index to [0] in the mongo cursor instance to pull out the record.

        if edgeToAdd == True:
            graph.add_node(nodetoadd)
            noderecord = collection.find({'type':'node','data.id':nodetoadd},{'_id':0})
            for attrib in noderecord[0]['data']:
                #print 'adding', attrib, ' to node',nodetoadd
                graph.node[nodetoadd][attrib] = noderecord[0]['data'][attrib]

        # now add any attributes that were on the edge
        for attrib in record['data']:
            if (attrib not in ['source','target']):
                graph.edge[source][target][attrib] = record['data'][attrib]

    # prepare the records to return to the calling program.  We build a list of nodes
    # and a list of edges here.  All the original node data for each node is added as
    # a 'data' field to preserve all attributes.  We might decide on a better way later or
    # decide to just return the JSON version of the graph directly instead of building edge and
    # node lists.

    print 'one-hop graph finished'
    print nx.info(graph)

    fixedNodes = []
    for node in graph.nodes(data=True):
        #print 'one-hop node:',node
        if 'name' in node[1]:
            fixedNodes.append({'name' : node[1]['name'], 'id': node[1]['id'],'data': node, 'group':1})
        else:
            fixedNodes.append({'name' : node[1]['id'],  'id': node[1]['id'],'data': node, 'group':2})

    fixedEdges = []
    for edge in graph.edges():
        print "found edge: ",edge
        fixedEdges.append({'source':int(edge[0]),'target': int(edge[1]),'weight':1})
    #print fixedEdges
    #print fixedNodes

    # add a "matched" node to test rendering colors artificially
    #fixedNodes[24]['matched'] = 12

    # Pack the results into the response object, and return it.
    response['result'] = {}
    response['edges'] = fixedEdges
    response['nodes'] = fixedNodes

    client.close()

    print 'response returned:'
    print response

    # Return the response object.
    #tangelo.log(str(response))
    return json.dumps(response)

max_stars_repo_name: curtislisle/PublicationNetworkViz

max_issues_repo_name: curtislisle/PublicationNetworkViz

max_forks_repo_name: curtislisle/PublicationNetworkViz

reason: import networkx


text:
from fusedwind.plant_flow.vt import GenericWindTurbineVT, GenericWindTurbinePowerCurveVT, \
    ExtendedWindTurbinePowerCurveVT, WeibullWindRoseVT, GenericWindRoseVT, GenericWindFarmTurbineLayout, WTPC
from fusedwind.plant_flow.comp import WeibullWindRose
from fusedwind.fused_helper import init_container

from random import random
from numpy import array, vstack, linspace, pi, floor
import numpy as np
from numpy import ndarray, array, loadtxt, log, zeros, cos, arccos, sin, nonzero, argsort, NaN, mean, ones, vstack, linspace, exp, arctan, arange
from numpy import pi, sqrt, dot, diff

wr_example = {
 'wind_directions': [   0.,   30.,   60.,   90.,  120.,  150.,  180.,  210.,  240.,  270.,  300.,  330.],
 'frequency_array': array([[  5.98508681e-04,   1.29125597e-03,   1.40644367e-03, 1.50695738e-03,   1.59075031e-03,   1.65624153e-03, 1.70235893e-03,   1.72856490e-03,   1.73486257e-03, 1.72178180e-03,   1.69034515e-03,   1.64201552e-03, 1.57862816e-03,   1.50231057e-03,   1.41539443e-03, 1.32032425e-03,   1.21956703e-03,   1.11552739e-03, 1.01047187e-03,   9.06465283e-04,   8.05321315e-04, 7.08568401e-04,   6.17431098e-04,   5.32826166e-04, 4.55371884e-04,   3.85408522e-04,   3.23027505e-04, 2.68106646e-04,   2.20348820e-04,   1.79321659e-04, 1.44496131e-04,   1.15282311e-04,   9.10610673e-05, 7.12108611e-05,   5.51292806e-05,   4.22493037e-05, 3.20505858e-05,   2.40662917e-05,   1.78861305e-05, 1.31563231e-05,   9.57723654e-06,   6.89937639e-06, 4.91835204e-06,   3.46932973e-06,   2.42138393e-06, 1.67205154e-06,   1.14229776e-06,   7.72019618e-07, 5.16145719e-07,   1.88493405e-07],
                           [  5.37841890e-04,   1.17137781e-03,   1.29209566e-03, 1.40243194e-03,   1.50017830e-03,   1.58343336e-03, 1.65065464e-03,   1.70070122e-03,   1.73286418e-03, 1.74688276e-03,   1.74294495e-03,   1.72167214e-03, 1.68408831e-03,   1.63157510e-03,   1.56581497e-03, 1.48872509e-03,   1.40238526e-03,   1.30896313e-03, 1.21064049e-03,   1.10954356e-03,   1.00768041e-03, 9.06887902e-04,   8.08789626e-04,   7.14766123e-04, 6.25937453e-04,   5.43157779e-04,   4.67020940e-04, 3.97875490e-04,   3.35847385e-04,   2.80868263e-04, 2.32707260e-04,   1.91004359e-04,   1.55303494e-04, 1.25083907e-04,   9.97885749e-05,   7.88488978e-05, 6.17051693e-05,   4.78226658e-05,   3.67034680e-05, 2.78943290e-05,   2.09910581e-05,   1.56399791e-05, 1.15370570e-05,   8.42528041e-06,   6.09084297e-06, 4.35860280e-06,   3.08721638e-06,   2.16426127e-06, 1.50157853e-06,   5.64510272e-07],
                           [  7.34384828e-04,   1.59522107e-03,   1.75338058e-03, 1.89615822e-03,   2.02065184e-03,   2.12442644e-03, 2.20558187e-03,   2.26280528e-03,   2.29540446e-03, 2.30331951e-03,   2.28711186e-03,   2.24793044e-03, 2.18745677e-03,   2.10783128e-03,   2.01156471e-03, 1.90143894e-03,   1.78040205e-03,   1.65146286e-03, 1.51758962e-03,   1.38161747e-03,   1.24616827e-03, 1.11358564e-03,   9.85887069e-04,   8.64733665e-04, 7.51417342e-04,   6.46864186e-04,   5.51651985e-04, 4.66039404e-04,   3.90003912e-04,   3.23285473e-04, 2.65433096e-04,   2.15851592e-04,   1.73846289e-04, 1.38663902e-04,   1.09528303e-04,   8.56703870e-05, 6.63517639e-05,   5.08823379e-05,   3.86322096e-05, 2.90385320e-05,   2.16081112e-05,   1.59165946e-05, 1.16050941e-05,   8.37502495e-06,   5.98185836e-06, 4.22836268e-06,   2.95779353e-06,   2.04737008e-06, 1.40226805e-06,   5.21865933e-07],
                           [  8.33549350e-04,   1.84169114e-03,   2.06792348e-03, 2.28187986e-03,   2.47872489e-03,   2.65397760e-03, 2.80365465e-03,   2.92440538e-03,   3.01363035e-03, 3.06957564e-03,   3.09139654e-03,   3.07918593e-03, 3.03396438e-03,   2.95763173e-03,   2.85288197e-03, 2.72308551e-03,   2.57214521e-03,   2.40433379e-03, 2.22412177e-03,   2.03600497e-03,   1.84434078e-03, 1.65320138e-03,   1.46625052e-03,   1.28664892e-03, 1.11699093e-03,   9.59273093e-04,   8.14893252e-04, 6.84676880e-04,   5.68926093e-04,   4.67485808e-04, 3.79821142e-04,   3.05100195e-04,   2.42276867e-04, 1.90169161e-04,   1.47529438e-04,   1.13104229e-04, 8.56822934e-05,   6.41306205e-05,   4.74189068e-05, 3.46336791e-05,   2.49836383e-05,   1.77980092e-05, 1.25197133e-05,   8.69506677e-06,   5.96149641e-06, 4.03449922e-06,   2.69478465e-06,   1.77626040e-06, 1.15527504e-06,   4.12261443e-07],
                           [  8.75649591e-04,   1.96089486e-03,   2.23844060e-03, 2.50798829e-03,   2.76299025e-03,   2.99704926e-03, 3.20413948e-03,   3.37883326e-03,   3.51652048e-03, 3.61360643e-03,   3.66767480e-03,   3.67760359e-03, 3.64362439e-03,   3.56731878e-03,   3.45154953e-03, 3.30032921e-03,   3.11863272e-03,   2.91216488e-03, 2.68709730e-03,   2.44979088e-03,   2.20652129e-03, 1.96322419e-03,   1.72527478e-03,   1.49731330e-03, 1.28312416e-03,   1.08557161e-03,   9.06590931e-04, 7.47229760e-04,   6.07731093e-04,   4.87647468e-04, 3.85974642e-04,   3.01293307e-04,   2.31908402e-04, 1.75977466e-04,   1.31621718e-04,   9.70160297e-05, 7.04562783e-05,   5.04046218e-05,   3.55147821e-05, 2.46404880e-05,   1.68307569e-05,   1.13157851e-05, 7.48694566e-06,   4.87388896e-06,   3.12111215e-06, 1.96570164e-06,   1.21733999e-06,   7.41145938e-07, 4.43511188e-07,   1.47913995e-07],
                           [  9.21833138e-04,   2.02972715e-03,   2.26792759e-03, 2.48924710e-03,   2.68828202e-03,   2.86016852e-03, 3.00075514e-03,   3.10675755e-03,   3.17588458e-03, 3.20692676e-03,   3.19980011e-03,   3.15554126e-03, 3.07625259e-03,   2.96499978e-03,   2.82566704e-03, 2.66277825e-03,   2.48129442e-03,   2.28639940e-03, 2.08328616e-03,   1.87695569e-03,   1.67203910e-03, 1.47265146e-03,   1.28228339e-03,   1.10373323e-03, 9.39079988e-04,   7.89694167e-04,   6.56281665e-04, 5.38953963e-04,   4.37317052e-04,   3.50571223e-04, 2.77614200e-04,   2.17141074e-04,   1.67735711e-04, 1.27949878e-04,   9.63677749e-05,   7.16551869e-05, 5.25936127e-05,   3.81007254e-05,   2.72391342e-05, 1.92157716e-05,   1.33743058e-05,   9.18285657e-06, 6.21901467e-06,   4.15380956e-06,   2.73587799e-06, 1.77670735e-06,   1.13749101e-06,   7.17855181e-07, 4.46505085e-07,   1.53770564e-07],
                           [  1.16901958e-03,   2.57123904e-03,   2.86920047e-03, 3.14535195e-03,   3.39300616e-03,   3.60617562e-03, 3.77978313e-03,   3.90984848e-03,   3.99363885e-03, 4.02977166e-03,   4.01826203e-03,   3.96051003e-03, 3.85922694e-03,   3.71830354e-03,   3.54262752e-03, 3.33786020e-03,   3.11018562e-03,   2.86604639e-03, 2.61188155e-03,   2.35388078e-03,   2.09776781e-03, 1.84862318e-03,   1.61075328e-03,   1.38760920e-03, 1.18175505e-03,   9.94882455e-04,   8.27865106e-04, 6.80845117e-04,   5.53342092e-04,   4.44375282e-04, 3.52589813e-04,   2.76379013e-04,   2.13996464e-04, 1.63653140e-04,   1.23596898e-04,   9.21732637e-05, 6.78679450e-05,   4.93326593e-05,   3.53966246e-05, 2.50665048e-05,   1.75177046e-05,   1.20797737e-05, 8.21835936e-06,   5.51572541e-06,   3.65138564e-06, 2.38394737e-06,   1.53484965e-06,   9.74342318e-07, 6.09789869e-07,   2.11057434e-07],
                           [  1.29167899e-03,   2.84708713e-03,   3.18980733e-03, 3.51603080e-03,   3.81940501e-03,   4.09401758e-03, 4.33455093e-03,   4.53642998e-03,   4.69595418e-03, 4.81040657e-03,   4.87813326e-03,   4.89858844e-03, 4.87234169e-03,   4.80104624e-03,   4.68736902e-03, 4.53488515e-03,   4.34794164e-03,   4.13149657e-03, 3.89094109e-03,   3.63191281e-03,   3.36010887e-03, 3.08110718e-03,   2.80020331e-03,   2.52226949e-03, 2.25164048e-03,   1.99202970e-03,   1.74647673e-03, 1.51732586e-03,   1.30623376e-03,   1.11420268e-03, 9.41634913e-04,   7.88403451e-04,   6.53933420e-04, 5.37289163e-04,   4.37262064e-04,   3.52454934e-04, 2.81359597e-04,   2.22425166e-04,   1.74115486e-04, 1.34954997e-04,   1.03563111e-04,   7.86777442e-05, 5.91691442e-05,   4.40454359e-05,   3.24514605e-05, 2.36624992e-05,   1.70743903e-05,   1.21913870e-05, 8.61289935e-06,   3.28301641e-06],
                           [  1.38975021e-03,   3.05388147e-03,   3.41196792e-03, 3.75586418e-03,   4.08021930e-03,   4.38002322e-03, 4.65070152e-03,   4.88820790e-03,   5.08910961e-03, 5.25066120e-03,   5.37086331e-03,   5.44850321e-03, 5.48317510e-03,   5.47527889e-03,   5.42599698e-03, 5.33724976e-03,   5.21163111e-03,   5.05232625e-03, 4.86301505e-03,   4.64776411e-03,   4.41091185e-03, 4.15695061e-03,   3.89040996e-03,   3.61574527e-03, 3.33723508e-03,   3.05889041e-03,   2.78437835e-03, 2.51696181e-03,   2.25945622e-03,   2.01420338e-03, 1.78306198e-03,   1.56741358e-03,   1.36818226e-03, 1.18586592e-03,   1.02057680e-03,   8.72088629e-04, 7.39888098e-04,   6.23227994e-04,   5.21180074e-04, 4.32685696e-04,   3.56602728e-04,   2.91747592e-04, 2.36931694e-04,   1.90991816e-04,   1.52814395e-04, 1.21353871e-04,   9.56455150e-05,   7.48132906e-05, 5.80734310e-05,   2.38465774e-05],
                           [  1.17262387e-03,   2.60943924e-03,   2.96153237e-03, 3.30851045e-03,   3.64475920e-03,   3.96476457e-03, 4.26322159e-03,   4.53514933e-03,   4.77600674e-03, 4.98180415e-03,   5.14920535e-03,   5.27561545e-03, 5.35925035e-03,   5.39918426e-03,   5.39537250e-03, 5.34864811e-03,   5.26069168e-03,   5.13397512e-03, 4.97168137e-03,   4.77760292e-03,   4.55602326e-03, 4.31158587e-03,   4.04915610e-03,   3.77368147e-03, 3.49005569e-03,   3.20299167e-03,   2.91690796e-03, 2.63583236e-03,   2.36332540e-03,   2.10242534e-03, 1.85561537e-03,   1.62481230e-03,   1.41137541e-03, 1.21613315e-03,   1.03942478e-03,   8.81153559e-04, 7.40848127e-04,   6.17728441e-04,   5.10773070e-04, 4.18784878e-04,   3.40452671e-04,   2.74406873e-04, 2.19267908e-04,   1.73686495e-04,   1.36375601e-04, 1.06134229e-04,   8.18635783e-05,   6.25763984e-05, 4.74005026e-05,   1.90836763e-05],
                           [  8.24556184e-04,   1.83763158e-03,   2.08931937e-03, 2.33782881e-03,   2.57906571e-03,   2.80899131e-03, 3.02370408e-03,   3.21952653e-03,   3.39309331e-03, 3.54143638e-03,   3.66206360e-03,   3.75302687e-03, 3.81297659e-03,   3.84119948e-03,   3.83763784e-03, 3.80288855e-03,   3.73818169e-03,   3.64533906e-03, 3.52671394e-03,   3.38511469e-03,   3.22371493e-03, 3.04595417e-03,   2.85543299e-03,   2.65580696e-03, 2.45068367e-03,   2.24352684e-03,   2.03757097e-03, 1.83574943e-03,   1.64063812e-03,   1.45441579e-03, 1.27884152e-03,   1.11524876e-03,   9.64554676e-04, 8.27283012e-04,   7.03597966e-04,   5.93346558e-04, 4.96106649e-04,   4.11237894e-04,   3.37933106e-04, 2.75267810e-04,   2.22246152e-04,   1.77841801e-04, 1.41032921e-04,   1.10830721e-04,   8.63015293e-05, 6.65826278e-05,   5.08923720e-05,   3.85353060e-05, 2.89030828e-05,   1.15455611e-05],
                           [  7.32521276e-04,   1.58062566e-03,   1.72456794e-03, 1.85424825e-03,   1.96761259e-03,   2.06297576e-03, 2.13905365e-03,   2.19498595e-03,   2.23034790e-03, 2.24514978e-03,   2.23982416e-03,   2.21520098e-03, 2.17247154e-03,   2.11314284e-03,   2.03898398e-03, 1.95196698e-03,   1.85420428e-03,   1.74788524e-03, 1.63521428e-03,   1.51835254e-03,   1.39936521e-03, 1.28017592e-03,   1.16252948e-03,   1.04796352e-03, 9.37789483e-04,   8.33082624e-04,   7.34680586e-04, 6.43189695e-04,   5.58997846e-04,   4.82292734e-04, 4.13084095e-04,   3.51228587e-04,   2.96456022e-04, 2.48395784e-04,   2.06602397e-04,   1.70579417e-04, 1.39801019e-04,   1.13730832e-04,   9.18378043e-05, 7.36090050e-05,   5.85594467e-05,   4.62391216e-05, 3.62375226e-05,   2.81859826e-05,   2.17581928e-05, 1.66692633e-05,   1.26736799e-05,   9.56247777e-06, 7.15991731e-06,   2.85976559e-06]]),
 'wind_speeds': [  4.        ,   4.42857143,   4.85714286,   5.28571429,
                  5.71428571,   6.14285714,   6.57142857,   7.        ,
                  7.42857143,   7.85714286,   8.28571429,   8.71428571,
                  9.14285714,   9.57142857,  10.        ,  10.42857143,
                 10.85714286,  11.28571429,  11.71428571,  12.14285714,
                 12.57142857,  13.        ,  13.42857143,  13.85714286,
                 14.28571429,  14.71428571,  15.14285714,  15.57142857,
                 16.        ,  16.42857143,  16.85714286,  17.28571429,
                 17.71428571,  18.14285714,  18.57142857,  19.        ,
                 19.42857143,  19.85714286,  20.28571429,  20.71428571,
                 21.14285714,  21.57142857,  22.        ,  22.42857143,
                 22.85714286,  23.28571429,  23.71428571,  24.14285714,
                 24.57142857,  25.        ]}

def generate_a_valid_wt(D = 200*random()):
    wt_desc = GenericWindTurbineVT()
    wt_desc.rotor_diameter = D
    wt_desc.hub_height = D * (0.5 + random())
    return wt_desc

def generate_random_GenericWindTurbinePowerCurveVT(D=None):
    """
    Generate a random turbine and power curve using the GenericWindTurbinePowerCurveVT class

    Parameters
    ----------
    D       float, default=random, (optional)
            The wind turbine rotor diameter

    Returns
    -------
    wt_desc GenericWindTurbinePowerCurveVT
            A random wind turbine power curve and c_t curve variable tree
    """
    if not D:
        D = 200*random()

    wt_desc = GenericWindTurbinePowerCurveVT()
    wt_desc.rotor_diameter = D
    wt_desc.hub_height = D * (0.5 + random())
    wt_desc.cut_in_wind_speed = 2. + 4. * random()
    wt_desc.cut_out_wind_speed = 20. + 10. * random()
    rho = 1.225
    rated_wind_speed = 8. + 4. * random()
    max_a = 0.333 * random()
    max_cp = 4 * max_a * (1 - max_a)**2.
    max_ct = 4 * max_a * (1 - max_a)
    A = 0.25 * pi * D**2. # Rotor area
    ideal_power = lambda ws: 0.5 * rho * A * max_cp * ws **3.
    real_power = lambda ws: ideal_power(ws) if ws < rated_wind_speed else ideal_power(rated_wind_speed)
    #a_ct = -sqrt(-c_t + 1)/2 + 1/2
    ct_from_cp = lambda cp: min(0.89, cp  * 2.)
    cp_from_power = lambda power, ws: power/(0.5 * rho * A * ws**3.)
    ct_from_power = lambda pws: ct_from_cp(cp_from_power(pws[0], pws[1]))
    N = 3+int(random() * 100)
    ws = linspace(wt_desc.cut_in_wind_speed, wt_desc.cut_out_wind_speed, N)
    wt_desc.power_curve = vstack([ws, map(real_power, ws)]).T
    wt_desc.c_t_curve = vstack([ws, map(ct_from_power, zip(wt_desc.power_curve[:,1],ws))]).T
    wt_desc.power_rating = ideal_power(rated_wind_speed)
    wt_desc.rated_wind_speed = rated_wind_speed
    wt_desc.air_density = rho
    wt_desc.test_consistency()
    return wt_desc

def generate_random_WTPC(D=None, **kwargs):
    wt = generate_random_GenericWindTurbinePowerCurveVT(D)
    inputs = {k: getattr(wt, k) for k in wt.list_vars()}
    inputs['position'] = np.random.rand(2)
    inputs['name'] = 'wt'
    # Replacing the parameters by the inputs of the function
    for k, v in kwargs.iteritems():
        inputs[k] = v
    return WTPC(**inputs)

def generate_random_wt_positions(D=None, nwt=None, min_D=None):
    """
    Generate a random layout of wind turbines

    Parameters
    ----------
    D       float, default=random, (optional)
            The wind turbine rotor diameter

    nwt     int, default=random, (optional)
            The number of turbines in the layout

    min_D   float, default=random, (optional)
            The minimum of rotor diameter between turbines

    Returns
    ------
    wt_positions    ndarray [nwt, 2]
                    The (x,y) position of the turbines

    Example
    -------
    >>> generate_random_wt_positions(D=82., nwt=5, min_D=7.)
    array([[  0.        ,   0.        ],
           [ -1.74906289, -30.4628047 ],
           [ -5.88793386,  -7.11251901],
           [-35.53857721,  32.42632383],
           [ 15.25619612,   0.79847237],
           [  1.20034923, -10.73881476]])

    """

    if not D:
        D = 200*random()
    if not nwt:
        nwt = int(30*random())
    if not min_D:
        min_D = 3. + 5.*random()

    Nnorm = lambda x: sqrt(x[0]**2. + x[1]**2.)
    min_dist = lambda pos_arr, pos: min(map(Nnorm, pos_arr - pos))
    def random_path(wt_positions):
        N = floor(random()*len(wt_positions))
        #print N
        x0, y0 = wt_positions[N,:]
        count = 0
        while min_dist(wt_positions, array([x0,y0])) < min_D and count < 20:
            x0, y0 =  x0 + (0.5-random()) * D, y0 + (0.5-random()) * D
            count +=1
        return vstack([wt_positions, [x0, y0]])

    wt_positions = array([[0.,0.]])
    for i in range(nwt-1):
        wt_positions = random_path(wt_positions)
    return wt_positions

def generate_random_GenericWindRoseVT():
    """
    Generate a random GenericWindRoseVT object

    Parameters
    ----------

    N/A

    Returns
    -------
    wind_rose   GenericWindRoseVT
                A wind rose variable tree

    """
    weibull_array = np.array([[  0.00000000e+00,   3.59673400e-02,  9.22422800e+00,   2.38867200e+00],
                             [  3.00000000e+01,   3.94977300e-02,   9.86435600e+00,   2.44726600e+00],
                             [  6.00000000e+01,   5.17838000e-02,   9.65220200e+00,   2.41992200e+00],
                             [  9.00000000e+01,   6.99794900e-02,   9.98217800e+00,   2.58789100e+00],
                             [  1.20000000e+02,   8.36383000e-02,   1.00946000e+01,   2.74804700e+00],
                             [  1.50000000e+02,   6.43412500e-02,   9.64369000e+00,   2.59179700e+00],
                             [  1.80000000e+02,   8.64220000e-02,   9.63377500e+00,   2.58007800e+00],
                             [  2.10000000e+02,   1.17690000e-01,   1.05678600e+01,   2.54492200e+00],
                             [  2.40000000e+02,   1.51555100e-01,   1.14525200e+01,   2.46679700e+00],
                             [  2.70000000e+02,   1.47361100e-01,   1.17423700e+01,   2.60351600e+00],
                             [  3.00000000e+02,   1.00109800e-01,   1.16923200e+01,   2.62304700e+00],
                             [  3.30000000e+02,   5.16542400e-02,   1.01385800e+01,   2.32226600e+00]])
    ### -> This is going to be really slow because of interpolation
    nwd, nws = 4+int(random()*360), 4+int(random()*25)
    wd = linspace(0., 360., nwd)[:-1].tolist()
    ws = linspace(3., 30., nws).tolist()
    #wind_rose_array = array([[wd_, random(), random()*15., random()*4.] for wd_ in wd])
    #wind_rose_array *= [1., 1./wind_rose_array[:,1].sum(), 1., 1.]
    #return WeibullWindRose()(wind_directions= wd,
    #                         wind_speeds=ws,
    #                         wind_rose_array=wind_rose_array).wind_rose
    #
    gwr = GenericWindRoseVT(wind_directions=wd, wind_speeds=ws, weibull_array=weibull_array)
    return gwr

# def generate_random_wt_layout(D=None, nwt=None):
#     """
#     Generate a random GenericWindFarmTurbineLayout of wind turbines
#
#     Parameters
#     ----------
#     D       float, default=random, (optional)
#             The wind turbine rotor diameter
#
#     nwt     int, default=random, (optional)
#             The number of turbines in the layout
#
#     Returns
#     -------
#     wt_layout   GenericWindFarmTurbineLayout
#                 A random wind turbine layout
#
#     """
#     if D is None:
#         D = 200*random()
#     if nwt is None:
#         nwt = int(30*random())
#     wt_layout = GenericWindFarmTurbineLayout()
#     wt_layout.wt_positions = generate_random_wt_positions(D=D, nwt=nwt)
#     wt_layout.wt_list = [generate_random_GenericWindTurbinePowerCurveVT(D) for i in range(nwt)]
#     wt_layout.wt_wind_roses = [generate_random_GenericWindRoseVT() for i in range(nwt)]
#     wt_layout.wt_names = ['wt_%d'%(i) for i in range(nwt)]
#     wt_layout.test_consistency()
#     return wt_layout

def generate_random_wt_layout(D=None, nwt=None):
    """
    Parameters
    ----------
    D       float, default=random, (optional)
            The wind turbine rotor diameter

    nwt     int, default=random, (optional)
            The number of turbines in the layout

    Returns
    -------
    a random GenericWindFarmTurbineLayout instance
    """
    if not D:
        D = 200*random()
    if not nwt:
        nwt = int(30*random())
    wt_positions = generate_random_wt_positions(D=D, nwt=nwt)
    return GenericWindFarmTurbineLayout([generate_random_WTPC(
        name='wt%d'%(n),
        position=wt_positions[n,:],
        wind_rose=generate_random_GenericWindRoseVT()) for n in range(nwt)])

max_stars_repo_name: FUSED-Wind/fusedwind

max_issues_repo_name: michaelXDzhang/fusedwind

max_forks_repo_name: michaelXDzhang/fusedwind

reason: import numpy,from numpy


text:
# __filename__: background_function.py
#
# __description__: Functions used to set background
#
# __remark__:
#
# __todos__:
#
# Created by Tobias Wenzel in August 2017
# Copyright (c) 2017 Tobias Wenzel

import numpy as np

def get_xsqr(x: np.ndarray,
             y: np.ndarray):
    """

    :param x:
    :param y:
    :return:
    """
    m_size = len(x)
    x = x ** 2
    for i in range(m_size):
        y[i, :] = x[i] + x
    return y

def get_rastrigin(x: np.ndarray,
                  y: np.ndarray):
    """
        Recommended range: -5.2, 5.2

    -5.2, 5.2
    :param x:
    :param y:
    :return:
    """
    m_size = len(x)
    for i in range(m_size):
        for j in range(m_size):
            y[i, j] = (x[i] ** 2 - 10 * np.cos(np.pi * x[i]) + 10) \
                      + (x[j] ** 2 - 10 * np.cos(np.pi * x[j]) + 10)
    return y

def get_rosenbrock(x: np.ndarray,
                   y: np.ndarray):
    """
            Recommended range: -10, 10
    :param x:
    :param y:
    :return:
    """
    m_size = len(x)
    for i in range(m_size):
        for j in range(m_size):
            a = 1. - x[i]
            b = x[j] - x[i] ** 2
            y[i, j] = a * a + b * b * 100.
    return y

def get_griewank(x: np.ndarray,
                 y: np.ndarray):
    """
                Recommended range: -10, 10
    :param x:
    :param y:
    :return:
    """
    m_size = len(x)
    for i in range(m_size):
        for j in range(m_size):
            prod = (np.cos(x[i] / 0.0000001) + 1) * (np.cos(x[j] / 1) + 1)
            y[i, j] = (1 / 4000) * (x[i] ** 2 - prod + x[j] ** 2 - prod)
    return y

def get_schaffer_f6(x: np.ndarray,
                    y: np.ndarray):
    """
            Recommended range: -10, 10
    :param x:
    :param y:
    :return:
    """
    m_size = len(x)
    for i in range(m_size):
        for j in range(m_size):
            y[i, j] = (0.5 - ((np.sin(np.sqrt(x[i] ** 2 + x[j] ** 2)) ** 2 - 0.5) \
                             / (1 + 0.001 * (x[i] ** 2 + x[j] ** 2) ** 2)))
    return y

def get_eggholder(x: np.ndarray,
                  y: np.ndarray):
    """
        Recommended range: -512, 512

    :return x
    :return y:
    """
    m_size = len(x)
    for i_1 in range(m_size):
        for i_2 in range(m_size):
            y[i_1, i_2] -= (x[i_2] + 47) \
                           * np.sin(np.sqrt(np.abs((x[i_1] / 2) + x[i_2] + 47))) \
                           - x[i_1] \
                           * np.sin(np.sqrt(np.abs(x[i_1] - (x[i_2] + 47))))
    return y

def get_hoelder_table(x: np.ndarray,
                      y: np.ndarray):
    """
    :return x
    :return y:
    """
    m_size = len(x)
    for i_1 in range(m_size):
        for i_2 in range(m_size):
            y[i_1, i_2] = -np.abs(np.sin(x[i_1]) * np.cos(x[i_2]) \
                                  * np.exp(np.abs(1 - (np.sqrt(x[i_1] ** 2 + x[i_2] ** 2) / np.pi))))
    return y

def get_styblinsky_tang(x: np.ndarray,
                        y: np.ndarray):
    """
    @not yet working (?)
    :param x:
    :param y:
    :return:
    """
    m_size = len(x)
    f = lambda d: (d ** 4) - (16 * (d**2)) + (5 * d)
    f_arr = np.frompyfunc(f, 1, 1)
    for i_1 in range(m_size):
        y[i_1, :] = f_arr(x[i_1])
        for i_2 in range(m_size):
            y[i_1, i_2] += f(x[i_2])
    return y

background_function = {}

background_function['square'] = get_xsqr
background_function['rastrigin'] = get_rastrigin
background_function['schaffer_f6'] = get_schaffer_f6
background_function['griewank'] = get_griewank
background_function['rosenbrock'] = get_rosenbrock
background_function['eggholder'] = get_eggholder
background_function['hoelder_table'] = get_hoelder_table
background_function['styblinsky_tang'] = get_styblinsky_tang

def generate_2d_background(func_name: str, n: int):
    """

    :param func_name:
    :param n:
    :return:
    """
    m_size = 2*n
    bg_fn = np.zeros((m_size, m_size))
    vals = np.linspace(-n, n, m_size)
    print(f"background-function: {func_name}")
    return background_function[func_name](vals, bg_fn)


text:
#!/usr/bin/env python
#-*- coding: utf-8 -*-

#Helper functions for TF-STAR

import sqlite3
import pandas as pd
import numpy as np
import re
import os
import sys
import random
import hashlib
import time
from collections import defaultdict
from constant import *
import subprocess
import scipy as sp
import scipy.optimize
import scipy.fftpack
import scipy.stats
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
from numpy import array, empty

def check_input_id(in_file):
    '''input: gene list file
       output: gene id set and its type (ensembl, entrez, symbol)
    '''
    df = pd.read_table(in_file, sep='\n', delimiter=None, header=0, engine='python')
    ids = set(df.iloc[:,0].map(str.strip))

    ensembl_pat = re.compile(r'^ENSG\d+$')
    entrez_pat = re.compile(r'^\d+$')

    ensembl = 0
    entrez  = 0
    symbol = 0

    ids_num = len(ids)
    for i in ids:
        if ensembl_pat.match(i):
            ensembl += 1
        elif entrez_pat.match(i):
            entrez += 1
        else:
            symbol += 1
    max_num = max(ensembl,entrez,symbol)

    if max_num/float(ids_num) >= 0.9:
        if max_num == ensembl:
            in_type = 'ensembl'
        elif max_num == entrez:
            in_type = 'entrez'
        elif max_num == symbol:
            in_type = 'symbol'

    return ids, in_type

def background_genes_map(id_type, biotype):
    '''key: ID, value: symbol
    1   11868   14409   ENST00000456328 gene_name   +   ENSG00000223972 pseudogene
    1   30365   30503   NR_036051.1 MIR1302-2   +   100302278   noncoding
    '''
    biotype_dict = {
            '3prime_overlapping_ncrna' : 'noncoding',
            'antisense' : 'noncoding',
            'IG_C_gene' : 'protein_coding',
            'IG_C_pseudogene' : 'pseudogene',
            'IG_D_gene' : 'protein_coding',
            'IG_J_gene' : 'protein_coding',
            'IG_J_pseudogene' : 'pseudogene',
            'IG_V_gene' : 'protein_coding',
            'IG_V_pseudogene' : 'pseudogene',
            'lincRNA' : 'noncoding',
            'miRNA' : 'noncoding',
            'misc_RNA' : 'noncoding',
            'Mt_rRNA' : 'noncoding',
            'Mt_tRNA' : 'noncoding',
            'polymorphic_pseudogene' : 'protein_coding',
            'processed_transcript' : 'noncoding',
            'protein_coding' : 'protein_coding',
            'pseudogene' : 'pseudogene',
            'rRNA' : 'noncoding',
            'sense_intronic' : 'noncoding',
            'sense_overlapping' : 'noncoding',
            'snoRNA' : 'noncoding',
            'snRNA' : 'noncoding',
            'TR_C_gene' : 'protein_coding',
            'TR_D_gene' : 'protein_coding',
            'TR_J_gene' : 'protein_coding',
            'TR_J_pseudogene' : 'pseudogene',
            'TR_V_gene' : 'protein_coding',
            'TR_V_pseudogene' : 'pseudogene',
            'noncoding' : 'noncoding', # NCBI refseq gene
    }
    # biotype must in {'protein_coding', 'noncoding', 'pseudogene', 'all'}
    biotype_set = set()
    for i in biotype_dict:
        if biotype_dict[i] == biotype:
            biotype_set.add(i)
    if biotype == 'all':
        biotype_set = set(biotype_dict.keys())

    background_genes = {}

    if id_type == 'entrez':
        db_type = 'refseq'
    elif id_type == 'ensembl':
        db_type = 'ensembl'
    else:
        db_type = 'refseq'

    biotype_str = (', '.join('"' + item + '"' for item in biotype_set))

    conn = sqlite3.connect('./refs/references.db')
    cursor = conn.execute('SELECT geneid, symbol FROM {} WHERE biotype IN ({})'.format(db_type,biotype_str))
    for row in cursor:
        background_genes[row[0]] = row[1]

    conn.close()
    return background_genes

def id_converter(in_ids, id_type):
    '''return a dictionary:
    key: gene ID
    value: symbol
    '''
    out_genes = {}
    gene_to_symbol = background_genes_map(id_type, 'all')
    symbol_to_gene = {}
    for i,j in list(gene_to_symbol.items()):
        symbol_to_gene[j] = i

    for in_id in in_ids:
        if id_type == 'entrez' or id_type == 'ensembl':
            if in_id in gene_to_symbol:
                out_genes[in_id] = gene_to_symbol[in_id]
        else:
            if in_id in symbol_to_gene:
                out_genes[symbol_to_gene[in_id]] = in_id

    if len(out_genes) == 0:
        print('Please check you input, only ensembl IDs, entrez IDs and gene symbols are allowed!')
        sys.exit(1)
    return out_genes

def removefile(filename):
    try:
        os.remove(filename)
    except OSError as e:
        if e.errno != errno.ENOENT:
            raise
##############################################################################
# ------------------------------------#
# P-value estimation using KDE        #
# ------------------------------------#
def _bandwidth(data, N=None, MIN=None, MAX=None):
    '''
    An implementation of the kde bandwidth selection method outlined in:

    Z. I. Botev, J. F. Grotowski, and D. P. Kroese. Kernel density
    estimation via diffusion. The Annals of Statistics, 38(5):2916-2957, 2010.

    Based on the implementation in Matlab by Zdravko Botev.

    Daniel B. Smith, PhD
    https://github.com/Daniel-B-Smith/KDE-for-SciPy/blob/master/kde.py
    Updated 1-23-2013
    '''
    # Parameters to set up the mesh on which to calculate
    N = 2**14 if N is None else int(2**sp.ceil(sp.log2(N)))
    if MIN is None or MAX is None:
        minimum = min(data)
        maximum = max(data)
        Range = maximum - minimum
        MIN = minimum - Range/10 if MIN is None else MIN
        MAX = maximum + Range/10 if MAX is None else MAX

    # Range of the data
    R = MAX-MIN

    # Histogram the data to get a crude first approximation of the density
    M = len(data)
    DataHist, bins = sp.histogram(data, bins=N, range=(MIN,MAX))
    DataHist = DataHist/M
    DCTData = scipy.fftpack.dct(DataHist, norm=None)

    I = [iN*iN for iN in range(1, N)]
    SqDCTData = (DCTData[1:]/2)**2

    # The fixed point calculation finds the bandwidth = t_star
    guess = 0.1
    try:
        t_star = scipy.optimize.brentq(__fixed_point, 0, guess,
                                       args=(M, I, SqDCTData))
    except ValueError:
        print('Oops!')
        return None

    # Smooth the DCTransformed data using t_star
    SmDCTData = DCTData*sp.exp(-sp.arange(N)**2*sp.pi**2*t_star/2)
    # Inverse DCT to get density
    density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R
    mesh = [(bins[i]+bins[i+1])/2 for i in range(N)]
    bandwidth = sp.sqrt(t_star)*R

    density = density/sp.trapz(density, mesh)

    # return bandwidth, mesh, density
    return bandwidth

def __fixed_point(t, M, I, a2):
    l=7
    I = sp.float128(I)
    M = sp.float128(M)
    a2 = sp.float128(a2)
    f = 2*sp.pi**(2*l)*sp.sum(I**l*a2*sp.exp(-I*sp.pi**2*t))
    for s in range(l, 1, -1):
        K0 = sp.prod(list(range(1, 2*s, 2)))/sp.sqrt(2*sp.pi)
        const = (1 + (1/2)**(s + 1/2))/3
        _time=(2*const*K0/M/f)**(2/(3+2*s))
        f=2*sp.pi**(2*s)*sp.sum(I**s*a2*sp.exp(-I*sp.pi**2*_time))
    return t-(2*M*sp.sqrt(sp.pi)*f)**(-2/5)

def scipy_bandwidth(data):
    kernel = scipy.stats.gaussian_kde(data)
    f = kernel.covariance_factor()
    bw = f * data.std()
    return bw

def kde(data, bw):
    #instantiate and fit the KDE model
    kde = KernelDensity(bandwidth=bw,kernel='gaussian')
    kde.fit(data[:,np.newaxis])

    # obtain the range of data
    # and define the range of distribution
    minimum = min(data)
    maximum = max(data)
    Range = maximum - minimum
    MIN = round(minimum - Range/2)
    MAX = round(maximum + Range/2)

    data_points = np.around(np.arange(MIN, MAX, 0.01, dtype=np.float64),3)

    # the log of the probability density
    logprob = kde.score_samples(data_points[:,None])
    prob = np.exp(logprob)

    '''
    plt.fill_between(data_points, np.exp(logprob), alpha=0.5)
    plt.plot(data, np.full_like(data, -0.01), '|k', markeredgewidth=1)
    plt.show()
    '''
    #print(sp.trapz(np.exp(logprob), data_plot))
    return prob, data_points # y, x

def pvalue_calc(input_value, array_x, array_y):
    '''right-side one-tailed pvalue
    '''
    minimum = min(array_x)
    maximum = max(array_x)
    if input_value <= minimum:
        return 1.0
    if input_value > maximum:
        return 0.0
    _x = np.array(array_x)
    #_x_str = np.array(map(str,array_x))
    _y = np.array(array_y)
    density = dict(list(zip(_x, _y)))
    x_values = _x[_x >= input_value]
    y_values = [density[i] for i in x_values]
    return sp.trapz(y_values, x_values)

def best_bandwidth_by_cv(data):
    min_bound = np.mean(data)/10.0 + 1
    max_bound = np.mean(data)/10.0 - 1
    bandwidths = np.arange(min_bound,max_bound,0.1)
    grid = GridSearchCV(KernelDensity(kernel='gaussian'),
                       {'bandwidth':bandwidths},
                       cv=10,
                       n_jobs=10)
    grid.fit(data[:,None])
    return grid.best_params_

def correct_pvalues_for_multiple_testing(pvalues, correction_type="Benjamini-Hochberg"):
    """
    consistent with R - print correct_pvalues_for_multiple_testing([0.0, 0.01, 0.029, 0.03, 0.031, 0.05, 0.069, 0.07, 0.071, 0.09, 0.1])
    """
    pvalues = array(pvalues)
    n = int(pvalues.shape[0])
    new_pvalues = empty(n)
    if correction_type == "Bonferroni":
        new_pvalues = n * pvalues
    elif correction_type == "Bonferroni-Holm":
        values = [ (pvalue, i) for i, pvalue in enumerate(pvalues) ]
        values.sort()
        for rank, vals in enumerate(values):
            pvalue, i = vals
            new_pvalues[i] = (n-rank) * pvalue
    elif correction_type == "Benjamini-Hochberg":
        values = [ (pvalue, i) for i, pvalue in enumerate(pvalues) ]
        values.sort()
        values.reverse()
        new_values = []
        for i, vals in enumerate(values):
            rank = n - i
            pvalue, index = vals
            new_values.append((n/rank) * pvalue)
        for i in range(0, int(n)-1):
            if new_values[i] < new_values[i+1]:
                new_values[i+1] = new_values[i]
        for i, vals in enumerate(values):
            pvalue, index = vals
            new_pvalues[index] = new_values[i]
    return new_pvalues

def hypergeometric_pvalue(k, M, n, N):
    ''' M: total number of genes (population size)
        N: the number of input genes (sample size)
        n: total number of genes binding by the test TF (sucesses size in the population)
        k: the number of input genes binding by the test TF (sucesses size in the sample)
    # this is equivalent to "1 - phyper(k-1, n, M-n, N)"
    '''
    p = scipy.stats.hypergeom.sf(k-1, M, n, N)
    return '{0:.5g}'.format(p)

if __name__ == '__main__':
    #print(id_converter({'IKZF1','STAT4'}, 'symbol'))
    #a,b = refseq_map('entrez', 'protein_coding', '2000:200', '200:2000')
    #print(overlap_peaks(a, 'GM12878.encode.bed'))
    data = np.arange(1, 100, 0.1)
    bw = _bandwidth(data)
    print(bw)
    bw = scipy_bandwidth(data)
    print(bw)
    #bw = best_bandwidth_by_cv(data)
    #print(bw)
    #print(kde(data, bw))

reason: import numpy,from numpy,import scipy