text:
# Copyright 2018-2020 Xanadu Quantum Technologies Inc.

# 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.
"""This module contains functions to construct many-body observables whose expectation
values can be used to simulate molecular properties.
"""
# pylint: disable=too-many-arguments, too-few-public-methods
import numpy as np
from openfermion.ops import FermionOperator
from openfermion.transforms import bravyi_kitaev, jordan_wigner

from . import structure

def _spin2_matrix_elements(sz):
    r"""Builds the table of matrix elements
    :math:`\langle \bm{\alpha}, \bm{\beta} \vert \hat{s}_1 \cdot \hat{s}_2 \vert
    \bm{\gamma}, \bm{\delta} \rangle` of the two-particle spin operator
    :math:`\hat{s}_1 \cdot \hat{s}_2`.

    The matrix elements are evaluated using the expression

    .. math::

        \langle ~ (\alpha, s_{z_\alpha});~ (\beta, s_{z_\beta}) ~ \vert \hat{s}_1 &&
        \cdot \hat{s}_2 \vert ~ (\gamma, s_{z_\gamma}); ~ (\delta, s_{z_\gamma}) ~ \rangle =
        \delta_{\alpha,\delta} \delta_{\beta,\gamma} \\
        && \times \left( \frac{1}{2} \delta_{s_{z_\alpha}, s_{z_\delta}+1}
        \delta_{s_{z_\beta}, s_{z_\gamma}-1} + \frac{1}{2} \delta_{s_{z_\alpha}, s_{z_\delta}-1}
        \delta_{s_{z_\beta}, s_{z_\gamma}+1} + s_{z_\alpha} s_{z_\beta}
        \delta_{s_{z_\alpha}, s_{z_\delta}} \delta_{s_{z_\beta}, s_{z_\gamma}} \right),

    where :math:`\alpha` and :math:`s_{z_\alpha}` refer to the quantum numbers of the spatial
    function and the spin projection, respectively, of the single-particle state
    :math:`\vert \bm{\alpha} \rangle \equiv \vert \alpha, s_{z_\alpha} \rangle`.

    Args:
        sz (array[float]): spin-projection of the single-particle states

    Returns:
        array: NumPy array with the table of matrix elements. The first four columns
        contain the indices :math:`\bm{\alpha}`, :math:`\bm{\beta}`, :math:`\bm{\gamma}`,
        :math:`\bm{\delta}` and the fifth column stores the computed matrix element.

    **Example**

    >>> sz = np.array([0.5, -0.5])
    >>> print(_spin2_matrix_elements(sz))
    [[ 0.    0.    0.    0.    0.25]
     [ 0.    1.    1.    0.   -0.25]
     [ 1.    0.    0.    1.   -0.25]
     [ 1.    1.    1.    1.    0.25]
     [ 0.    1.    0.    1.    0.5 ]
     [ 1.    0.    1.    0.    0.5 ]]
    """

    n = np.arange(sz.size)

    alpha = n.reshape(-1, 1, 1, 1)
    beta = n.reshape(1, -1, 1, 1)
    gamma = n.reshape(1, 1, -1, 1)
    delta = n.reshape(1, 1, 1, -1)

    # we only care about indices satisfying the following boolean mask
    mask = np.logical_and(alpha // 2 == delta // 2, beta // 2 == gamma // 2)

    # diagonal elements
    diag_mask = np.logical_and(sz[alpha] == sz[delta], sz[beta] == sz[gamma])
    diag_indices = np.argwhere(np.logical_and(mask, diag_mask))
    diag_values = (sz[alpha] * sz[beta]).flatten()

    diag = np.vstack([diag_indices.T, diag_values]).T

    # off-diagonal elements
    m1 = np.logical_and(sz[alpha] == sz[delta] + 1, sz[beta] == sz[gamma] - 1)
    m2 = np.logical_and(sz[alpha] == sz[delta] - 1, sz[beta] == sz[gamma] + 1)

    off_diag_mask = np.logical_and(mask, np.logical_or(m1, m2))
    off_diag_indices = np.argwhere(off_diag_mask)
    off_diag_values = np.full([len(off_diag_indices)], 0.5)

    off_diag = np.vstack([off_diag_indices.T, off_diag_values]).T

    # combine the off diagonal and diagonal tables into a single table
    return np.vstack([diag, off_diag])

def spin2(electrons, orbitals, mapping="jordan_wigner", wires=None):
    r"""Computes the total spin operator :math:`\hat{S}^2`.

    The total spin operator :math:`\hat{S}^2` is given by

    .. math::

        \hat{S}^2 = \frac{3}{4}N + \sum_{ \bm{\alpha}, \bm{\beta}, \bm{\gamma}, \bm{\delta} }
        \langle \bm{\alpha}, \bm{\beta} \vert \hat{s}_1 \cdot \hat{s}_2
        \vert \bm{\gamma}, \bm{\delta} \rangle ~
        \hat{c}_\bm{\alpha}^\dagger \hat{c}_\bm{\beta}^\dagger
        \hat{c}_\bm{\gamma} \hat{c}_\bm{\delta},

    where the two-particle matrix elements are computed as,

    .. math::

        \langle \bm{\alpha}, \bm{\beta} \vert \hat{s}_1 \cdot \hat{s}_2
        \vert \bm{\gamma}, \bm{\delta} \rangle = && \delta_{\alpha,\delta} \delta_{\beta,\gamma} \\
        && \times \left( \frac{1}{2} \delta_{s_{z_\alpha}, s_{z_\delta}+1}
        \delta_{s_{z_\beta}, s_{z_\gamma}-1} + \frac{1}{2} \delta_{s_{z_\alpha}, s_{z_\delta}-1}
        \delta_{s_{z_\beta}, s_{z_\gamma}+1} + s_{z_\alpha} s_{z_\beta}
        \delta_{s_{z_\alpha}, s_{z_\delta}} \delta_{s_{z_\beta}, s_{z_\gamma}} \right).

    In the equations above :math:`N` is the number of electrons, :math:`\alpha` refer to the
    quantum numbers of the spatial wave function and :math:`s_{z_\alpha}` is
    the spin projection of the single-particle state
    :math:`\vert \bm{\alpha} \rangle \equiv \vert \alpha, s_{z_\alpha} \rangle`.
    The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation
    and annihilation operators, respectively.

    Args:
        electrons (int): Number of electrons. If an active space is defined, this is
            the number of active electrons.
        orbitals (int): Number of *spin* orbitals. If an active space is defined,  this is
            the number of active spin-orbitals.
        mapping (str): Specifies the transformation to map the fermionic operator to the
            Pauli basis. Input values can be ``'jordan_wigner'`` or ``'bravyi_kitaev'``.
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator
            to an observable measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted.
            If None, will use identity map (e.g. 0->0, 1->1, ...).

    Returns:
        pennylane.Hamiltonian: the total spin observable :math:`\hat{S}^2`

    **Example**

    >>> electrons = 2
    >>> orbitals = 4
    >>> S2 = spin2(electrons, orbitals, mapping="jordan_wigner")
    >>> print(S2)
    (0.75) [I0]
    + (0.375) [Z1]
    + (-0.375) [Z0 Z1]
    + (0.125) [Z0 Z2]
    + (0.375) [Z0]
    + (-0.125) [Z0 Z3]
    + (-0.125) [Z1 Z2]
    + (0.125) [Z1 Z3]
    + (0.375) [Z2]
    + (0.375) [Z3]
    + (-0.375) [Z2 Z3]
    + (0.125) [Y0 X1 Y2 X3]
    + (0.125) [Y0 Y1 X2 X3]
    + (0.125) [Y0 Y1 Y2 Y3]
    + (-0.125) [Y0 X1 X2 Y3]
    + (-0.125) [X0 Y1 Y2 X3]
    + (0.125) [X0 X1 X2 X3]
    + (0.125) [X0 X1 Y2 Y3]
    + (0.125) [X0 Y1 X2 Y3]

    >>> S2 = spin2(electrons, orbitals, mapping="jordan_wigner", wires=['w0','w1','w2','w3'])
    >>> print(S2)
    (0.75) [Iw0]
    + (0.375) [Zw1]
    + (-0.375) [Zw0 Zw1]
    + (0.125) [Zw0 Zw2]
    + (0.375) [Zw0]
    + (-0.125) [Zw0 Zw3]
    + (-0.125) [Zw1 Zw2]
    + (0.125) [Zw1 Zw3]
    + (0.375) [Zw2]
    + (0.375) [Zw3]
    + (-0.375) [Zw2 Zw3]
    + (0.125) [Yw0 Xw1 Yw2 Xw3]
    + (0.125) [Yw0 Yw1 Xw2 Xw3]
    + (0.125) [Yw0 Yw1 Yw2 Yw3]
    + (-0.125) [Yw0 Xw1 Xw2 Yw3]
    + (-0.125) [Xw0 Yw1 Yw2 Xw3]
    + (0.125) [Xw0 Xw1 Xw2 Xw3]
    + (0.125) [Xw0 Xw1 Yw2 Yw3]
    + (0.125) [Xw0 Yw1 Xw2 Yw3]
    """

    if electrons <= 0:
        raise ValueError(
            "'electrons' must be greater than 0; got for 'electrons' {}".format(electrons)
        )

    if orbitals <= 0:
        raise ValueError(
            "'orbitals' must be greater than 0; got for 'orbitals' {}".format(orbitals)
        )

    sz = np.where(np.arange(orbitals) % 2 == 0, 0.5, -0.5)

    table = _spin2_matrix_elements(sz)
    # To be consistent with the 'observable' function
    table[:, 4] = 2 * table[:, 4]

    return observable([table], init_term=3 / 4 * electrons, mapping=mapping, wires=wires)

def observable(matrix_elements, init_term=0, mapping="jordan_wigner", wires=None):

    r"""Builds the many-body observable whose expectation value can be
    measured in PennyLane.

    This function can be used to build second-quantized operators in the basis
    of single-particle states (e.g., HF states) and to transform them into
    PennyLane observables. In general, the many-body observable :math:`\hat{O}` can combine
    one-particle and two-particle operators, as in the case for electronic Hamiltonians:

    .. math::

        \hat{O} = \sum_{\alpha, \beta} \langle \alpha \vert \hat{t}^{(1)} +
        \cdots + \hat{t}^{(n)} \vert \beta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta
        + \frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta}
        \langle \alpha, \beta \vert \hat{v}^{(1)} + \cdots + \hat{v}^{(n)}
        \vert \gamma, \delta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta^\dagger
        \hat{c}_\gamma \hat{c}_\delta

    In the latter equations the indices :math:`\alpha, \beta, \gamma, \delta` run over the
    basis of single-particle states. The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}`
    are the particle creation and annihilation operators, respectively.
    :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` denotes the matrix element of
    the single-particle operator :math:`\hat{t}` entering the observable. For example,
    in electronic structure calculations, this is the case for: the kinetic energy operator,
    the nuclei Coulomb potential, or any other external fields included in the model Hamiltonian.
    On the other hand, :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle`
    denotes the matrix element of the two-particle operator :math:`\hat{v}`, for example, the
    Coulomb interaction between the electrons.

    If an `active space <https://en.wikipedia.org/wiki/Complete_active_space>`_
    (see :func:`~.active_space`) is defined, the observable is expanded over the truncated
    basis of active orbitals. The contribution of core orbitals can be passed to the
    function using the keyword argument ``init_term``.

    The function uses tools of `OpenFermion <https://github.com/quantumlib/OpenFermion>`_
    to build the second-quantized operator and map it to basis of Pauli matrices via the
    Jordan-Wigner or Bravyi-Kitaev transformation. Finally, the qubit operator is
    converted to a a PennyLane observable by the function :func:`~.convert_observable`.

    Args:
        matrix_elements (list(array[float])): list containing 2D numpy arrays with the matrix
            elements of the operators :math:`\hat{t}` and :math:`\hat{v}`.
            For single-particle operators the :math:`i`-th array in the list will have shape
            ``(matrix_elements[i].shape[0], 3)`` with each row containing the indices
            :math:`\alpha`, :math:`\beta` and the matrix element
            :math:`\langle \alpha \vert \hat{t}^{(i)}\vert \beta \rangle`.
            For two-particle operators the :math:`j`-th array in the list
            will have shape ``(matrix_elements[j].shape[0], 5)`` with each row containing
            the indices :math:`\alpha`, :math:`\beta`, :math:`\gamma`, :math:`\delta` and
            the matrix element
            :math:`\langle \alpha, \beta \vert \hat{v}^{(j)}\vert \gamma, \delta \rangle`.
        init_term (float): the contribution of core orbitals, if any, or other quantity
            required to initialize the many-body observable.
        mapping (str): specifies the fermion-to-qubit mapping. Input values can
            be ``'jordan_wigner'`` or ``'bravyi_kitaev'``.
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator
            to an observable measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted.
            If None, will use identity map (e.g. 0->0, 1->1, ...).

    Returns:
        pennylane.Hamiltonian: the fermionic-to-qubit transformed observable

    **Example**

    >>> t = np.array([[0., 0., 0.5], [1.0, 1.0, -0.5], [1.0, 0., 0.]])
    >>> v = np.array([[ 0., 0., 0., 0., 0.25], [ 0., 1., 1., 0., -0.25], [ 1., 0., 0., 1., -0.5]])
    >>> matrix_elements = [t, v]
    >>> print(observable(matrix_elements, init_term=1/4, mapping="bravyi_kitaev"))
    (0.15625) [I0]
    + (-0.15625) [Z0]
    + (0.34375) [Z0 Z1]
    + (-0.09375) [Z1]
    """

    if mapping.strip().lower() not in ("jordan_wigner", "bravyi_kitaev"):
        raise TypeError(
            "The '{}' transformation is not available. \n "
            "Please set 'mapping' to 'jordan_wigner' or 'bravyi_kitaev'.".format(mapping)
        )

    # Initialize the FermionOperator
    mb_obs = FermionOperator() + FermionOperator("") * init_term

    for table in matrix_elements:

        if len(table.shape) != 2:
            raise ValueError(
                "Expected dimension for arrays in 'matrix_elements' is 2; got {}".format(
                    table.shape
                )
            )

        if table.shape[1] not in (3, 5):
            raise ValueError(
                "Expected entries of matrix element tables to be of shape (3,) or (5,); got {}".format(
                    table.shape[1]
                )
            )

        if table.shape[1] == 5:
            # two-particle operator
            for i in table:
                mb_obs += FermionOperator(
                    ((int(i[0]), 1), (int(i[1]), 1), (int(i[2]), 0), (int(i[3]), 0)), i[4] / 2
                )
        else:
            # single-particle operator
            for i in table:
                mb_obs += FermionOperator(((int(i[0]), 1), (int(i[1]), 0)), i[2])

    # Map the fermionic operator to a qubit operator
    if mapping.strip().lower() == "bravyi_kitaev":
        return structure.convert_observable(bravyi_kitaev(mb_obs), wires=wires)

    return structure.convert_observable(jordan_wigner(mb_obs), wires=wires)

def spin_z(orbitals, mapping="jordan_wigner", wires=None):
    r"""Computes the total spin projection operator :math:`\hat{S}_z` in the Pauli basis.

    The total spin projection operator :math:`\hat{S}_z` is given by

    .. math::

        \hat{S}_z = \sum_{\alpha, \beta} \langle \alpha \vert \hat{s}_z \vert \beta \rangle
        ~ \hat{c}_\alpha^\dagger \hat{c}_\beta, ~~ \langle \alpha \vert \hat{s}_z
        \vert \beta \rangle = s_{z_\alpha} \delta_{\alpha,\beta},

    where :math:`s_{z_\alpha} = \pm 1/2` is the spin-projection of the single-particle state
    :math:`\vert \alpha \rangle`. The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}`
    are the particle creation and annihilation operators, respectively.

    Args:
        orbitals (str): Number of *spin* orbitals. If an active space is defined, this is
            the number of active spin-orbitals.
        mapping (str): Specifies the transformation to map the fermionic operator to the
            Pauli basis. Input values can be ``'jordan_wigner'`` or ``'bravyi_kitaev'``.
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator
            to an observable measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted.
            If None, will use identity map (e.g. 0->0, 1->1, ...).

    Returns:
        pennylane.Hamiltonian: the total spin projection observable :math:`\hat{S}_z`

    **Example**

    >>> orbitals = 4
    >>> Sz = spin_z(orbitals, mapping="jordan_wigner")
    >>> print(Sz)
    (-0.25) [Z0]
    + (0.25) [Z1]
    + (-0.25) [Z2]
    + (0.25) [Z3]
    """

    if orbitals <= 0:
        raise ValueError(
            "'orbitals' must be greater than 0; got for 'orbitals' {}".format(orbitals)
        )

    r = np.arange(orbitals)
    sz_orb = np.where(np.arange(orbitals) % 2 == 0, 0.5, -0.5)
    table = np.vstack([r, r, sz_orb]).T

    return observable([table], mapping=mapping, wires=wires)

def particle_number(orbitals, mapping="jordan_wigner", wires=None):
    r"""Computes the particle number operator :math:`\hat{N}=\sum_\alpha \hat{n}_\alpha`
    in the Pauli basis.

    The particle number operator is given by

    .. math::

        \hat{N} = \sum_\alpha \hat{c}_\alpha^\dagger \hat{c}_\alpha,

    where the index :math:`\alpha` runs over the basis of single-particle states
    :math:`\vert \alpha \rangle`, and the operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are
    the particle creation and annihilation operators, respectively.

    Args:
        orbitals (int): Number of *spin* orbitals. If an active space is defined, this is
            the number of active spin-orbitals.
        mapping (str): Specifies the transformation to map the fermionic operator to the
            Pauli basis. Input values can be ``'jordan_wigner'`` or ``'bravyi_kitaev'``.
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator
            to an observable measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted.
            If None, will use identity map (e.g. 0->0, 1->1, ...).

    Returns:
        pennylane.Hamiltonian: the fermionic-to-qubit transformed observable

    **Example**

    >>> orbitals = 4
    >>> N = particle_number(orbitals, mapping="jordan_wigner")
    >>> print(N)
    (2.0) [I0]
    + (-0.5) [Z0]
    + (-0.5) [Z1]
    + (-0.5) [Z2]
    + (-0.5) [Z3]
    >>> N = particle_number(orbitals, mapping="jordan_wigner", wires=['w0','w1','w2','w3'])
    >>> print(N)
    (2.0) [Iw0]
    + (-0.5) [Zw0]
    + (-0.5) [Zw1]
    + (-0.5) [Zw2]
    + (-0.5) [Zw3]
    """

    if orbitals <= 0:
        raise ValueError(
            "'orbitals' must be greater than 0; got for 'orbitals' {}".format(orbitals)
        )

    r = np.arange(orbitals)
    table = np.vstack([r, r, np.ones([orbitals])]).T

    return observable([table], mapping=mapping, wires=wires)

def one_particle(matrix_elements, core=None, active=None, cutoff=1.0e-12):
    r"""Generates the table of matrix elements of a given one-particle operator
    required to build many-body qubit observables.

    Second quantized one-particle operators are expanded in the basis of single-particle
    states as

    .. math::

        \hat{T} = \sum_{\alpha, \beta} \langle \alpha \vert \hat{t} \vert \beta \rangle
        [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} +
        \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}].

    In the equation above the indices :math:`\alpha, \beta` run over the basis of spatial
    orbitals :math:`\phi_\alpha(r)`. Since the operator :math:`\hat{t}` acts only on the
    spatial coordinates, the spin quantum numbers are indicated explicitly with the up/down arrows.
    The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation
    and annihilation operators, respectively, and
    :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` denotes the matrix elements of
    the operator :math:`\hat{t}`

    .. math::

        \langle \alpha \vert \hat{t} \vert \beta \rangle = \int dr ~ \phi_\alpha^*(r)
        \hat{t}(r) \phi_\beta(r).

    If an active space is defined (see :func:`~.active_space`), the summation indices
    run over the active orbitals and the contribution due to core orbitals is computed as
    :math:`T_\mathrm{core} = 2 \sum_{\alpha\in \mathrm{core}}
    \langle \alpha \vert \hat{t} \vert \alpha \rangle`.

    Args:
        matrix_elements (array[float]): 2D NumPy array with the matrix elements
            :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle`
        core (list): indices of core orbitals, i.e., the orbitals that are
            not correlated in the many-body wave function
        active (list): indices of active orbitals, i.e., the orbitals used to
            build the correlated many-body wave function
        cutoff (float): Cutoff value for including matrix elements. The
            matrix elements with absolute value less than ``cutoff`` are neglected.

    Returns:
        tuple: Table of indices and matrix elements of the one-particle operator
        and the contribution due to core orbitals. The returned table is a 2D Numpy
        array where each row contains three elements, the *spin*-orbital indices
        :math:`\alpha, \beta` and the matrix element
        :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle`.

    **Example**

    >>> matrix_elements = np.array(
    ...     [[-4.72739313e+00, -1.05499666e-01, -1.66961416e-01,
    ...       6.18014041e-16,  2.86964662e-16, -3.46772026e-02],
    ...      [-1.05499666e-01, -1.49264622e+00, 3.28928073e-02,
    ...       -2.20398308e-16,  1.93277291e-16, 5.27078882e-02],
    ...      [-1.66961416e-01,  3.28928073e-02, -1.12554473e+00,
    ...       -2.82912389e-17,  2.55224784e-16, -3.04455743e-02],
    ...      [ 6.18014041e-16, -2.20398308e-16, -2.82912389e-17,
    ...       -1.13579985e+00, -1.94289029e-16, -2.36158697e-16],
    ...      [ 2.86964662e-16,  1.93277291e-16,  2.55224784e-16,
    ...       -2.77555756e-16, -1.13579985e+00,  2.06665432e-16],
    ...      [-3.46772026e-02,  5.27078882e-02, -3.04455743e-02,
    ...       -2.36158697e-16,  2.06665432e-16, -9.50966595e-01]]
    ... )
    >>> table, t_core = one_particle(matrix_elements, core=[0], active=[1, 2])
    >>> print(table)
    [[ 0.          0.         -1.49264622]
     [ 1.          1.         -1.49264622]
     [ 0.          2.          0.03289281]
     [ 1.          3.          0.03289281]
     [ 2.          0.          0.03289281]
     [ 3.          1.          0.03289281]
     [ 2.          2.         -1.12554473]
     [ 3.          3.         -1.12554473]]
    >>> print(t_core)
    -9.45478626
    """

    orbitals = matrix_elements.shape[0]

    if matrix_elements.ndim != 2:
        raise ValueError(
            "'matrix_elements' must be a 2D array; got matrix_elements.ndim = {}".format(
                matrix_elements.ndim
            )
        )

    if core is None:
        t_core = 0
    else:
        if any([i > orbitals - 1 or i < 0 for i in core]):
            raise ValueError(
                "Indices of core orbitals must be between 0 and {}; got core = {}".format(
                    orbitals - 1, core
                )
            )

        # Compute contribution due to core orbitals
        t_core = 2 * np.sum(matrix_elements[np.array(core), np.array(core)])

    if active is None:
        if core is None:
            active = list(range(orbitals))
        else:
            active = [i for i in range(orbitals) if i not in core]

    if any([i > orbitals - 1 or i < 0 for i in active]):
        raise ValueError(
            "Indices of active orbitals must be between 0 and {}; got active = {}".format(
                orbitals - 1, active
            )
        )

    # Indices of the matrix elements with absolute values >= cutoff
    indices = np.nonzero(np.abs(matrix_elements) >= cutoff)

    # Single out the indices of active orbitals
    num_indices = len(indices[0])
    pairs = [
        [indices[0][i], indices[1][i]]
        for i in range(num_indices)
        if all(indices[j][i] in active for j in range(len(indices)))
    ]

    # Building the table of indices and matrix elements
    table = np.zeros((2 * len(pairs), 3))

    for i, pair in enumerate(pairs):
        alpha, beta = pair
        element = matrix_elements[alpha, beta]

        # spin-up term
        table[2 * i, 0] = 2 * active.index(alpha)
        table[2 * i, 1] = 2 * active.index(beta)
        table[2 * i, 2] = element

        # spin-down term
        table[2 * i + 1, 0] = 2 * active.index(alpha) + 1
        table[2 * i + 1, 1] = 2 * active.index(beta) + 1
        table[2 * i + 1, 2] = element

    return table, t_core

def two_particle(matrix_elements, core=None, active=None, cutoff=1.0e-12):
    r"""Generates the table of matrix elements of a given two-particle operator
    required to build many-body qubit observables.

    Second quantized two-particle operators are expanded in the basis of single-particle
    states as

    .. math::

        \hat{V} = \frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta}
        \langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle
        ~ &[& \hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow}^\dagger
        \hat{c}_{\gamma\uparrow} \hat{c}_{\delta\uparrow} + \hat{c}_{\alpha\uparrow}^\dagger
        \hat{c}_{\beta\downarrow}^\dagger \hat{c}_{\gamma\downarrow} \hat{c}_{\delta\uparrow} \\
        &+& \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\uparrow}^\dagger
        \hat{c}_{\gamma\uparrow} \hat{c}_{\delta\downarrow} + \hat{c}_{\alpha\downarrow}^\dagger
        \hat{c}_{\beta\downarrow}^\dagger \hat{c}_{\gamma\downarrow} \hat{c}_{\delta\downarrow}~].

    In the equation above the indices :math:`\alpha, \beta, \gamma, \delta` run over the basis
    of spatial orbitals :math:`\phi_\alpha(r)`. Since the operator :math:`v` acts only on the
    spatial coordinates the spin quantum numbers are indicated explicitly with the up/down arrows.
    The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation and
    annihilation operators, respectively, and
    :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` denotes the
    matrix elements of the operator :math:`\hat{v}`

    .. math::

        \langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle =
        \int dr_1 \int dr_2 ~ \phi_\alpha^*(r_1) \phi_\beta^*(r_2) ~\hat{v}(r_1, r_2)~
        \phi_\gamma(r_2) \phi_\delta(r_1).

    If an active space is defined (see :func:`~.active_space`), the summation indices
    run over the active orbitals and the contribution due to core orbitals, is computed as
    :math:`V_\mathrm{core} = \sum_{\alpha,\beta \in \mathrm{core}}
    [2 \langle \alpha, \beta \vert \hat{v} \vert \beta, \alpha \rangle
    - \langle \alpha, \beta \vert \hat{v} \vert \alpha, \beta \rangle ]`.

    Args:
        matrix_elements (array[float]): 4D NumPy array with the matrix elements
            :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle`
        core (list): indices of core orbitals, i.e., the orbitals that are
            not correlated in the many-body wave function
        active (list): indices of active orbitals, i.e., the orbitals used to
            build the correlated many-body wave function
        cutoff (float): Cutoff value for including matrix elements. The
            matrix elements with absolute value less than ``cutoff`` are neglected.

    Returns:
        tuple: Table with indices and matrix elements of the two-particle operator
        and the contribution due to core orbitals. The returned table is a 2D Numpy
        array where each row contains five elements, the *spin*-orbital indices
        :math:`\alpha, \beta, \gamma, \delta` and the matrix element
        :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle`.

    **Example**

    >>> matrix_elements = np.array([[[[ 6.82389533e-01, -1.45716772e-16],
    ...                               [-2.77555756e-17,  1.79000576e-01]],
    ...                              [[-2.77555756e-17,  1.79000576e-16],
    ...                               [ 6.70732778e-01, 0.00000000e+00]]],
    ...                             [[[-1.45716772e-16,  6.70732778e-16],
    ...                               [ 1.79000576e-01, -8.32667268e-17]],
    ...                              [[ 1.79000576e-16, -8.32667268e-17],
    ...                               [ 0.00000000e+00,  7.05105632e-01]]]])
    >>> table, v_core = two_particle(matrix_elements, core=[0], active=[1])
    >>> print(table)
    [[0.         0.         0.         0.         0.70510563]
     [0.         1.         1.         0.         0.70510563]
     [1.         0.         0.         1.         0.70510563]
     [1.         1.         1.         1.         0.70510563]]
    >>> print(v_core)
    0.682389533
    """

    orbitals = matrix_elements.shape[0]

    if matrix_elements.ndim != 4:
        raise ValueError(
            "'matrix_elements' must be a 4D array; got 'matrix_elements.ndim = ' {}".format(
                matrix_elements.ndim
            )
        )

    if core is None:
        v_core = 0
    else:
        if any([i > orbitals - 1 or i < 0 for i in core]):
            raise ValueError(
                "Indices of core orbitals must be between 0 and {}; got core = {}".format(
                    orbitals - 1, core
                )
            )

        # Compute the contribution of core orbitals
        v_core = sum(
            [
                2 * matrix_elements[alpha, beta, beta, alpha]
                - matrix_elements[alpha, beta, alpha, beta]
                for alpha in core
                for beta in core
            ]
        )

    if active is None:
        if core is None:
            active = list(range(orbitals))
        else:
            active = [i for i in range(orbitals) if i not in core]

    if any([i > orbitals - 1 or i < 0 for i in active]):
        raise ValueError(
            "Indices of active orbitals must be between 0 and {}; got active = {}".format(
                orbitals - 1, active
            )
        )

    # Indices of the matrix elements with absolute values >= cutoff
    indices = np.nonzero(np.abs(matrix_elements) >= cutoff)

    # Single out the indices of active orbitals
    num_indices = len(indices[0])
    quads = [
        [indices[0][i], indices[1][i], indices[2][i], indices[3][i]]
        for i in range(num_indices)
        if all(indices[j][i] in active for j in range(len(indices)))
    ]

    # Building the table of matrix elements
    table = np.zeros((4 * len(quads), 5))

    for i, quad in enumerate(quads):
        alpha, beta, gamma, delta = quad
        element = matrix_elements[alpha, beta, gamma, delta]

        # up-up-up-up term
        table[4 * i, 0] = 2 * active.index(alpha)
        table[4 * i, 1] = 2 * active.index(beta)
        table[4 * i, 2] = 2 * active.index(gamma)
        table[4 * i, 3] = 2 * active.index(delta)
        table[4 * i, 4] = element

        # up-down-down-up term
        table[4 * i + 1, 0] = 2 * active.index(alpha)
        table[4 * i + 1, 1] = 2 * active.index(beta) + 1
        table[4 * i + 1, 2] = 2 * active.index(gamma) + 1
        table[4 * i + 1, 3] = 2 * active.index(delta)
        table[4 * i + 1, 4] = element

        # down-up-up-down term
        table[4 * i + 2, 0] = 2 * active.index(alpha) + 1
        table[4 * i + 2, 1] = 2 * active.index(beta)
        table[4 * i + 2, 2] = 2 * active.index(gamma)
        table[4 * i + 2, 3] = 2 * active.index(delta) + 1
        table[4 * i + 2, 4] = element

        # down-down-down-down term
        table[4 * i + 3, 0] = 2 * active.index(alpha) + 1
        table[4 * i + 3, 1] = 2 * active.index(beta) + 1
        table[4 * i + 3, 2] = 2 * active.index(gamma) + 1
        table[4 * i + 3, 3] = 2 * active.index(delta) + 1
        table[4 * i + 3, 4] = element

    return table, v_core

__all__ = [
    "observable",
    "particle_number",
    "spin_z",
    "spin2",
    "one_particle",
    "two_particle",
    "_spin2_matrix_elements",
]

max_stars_repo_name: theRoughCode/pennylane

max_issues_repo_name: theRoughCode/pennylane

max_forks_repo_name: theRoughCode/pennylane


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

def read_meta_file(path):
    with open(path, 'r') as f:
        data = json.load(f)
    return data

def parse_day_night(data, DEBUG=False):
    if DEBUG:
        print(data['daytime'])
    if data is None:
        return True
    day = data['daytime']['day']
    night = data['daytime']['night']
    assert day != night
    return day

def read_tiff_image(total_path):
    """
    reads raw 16bit tiff images
    :param path: folder Path
    :param name: image name
    :return: dtype= unit16 numpy array containing raw images
    """
    #return cv2.cvtColor(cv2.imread(os.path.join(path, name), -1), cv2.__getattribute__("COLOR_BAYER_" + "GR2RGB"))
    return cv2.imread(total_path, -1)

def save_tiff_image(path, name, image):
    """
    :param path: destination folder
    :param name: image name .png, .tiff denotes saving operation
    :param image: input image which has to be saved
    :return: None
    """
    cv2.imwrite(os.path.join(path, name), image)

def read_image_intrinsic(path_dataset, name_camera_calib):

    with open(os.path.join(path_dataset, name_camera_calib), 'r') as f:
        data_camera = json.load(f)

    return data_camera

def check_image(image, name='image'):

    print('#### entered %s check loop ####'%name)
    print('max', np.max(image))
    print('min', np.min(image))
    print('shape', np.shape(image))
    print(image.dtype)

BGR2LAB=[[0.1804375, 0.3575761, 0.4124564]/np.sum([0.1804375, 0.3575761, 0.4124564]),
         [0.0721750, 0.7151522, 0.2126729]/np.sum([0.0721750, 0.7151522, 0.2126729]),
         [0.0193339, 0.1191920, 0.9503041]/np.sum([[0.0193339, 0.1191920, 0.9503041]])
]

a=np.asarray(BGR2LAB)

LAB2BGR=np.linalg.inv(BGR2LAB)

def proj_BGR2LAB(image):
    img=np.copy(image).astype(np.float64)
    img2=np.copy(img)
    for i in range(img.shape[0]):
        a=np.matmul(BGR2LAB,img[i,:,:].transpose([1,0])).transpose([1,0])
        img2[i, :,:]=a

    return np.clip(img2, 0, 2**16).astype(np.uint16)

def proj_LAB2BGR(image):
    img=np.copy(image).astype(np.float64)
    img2=np.copy(img)
    for i in range(img.shape[0]):
        a=np.matmul(LAB2BGR,img[i,:,:].transpose([1,0])).transpose([1,0])
        img2[i, :,:]=a

    return np.clip(img2, 0, 2 ** 16).astype(np.uint16)

def apply_clahe_16bit(image):
    image = (2**-8*image).astype(np.uint8)
    clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
    image = cv2.cvtColor(image, cv2.COLOR_BGR2LAB)
    image[:, :, 0] = clahe.apply(image[:,:,0])
    clahe2 = cv2.createCLAHE(clipLimit=0.5, tileGridSize=(8, 8))
    image[:, :, 1] = clahe2.apply(image[:,:,1])
    image[:, :, 2] = clahe2.apply(image[:,:,2])
    image = cv2.cvtColor(image, cv2.COLOR_LAB2BGR)

    return image

def apply_clahe_8bit(image):
    clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
    image = cv2.cvtColor(image, cv2.COLOR_BGR2LAB)
    image[:, :, 0] = clahe.apply(image[:,:,0])
    clahe2 = cv2.createCLAHE(clipLimit=0.5, tileGridSize=(8, 8))
    image[:, :, 1] = clahe2.apply(image[:,:,1])
    image[:, :, 2] = clahe2.apply(image[:,:,2])
    image = cv2.cvtColor(image, cv2.COLOR_LAB2BGR)

    return image

ext: py

max_stars_repo_name: ladt/SeeingThroughFog

max_issues_repo_name: ladt/SeeingThroughFog

max_forks_repo_name: ladt/SeeingThroughFog

reason: import numpy


text:
import asyncio
from datetime import datetime, timezone
import logging
import textwrap
from typing import Any, Dict, Iterable, List, Optional, Sequence, Set, Tuple, Type, Union

import numpy as np
import pandas as pd
from pandas.core.dtypes.cast import OutOfBoundsDatetime, tslib
from pandas.core.internals import make_block
from pandas.core.internals.managers import BlockManager, form_blocks
import sentry_sdk
from sqlalchemy import Boolean, Column, DateTime, Integer
from sqlalchemy.orm.attributes import InstrumentedAttribute
from sqlalchemy.sql import ClauseElement
from sqlalchemy.sql.elements import Label

from athenian.api import metadata
from athenian.api.db import Database, DatabaseLike
from athenian.api.models.metadata.github import Base as MetadataBase
from athenian.api.models.persistentdata.models import Base as PerdataBase
from athenian.api.models.precomputed.models import GitHubBase as PrecomputedBase
from athenian.api.models.state.models import Base as StateBase
from athenian.api.to_object_arrays import to_object_arrays_split
from athenian.api.tracing import MAX_SENTRY_STRING_LENGTH

async def read_sql_query(sql: ClauseElement,
                         con: DatabaseLike,
                         columns: Union[Sequence[str], Sequence[InstrumentedAttribute],
                                        MetadataBase, PerdataBase, PrecomputedBase, StateBase],
                         index: Optional[Union[str, Sequence[str]]] = None,
                         soft_limit: Optional[int] = None,
                         ) -> pd.DataFrame:
    """Read SQL query into a DataFrame.

    Returns a DataFrame corresponding to the result set of the query.
    Optionally provide an `index_col` parameter to use one of the
    columns as the index, otherwise default integer index will be used.

    Parameters
    ----------
    sql     : SQLAlchemy query object to be executed.
    con     : async SQLAlchemy database engine.
    columns : list of the resulting columns names, column objects or the model if SELECT *
    index   : Name(s) of the index column(s).
    soft_limit
            : Load this number of rows at maximum.

    Returns
    -------
    DataFrame

    Notes
    -----
    Any datetime values with time zone information parsed via the `parse_dates`
    parameter will be converted to UTC.
    """
    try:
        data = await con.fetch_all(query=sql)
    except Exception as e:
        try:
            sql = str(sql)
        except Exception:
            sql = repr(sql)
        sql = textwrap.shorten(sql, MAX_SENTRY_STRING_LENGTH - 500)
        logging.getLogger("%s.read_sql_query" % metadata.__package__).error(
            "%s: %s; %s", type(e).__name__, e, sql)
        raise e from None
    if soft_limit is not None and len(data) > soft_limit:
        data = data[:soft_limit]
    return wrap_sql_query(data, columns, index)

def _create_block_manager_from_arrays(
    arrays_typed: List[np.ndarray],
    arrays_obj: np.ndarray,
    names_typed: List[str],
    names_obj: List[str],
    size: int,
) -> BlockManager:
    assert len(arrays_typed) == len(names_typed)
    assert len(arrays_obj) == len(names_obj)
    range_index = pd.RangeIndex(stop=size)
    typed_index = pd.Index(names_typed)
    blocks = form_blocks(arrays_typed, typed_index, [typed_index, range_index])
    blocks.append(make_block(arrays_obj, placement=np.arange(len(arrays_obj)) + len(arrays_typed)))
    return BlockManager(blocks, [pd.Index(names_typed + names_obj), range_index])

def wrap_sql_query(data: List[Sequence[Any]],
                   columns: Union[Sequence[str], Sequence[InstrumentedAttribute],
                                  MetadataBase, StateBase],
                   index: Optional[Union[str, Sequence[str]]] = None,
                   ) -> pd.DataFrame:
    """Turn the fetched DB records to a pandas DataFrame."""
    try:
        columns[0]
    except TypeError:
        dt_columns = _extract_datetime_columns(columns.__table__.columns)
        int_columns = _extract_integer_columns(columns.__table__.columns)
        bool_columns = _extract_boolean_columns(columns.__table__.columns)
        columns = [c.name for c in columns.__table__.columns]
    else:
        dt_columns = _extract_datetime_columns(columns)
        int_columns = _extract_integer_columns(columns)
        bool_columns = _extract_boolean_columns(columns)
        columns = [(c.name if not isinstance(c, str) else c) for c in columns]
    typed_cols_indexes = []
    typed_cols_names = []
    obj_cols_indexes = []
    obj_cols_names = []
    for i, column in enumerate(columns):
        if column in dt_columns or column in int_columns or column in bool_columns:
            cols_indexes = typed_cols_indexes
            cols_names = typed_cols_names
        else:
            cols_indexes = obj_cols_indexes
            cols_names = obj_cols_names
        cols_indexes.append(i)
        cols_names.append(column)
    log = logging.getLogger(f"{metadata.__package__}.wrap_sql_query")
    # we used to have pd.DataFrame.from_records + bunch of convert_*() in relevant columns
    # the current approach is faster for several reasons:
    # 1. avoid an expensive copy of the object dtype columns in the BlockManager construction
    # 2. call tslib.array_to_datetime directly without surrounding Pandas bloat
    # 3. convert to int in the numpy domain and thus do not have to mess with indexes
    #
    # an ideal conversion would be loading columns directly from asyncpg but that requires
    # quite some changes in their internals
    with sentry_sdk.start_span(op="wrap_sql_query/convert", description=str(size := len(data))):
        data_typed, data_obj = to_object_arrays_split(data, typed_cols_indexes, obj_cols_indexes)
        converted_typed = []
        discard_mask = None
        for column, values in zip(typed_cols_names, data_typed):
            if column in dt_columns:
                converted_typed.append(_convert_datetime(values))
            elif column in int_columns:
                values, discarded = _convert_integer(values, column, int_columns[column], log)
                converted_typed.append(values)
                if discarded is not None:
                    if discard_mask is None:
                        discard_mask = np.zeros(len(data), dtype=bool)
                    discard_mask[discarded] = True
            elif column in bool_columns:
                converted_typed.append(values.astype(bool))
            else:
                raise AssertionError("impossible: typed columns are either dt or int")
        if discard_mask is not None:
            left = ~discard_mask
            size = left.sum()
            converted_typed = [arr[left] for arr in converted_typed]
            data_obj = data_obj[:, left]
    with sentry_sdk.start_span(op="wrap_sql_query/pd.DataFrame()", description=str(size)):
        block_mgr = _create_block_manager_from_arrays(
            converted_typed, data_obj, typed_cols_names, obj_cols_names, size)
        frame = pd.DataFrame(block_mgr, columns=typed_cols_names + obj_cols_names, copy=False)
        for column in dt_columns:
            try:
                frame[column] = frame[column].dt.tz_localize(timezone.utc)
            except (AttributeError, TypeError):
                continue
        if index is not None:
            frame.set_index(index, inplace=True)
    return frame

def _extract_datetime_columns(columns: Iterable[Union[Column, str]]) -> Set[str]:
    return {
        c.name for c in columns
        if not isinstance(c, str) and (
            isinstance(c.type, DateTime) or
            (isinstance(c.type, type) and issubclass(c.type, DateTime))
        )
    }

def _extract_boolean_columns(columns: Iterable[Union[Column, str]]) -> Set[str]:
    return {
        c.name for c in columns
        if not isinstance(c, str) and (
            isinstance(c.type, Boolean) or
            (isinstance(c.type, type) and issubclass(c.type, Boolean))
        )
    }

def _extract_integer_columns(columns: Iterable[Union[Column, str]],
                             ) -> Dict[str, bool]:
    return {
        c.name: getattr(
            c, "info", {} if not isinstance(c, Label) else getattr(c.element, "info", {}),
        ).get("erase_nulls", False)
        for c in columns
        if not isinstance(c, str) and (
            isinstance(c.type, Integer) or
            (isinstance(c.type, type) and issubclass(c.type, Integer))
        )
        and not getattr(c, "nullable", False)
        and (not isinstance(c, Label) or (
            (not getattr(c.element, "nullable", False))
            and (not getattr(c, "nullable", False))
        ))
    }

def _convert_datetime(arr: np.ndarray) -> np.ndarray:
    # None converts to NaT
    try:
        ts, offset = tslib.array_to_datetime(arr, utc=True, errors="raise")
        assert offset is None
    except OutOfBoundsDatetime:
        # TODO(vmarkovtsev): copy the function and set OOB values to NaT
        # this comparison is very slow but still faster than removing tzinfo and taking np.array()
        arr[arr == datetime(1, 1, 1)] = None
        arr[arr == datetime(1, 1, 1, tzinfo=timezone.utc)] = None
        try:
            return _convert_datetime(arr)
        except OutOfBoundsDatetime as e:
            raise e from None
    # 0 converts to 1970-01-01T00:00:00
    ts[ts == np.zeros(1, ts.dtype)[0]] = None
    return ts

def postprocess_datetime(frame: pd.DataFrame,
                         columns: Optional[Iterable[str]] = None,
                         ) -> pd.DataFrame:
    """Ensure *inplace* that all the timestamps inside the dataframe are valid UTC or NaT.

    :return: Fixed dataframe - the same instance as `frame`.
    """
    utc_dt1 = datetime(1, 1, 1, tzinfo=timezone.utc)
    dt1 = datetime(1, 1, 1)
    if columns is not None:
        obj_cols = dt_cols = columns
    else:
        obj_cols = frame.select_dtypes(include=[object])
        dt_cols = frame.select_dtypes(include=["datetime"])
    for col in obj_cols:
        fc = frame[col]
        if utc_dt1 in fc:
            fc.replace(utc_dt1, pd.NaT, inplace=True)
        if dt1 in fc:
            fc.replace(dt1, pd.NaT, inplace=True)
    for col in dt_cols:
        fc = frame[col]
        if 0 in fc:
            fc.replace(0, pd.NaT, inplace=True)
        try:
            frame[col] = fc.dt.tz_localize(timezone.utc)
        except (AttributeError, TypeError):
            continue
    return frame

def _convert_integer(arr: np.ndarray,
                     name: str,
                     erase_null: bool,
                     log: logging.Logger,
                     ) -> Tuple[np.ndarray, Optional[np.ndarray]]:
    nulls = None
    while True:
        try:
            return arr.astype(int), nulls
        except TypeError as e:
            nulls = np.equal(arr, None)
            if not nulls.any() or not erase_null:
                raise ValueError(f"Column {name} is not all-integer") from e
            log.error("fetched nulls instead of integers in %s", name)
            arr[nulls] = 0

def postprocess_integer(frame: pd.DataFrame, columns: Iterable[Tuple[str, int]]) -> pd.DataFrame:
    """Ensure *inplace* that all the integers inside the dataframe are not objects.

    :return: Fixed dataframe, a potentially different instance.
    """
    dirty_index = False
    log = None
    for col, erase_null in columns:
        while True:
            try:
                frame[col] = frame[col].astype(int, copy=False)
                break
            except TypeError as e:
                nulls = frame[col].isnull().values
                if not nulls.any():
                    raise ValueError(f"Column {col} is not all-integer") from e
                if not erase_null:
                    raise ValueError(f"Column {col} is not all-integer\n"
                                     f"{frame.loc[nulls].to_dict('records')}") from e
                if log is None:
                    log = logging.getLogger(f"{metadata.__package__}.read_sql_query")
                log.error("fetched nulls instead of integers in %s: %s",
                          col, frame.loc[nulls].to_dict("records"))
                frame = frame.take(np.flatnonzero(~nulls))
                dirty_index = True
    if dirty_index:
        frame.reset_index(drop=True, inplace=True)
    return frame

async def gather(*coros_or_futures,
                 op: Optional[str] = None,
                 description: Optional[str] = None,
                 catch: Type[BaseException] = Exception,
                 ) -> Tuple[Any, ...]:
    """Return a future aggregating results/exceptions from the given coroutines/futures.

    This is equivalent to `asyncio.gather(*coros_or_futures, return_exceptions=True)` with
    subsequent exception forwarding.

    :param op: Wrap the execution in a Sentry span with this `op`.
    :param description: Sentry span description.
    :param catch: Forward exceptions of this type.
    """
    async def body():
        if len(coros_or_futures) == 0:
            return tuple()
        if len(coros_or_futures) == 1:
            return (await coros_or_futures[0],)
        results = await asyncio.gather(*coros_or_futures, return_exceptions=True)
        for r in results:
            if isinstance(r, catch):
                raise r from None
        return results

    if op is not None:
        with sentry_sdk.start_span(op=op, description=description):
            return await body()
    return await body()

async def read_sql_query_with_join_collapse(
        query: ClauseElement,
        db: Database,
        columns: Union[Sequence[str], Sequence[InstrumentedAttribute],
                       MetadataBase, PerdataBase, PrecomputedBase, StateBase],
        index: Optional[Union[str, Sequence[str]]] = None,
        soft_limit: Optional[int] = None,
) -> pd.DataFrame:
    """Enforce the predefined JOIN order in read_sql_query()."""
    query = query.with_statement_hint("Set(join_collapse_limit 1)")
    return await read_sql_query(query, db, columns=columns, index=index, soft_limit=soft_limit)

# Allow other coroutines to execute every Nth iteration in long loops
COROUTINE_YIELD_EVERY_ITER = 250

async def list_with_yield(iterable: Iterable[Any], sentry_op: str) -> List[Any]:
    """Drain an iterable to a list, tracing the loop in Sentry and respecting other coroutines."""
    with sentry_sdk.start_span(op=sentry_op) as span:
        things = []
        for i, thing in enumerate(iterable):
            if (i + 1) % COROUTINE_YIELD_EVERY_ITER == 0:
                await asyncio.sleep(0)
            things.append(thing)
        try:
            span.description = str(i)
        except UnboundLocalError:
            pass
    return things

max_stars_repo_name: athenianco/athenian-api-open

max_issues_repo_name: athenianco/athenian-api-open

max_forks_repo_name: athenianco/athenian-api-open


text:
# output.py module dealing with outputs of DEMENTpy.
# Bin Wang, January, 2020

import numpy as np
import pandas as pd

class Output():
    """
    This class deals with outputs.

    Accepts data derived from the initialization.py and grid.py modules, and have two methods:
        output():             stores all time series
        microbes_abundance(): a special method
        microbes_tradeoff():  another special method tracking intracellular carbon allocation
    """

    def __init__(self,runtime,data_init):
        """
        The constructor of Output class.

        Parameters:
            runtime:   user-specified parameters when running the model
            data_init: data dictionary;all inititialized data from the 'initialization.py' module
        Returns:
            Initialization:   all data initialized preceding the execution of grid.py: dictionary
            Microbial_traits: microbial traits only pulled out from Initialization: dataframe
            SubstratesSeries: substrate-specific total mass over the grid: Substrate * C
            Substrates_Sum:   total substrates over the grid: day * C
            MonomersSeries:   monomoer-specific total mass over the grid: Monomer * C
            Monomers_Sum:     total monomers over the grid: day * C
            MicrobesSeries:   taxon-specific total biomass over the grid: Taxon * (C,N,P)
            Microbes_Sum:     total biomass over the grid: day * (C,N,P)
            TransporterSeries: taxon-specific total transporter production over the grid
            EnzymesSeries:    enzyme-specific total mass over the grid: Enzyme * C
            Enzymes_Sum:      total enzyme summed up over the grid: list
            OsmolyteSeries:   taxon-specific total osmolyte production over the grid
            RespSeries:       total respiration over the grid
            CUE:
        Notes:
            Variables that are not tracked but instead by DEMENT(R version):
                NH4Series:
                PO4Series:
        """

        # A couple of vars used in processing outputs
        n_taxa    = int(runtime.loc['n_taxa',1])                  # number of taxa
        Mic_index = ["Tax" + str(i) for i in range(1,n_taxa + 1)] # microbial taxa index
        # number of time steps in each pulse
        self.cycle = int(runtime.loc['end_time',1])
        # Pass all runtime parameters to Runtime
        self.Runtime = runtime
        # Pass all initialized data (a dictionary) to 'Initialization'
        self.Initialization = data_init

        # Pull out microbial traits data and put them in a dataframe:Microbial_traits
        data = np.concatenate(
            (data_init['fb'][0:n_taxa],
            data_init['Microbes_pp']['C'][0:n_taxa],
            data_init['Microbes_pp']['N'][0:n_taxa],
            data_init['Microbes_pp']['P'][0:n_taxa],
            data_init['UptakeGenes'].sum(axis=1)[0:n_taxa],
            data_init['EnzGenes'].sum(axis=1)[0:n_taxa],
            data_init['OsmoGenes'].sum(axis=1)[0:n_taxa],
            data_init['UptakeGenes_trait'][0:n_taxa],
            data_init['EnzProdConsti_trait'][0:n_taxa],
            data_init['EnzProdInduci_trait'][0:n_taxa],
            data_init['OsmoProdConsti_trait'][0:n_taxa],
            data_init['OsmoProdInduci_trait'][0:n_taxa],
            data_init['TaxDroughtTol'][0:n_taxa]),
            axis=0
        )
        columns = [
            'F/B','C','N','P','Uptake_Gene','Enz_Gene','Osmo_Gene','Uptake_Cost','Enz_Consti_Cost',
            'Enz_Induci_Cost','Osmo_Consti_Cost','Osmo_Induci_Cost','Drought_tolerance'
        ]
        self.Microbial_traits = pd.DataFrame(data=data.reshape(13,n_taxa).T, index=Mic_index, columns=columns, dtype='float32')

        # Account for inputs in mass balance
        # self.Cum_Substrate      = (data_init['Substrates'].groupby(level=0,sort=False).sum()).sum(axis=0)
        # self.Cum_Monomer        = (data_init['Monomers'].groupby(level=0,sort=False).sum()).sum(axis=0)
        # self.Cum_Monomer_ratios = ecosystem.Cum_Monomer_ratios

        # Degradation rates
        #DecayRates_list       = [0] * n_substrates * gridsize
        #DecayRates_array      = np.array(DecayRates_list).reshape(n_substrates * gridsize,1)
        #DecayRates_df         = pd.DataFrame(data = DecayRates_array, index= data_init['Substrates'].index, columns= ['C'])
        #self.DecayRatesSeries = DecayRates_df.groupby(level=0,sort=False).sum()

        # Leaching
        #self.Cum_Leaching_N   = pd.Series([0],index=[0])
        #self.Cum_Leaching_P   = pd.Series([0],index=[0])

        # Substrates
        Substrates_grid           = data_init['Substrates'].groupby(level=0,sort=False).sum()
        Substrates_grid['C'].name = 0  # set the series name to 0
        self.SubstratesSeries     = Substrates_grid['C']
        #self.Substrates_Sum      = pd.Series([Substrates_grid['C'].sum()],index=[0])

        # Monomers
        Monomers_grid           = data_init['Monomers'].groupby(level=0,sort=False).sum()
        Monomers_grid['C'].name = 0
        self.MonomersSeries     = Monomers_grid['C']
        #self.Monomers_Sum      = pd.Series([sum(Monomers_grid["C"])],index=[0])
        #self.NH4Series         = pd.Series([Monomers_grid.loc["NH4","N"]],index=[0])
        #self.PO4Series         = pd.Series([Monomers_grid.loc["PO4","P"]],index=[0])

        # Microbes
        Microbes_grid             = data_init['Microbes'].groupby(level=0,sort=False).sum()
        Microbes_grid['C'].name   = 0                   # rename the series name to 0
        self.MicrobesSeries       = Microbes_grid['C']  # taxon-specific total biomass summed over the grid
        self.MicrobesSeries_repop = Microbes_grid['C']  # created for reinitialization
        # define a new var placing biomass (in C) by taxon grid-specifically
        self.microbes_grid_taxa   = data_init['Microbes']['C']

        # Number of individuals of Taxon count
        Taxon_index            = data_init['Microbes']['C'] > 0
        Taxon_index.name       = 0
        self.Taxon_count       = Taxon_index.groupby(level=0,sort=False).sum().astype('uint32')
        #self.Taxon_count_repop = Taxon_index.groupby(level=0,sort=False).sum().astype('uint32')

        # Taxon-specific CUE
        #self.CUE_TaxonSeries = pd.Series(data=[0]*n_taxa,index=Mic_index)

        # Enzymes
        Enzymes_grid       = data_init['Enzymes'].groupby(level=0,sort=False).sum() # total of each enzyme summed over the spatial grid
        Enzymes_grid.name  = 0
        self.EnzymesSeries = Enzymes_grid
        #self.Enzymes_Sum  = pd.Series([Enzymes_grid.sum()],index=[0])

        # Emergent properties over the grid
        self.RespSeries = pd.Series([0],index=[0], dtype='float32')  # respiration
        self.CUE_system = pd.Series([0],index=[0], dtype='float32')  # emergent CUE
        self.Kill       = pd.Series([0],index=[0], dtype='uint32')   # stochastic death toll

        # Transporters: taxon-specific production summed over the grid by taxon
        #self.TransporterSeries = pd.Series(data=[0]*n_taxa,index=Mic_index)

        # Enyzmes: taxon-specific production summed over the grid by taxon
        #self.EnzymeConSeries    = pd.Series(data=[0]*n_taxa,index=Mic_index)
        #self.EnzymeIndSeries    = pd.Series(data=[0]*n_taxa,index=Mic_index)
        self.Enzyme_TaxonSeries = pd.Series(data=[0]*n_taxa, index=Mic_index, dtype='float32')

        # Osmolytes: taxon-specific production summed over the grid by taxon
        #self.OsmolyteConSeries    = pd.Series(data=[0]*n_taxa,index=Mic_index)
        #self.OsmolyteIndSeries    = pd.Series(data=[0]*n_taxa,index=Mic_index)
        self.Osmolyte_TaxonSeries = pd.Series(data=[0]*n_taxa, index=Mic_index, dtype='float32')

        # Growth yield
        self.Growth_yield = pd.Series(data=[0]*n_taxa, index=Mic_index, dtype='float32')

    def output(self,ecosystem,year,day):
        """
        Records outputs in various variables of each iteration.

        Parameters:
            ecosystem: object from the grid.py module
            day:       the day to record to outputs
        Returns:
        """

        # Account for inputs in mass balance
        # self.Cum_Substrate = self.Cum_Substrate + ecosystem.SubInput.sum(axis = 0)
        # self.Cum_Monomer = self.Cum_Monomer + (ecosystem.MonInput.mul(ecosystem.Cum_Monomer_ratios,axis=0)).sum(axis=0)

        # DecayRates
        # DecayRates_grid = ecosystem.DecayRates.groupby(level=0,sort=False).sum()
        # self.DecayRatesSeries = pd.concat([self.DecayRatesSeries,DecayRates_grid],axis=1,sort=False)

        # Substrates
        Substrates_grid           = ecosystem.Substrates.groupby(level=0,sort=False).sum()
        Substrates_grid['C'].name = self.cycle*year + (day+1) # index the output continuously by day
        self.SubstratesSeries     = pd.concat([self.SubstratesSeries,Substrates_grid['C']], axis=1, sort=False)
        #self.Substrates_Sum   = pd.concat([self.Substrates_Sum,pd.Series([Substrates_grid['C'].sum()],index=[day+1])],axis=0,sort=False)

        # Monomers
        Monomers_grid           = ecosystem.Monomers.groupby(level=0,sort=False).sum()
        Monomers_grid['C'].name = self.cycle*year + (day+1)
        self.MonomersSeries     = pd.concat([self.MonomersSeries,Monomers_grid['C']], axis=1, sort=False)
        #self.Monomers_Sum   = pd.concat([self.Monomers_Sum, pd.Series([sum(Monomers_grid["C"])],  index=[day+1])], axis=0, sort=False)
        #self.NH4Series      = pd.concat([self.NH4Series, pd.Series([Monomers_grid.loc["NH4","N"]],index=[day+1])], axis=0, sort=False)
        #self.PO4Series      = pd.concat([self.PO4Series, pd.Series([Monomers_grid.loc["PO4","P"]],index=[day+1])], axis=0, sort=False)

        # Interim Microbes
        #Microbes_interim_grid = ecosystem.Microbes_interim.groupby(level=0,sort=False).sum()
        #self.Microbes_Interim = pd.concat([self.MicrobesSeries,Microbes_interim_grid['C']],axis=1,sort=False)
        # Count taxon for averaging taxon CUE
        #Taxon_index = (ecosystem.Microbes_interim)['C'] > 0
        #Taxon_index.name = day + 1
        #taxon_count = Taxon_index.groupby(level=0,sort=False).sum()
        #self.Taxon_count = pd.concat([self.Taxon_count,taxon_count],axis=1,sort=False)

        # Microbe
        ## Taxon abundance
        Taxon_index      = ecosystem.Microbes['C'] > 0
        Taxon_index.name = self.cycle*year + (day+1)
        self.Taxon_count = pd.concat([self.Taxon_count, Taxon_index.groupby(level=0,sort=False).sum().astype('uint32')], axis=1, sort=False)
        ## Taxon biomass
        Microbes_grid           = ecosystem.Microbes.groupby(level=0,sort=False).sum()
        Microbes_grid['C'].name = self.cycle*year + (day+1)
        self.MicrobesSeries     = pd.concat([self.MicrobesSeries, Microbes_grid['C']], axis=1, sort=False)
        #self.Microbes_Sum   = pd.concat([self.Microbes_Sum,pd.Series([Microbes_grid['C'].sum()],index=[day+1])],axis=0,sort=False)
        self.microbes_grid_taxa  = pd.concat([self.microbes_grid_taxa,ecosystem.Microbes['C']],axis=1,sort=False)

        # Microbe-Taxon-specific CUE
        #CUE_Taxon_grid = ecosystem.CUE_Taxon.groupby(level=0,sort=False).sum()
        #CUE_Taxon_grid = CUE_Taxon_grid/taxon_count
        #CUE_Taxon_grid.name = day + 1
        #self.CUE_TaxonSeries = pd.concat([self.CUE_TaxonSeries,CUE_Taxon_grid],axis=1,sort=False)

        # Enzymes over the grid
        # Derive the enzyme-specific enzyme production summed over the grid by enzyme
        Enzymes_grid       = ecosystem.Enzymes.groupby(level=0,sort=False).sum()
        Enzymes_grid.name  = self.cycle*year + (day+1)
        self.EnzymesSeries = pd.concat([self.EnzymesSeries,Enzymes_grid], axis=1, sort=False)

        # Respiration
        self.RespSeries = pd.concat([self.RespSeries, pd.Series([ecosystem.Respiration],index=[day+1],dtype='float32')], axis=0, sort=False)
        # Carbon use efficiency
        self.CUE_system = pd.concat([self.CUE_system, pd.Series([ecosystem.CUE_system], index=[day+1],dtype='float32')], axis=0, sort=False)
        # Death toll of stochasticity origin
        self.Kill       = pd.concat([self.Kill, pd.Series([ecosystem.Kill],index=[day+1],dtype='uint32')], axis=0, sort=False)

    def microbes_abundance(self,ecosystem,year, day):
        """
        Seperately output Microbes from each time step and put them in a dataframe.

        Aims to eventually deal with reinitializing microbial community on the
        grid in a new pulse via calculating the cumulative frequency of different
        taxa over each cycle/pulse. The reason to have this seperate method instead of using the method
        above is b/c of the need to track outputs in every iteration, whereas the
        method above tracks outputs with a certain time interval (unless the interval set to 1)

        Parameters:
            ecosystem: an instance of the Grid object, in which only the Microbes is used.
            day:       the iteration index
        Returns:
            MicrobesSeries_repop: dataframe; tracking total biomass of differing taxa
            Taxon_count_repop:    dataframe; tracking abundances of differing taxa
        """

        # Track biomass of every taxon
        Microbes_grid             = ecosystem.Microbes.groupby(level=0,sort=False).sum()
        Microbes_grid['C'].name   = self.cycle*year + (day+1)
        self.MicrobesSeries_repop = pd.concat([self.MicrobesSeries_repop,Microbes_grid['C']],axis=1,sort=False)

        # Track abundance of every taxon
        #Taxon_index      = ecosystem.Microbes['C'] > 0
        #Taxon_index.name = day + 1
        #taxon_count      = Taxon_index.groupby(level=0, sort=False).sum().astype('uint32')
        #self.Taxon_count_repop = pd.concat([self.Taxon_count_repop,taxon_count], axis=1, sort=False)

    def microbes_tradeoff(self, ecosystem, year, day):
        """
        Track allocation of C uptake to different components.

        These components include: osmolytes, enzymes, and carbon gain
        """

        # Microbe-Transporters
        #Transporter_grid = ecosystem.Transporters.groupby(level=0,sort=False).sum() #taxon-specific transpotors summed over the grid by taxon
        #Transporter_grid.name = day + 1
        #self.TransporterSeries = pd.concat([self.TransporterSeries,Transporter_grid],axis=1,sort=False)

        # Microbe-Enzymes
        # Constitutive
        #Enzyme_Con_grid = ecosystem.Enzyme_Con.groupby(level=0,sort=False).sum() # Taxon-specific osmolytes summed over the grid by taxon
        #Enzyme_Con_grid.name = day + 1
        #self.EnzymeConSeries = pd.concat([self.EnzymeConSeries,Enzyme_Con_grid],axis=1,sort=False)
        # Inducible
        Enzyme_Ind_grid = ecosystem.Taxon_Enzyme_Cost_C.groupby(level=0,sort=False).sum()
        Enzyme_Ind_grid.name = self.cycle*year + (day+1)
        #self.EnzymeIndSeries = pd.concat([self.EnzymeIndSeries,Enzyme_Ind_grid],axis=1,sort=False)
        # Total
        self.Enzyme_TaxonSeries = pd.concat([self.Enzyme_TaxonSeries, Enzyme_Ind_grid],axis=1,sort=False)

        # Microbe-Osmolytes
        # Constitutive
        #Osmolyte_Con_grid = ecosystem.Osmolyte_Con.groupby(level=0,sort=False).sum()
        #Osmolyte_Con_grid.name = day + 1
        #self.OsmolyteConSeries = pd.concat([self.OsmolyteConSeries,Osmolyte_Con_grid],axis=1,sort=False)
        # Inducible
        Osmolyte_Ind_grid = ecosystem.Taxon_Osmo_Cost_C.groupby(level=0,sort=False).sum()
        Osmolyte_Ind_grid.name = self.cycle*year + (day+1)
        #self.OsmolyteIndSeries = pd.concat([self.OsmolyteIndSeries,Osmolyte_Ind_grid],axis=1,sort=False)
        # Total
        self.Osmolyte_TaxonSeries = pd.concat([self.Osmolyte_TaxonSeries, Osmolyte_Ind_grid],axis=1,sort=False)

        # Growth yield by Taxon
        GY_grid = ecosystem.Microbe_C_Gain.groupby(level=0,sort=False).sum()
        GY_grid.name = self.cycle*year + (day+1)
        self.Growth_yield = pd.concat([self.Growth_yield,GY_grid],axis=1,sort=False)

max_stars_repo_name: bioatmosphere/DEMENTpy

max_issues_repo_name: bioatmosphere/DEMENTpy

max_forks_repo_name: bioatmosphere/DEMENTpy


text:
from sympy import partition as p
from time import time

t = time()
i = 10
while str(p(i))[-6:] != "000000":
    i += 1
print(f"{p(i)}\n{i} \n{time() - t}")

reason: from sympy


text:
import numpy as np
from Jabc import Jabc, Jabc2
import matplotlib.pyplot as plt

omega = np.exp(np.pi * 1j * 2 / 3)
r = (np.sqrt(3)+1) / (np.sqrt(3)-1)

# 4 qubits
zs = [1, omega, omega**2, 0]
zs_A = [zs[0]]
zs_B = [zs[1]]
zs_C = [zs[2]]
zs_D = [zs[3]]

X = [z.real for z in zs]
Y = [z.imag for z in zs]
color= ["red", "green", "blue", "black"]
plt.scatter(X,Y, c=color)
plt.show()

print("4 qubit case. J(A,B,C)={}".format(Jabc(zs_A, zs_B, zs_C, zs_D, verbose=False)))
print("4 qubit case. (compressed) J(A,B,C)={}".format(Jabc2(zs_A, zs_B, zs_C, zs_D, verbose=False)))

# 6 qubits
zs = [1, omega, omega**2, r * np.sqrt(omega), r * np.sqrt(omega) * omega, r * np.sqrt(omega) * omega * omega]
zs_A = [zs[0]]
zs_B = [zs[1]]
zs_C = [zs[2]]
zs_D = zs[3:]

X = [z.real for z in zs]
Y = [z.imag for z in zs]
color= ["red", "green", "blue", "black", "black", "black"]
plt.scatter(X,Y, c=color)
plt.show()
print("6 qubit case. J(A,B,C)={}".format(Jabc(zs_A, zs_B, zs_C, zs_D, verbose=False)))
print("6 qubit case. (compressed) J(A,B,C)={}".format(Jabc2(zs_A, zs_B, zs_C, zs_D, verbose=False)))

# Fibonacci
N=21
zs = []
phi= (1+np.sqrt(5))/2
for i in range(N):
    t1 = i/N
    t2 = i/phi
    zs.append(np.sqrt(t1) * np.exp(2*np.pi * 1j * t2))

X = [z.real for z in zs]
Y = [z.imag for z in zs]
plt.scatter(X,Y)
plt.show()

max_stars_repo_name: ikim-quantum/Jabc

max_issues_repo_name: ikim-quantum/Jabc

max_forks_repo_name: ikim-quantum/Jabc


text:
"""
Copyright 2013 Steven Diamond

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.
"""
from typing import List, Tuple, Union

import numpy as np
import scipy.sparse as sp

from cvxpy.atoms.axis_atom import AxisAtom
from cvxpy.atoms.norm1 import norm1
from cvxpy.atoms.norm_inf import norm_inf
from cvxpy.constraints.constraint import Constraint
from cvxpy.utilities.power_tools import pow_high, pow_mid, pow_neg

def pnorm(x, p: Union[int, str] = 2, axis=None, keepdims: bool = False, max_denom: int = 1024):
    """Factory function for a mathematical p-norm.

    Parameters
    ----------
    p : numeric type or string
       The type of norm to construct; set this to np.inf or 'inf' to
       construct an infinity norm.

    Returns
    -------
    Atom
       A norm1, norm_inf, or Pnorm object.
    """
    if p == 1:
        return norm1(x, axis=axis, keepdims=keepdims)
    elif p in [np.inf, 'inf', 'Inf']:
        return norm_inf(x, axis=axis, keepdims=keepdims)
    else:
        return Pnorm(x, p=p, axis=axis, keepdims=keepdims, max_denom=max_denom)

class Pnorm(AxisAtom):
    r"""The vector p-norm, for p not equal to 1 or infinity.

    If given a matrix variable, ``pnorm`` will treat it as a vector, and
    compute the p-norm of the concatenated columns. Only accepts p values
    that are not equal to 1 or infinity; the norm1 and norm_inf classes
    handle those norms.

    For :math:`p > 1`, the p-norm is given by

    .. math::

        \|x\|_p = \left(\sum_i |x_i|^p \right)^{1/p},

    with domain :math:`x \in \mathbf{R}^n`.

    For :math:`p < 1,\ p \neq 0`, the p-norm is given by

    .. math::

        \|x\|_p = \left(\sum_i x_i^p \right)^{1/p},

    with domain :math:`x \in \mathbf{R}^n_+`.

    - Note that the "p-norm" is actually a **norm** only when
      :math:`p > 1`. For these cases, it is convex.
    - The expression is not defined when :math:`p = 0`.
    - Otherwise, when :math:`p < 1`, the expression is
      concave, but it is not a true norm.

    .. note::

        Generally, ``p`` cannot be represented exactly, so a rational,
        i.e., fractional, **approximation** must be made.

        Internally, ``pnorm`` computes a rational approximation
        to the reciprocal :math:`1/p` with a denominator up to ``max_denom``.
        The resulting
        approximation can be found through the attribute ``pnorm.p``.
        The approximation error is given by the attribute ``pnorm.approx_error``.
        Increasing ``max_denom`` can give better approximations.

        When ``p`` is an ``int`` or ``Fraction`` object, the approximation
        is usually **exact**.

    Parameters
    ----------
    x : cvxpy.Variable
        The value to take the norm of.

    p : int, float, or Fraction
        We require that :math:`p > 1`, but :math:`p \neq \infty`. See the
        norm1 and norm_inf classes for these norms, or use the pnorm
        function wrapper to instantiate them.

    max_denom : int
        The maximum denominator considered in forming a rational approximation
        for ``p``.

    axis : 0 or 1
           The axis to apply the norm to.

    Returns
    -------
    Expression
        An Expression representing the norm.
    """
    _allow_complex = True

    def __init__(self, x, p: int = 2, axis=None,
                 keepdims: bool = False, max_denom: int = 1024) -> None:
        if p < 0:
            # TODO(akshayka): Why do we accept p < 0?
            self.p, _ = pow_neg(p, max_denom)
        elif 0 < p < 1:
            self.p, _ = pow_mid(p, max_denom)
        elif p > 1:
            self.p, _ = pow_high(p, max_denom)
        elif p == 1:
            raise ValueError('Use the norm1 class to instantiate a one norm.')
        elif p == 'inf' or p == 'Inf' or p == np.inf:
            raise ValueError('Use the norm_inf class to instantiate an '
                             'infinity norm.')
        else:
            raise ValueError('Invalid p: {}'.format(p))
        self.approx_error = float(abs(self.p - p))
        self.original_p = p
        super(Pnorm, self).__init__(x, axis=axis, keepdims=keepdims)

    def numeric(self, values):
        """Returns the p-norm of x.
        """

        if self.axis is None:
            values = np.array(values[0]).flatten()
        else:
            values = np.array(values[0])

        if self.p < 1 and np.any(values < 0):
            return -np.inf
        if self.p < 0 and np.any(values == 0):
            return 0.0

        return np.linalg.norm(values, float(self.p), axis=self.axis,
                              keepdims=self.keepdims)

    def validate_arguments(self) -> None:
        super(Pnorm, self).validate_arguments()
        # TODO(akshayka): Why is axis not supported for other norms?
        if self.axis is not None and self.p != 2:
            raise ValueError(
                "The axis parameter is only supported for p=2.")
        if self.p < 1 and self.args[0].is_complex():
            raise ValueError("pnorm(x, p) cannot have x complex for p < 1.")

    def sign_from_args(self) -> Tuple[bool, bool]:
        """Returns sign (is positive, is negative) of the expression.
        """
        # Always positive.
        return (True, False)

    def is_atom_convex(self) -> bool:
        """Is the atom convex?
        """
        return self.p > 1

    def is_atom_concave(self) -> bool:
        """Is the atom concave?
        """
        return self.p < 1

    def is_atom_log_log_convex(self) -> bool:
        """Is the atom log-log convex?
        """
        return True

    def is_atom_log_log_concave(self) -> bool:
        """Is the atom log-log concave?
        """
        return False

    def is_incr(self, idx) -> bool:
        """Is the composition non-decreasing in argument idx?
        """
        return self.p < 1 or (self.p > 1 and self.args[0].is_nonneg())

    def is_decr(self, idx) -> bool:
        """Is the composition non-increasing in argument idx?
        """
        return self.p > 1 and self.args[0].is_nonpos()

    def is_pwl(self) -> bool:
        """Is the atom piecewise linear?
        """
        return False

    def get_data(self):
        return [self.p, self.axis]

    def name(self) -> str:
        return "%s(%s, %s)" % (self.__class__.__name__,
                               self.args[0].name(),
                               self.p)

    def _domain(self) -> List[Constraint]:
        """Returns constraints describing the domain of the node.
        """
        if self.p < 1 and self.p != 0:
            return [self.args[0] >= 0]
        else:
            return []

    def _grad(self, values):
        """Gives the (sub/super)gradient of the atom w.r.t. each argument.

        Matrix expressions are vectorized, so the gradient is a matrix.

        Args:
            values: A list of numeric values for the arguments.

        Returns:
            A list of SciPy CSC sparse matrices or None.
        """
        return self._axis_grad(values)

    def _column_grad(self, value):
        """Gives the (sub/super)gradient of the atom w.r.t. a column argument.

        Matrix expressions are vectorized, so the gradient is a matrix.

        Args:
            value: A numeric value for a column.

        Returns:
            A NumPy ndarray or None.
        """
        rows = value.size
        # Outside domain.
        if self.p < 1 and np.any(value <= 0):
            return None
        D_null = sp.csc_matrix((rows, 1), dtype='float64')
        denominator = np.linalg.norm(value, float(self.p))
        denominator = np.power(denominator, self.p - 1)
        # Subgrad is 0 when denom is 0 (or undefined).
        if denominator == 0:
            if self.p > 1:
                return D_null.todense()
            else:
                return None
        else:
            nominator = np.power(value, self.p - 1)
            frac = np.divide(nominator, denominator)
            return np.reshape(frac, (frac.size, 1))

max_stars_repo_name: QiuWJX/cvxpy

max_issues_repo_name: QiuWJX/cvxpy

max_forks_repo_name: QiuWJX/cvxpy

reason: import numpy,import scipy,from cvxpy


text:
import numpy as np
import pandas as pd

aminos_df = pd.read_csv("./aminos_table.csv", index_col=0)

def prot_to_vector(seq: str) -> np.ndarray:
    """Concatenate the amino acid features for each position of the sequence.
    Args:
        seq: A string representing an amino acid sequence.

    Returns:
        A numpy array of features, shape (len(seq), features)"""
    try:
        chain = [aminos_df.loc[residue].values for residue in seq]
    except KeyError:
        raise FeaturesError('Invalid residue encountered in fixed_prot')
    return np.concatenate(chain, axis=0).reshape(len(seq), -1)