text:
# This code is part of Qiskit.
#
# (C) Copyright IBM 2020.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

""" ListOp Operator Class """

from functools import reduce
from numbers import Number
from typing import Callable, Dict, Iterator, List, Optional, Set, Sequence, Union, cast

import numpy as np
from scipy.sparse import spmatrix

from qiskit.circuit import ParameterExpression, QuantumCircuit
from qiskit.opflow.exceptions import OpflowError
from qiskit.opflow.operator_base import OperatorBase
from qiskit.quantum_info import Statevector
from qiskit.utils import arithmetic

class ListOp(OperatorBase):
    """
    A Class for manipulating List Operators, and parent class to ``SummedOp``, ``ComposedOp``,
    and ``TensoredOp``.

    List Operators are classes for storing and manipulating lists of Operators, State functions,
    or Measurements, and include some rule or ``combo_fn`` defining how the Operator functions
    of the list constituents should be combined to form to cumulative Operator function of the
    ``ListOp``. For example, a ``SummedOp`` has an addition-based ``combo_fn``, so once the
    Operators in its list are evaluated against some bitstring to produce a list of results,
    we know to add up those results to produce the final result of the ``SummedOp``'s
    evaluation. In theory, this ``combo_fn`` can be any function over classical complex values,
    but for convenience we've chosen for them to be defined over NumPy arrays and values. This way,
    large numbers of evaluations, such as after calling ``to_matrix`` on the list constituents,
    can be efficiently combined. While the combination function is defined over classical
    values, it should be understood as the operation by which each Operators' underlying
    function is combined to form the underlying Operator function of the ``ListOp``. In this
    way, the ``ListOps`` are the basis for constructing large and sophisticated Operators,
    State Functions, and Measurements.

    The base ``ListOp`` class is particularly interesting, as its ``combo_fn`` is "the identity
    list Operation". Meaning, if we understand the ``combo_fn`` as a function from a list of
    complex values to some output, one such function is returning the list as-is. This is
    powerful for constructing compact hierarchical Operators which return many measurements in
    multiple dimensional lists.
    """

    def __init__(
        self,
        oplist: Sequence[OperatorBase],
        combo_fn: Callable = lambda x: x,
        coeff: Union[complex, ParameterExpression] = 1.0,
        abelian: bool = False,
        grad_combo_fn: Optional[Callable] = None,
    ) -> None:
        """
        Args:
            oplist: The list of ``OperatorBases`` defining this Operator's underlying function.
            combo_fn (callable): The recombination function to combine classical results of the
                ``oplist`` Operators' eval functions (e.g. sum).
            coeff: A coefficient multiplying the operator
            abelian: Indicates whether the Operators in ``oplist`` are known to mutually commute.
            grad_combo_fn: The gradient of recombination function. If None, the gradient will
                be computed automatically.
            Note that the default "recombination function" lambda above is essentially the
            identity - it accepts the list of values, and returns them in a list.
        """
        super().__init__()
        self._oplist = self._check_input_types(oplist)
        self._combo_fn = combo_fn
        self._coeff = coeff
        self._abelian = abelian
        self._grad_combo_fn = grad_combo_fn

    def _check_input_types(self, oplist):
        if all(isinstance(x, OperatorBase) for x in oplist):
            return list(oplist)
        else:
            badval = next(x for x in oplist if not isinstance(x, OperatorBase))
            raise TypeError(f"ListOp expecting objects of type OperatorBase, got {badval}")

    def _state(
        self,
        coeff: Optional[Union[complex, ParameterExpression]] = None,
        combo_fn: Optional[Callable] = None,
        abelian: Optional[bool] = None,
        grad_combo_fn: Optional[Callable] = None,
    ) -> Dict:
        return {
            "coeff": coeff if coeff is not None else self.coeff,
            "combo_fn": combo_fn if combo_fn is not None else self.combo_fn,
            "abelian": abelian if abelian is not None else self.abelian,
            "grad_combo_fn": grad_combo_fn if grad_combo_fn is not None else self.grad_combo_fn,
        }

    @property
    def oplist(self) -> List[OperatorBase]:
        """The list of ``OperatorBases`` defining the underlying function of this
        Operator.

        Returns:
            The Operators defining the ListOp
        """
        return self._oplist

    @property
    def combo_fn(self) -> Callable:
        """The function defining how to combine ``oplist`` (or Numbers, or NumPy arrays) to
        produce the Operator's underlying function. For example, SummedOp's combination function
        is to add all of the Operators in ``oplist``.

        Returns:
            The combination function.
        """
        return self._combo_fn

    @property
    def grad_combo_fn(self) -> Optional[Callable]:
        """The gradient of ``combo_fn``."""
        return self._grad_combo_fn

    @property
    def abelian(self) -> bool:
        """Whether the Operators in ``oplist`` are known to commute with one another.

        Returns:
            A bool indicating whether the ``oplist`` is Abelian.
        """
        return self._abelian

    @property
    def distributive(self) -> bool:
        """Indicates whether the ListOp or subclass is distributive under composition.
        ListOp and SummedOp are, meaning that (opv @ op) = (opv[0] @ op + opv[1] @ op)
        (using plus for SummedOp, list for ListOp, etc.), while ComposedOp and TensoredOp
        do not behave this way.

        Returns:
            A bool indicating whether the ListOp is distributive under composition.
        """
        return True

    @property
    def coeff(self) -> Union[complex, ParameterExpression]:
        """The scalar coefficient multiplying the Operator.

        Returns:
            The coefficient.
        """
        return self._coeff

    def primitive_strings(self) -> Set[str]:
        return reduce(set.union, [op.primitive_strings() for op in self.oplist])

    @property
    def num_qubits(self) -> int:
        num_qubits0 = self.oplist[0].num_qubits
        if not all(num_qubits0 == op.num_qubits for op in self.oplist):
            raise ValueError("Operators in ListOp have differing numbers of qubits.")
        return num_qubits0

    def add(self, other: OperatorBase) -> "ListOp":
        if self == other:
            return self.mul(2.0)

        # Avoid circular dependency
        # pylint: disable=cyclic-import
        from .summed_op import SummedOp

        return SummedOp([self, other])

    def adjoint(self) -> "ListOp":
        # TODO do this lazily? Basically rebuilds the entire tree, and ops and adjoints almost
        #  always come in pairs, so an AdjointOp holding a reference could save copying.
        if self.__class__ == ListOp:
            return ListOp(
                [op.adjoint() for op in self.oplist], **self._state(coeff=self.coeff.conjugate())
            )  # coeff is conjugated
        return self.__class__(
            [op.adjoint() for op in self.oplist], coeff=self.coeff.conjugate(), abelian=self.abelian
        )

    def traverse(
        self, convert_fn: Callable, coeff: Optional[Union[complex, ParameterExpression]] = None
    ) -> "ListOp":
        """Apply the convert_fn to each node in the oplist.

        Args:
            convert_fn: The function to apply to the internal OperatorBase.
            coeff: A coefficient to multiply by after applying convert_fn.
                If it is None, self.coeff is used instead.

        Returns:
            The converted ListOp.
        """
        if coeff is None:
            coeff = self.coeff

        if self.__class__ == ListOp:
            return ListOp([convert_fn(op) for op in self.oplist], **self._state(coeff=coeff))
        return self.__class__(
            [convert_fn(op) for op in self.oplist], coeff=coeff, abelian=self.abelian
        )

    def equals(self, other: OperatorBase) -> bool:
        if not isinstance(other, type(self)) or not len(self.oplist) == len(other.oplist):
            return False
        # Note, ordering matters here (i.e. different list orders will return False)
        return self.coeff == other.coeff and all(
            op1 == op2 for op1, op2 in zip(self.oplist, other.oplist)
        )

    # We need to do this because otherwise Numpy takes over scalar multiplication and wrecks it if
    # isinstance(scalar, np.number) - this started happening when we added __get_item__().
    __array_priority__ = 10000

    def mul(self, scalar: Union[complex, ParameterExpression]) -> "ListOp":
        if not isinstance(scalar, (int, float, complex, ParameterExpression)):
            raise ValueError(
                "Operators can only be scalar multiplied by float or complex, not "
                "{} of type {}.".format(scalar, type(scalar))
            )
        if self.__class__ == ListOp:
            return ListOp(self.oplist, **self._state(coeff=scalar * self.coeff))
        return self.__class__(self.oplist, coeff=scalar * self.coeff, abelian=self.abelian)

    def tensor(self, other: OperatorBase) -> OperatorBase:
        # Avoid circular dependency
        # pylint: disable=cyclic-import
        from .tensored_op import TensoredOp

        return TensoredOp([self, other])

    def tensorpower(self, other: int) -> Union[OperatorBase, int]:
        # Hack to make op1^(op2^0) work as intended.
        if other == 0:
            return 1
        if not isinstance(other, int) or other <= 0:
            raise TypeError("Tensorpower can only take positive int arguments")

        # Avoid circular dependency
        # pylint: disable=cyclic-import
        from .tensored_op import TensoredOp

        return TensoredOp([self] * other)

    def _expand_dim(self, num_qubits: int) -> "ListOp":
        oplist = [
            op._expand_dim(num_qubits + self.num_qubits - op.num_qubits) for op in self.oplist
        ]
        return ListOp(oplist, **self._state())

    def permute(self, permutation: List[int]) -> "OperatorBase":
        """Permute the qubits of the operator.

        Args:
            permutation: A list defining where each qubit should be permuted. The qubit at index
                j should be permuted to position permutation[j].

        Returns:
            A new ListOp representing the permuted operator.

        Raises:
            OpflowError: if indices do not define a new index for each qubit.
        """
        new_self = self
        circuit_size = max(permutation) + 1

        try:
            if self.num_qubits != len(permutation):
                raise OpflowError("New index must be defined for each qubit of the operator.")
        except ValueError:
            raise OpflowError(
                "Permute is only possible if all operators in the ListOp have the "
                "same number of qubits."
            ) from ValueError
        if self.num_qubits < circuit_size:
            # pad the operator with identities
            new_self = self._expand_dim(circuit_size - self.num_qubits)
        qc = QuantumCircuit(circuit_size)
        # extend the indices to match the size of the circuit
        permutation = (
            list(filter(lambda x: x not in permutation, range(circuit_size))) + permutation
        )

        # decompose permutation into sequence of transpositions
        transpositions = arithmetic.transpositions(permutation)
        for trans in transpositions:
            qc.swap(trans[0], trans[1])

        # pylint: disable=cyclic-import
        from ..primitive_ops.circuit_op import CircuitOp

        return CircuitOp(qc.reverse_ops()) @ new_self @ CircuitOp(qc)

    def compose(
        self, other: OperatorBase, permutation: Optional[List[int]] = None, front: bool = False
    ) -> OperatorBase:

        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
        new_self = cast(ListOp, new_self)

        if front:
            return other.compose(new_self)
        # Avoid circular dependency
        # pylint: disable=cyclic-import
        from .composed_op import ComposedOp

        return ComposedOp([new_self, other])

    def power(self, exponent: int) -> OperatorBase:
        if not isinstance(exponent, int) or exponent <= 0:
            raise TypeError("power can only take positive int arguments")

        # Avoid circular dependency
        # pylint: disable=cyclic-import
        from .composed_op import ComposedOp

        return ComposedOp([self] * exponent)

    def to_matrix(self, massive: bool = False) -> np.ndarray:
        OperatorBase._check_massive("to_matrix", True, self.num_qubits, massive)

        # Combination function must be able to handle classical values.
        # Note: this can end up, when we have list operators containing other list operators, as a
        #       ragged array and numpy 1.19 raises a deprecation warning unless this is explicitly
        #       done as object type now - was implicit before.
        mat = self.combo_fn(
            np.asarray(
                [op.to_matrix(massive=massive) * self.coeff for op in self.oplist], dtype=object
            )
        )
        # Note: As ComposedOp has a combo function of inner product we can end up here not with
        # a matrix (array) but a scalar. In which case we make a single element array of it.
        if isinstance(mat, Number):
            mat = [mat]
        return np.asarray(mat, dtype=complex)

    def to_spmatrix(self) -> Union[spmatrix, List[spmatrix]]:
        """Returns SciPy sparse matrix representation of the Operator.

        Returns:
            CSR sparse matrix representation of the Operator, or List thereof.
        """

        # Combination function must be able to handle classical values
        return self.combo_fn([op.to_spmatrix() for op in self.oplist]) * self.coeff

    def eval(
        self,
        front: Optional[
            Union[str, Dict[str, complex], np.ndarray, OperatorBase, Statevector]
        ] = None,
    ) -> Union[OperatorBase, complex]:
        """
        Evaluate the Operator's underlying function, either on a binary string or another Operator.
        A square binary Operator can be defined as a function taking a binary function to another
        binary function. This method returns the value of that function for a given StateFn or
        binary string. For example, ``op.eval('0110').eval('1110')`` can be seen as querying the
        Operator's matrix representation by row 6 and column 14, and will return the complex
        value at those "indices." Similarly for a StateFn, ``op.eval('1011')`` will return the
        complex value at row 11 of the vector representation of the StateFn, as all StateFns are
        defined to be evaluated from Zero implicitly (i.e. it is as if ``.eval('0000')`` is already
        called implicitly to always "indexing" from column 0).

        ListOp's eval recursively evaluates each Operator in ``oplist``,
        and combines the results using the recombination function ``combo_fn``.

        Args:
            front: The bitstring, dict of bitstrings (with values being coefficients), or
                StateFn to evaluated by the Operator's underlying function.

        Returns:
            The output of the ``oplist`` Operators' evaluation function, combined with the
            ``combo_fn``. If either self or front contain proper ``ListOps`` (not ListOp
            subclasses), the result is an n-dimensional list of complex or StateFn results,
            resulting from the recursive evaluation by each OperatorBase in the ListOps.

        Raises:
            NotImplementedError: Raised if called for a subclass which is not distributive.
            TypeError: Operators with mixed hierarchies, such as a ListOp containing both
                PrimitiveOps and ListOps, are not supported.
            NotImplementedError: Attempting to call ListOp's eval from a non-distributive subclass.

        """
        # pylint: disable=cyclic-import
        from ..state_fns.dict_state_fn import DictStateFn
        from ..state_fns.vector_state_fn import VectorStateFn
        from ..state_fns.sparse_vector_state_fn import SparseVectorStateFn

        # The below code only works for distributive ListOps, e.g. ListOp and SummedOp
        if not self.distributive:
            raise NotImplementedError(
                "ListOp's eval function is only defined for distributive " "ListOps."
            )

        evals = [op.eval(front) for op in self.oplist]

        # Handle application of combo_fn for DictStateFn resp VectorStateFn operators
        if self._combo_fn != ListOp([])._combo_fn:
            if (
                all(isinstance(op, DictStateFn) for op in evals)
                or all(isinstance(op, VectorStateFn) for op in evals)
                or all(isinstance(op, SparseVectorStateFn) for op in evals)
            ):
                if not all(
                    op.is_measurement == evals[0].is_measurement for op in evals  # type: ignore
                ):
                    raise NotImplementedError(
                        "Combo_fn not yet supported for mixed measurement "
                        "and non-measurement StateFns"
                    )
                result = self.combo_fn(evals)
                if isinstance(result, list):
                    multiplied = self.coeff * np.array(result)
                    return multiplied.tolist()
                return self.coeff * result

        if all(isinstance(op, OperatorBase) for op in evals):
            return self.__class__(evals)  # type: ignore
        elif any(isinstance(op, OperatorBase) for op in evals):
            raise TypeError("Cannot handle mixed scalar and Operator eval results.")
        else:
            result = self.combo_fn(evals)
            if isinstance(result, list):
                multiplied = self.coeff * np.array(result)
                return multiplied.tolist()
            return self.coeff * result

    def exp_i(self) -> OperatorBase:
        """Return an ``OperatorBase`` equivalent to an exponentiation of self * -i, e^(-i*op)."""
        # pylint: disable=unidiomatic-typecheck
        if type(self) == ListOp:
            return ListOp(
                [op.exp_i() for op in self.oplist], **self._state(abelian=False)  # type: ignore
            )

        # pylint: disable=cyclic-import
        from ..evolutions.evolved_op import EvolvedOp

        return EvolvedOp(self)

    def log_i(self, massive: bool = False) -> OperatorBase:
        """Return a ``MatrixOp`` equivalent to log(H)/-i for this operator H. This
        function is the effective inverse of exp_i, equivalent to finding the Hermitian
        Operator which produces self when exponentiated. For proper ListOps, applies ``log_i``
        to all ops in oplist.
        """
        if self.__class__.__name__ == ListOp.__name__:
            return ListOp(
                [op.log_i(massive=massive) for op in self.oplist],  # type: ignore
                **self._state(abelian=False),
            )

        return self.to_matrix_op(massive=massive).log_i(massive=massive)

    def __str__(self) -> str:
        content_string = ",\n".join([str(op) for op in self.oplist])
        main_string = "{}([\n{}\n])".format(
            self.__class__.__name__, self._indent(content_string, indentation=self.INDENTATION)
        )
        if self.abelian:
            main_string = "Abelian" + main_string
        if self.coeff != 1.0:
            main_string = f"{self.coeff} * " + main_string
        return main_string

    def __repr__(self) -> str:
        return "{}({}, coeff={}, abelian={})".format(
            self.__class__.__name__, repr(self.oplist), self.coeff, self.abelian
        )

    @property
    def parameters(self):
        params = set()
        for op in self.oplist:
            params.update(op.parameters)
        if isinstance(self.coeff, ParameterExpression):
            params.update(self.coeff.parameters)
        return params

    def assign_parameters(self, param_dict: dict) -> OperatorBase:
        param_value = self.coeff
        if isinstance(self.coeff, ParameterExpression):
            unrolled_dict = self._unroll_param_dict(param_dict)
            if isinstance(unrolled_dict, list):
                return ListOp([self.assign_parameters(param_dict) for param_dict in unrolled_dict])
            if self.coeff.parameters <= set(unrolled_dict.keys()):
                binds = {param: unrolled_dict[param] for param in self.coeff.parameters}
                param_value = float(self.coeff.bind(binds))
        return self.traverse(lambda x: x.assign_parameters(param_dict), coeff=param_value)

    def reduce(self) -> OperatorBase:
        reduced_ops = [op.reduce() for op in self.oplist]
        if self.__class__ == ListOp:
            return ListOp(reduced_ops, **self._state())
        return self.__class__(reduced_ops, coeff=self.coeff, abelian=self.abelian)

    def to_matrix_op(self, massive: bool = False) -> "ListOp":
        """Returns an equivalent Operator composed of only NumPy-based primitives, such as
        ``MatrixOp`` and ``VectorStateFn``."""
        if self.__class__ == ListOp:
            return cast(
                ListOp,
                ListOp(
                    [op.to_matrix_op(massive=massive) for op in self.oplist], **self._state()
                ).reduce(),
            )
        return cast(
            ListOp,
            self.__class__(
                [op.to_matrix_op(massive=massive) for op in self.oplist],
                coeff=self.coeff,
                abelian=self.abelian,
            ).reduce(),
        )

    def to_circuit_op(self) -> OperatorBase:
        """Returns an equivalent Operator composed of only QuantumCircuit-based primitives,
        such as ``CircuitOp`` and ``CircuitStateFn``."""
        # pylint: disable=cyclic-import
        from ..state_fns.operator_state_fn import OperatorStateFn

        if self.__class__ == ListOp:
            return ListOp(
                [
                    op.to_circuit_op() if not isinstance(op, OperatorStateFn) else op
                    for op in self.oplist
                ],
                **self._state(),
            ).reduce()
        return self.__class__(
            [
                op.to_circuit_op() if not isinstance(op, OperatorStateFn) else op
                for op in self.oplist
            ],
            coeff=self.coeff,
            abelian=self.abelian,
        ).reduce()

    def to_pauli_op(self, massive: bool = False) -> "ListOp":
        """Returns an equivalent Operator composed of only Pauli-based primitives,
        such as ``PauliOp``."""
        # pylint: disable=cyclic-import
        from ..state_fns.state_fn import StateFn

        if self.__class__ == ListOp:
            return ListOp(
                [
                    op.to_pauli_op(massive=massive)  # type: ignore
                    if not isinstance(op, StateFn)
                    else op
                    for op in self.oplist
                ],
                **self._state(),
            ).reduce()
        return self.__class__(
            [
                op.to_pauli_op(massive=massive)  # type: ignore
                if not isinstance(op, StateFn)
                else op
                for op in self.oplist
            ],
            coeff=self.coeff,
            abelian=self.abelian,
        ).reduce()

    def _is_empty(self):
        return len(self.oplist) == 0

    # Array operations:

    def __getitem__(self, offset: Union[int, slice]) -> OperatorBase:
        """Allows array-indexing style access to the Operators in ``oplist``.

        Args:
            offset: The index of ``oplist`` desired.

        Returns:
            The ``OperatorBase`` at index ``offset`` of ``oplist``,
            or another ListOp with the same properties as this one if offset is a slice.
        """
        if isinstance(offset, int):
            return self.oplist[offset]

        if self.__class__ == ListOp:
            return ListOp(oplist=self._oplist[offset], **self._state())

        return self.__class__(oplist=self._oplist[offset], coeff=self._coeff, abelian=self._abelian)

    def __iter__(self) -> Iterator:
        """Returns an iterator over the operators in ``oplist``.

        Returns:
            An iterator over the operators in ``oplist``
        """
        return iter(self.oplist)

    def __len__(self) -> int:
        """Length of ``oplist``.

        Returns:
            An int equal to the length of ``oplist``.
        """
        return len(self.oplist)

max_stars_repo_name: artmenlope/qiskit-terra

max_issues_repo_name: artmenlope/qiskit-terra

max_forks_repo_name: artmenlope/qiskit-terra


text:
# ------------------------------------------------------------------------
# Deformable DETR
# Copyright (c) 2020 SenseTime. All Rights Reserved.
# Licensed under the Apache License, Version 2.0 [see LICENSE for details]
# ------------------------------------------------------------------------
# Modified from DETR (https://github.com/facebookresearch/detr)
# Copyright (c) Facebook, Inc. and its affiliates. All Rights Reserved
# ------------------------------------------------------------------------

"""
This file provides the definition of the convolutional heads used to predict masks, as well as the losses
"""
from __future__ import division
import numpy as np
import h5py
from .smpl_utils import rot6d_to_rotmat, batch_rodrigues
import torch
import torch.nn as nn
import torch.nn.functional as F
from util.misc import NestedTensor, interpolate, nested_tensor_from_tensor_list

try:
    from panopticapi.utils import id2rgb, rgb2id
except ImportError:
    pass

from .deformable_detr import inverse_sigmoid

class DETRsmpl(nn.Module):
    def __init__(self, detr, freeze_detr=False, head_type='hmr'):
        super().__init__()
        self.detr = detr

        if freeze_detr:
            for p in self.parameters():
                p.requires_grad_(False)

        hidden_dim, nheads = detr.transformer.d_model, detr.transformer.nhead
        # self.bbox_attention = MHAttentionMap(hidden_dim, hidden_dim, nheads, dropout=0)
        # self.mask_head = MaskHeadSmallConv(hidden_dim + nheads, [1024, 512, 256], hidden_dim)
        self.head_type = head_type
        if self.head_type == 'hmr':
            self.smpl_head = HMRHead(hidden_dim)
        elif self.head_type == 'cmr':
            self.smpl_head = CMRHead(hidden_dim)
        else:
            raise KeyError('Unknown smpl head type: {}'.format(self.head_type))

    def forward(self, samples: NestedTensor):
        """ The forward expects a NestedTensor, which consists of:
               - samples.tensor: batched images, of shape [batch_size x 3 x H x W]
               - samples.mask: a binary mask of shape [batch_size x H x W], containing 1 on padded pixels

            It returns a dict with the following elements:
               - "pred_logits": the classification logits (including no-object) for all queries.
                                Shape= [batch_size x num_queries x (num_classes + 1)]
               - "pred_boxes": The normalized boxes coordinates for all queries, represented as
                               (center_x, center_y, height, width). These values are normalized in [0, 1],
                               relative to the size of each individual image (disregarding possible padding).
                               See PostProcess for information on how to retrieve the unnormalized bounding box.
               - "aux_outputs": Optional, only returned when auxilary losses are activated. It is a list of
                                dictionnaries containing the two above keys for each decoder layer.
        """
        if not isinstance(samples, NestedTensor):
            samples = nested_tensor_from_tensor_list(samples)
        features, pos = self.detr.backbone(samples)

        srcs = []
        masks = []
        for l, feat in enumerate(features):
            src, mask = feat.decompose()
            srcs.append(self.detr.input_proj[l](src))
            masks.append(mask)
            assert mask is not None
        if self.detr.num_feature_levels > len(srcs):
            _len_srcs = len(srcs)
            for l in range(_len_srcs, self.detr.num_feature_levels):
                if l == _len_srcs:
                    src = self.detr.input_proj[l](features[-1].tensors)
                else:
                    src = self.detr.input_proj[l](srcs[-1])
                m = samples.mask
                mask = F.interpolate(m[None].float(), size=src.shape[-2:]).to(torch.bool)[0]
                pos_l = self.detr.backbone[1](NestedTensor(src, mask)).to(src.dtype)
                srcs.append(src)
                masks.append(mask)
                pos.append(pos_l)

        query_embeds = None
        if not self.detr.two_stage:
            query_embeds = self.detr.query_embed.weight
        hs, init_reference, inter_references, enc_outputs_class, enc_outputs_coord_unact = self.detr.transformer(srcs, masks,
                                                                                pos,
                                                                                query_embeds)
        outputs_classes = []
        outputs_coords = []
        for lvl in range(hs.shape[0]):
            if lvl == 0:
                reference = init_reference
            else:
                reference = inter_references[lvl - 1]
            reference = inverse_sigmoid(reference)
            outputs_class = self.detr.class_embed[lvl](hs[lvl])
            tmp = self.detr.bbox_embed[lvl](hs[lvl])
            if reference.shape[-1] == 4:
                tmp += reference
            else:
                assert reference.shape[-1] == 2
                tmp[..., :2] += reference
            outputs_coord = tmp.sigmoid()
            outputs_classes.append(outputs_class)
            outputs_coords.append(outputs_coord)
        outputs_class = torch.stack(outputs_classes)
        outputs_coord = torch.stack(outputs_coords)

        out = {'pred_logits': outputs_class[-1], 'pred_boxes': outputs_coord[-1]}
        if self.detr.aux_loss:
            out['aux_outputs'] = self.detr._set_aux_loss(outputs_class, outputs_coord)

        if self.detr.two_stage:
            enc_outputs_coord = enc_outputs_coord_unact.sigmoid()
            out['enc_outputs'] = {'pred_logits': enc_outputs_class, 'pred_boxes': enc_outputs_coord}

        # SMPL paras
        ##FIXME There is no intermediate loss for SMPL paras yet
        if self.detr.aux_loss and False:
            smpl_para = self.smpl_head(hs)
            out["pred_smpl_pose"] = smpl_para[0][-1]
            out["pred_smpl_shape"] = smpl_para[1][-1]
            out["pred_camera"] = smpl_para[2][-1]

            for lvl in range(len(out['aux_outputs'])):
                out['aux_outputs'][lvl]["pred_smpl_pose"] = smpl_para[0][lvl]
                out['aux_outputs'][lvl]["pred_smpl_shape"] = smpl_para[1][lvl]
                out['aux_outputs'][lvl]["pred_camera"] = smpl_para[2][lvl]
        else:
            smpl_para = self.smpl_head(hs[-1].unsqueeze(0), outputs_class[-1].unsqueeze(0))
            out["pred_smpl_pose"] = smpl_para[0][-1]
            out["pred_smpl_shape"] = smpl_para[1][-1]
            out["pred_camera"] = smpl_para[2][-1]

        return out

"""
This file contains definitions of layers used as building blocks in SMPLParamRegressor
"""

class FCBlock(nn.Module):
    """Wrapper around nn.Linear that includes batch normalization and activation functions."""

    def __init__(self, in_size, out_size, batchnorm=True, activation=nn.ReLU(inplace=True), dropout=False):
        super(FCBlock, self).__init__()
        module_list = [nn.Linear(in_size, out_size)]
        if batchnorm:
            module_list.append(nn.BatchNorm1d(out_size))
        if activation is not None:
            module_list.append(activation)
        if dropout:
            module_list.append(dropout)
        self.fc_block = nn.Sequential(*module_list)

    def forward(self, x):
        return self.fc_block(x)

class FCResBlock(nn.Module):
    """Residual block using fully-connected layers."""

    def __init__(self, in_size, out_size, batchnorm=True, activation=nn.ReLU(inplace=True), dropout=False):
        super(FCResBlock, self).__init__()
        self.fc_block = nn.Sequential(nn.Linear(in_size, out_size),
                                      nn.BatchNorm1d(out_size),
                                      nn.ReLU(inplace=True),
                                      nn.Linear(out_size, out_size),
                                      nn.BatchNorm1d(out_size))

    def forward(self, x):
        return F.relu(x + self.fc_block(x))

"""
Definition of SMPL Parameter Regressor used for regressing the SMPL parameters from the 3D shape
"""

class BaseSMPLHead(nn.Module):
    def __init__(self):
        super().__init__()
        self.max_num = 10

    def forward(self, x, pred_class):

        stage, bs, num_query, channel = x.shape

        # only the first K bbox for saving memory
        topk_value, topk_idx = torch.topk(pred_class[..., 1], self.max_num, 2)
        idx_stage = torch.arange(stage, device=x.device, dtype=torch.long)[:, None, None]
        idx_bs = torch.arange(bs, device=x.device, dtype=torch.long)[None, :, None]
        idx = (idx_stage * bs * num_query + idx_bs * num_query + topk_idx)
        valid_top = torch.zeros([stage * bs * num_query], dtype=torch.bool, device=x.device)
        valid_top[idx] = 1

        # only foreground bbox
        x = x.flatten(start_dim=0, end_dim=-2)
        pred_class = pred_class.flatten(start_dim=0, end_dim=-2)
        valid_class = pred_class[..., 1] > pred_class[..., 0]

        valid = valid_class * valid_top

        valid[:] = True     # for debug

        num_all = x.shape[0]
        # rotmat = x.new_zeros([num_all, 24, 3, 3])
        rotmat = torch.eye(3, dtype=x.dtype, device=x.device)[None, None, :, :].repeat([num_all, 24, 1, 1])
        betas = x.new_zeros([num_all, 10])
        camera = x.new_zeros([num_all, 3])
        camera[:, 0] = 1

        # forward the head
        rotmat[valid], betas[valid], camera[valid] = self.head_forward(x[valid])

        rotmat = rotmat.view(stage, bs, num_query, 24, 3, 3)
        betas = betas.view(stage, bs, num_query, 10)
        camera = camera.view(stage, bs, num_query, 3)
        return rotmat, betas, camera

class CMRHead(BaseSMPLHead):

    def __init__(self, in_channels, use_cpu_svd=True):
        super().__init__()
        # self.layers = nn.Sequential(FCBlock(in_channels, 1024),
        #                             # FCResBlock(1024, 1024),
        #                             # FCResBlock(1024, 1024),
        #                             nn.Linear(1024, 24 * 3 * 3 + 10))

        self.layers = nn.Sequential(FCBlock(in_channels, in_channels),
                                    FCResBlock(in_channels, in_channels),
                                    nn.Linear(in_channels, 24 * 3 * 3 + 10 + 3))

        self.use_cpu_svd = use_cpu_svd

    def head_forward(self, x):
        batch_size = x.shape[0]
        x = x.view(batch_size, -1)
        x = self.layers(x)
        rotmat = x[:, :24*3*3].view(-1, 24, 3, 3).contiguous()
        betas = x[:, 24*3*3:-3].contiguous()
        camera = x[:, -3:].contiguous()

        rotmat = rotmat.view(-1, 3, 3).contiguous()
        orig_device = rotmat.device
        if self.use_cpu_svd:
            rotmat = rotmat.cpu()
        # U, S, V = batch_svd(rotmat)
        U, S, V = torch.svd(rotmat)

        rotmat = torch.matmul(U, V.transpose(1,2))

        with torch.no_grad():
            det = torch.det(rotmat)
            det = det[:, None, None]
        # det = torch.zeros(rotmat.shape[0], 1, 1).to(rotmat.device)
        # with torch.no_grad():
        #     for i in range(rotmat.shape[0]):
        #         det[i] = torch.det(rotmat[i])

        rotmat = rotmat * det
        rotmat = rotmat.view(batch_size, 24, 3, 3)
        rotmat = rotmat.to(orig_device)
        return rotmat, betas, camera

def batch_svd(A):
    """Wrapper around torch.svd that works when the input is a batch of matrices."""
    U_list = []
    S_list = []
    V_list = []
    for i in range(A.shape[0]):
        U, S, V = torch.svd(A[i])
        U_list.append(U)
        S_list.append(S)
        V_list.append(V)
    U = torch.stack(U_list, dim=0)
    S = torch.stack(S_list, dim=0)
    V = torch.stack(V_list, dim=0)
    return U, S, V

class HMRHead(BaseSMPLHead):

    def __init__(self, in_channels,
                 init_param_file='data/neutral_smpl_mean_params.h5',
                 implicity_size=-1, iteraration=3):
        super().__init__()

        # Load SMPL mean parameters
        f = h5py.File(init_param_file, 'r')
        init_grot = np.array([np.pi, 0., 0.])
        init_pose = np.hstack([init_grot, f['pose'][3:]])
        init_pose = torch.tensor(init_pose.astype('float32'))
        init_rotmat = batch_rodrigues(init_pose.contiguous().view(-1, 3))
        init_contrep = init_rotmat.view(-1, 3, 3)[:, :, ::2].contiguous().view(-1)

        init_shape = torch.tensor(f['shape'][:].astype('float32'))
        init_cam = torch.tensor([0.9, 0, 0])

        self.register_buffer('init_contrep', init_contrep)
        self.register_buffer('init_shape', init_shape)
        self.register_buffer('init_cam', init_cam)
        self.npose = init_rotmat.shape[0] * 6
        # Multiply by 6 as we need to estimate two matrixes on each joins

        self.in_channels = in_channels
        if implicity_size < 0:
            implicity_size = in_channels
        self.implicity_size = implicity_size
        self.iteration = iteraration

        self.fc_blocks = nn.Sequential(
            nn.Linear(self.in_channels + 2 * 72 + 10 + 3, self.implicity_size),
            nn.ReLU(),
            nn.Dropout(),
            nn.Linear(self.implicity_size, self.implicity_size),
            nn.Dropout(),
            nn.Linear(self.implicity_size, self.npose + 10 + 3)
        )
        self.init_weights()

    def head_forward(self, x):
        batch_size = x.shape[0]

        init_pose = self.init_contrep.view(1, -1).expand(batch_size, -1)
        init_shape = self.init_shape.view(1, -1).expand(batch_size, -1)
        init_cam = self.init_cam.view(1, -1).expand(batch_size, -1)

        theta = torch.cat([init_pose, init_shape, init_cam], 1)
        thetas = []
        for _ in range(self.iteration):
            total_inputs = torch.cat([x, theta], 1)
            theta = theta + self.fc_blocks(total_inputs)
            thetas.append(theta)

        pred_pose = theta[:, :self.npose].contiguous()
        pred_betas = theta[:, self.npose:self.npose + 10].contiguous()
        pred_camera = theta[:, self.npose + 10:].contiguous()
        pred_rotmat = rot6d_to_rotmat(pred_pose).view(batch_size, 24, 3, 3)
        return pred_rotmat, pred_betas, pred_camera

    def init_weights(self):
        for n, p in self.named_parameters():
            if 'weight' in n:
                nn.init.xavier_uniform_(p)
            if 'bias' in n:
                nn.init.zeros_(p)

ext: py

reason: import numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Aug  6 10:24:07 2020

@author: TeSolva
"""

import numpy as np

import pandas as pd

import matplotlib.dates as mdates

import matplotlib.pyplot as plt

import os

import re

#%%

def my_int(my_string):

    try:

        out = int(my_string)

    except:

        out = -9999

    return out

#%%
def func_nummer(x):

    temp = re.findall(r"[\w']+", x)

    my_list=list(map(my_int, temp))

    temp = np.array(my_list)[np.array(my_list)>-9999]

    if len(temp) > 1:

        out = temp[0]

    else:

        out = False

    return out
#%%

def func_id(x):

    temp = re.findall(r"[\w']+", x)

    my_list=list(map(my_int, temp))

    temp = np.array(my_list)[np.array(my_list)>-9999]

    if len(temp) > 1:

        out = temp[1]

    else:

        out = False

    return out
#%%

def outlier_1d_mad_based(sample, thresh=3.5):

    """

    outlier_1d_mad_based(sample, thresh=3.5)

    routine to analyse a given 1d data sample to check for outliers.

    see reference for more details on the background of the used algorithm.

    the function returns a boolean array with True if a value in the sample

    is an outliers and False otherwise.

    Parameters:

    -----------

        sample : array_like

            An numobservations by numdimensions array of observations

        thresh : float

            The modified z-score to use as a threshold. Observations with

            a modified z-score (based on the median absolute deviation) greater

            than this value will be classified as outliers.

    Returns:

    --------

        A numobservations-length boolean array.

    Examples

    --------

        # Generate some data

        sample = np.random.normal(0, 0.5, 50)

        # Add three outliers...

        sample = np.r_[sample, -3, -10, 12]

        # call function and check for outliers

        out = outlier_1d_mad_based(sample)

    References:

    ----------

        Boris Iglewicz and David Hoaglin (1993), "Volume 16: How to Detect and

        Handle Outliers", The ASQC Basic References in Quality Control:

        Statistical Techniques, Edward F. Mykytka, Ph.D., Editor.

    """

    if len(sample.shape) == 1:

        sample = sample[:, None]

    median = np.median(sample, axis=0)

    diff = np.sum((sample - median) ** 2, axis=-1)

    diff = np.sqrt(diff)

    med_abs_deviation = np.median(diff)

    modified_z_score = 0.6745 * diff / med_abs_deviation

    return modified_z_score > thresh

#%%

# #path = "PATH\"

# # define all filenames in current directory

# path = os.getcwd() # path to this file here

# list_dir = os.listdir(os.getcwd()) # all filenames in the directory

# file_set = []

# for i in list_dir:

#     if i.endswith(".csv"):

#         file_set.append(i)
#%%

# define file name

file1 = "chronicreplay-audiothek-appapi.tesolva.dev_2020-08-08_07-29-08.csv"

file2 = "chronicreplay-audiothek-appapi.solr.tesolva.dev_2020-08-11_08-47-32.csv"

#

file_set = [file1, file2]

###############################################################################

#%% pre-process
X_file = file_set[1]
filename = os.path.splitext(X_file)[0]

path = os.getcwd()
outputPath = path + '/' + filename + '/' + filename
os.makedirs(filename, exist_ok=True)

df = pd.read_csv(X_file,sep='\t')

df.StartTime = pd.to_datetime(df.StartTime)

# create new columns

df['Request_No']=df.Request

df['Request_ID']=df.Request

# add values (No and ID) from URL to new columns

df.Request_No = df.Request_No.apply(func_nummer)

df.Request_ID = df.Request_ID.apply(func_id)

# Requests without /health and /metric

df = df[df['Request_No'] != False]

#%%
# Threshold to find URLs with high values of request time

THLD = 1000 # Threshold

df['Duration_greaterThan_1000'] = (df.Duration>THLD).astype(int)

#%%
plt.figure(figsize = (8.5,11))
myFmt = mdates.DateFormatter('%H:%M')
plt.gca().xaxis.set_major_formatter(myFmt)

plt.scatter(df.StartTime[df.Duration<THLD].values, df.Duration[df.Duration<THLD].values, marker='^', label='Duration < THLD [ms]',

               alpha=0.3, edgecolors='none')

plt.scatter(df.StartTime[df.Duration>THLD].values, df.Duration[df.Duration>THLD].values, marker='o', label='Duration > THLD [ms]',

               alpha=0.3, edgecolors='none')

# Show the boundary between the regions:
plt.plot(df.StartTime.values, np.ones(len(df))*THLD,'-.k')

plt.legend()
plt.ylabel("Duration of request time (ms)")
plt.xlabel("Time")
plt.title("file:" + X_file)
plt.axis([None, None, 0, THLD*2])

plt.savefig(outputPath + '-Duration.png')

#%%
plt.figure(figsize = (8.5,11))
myFmt = mdates.DateFormatter('%H:%M')
plt.gca().xaxis.set_major_formatter(myFmt)

plt.scatter(df.StartTime[df.Duration<THLD].values, df.Difference[df.Duration<THLD].values, marker='^', label='Duration < THLD [ms]',

               alpha=0.3, edgecolors='none')

plt.scatter(df.StartTime[df.Duration>THLD].values, df.Difference[df.Duration>THLD].values, marker='o', label='Duration > THLD [ms]',

               alpha=0.3, edgecolors='none')

# Show the boundary between the regions:
plt.plot(df.StartTime.values, np.ones(len(df))*THLD,'-.k')

plt.legend()
plt.ylabel("Difference of request time (ms)")
plt.xlabel("Time")
plt.title("file:" + X_file)
plt.axis([None, None, 0, THLD*2])

plt.savefig(outputPath + '-Difference.png')

#%%
# only interested in URLs with request time greater than THLD

df_Upper = df[df.Duration>THLD]

df_Upper['Percent'] = len(df_Upper)/len(df.index)*100
df_Upper['Mean'] = df_Upper["Duration"].mean()

# export csv files
df_Upper[['Duration', 'Request_No', 'Request_ID', 'Percent','Mean', 'Request']].to_csv(outputPath + '-URLs-request-time-Upper.csv', index=False)

#%%
# count Request_No for Requests time > THLD

df_temp = (df_Upper.drop_duplicates().Request_No.value_counts())

df_Counts = pd.DataFrame({'Request_No':df_temp.index, 'Counts':df_temp.values})

df_Counts.to_csv(outputPath + '-RequestsNo-Counts-Upper.csv', index=False)

#%%
# values between 500 and 1000 ms

df_Between = df[df['Duration'].between(500, 1000, inclusive=True)]

# export csv files

df_Between['Percent'] = len(df_Between)/len(df.index)*100

df_Between['Mean'] = df_Between["Duration"].mean()

df_Between[['Duration', 'Request_No', 'Request_ID', 'Percent','Mean', 'Request']].to_csv(outputPath + '-URLs-request-time-Between.csv', index=False)

# count Request_No for Requests time  500 < time < 1000

df_temp2 = (df_Between.drop_duplicates().Request_No.value_counts())

df_CountsBetween = pd.DataFrame({'RequestNo':df_temp2.index, 'Counts':df_temp2.values})

df_CountsBetween.to_csv(outputPath + '-RequestsNo-Counts-Between.csv', index=False)

#%%
# values lower than 500 ms
df_Lower = df[df['Duration'].between(0, 500, inclusive=False)]

# export csv files

df_Lower['Percent'] = len(df_Lower)/len(df.index)*100

df_Lower['Mean'] = df_Lower["Duration"].mean()

df_Lower[['Duration', 'Request_No', 'Request_ID', 'Percent','Mean', 'Request']].to_csv(outputPath + '-URLs-request-time-Lower.csv', index=False)

# count Request_No for Requests time  < 500

df_temp2 = (df_Lower.drop_duplicates().Request_No.value_counts())

df_CountsLower = pd.DataFrame({'RequestNo':df_temp2.index, 'Counts':df_temp2.values})

df_CountsLower.to_csv(outputPath + '-RequestsNo-Counts-Lower.csv', index=False)

#%%

# df_notSameStatus = df[df.SameStatus==False]

# # export csv files
# df_notSameStatus.to_csv(outputPath + '-URLs-wo-same-status.csv', index=False)

#%%
del(X_file,file1, file2, file_set) # delete the variable from workspace
#%%

max_stars_repo_name: paulwellnerbou/chronicreplay

max_issues_repo_name: paulwellnerbou/chronicreplay

max_forks_repo_name: paulwellnerbou/chronicreplay


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 18 11:15:00 2017

@author: Prashant Kumar, RCCIIT

@mail: kr.prashant94@gmail.com

@supervision: Dilip Kumar , NIT Jsr

"""

import matplotlib.patches as patches
import matplotlib.pyplot as plt
import csv
import time
import datetime
import random
import numpy as np
import matplotlib.animation as animation
import threading
import pandas
import os
import serial
import io

'''
Global declaration
-------
- dataColumns = ['Time','Data','Average','SD']
    column in pandas database.
- dataALog = pandas.DataFrame(columns = dataColumns)
    Pandas DataFrame form input A.
- dataBLog = pandas.DataFrame(columns = dataColumns)
    Pandas DataFrame form input A.

- fig = plt.figure()
    Figure in ploat.
- fig.subplots_adjust(hspace = 0.8, wspace = 0.4)
    Global subplot adjustment.
- plt.rc('font', size=8)
    Default text sizes in subplot.

- ax1 = fig.add_subplot(2,3,1)
    Axis 1 for A vs B graph.
- ax2 = fig.add_subplot(2,3,2)
    Axis 2 for A vs Average of A graph.
- ax3 = fig.add_subplot(2,3,3)
    Axis 3 for B vs Average of B graph.
- ax4 = fig.add_subplot(2,3,4)
    Axis 4 for A vs SD of A graph.
- ax5 = fig.add_subplot(2,3,5)
    Axis 5 for B vs SD of B graph.
- ax6 = fig.add_subplot(2,3,6)
    Axis 6 for showing current status.

- ax6.get_xaxis().set_visible(False), ax6.get_yaxis().set_visible(False)
    Hiding axis visibility *import matplotlib.pyplot as plt*.
- plt.axis('off')
    Off plt.axis('off'),
- xdata, dataA, dataB, dataAmean, dataBmean, dataAsd, dataBsd = [0], [0], [0], [0], [0], [0], [0]
    Initializing data for graph and log.
- labels = ['0:0']
    Initializing label for graph.
- ln1, = ax1.plot([0], [0], 'red', linewidth = 1 )
    Line 1 for plot of A in Axis 1.
- ln2, = ax1.plot([0], [0], 'green', linewidth = 1 )
    Line 2 for plot of B in Axis 1.
- ln3, = ax2.plot([0], [0], 'red', linewidth = 1 )
    Line 3 for plot of A in Axis 2.
- ln4, = ax2.plot([0], [0], 'blue', linestyle = '--', linewidth = 1)
    Line 4 for plot of Average of A in Axis 2.
- ln5, = ax3.plot([0], [0], 'green', linewidth = 1 )
    Line 5 for plot of B in Axis 3.
- ln6, = ax3.plot([0], [0], 'blue', linestyle = '--', linewidth = 1)
    Line 4 for plot of Average of B in Axis 3.
- ln7, = ax4.plot([0], [0], 'red', linewidth = 1 )
    Line 5 for plot of A in Axis 4.
- ln8, = ax4.plot([0], [0], 'black', linestyle = '--', linewidth = 1)
    Line 5 for plot SD of A in Axis 4.
- ln9, = ax5.plot([0], [0], 'green', linewidth = 1 )
    Line 5 for plot of B in Axis 5.
- ln10, = ax5.plot([0], [0], 'black', linestyle = '--', linewidth = 1)
    Line 5 for plot SD of B in Axis 5.

- red_patch = patches.Patch(color='red', label='Data A')
    Red colour patch information for Data A.
- green_patch = patches.Patch(color='green', label='Data B')
    Green colour patch information for Data B.
- avg_patch = patches.Patch(color='blue', label='Avgerage')
     Blue colour patch information for Average of Data.
- sd_patch = patches.Patch(color='black', label='SD')
        Black colour patch information for SD of Data.

- initTime = time.time()
    Initial time when program started.
- threads = []
    List of threads started by program. Thread used to save data in .CSV database.
- yLimits = 100
    Maximum height of y axis.

'''

class pyAnalysis:
    dataColumns = ['Time','Data','Average','SD']
    dataALog = pandas.DataFrame(columns = dataColumns)
    dataBLog = pandas.DataFrame(columns = dataColumns)

    fig = plt.figure()
    fig.subplots_adjust(hspace = 0.8, wspace = 0.4)
    plt.rc('font', size=8)          # controls default text sizes

    ax1 = fig.add_subplot(2,3,1)
    ax2 = fig.add_subplot(2,3,2)
    ax3 = fig.add_subplot(2,3,3)
    ax4 = fig.add_subplot(2,3,4)
    ax5 = fig.add_subplot(2,3,5)
    ax6 = fig.add_subplot(2,3,6)

    ax6.get_xaxis().set_visible(False)
    ax6.get_yaxis().set_visible(False)
    plt.axis('off')

    xdata, dataA, dataB, dataAmean, dataBmean, dataAsd, dataBsd = [0], [0], [0], [0], [0], [0], [0]
    labels = ['0:0']
    ln1, = ax1.plot([0], [0], 'red', linewidth = 1 )
    ln2, = ax1.plot([0], [0], 'green', linewidth = 1 )
    ln3, = ax2.plot([0], [0], 'red', linewidth = 1 )
    ln4, = ax2.plot([0], [0], 'blue', linestyle = '--', linewidth = 1)
    ln5, = ax3.plot([0], [0], 'green', linewidth = 1 )
    ln6, = ax3.plot([0], [0], 'blue', linestyle = '--', linewidth = 1)
    ln7, = ax4.plot([0], [0], 'red', linewidth = 1 )
    ln8, = ax4.plot([0], [0], 'black', linestyle = '--', linewidth = 1)
    ln9, = ax5.plot([0], [0], 'green', linewidth = 1 )
    ln10, = ax5.plot([0], [0], 'black', linestyle = '--', linewidth = 1)

    red_patch = patches.Patch(color='red', label='Pollutent (ppm)')
    green_patch = patches.Patch(color='green', label='CO (ppm)')
    avg_patch = patches.Patch(color='blue', label='Avgerage')
    sd_patch = patches.Patch(color='black', label='SD')

    initTime = time.time()
    threads = []

    yLimits = 100

    def __init__(self, port):
        #Initilize the dataBase.csv file to Start writing
        self.ser = serial.Serial()
        self.ser.port = port
        if(os.path.isfile(os.path.dirname(__file__)+"/dataBaseA.csv")):
            self.fetchDatabase()
            self.dataA[0] = dataALog.iloc[-1]['Data']
            self.dataB[0] = dataBLog.iloc[-1]['Data']

            self.dataAmean[0] = dataALog.iloc[-1]['Average']
            self.dataBmean[0] = dataBLog.iloc[-1]['Average']

            self.dataAsd[0] = dataALog.iloc[-1]['SD']
            self.dataBsd[0] = dataBLog.iloc[-1]['SD']
        else:
            self.createCSV()
        return

    #Data function

    def getNewData(self, number, error):
        """
        Description:
        -------
		NA

        Parameters
        -------
        **number** : *(int or float)* :
            The base number

        **error** : *(int or float)* :
            Percentage fluctuation or error in new number.

        Return
        -------
        *float* :
            A new random float within error limit and round up to 2 digits.

        Examples
        -------
        >>> self.getNewData(100, 5)
        102.22
        >>> self.getNewData(100, 5)
        97.70
        >>> self.getNewData(100, 5)
        104.35
        """
        self.ser.open()
        tmp = self.ser.readline().decode().strip()
        data = tmp.split(',')
        self.ser.close()
        return int(data[0].strip()) , int(data[1].strip())

    def average(self, data, log, n):
        """
        Description
        --------
        **average()** return the mathematical mean of the last *n-1* data form *log* and *data*

        Parameters
        --------
        **data** : *(float)*
            The input number.
        **log** : *(pandas.DataFrame)*
            last n-1 data log (at least) for getting mean.
        **n** : *(int)*
        (N-1)th number form log to consider mean.

        Return
        --------
        *float* :
            Return the mathematical mean of the last *n-1* data form *log* and *data*

        Examples
        --------
        >>> self.average(inputData, dataLog, n)
        """
        if(len(log) > n-1):
            dataList = list(log.iloc[-n:-1]['Data'])
            dataList.append(data)
            return round(np.mean(dataList, axis=0), 2)
        else:
            return 0

    def sd(self, data, log, n):
        """
        Description
        --------
        **sd()** return the mathematical Standard Deviation(SD) of the last *n-1* data form *log* and *data*

        Parameters
        --------
        **data** : *(float)*
            The input number.
        **log** : *(pandas.DataFrame)*
            last n-1 data log(at least) for getting SD.
        **n** : *(int)*
            (N-1)th number form log to consider SD.

        Return
        --------
        *float* :
            Return the mathematical SD of the last *n-1* data form *log* and *data*

        Examples
        --------
        >>> self.sd(inputData, dataLog, n)
        """
        if(len(log) > n-1):
            dataList = list(log.iloc[-n:-1]['Data'])
            dataList.append(data)
            return round(np.std(dataList, axis=0),2)
        else:
            return 0

    def Save(self, inputDataA, inputDataB):
        """
        Description
        --------
        **Save()** save the input data into .CSV database

        Parameters
        --------
        **inputDataA** : *(pandas.DataFrame)*
            The input Data A with avg. and SD.
        **inputDataB** : *(pandas.DataFrame)*
            The input Data B with avg. and SDs.

        Return
        --------
        *void*

        Examples
        --------
        >>> self.Save(inputDataA, inputDataB)
        """
        #os.system('cls')
        #print(pandas.concat([self.dataALog, self.dataBLog], axis=1))
        with open('dataBaseA.csv', 'a') as csvfile:
            writer = csv.writer(csvfile, delimiter=',', quotechar=' ', lineterminator='\n', quoting=csv.QUOTE_MINIMAL)
            writer.writerow(inputDataA)

        with open('dataBaseB.csv', 'a') as csvfile:
            writer = csv.writer(csvfile, delimiter=',', quotechar=' ', lineterminator='\n', quoting=csv.QUOTE_MINIMAL)
            writer.writerow(inputDataB)

    def createCSV(self):
        """
        Description
        --------
        **createCSV()** create two .CSV file dataBaseA.csv and dataBaseB.csv for data logging

        Parameters
        --------
        *void*

        Return
        --------
        *void*

        Examples
        --------
        >>> self.createCSV()
        """
        with open('dataBaseA.csv', 'w') as csvfile:
            writer = csv.writer(csvfile, delimiter=',', quotechar=' ', lineterminator='\n', quoting=csv.QUOTE_MINIMAL)
            writer.writerow(self.dataColumns)
            writer.writerow([''])
        with open('dataBaseB.csv', 'w') as csvfile:
            writer = csv.writer(csvfile, delimiter=',', quotechar=' ',lineterminator='\n', quoting=csv.QUOTE_MINIMAL)
            writer.writerow(self.dataColumns)
            writer.writerow([''])
    def fetchDatabase(self):
        """
        Description
        --------
        **fetchDatabase()** fetch the data in pandas.DataFrame format and save it to global variable dataALog and dataBLog

        Parameters
        --------
        *void*

        Return
        --------
        *void*

        Examples
        --------
        >>> self.fetchDatabase()
        """
        global dataALog, dataBLog
        dataALog = pandas.read_csv(os.path.dirname(__file__)+"/dataBaseA.csv")
        dataBLog = pandas.read_csv(os.path.dirname(__file__)+"/dataBaseB.csv")

    '''
        Graph Function
    '''
    def init(self):
        """
        Description
        --------
        **init()** is used for initialize the graph.

        Parameters
        --------
        *void*

        Returns
        --------
        *void*
        """
        self.ax1.set_title('All Pollutent vs CO')
        self.ax2.set_title('All Pollutent vs Average')
        self.ax3.set_title('CO vs Average')
        self.ax4.set_title('All Pollutent vs SD')
        self.ax5.set_title('CO vs SD')

        self.setAxisLimits(self.ax1)
        self.setAxisLimits(self.ax2)
        self.setAxisLimits(self.ax3)
        self.setAxisLimits(self.ax4)
        self.setAxisLimits(self.ax5)

        self.ax1.legend(handles=[self.red_patch, self.green_patch])

        self.ax2.legend(handles=[self.red_patch, self.avg_patch])
        self.ax3.legend(handles=[self.green_patch, self.avg_patch])

        self.ax4.legend(handles=[self.red_patch, self.sd_patch])
        self.ax5.legend(handles=[self.green_patch, self.sd_patch])

        return

    def graph(self, frame):
        """
        Description
        --------
        **graph()** plot and update graph in every unit time

        Parameters
        --------
        *frame*

        Return
        --------
        *ln1, ln2, ln3, ln4, ln5, ln6, ln7, ln8, ln9, ln10*
            All line in plot
        """
        global yLimits, dataALog, dataBLog

        #Time Sync.
        dateTime = time.time()
        time.sleep(1 - (dateTime-self.initTime - int(dateTime-self.initTime)))

        #Inputs
        inputDataA, inputDataB = self.getNewData(self.dataA[-1], 5)

        dateTime = time.time() #Timestamp
        dateTimeNow = str(datetime.datetime.fromtimestamp(dateTime).strftime('%H:%M:%S ')) # Time in Calander format
        timeWhole = int(dateTime-self.initTime)

        #change Y-axis
        if(inputDataA >= self.yLimits or inputDataB >= self.yLimits):
            self.yLimits = inputDataA+5 if (inputDataA > inputDataB) else inputDataB+5
            self.ax1.set_ylim(0, self.yLimits)
            self.ax2.set_ylim(0, self.yLimits)
            self.ax3.set_ylim(0, self.yLimits)
            self.ax4.set_ylim(0, self.yLimits)
            self.ax5.set_ylim(0, self.yLimits)

        #Change X-axis
        if timeWhole > 11:
            self.labels.pop(0)
            self.xdata.pop(0)
            self.dataA.pop(0)
            self.dataB.pop(0)
            self.dataAmean.pop(0)
            self.dataBmean.pop(0)
            self.dataAsd.pop(0)
            self.dataBsd.pop(0)
            self.ax1.set_xlim(timeWhole-10, timeWhole)
            self.ax2.set_xlim(timeWhole-10, timeWhole)
            self.ax3.set_xlim(timeWhole-10, timeWhole)
            self.ax4.set_xlim(timeWhole-10, timeWhole)
            self.ax5.set_xlim(timeWhole-10, timeWhole)

        dataALog = dataALog.append({'Time':dateTimeNow,'Data':round(inputDataA,2),'Average':self.average(inputDataA, dataALog, 10),'SD':self.sd(inputDataA, dataALog, 10)}, ignore_index=True)
        dataBLog = dataBLog.append({'Time':dateTimeNow,'Data':inputDataB,'Average':self.average(inputDataB, dataBLog, 10),'SD':self.sd(inputDataB, dataBLog, 10)}, ignore_index=True)

        self.ax6.clear()
        self.ax6.text(0.1, 0, '$Status$ \nTime    : '+dateTimeNow+' \nData A : '+str(inputDataA)+' \nData B : '+str(inputDataB)+' \nAvg A  : '+str(dataALog.iloc[-1]['Average'])+' \nAvg B  : '+str(dataBLog.iloc[-1]['Average'])+' \nSD A    : '+str(dataALog.iloc[-1]['SD'])+' \nSD B    : '+str(dataBLog.iloc[-1]['SD'])+' ', fontsize=10)
        plt.axis('off')

        self.xdata.append(timeWhole)
        self.dataA.append(inputDataA)
        self.dataB.append(inputDataB)
        self.dataAmean.append(dataALog.iloc[-1]['Average'])
        self.dataBmean.append(dataBLog.iloc[-1]['Average'])
        self.dataAsd.append(dataALog.iloc[-1]['SD'])
        self.dataBsd.append(dataBLog.iloc[-1]['SD'])
        self.labels.append(dateTimeNow)

        plt.setp(self.ax1, xticks=self.xdata, xticklabels=self.labels)
        plt.setp(self.ax2, xticks=self.xdata, xticklabels=self.labels)
        plt.setp(self.ax3, xticks=self.xdata, xticklabels=self.labels)
        plt.setp(self.ax4, xticks=self.xdata, xticklabels=self.labels)
        plt.setp(self.ax5, xticks=self.xdata, xticklabels=self.labels)

        self.ln1.set_data(self.xdata, self.dataA)
        self.ln2.set_data(self.xdata, self.dataB)
        self.ln3.set_data(self.xdata, self.dataA)
        self.ln4.set_data(self.xdata, self.dataAmean)
        self.ln5.set_data(self.xdata, self.dataB)
        self.ln6.set_data(self.xdata, self.dataBmean)
        self.ln7.set_data(self.xdata, self.dataA)
        self.ln8.set_data(self.xdata, self.dataAsd)
        self.ln9.set_data(self.xdata, self.dataB)
        self.ln10.set_data(self.xdata, self.dataBsd)

        #Threaded Work
        t = threading.Thread(target=self.Save, args=(list(dataALog.iloc[-1]), list(dataBLog.iloc[-1]),))
        self.threads.append(t)
        t.start()

        return

    def setAxisLimits(self, ax):
        """
        Description
        --------
        **setAxisLimits()** set the rang of x-axis and y-axis according to available data.

        Parameters
        --------
        **ax** : *(fig.subplot)*
            The axes of the subplot.

        Return
        --------
            *void*

        Examples
        --------
        >>> self.setAxisLimits(ax1)
        >>> self.setAxisLimits(ax2)
        >>> self.setAxisLimits(ax3)
        """
        ax.set_xlim(0, 11)
        ax.set_ylim(0, self.yLimits)
        self.rotateLabel(ax)

    def rotateLabel(self, ax):
        """
        Description
        -------
        **rotateLabel()** rotates the labels of the graph by 90deg.

        Parameters
        ----------
        **ax** : *(fig.add_subplot)*
            ax is subplot  graph in the figure. must be *matplotlib subplot*

        Return
        -------
        *void*

        Examples
        --------
        >>> self.rotateLabel(ax1)
        """
        for label in ax.get_xmajorticklabels() + ax.get_xmajorticklabels():
            label.set_rotation(90)
            label.set_horizontalalignment("right")

    def genrateLineAnalysis(self):
        ani = animation.FuncAnimation(self.fig, self.graph, init_func=self.init,
                                      interval=500, blit=False)
        plt.show()


text:
#!/usr/bin/env python
"""
strike_hist.py (python3)
Zhiang Chen, Sep 2019

Copyright (c) 2019 Distributed Robotic Exploration and Mapping Systems Laboratory, ASU
"""

import os
import cv2
import gdal
import pickle
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
from math import sqrt

class Strike_analysis(object):
    def __init__(self, skeleton_path, orth_path, instance_path, ):
        assert os.path.isfile(skeleton_path)
        assert os.path.isfile(orth_path)
        assert os.path.isfile(instance_path)

        self.skeleton = cv2.imread(skeleton_path)[:, :, 0]
        sk = np.max(self.skeleton, axis=1)
        sk = np.where(sk > 10)[0]
        self.sk_up, self.sk_down = sk[0], sk[-1]

        self.orth = gdal.Open(orth_path)
        coord = self.orth.GetGeoTransform()
        self.x, self.x_size, _, self.y, _, self.y_size = coord
        self.h, self.w = self.orth.GetRasterBand(1).ReadAsArray().shape
        self.X = self.x + self.x_size * self.w
        self.Y = self.y + self.y_size * self.h

        with open(instance_path, 'rb') as pk:
            self.instances = pickle.load(pk)

    def create_heatmap(self):
        B = self.orth.GetRasterBand(1).ReadAsArray()
        G = self.orth.GetRasterBand(2).ReadAsArray()
        R = self.orth.GetRasterBand(3).ReadAsArray()
        self.histmap = np.zeros((self.h, self.w, 3))
        self.histmap[:, :, 0], self.histmap[:, :, 1], self.histmap[:, :, 2] = R, G, B  # B, G, R

    def get_normal_vector(self, latitude_up, latitude_down):
        """
        get normal vector of the skeleton given latitude bound
        :param latitude_up:
        :param latitude_down:
        :return: intersection, norm vector (both in pixel coordinates)
        """
        pixel_up = self.__lad2pixel(latitude_up)
        pixel_down = self.__lad2pixel(latitude_down)
        assert self.sk_up < pixel_up < pixel_down < self.sk_down
        latitude = (latitude_up + latitude_down)/2
        intersection = self.__get_intersection(latitude)

        section = self.skeleton[pixel_up:pixel_down, :]
        points = np.where(section > 10)
        points = np.array(points)
        points = np.swapaxes(points, axis1=0, axis2=1)
        points[:, 0] = points[:, 0] + pixel_up
        a, b = np.polyfit(points[:, 1], points[:, 0], deg=1)  # y = a*x + b
        # Note the coordinates are the pixel coordinates
        norm_v = np.array((-1/a, 1.0))
        norm_v = norm_v / np.linalg.norm(norm_v)  # (y, x)

        return intersection, norm_v

    def generate_one_box(self, center, normal, offset, h, w):
        """
        generate one box that is 'offset' off the 'center' along 'normal' vector with height 'h' and width 'w'
        :param center: in pixel coordinates
        :param normal: in pixel coordinates
        :param offset: in pixel
        :param h: in pixel
        :param w: in pixel
        :return: box coordinates in pixel coordinates, top_left(y, x), top_right(y, x), bottom_right(y, x), bottom_left(y, x)
        """
        p1 = (-h / 2.0, -w / 2.0)
        p2 = (-h / 2.0, w / 2.0)
        p3 = (h / 2.0, w / 2.0)
        p4 = (h / 2.0, -w / 2.0)
        points = np.array((p1, p2, p3, p4))
        s, c = normal
        s, c = -c, s
        r = np.array(((c, -s), (s, c)))
        points_ = np.matmul(points, r)

        """
        print(normal)
        print(points)
        print(points_)
        plt.plot(points[:, 1], points[:, 0], 'x')
        plt.plot(points_[:, 1], points_[:, 0], 'ro')
        plt.plot((0, normal[1]*10), (0, normal[0]*10), 'o')
        plt.ylim(-100, 100)
        plt.xlim(-100, 100)
        plt.show()
        """
        points_ = points_ + center + offset*normal
        return points_.astype(int)

    def generate_boxes(self, center, normal, h, w, step, num, symmetric=True):
        """
        generate boxes
        :param center: pixel coordinates
        :param normal: pixel coordinates
        :param h: in pixel
        :param w: in pixel
        :param step: in pixel
        :param num: box number on one side
        :param symmetric: bool
        :return: boxes in pixel coordinates
        """
        offsets = [step*i for i in range(num+1)[1:]]

        if symmetric:
            offsets = [step * i for i in range(-num, num + 1)]

        boxes = [self.generate_one_box(center, normal, offset, h, w) for offset in offsets]
        return boxes

    def get_instances_in_box(self, box):
        """
        get instances in a box defined by points
        :param box:
        :return: instances in the box, list
        """
        polygon = Polygon(box)
        boxed_instances = []
        for instance in self.instances:
            y, x, Y, X = instance['bb']
            point = Point((y+Y)/2.0, (x+X)/2.0,)
            if polygon.contains(point):
                boxed_instances.append(instance)
        #print(len(boxed_instances))
        return boxed_instances

    def draw_box_in_skeleton(self, points, sk, color=255):
        """
        example:
            intersection, norm_v = sa.get_normal_vector(4146290, 4146285)
            points = sa.generate_one_box(intersection, norm_v, 0, 30, 100)
            box_sk = sa.draw_box_in_skeleton(points, sa.skeleton)
            cv2.imwrite("box_sk.jpg", box_sk)
        :param points:
        :param sk:
        :return:
        """
        n, _ = points.shape
        points = points.reshape((1, n, 2))
        box = points.copy()
        box[:, :, 0], box[:, :, 1] = points[:, :, 1], points[:, :, 0]
        sk_box = cv2.fillPoly(sk, box, color)
        return sk_box

    def draw_instances_in_skeleton(self, instances, sk):
        for instance in instances:
            mask = instance['mask']
            for y,x in mask:
                sk[y, x] = 255
        return sk

    def __get_intersection(self, latitude):
        """
        get an intersection with the skeleton given a latitude
        :param latitude: latitude in meters
        :return: intersection
        """
        pixel_y = self.__lad2pixel(latitude)
        assert (pixel_y >= 0) & (pixel_y < self.h)
        x = self.skeleton[pixel_y, :]
        pixel_x = np.where(x > 10)
        assert len(pixel_x[0].tolist()) > 0
        pixel_x = int(np.average(pixel_x))
        intersection = np.array((pixel_y, pixel_x)).astype(int)
        return intersection

    def __utm2pixel(self, utm_y, utm_x):
        """
        convert utm to pixel coordinates
        :param utm_x:
        :param utm_y:
        :return:
        """
        y = (utm_y - self.y) / self.y_size
        x = (utm_x - self.x) / self.x_size
        pixel_coord = np.array((y, x)).astype(int)
        return pixel_coord

    def __lad2pixel(self, utm_y):
        """
        convert utm latitude to pixel coordinate
        :param utm_y:
        :return:
        """
        y = (utm_y - self.y) / self.y_size
        return int(y)

    def __pixel2lad(self, pixel):
        utm_y = self.y + pixel * self.y_size
        return utm_y

    def __create_bool_mask(self, mask, size):
        """
        maybe there are more efficient ways to do this...
        :param mask: mask by index
        :param size: size of image
        :return: bool mask
        """
        mask_ = np.zeros(size)
        for y, x in mask:
            if (y < size[0]) & (x < size[1]):
                mask_[y, x] = 1
        return mask_

    def __add_hist2colormap(self, hist, box, side_box, box_id = None, box_num=None, width_vector=None, brightness=None, colormap_max=False):
        scale = 0.5 # the scale of side box
        num = len(hist)
        p1 = box[0]
        p4 = box[3]
        direction = box[1] - box[0]
        direction_step = direction / num
        if not side_box:
            for i in range(num):
                b1 = p1 + direction_step * i
                b2 = b1 + direction_step
                b4 = p4 + direction_step * i
                b3 = b4 + direction_step
                subbox = np.array((b1, b2, b3, b4), dtype=int)
                subbox = subbox.reshape((1, 4, 2))
                points = subbox.copy()
                subbox[:, :, 0], subbox[:, :, 1] = points[:, :, 1], points[:, :, 0]
                if not colormap_max:
                    if brightness == None:
                        self.histmap = cv2.fillPoly(self.histmap, subbox, (int(hist[i]), 0, int(hist[i])))
                    else:
                        if hist[i] > 1:
                            self.histmap = cv2.fillPoly(self.histmap, subbox, (brightness, 0, brightness))  # only for mean (one value) plotting
                        else:
                            self.histmap = cv2.fillPoly(self.histmap, subbox, (int(hist[i]), 0, int(hist[i])))
                else:
                    cmap = plt.cm.get_cmap('jet')
                    value = hist.max()
                    value = int(255.0*value/colormap_max)
                    color = np.array(cmap(value)[:3])*255
                    #color = cv2.applyColorMap(, cv2.COLORMAP_JET)  # it somehow doesn't work
                    self.histmap = cv2.fillPoly(self.histmap, subbox, (int(color[2]), int(color[1]), int(color[0])))

                #self.histmap = cv2.fillPoly(self.histmap, subbox, int(i*255/num) )
        elif side_box:
            direction_step = direction_step * scale

            p1 = p1 + width_vector * (box_num - box_id) * scale + width_vector * box_id             + width_vector*0.1
            p4 = p4 + width_vector * (box_id + 1) + width_vector * (box_num - box_id - 1) * scale   + width_vector*0.1

            for i in range(num):
                b1 = p1 + direction_step * i
                b2 = b1 + direction_step
                b4 = p4 + direction_step * i
                b3 = b4 + direction_step
                subbox = np.array((b1, b2, b3, b4), dtype=int)
                subbox = subbox.reshape((1, 4, 2))
                points = subbox.copy()
                subbox[:, :, 0], subbox[:, :, 1] = points[:, :, 1], points[:, :, 0]
                self.histmap = cv2.fillPoly(self.histmap, subbox, (int(hist[i]), 0, int(hist[i])))

    def __apply_mask(self, image, mask, alpha=0.5):
        """
        Apply the given mask to the image.
        :param image:
        :param mask:
        :param color: (0,0,1): blue
        :param alpha:
        :return:
        """
        for c in range(3):
            image[:, :, c] = image[:, :, c] * (1 - alpha) + mask[:, :, c] * alpha
        return image

    def rock_hist(self, boxes, size_max=4000, size_step=200, mode="size", side_box=False, statistics=False, colormap=False):
        """
        get the histograms of rocks that are in the boxes
        :param boxes:
        :return:
        """
        """ analysis"""
        rock_statistics = []
        for box in boxes:
            rock_box = dict()
            rock_box["box"] = box
            rocks = self.get_instances_in_box(box)
            rock_sizes = []
            orientations = []
            major_lengths = []
            for rock in rocks:
                bb = rock["bb"]
                mask = rock["mask"]
                top_left = bb[:2]
                mask = mask - top_left
                mask_shape = (bb[2] -bb[0], bb[3] - bb[1])
                mask = self.__create_bool_mask(mask, mask_shape)
                #_, contours, _ = cv2.findContours(mask.astype(np.uint8).copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)  # cv2.findContours has been changed in Ubuntu 18
                contours, _ = cv2.findContours(mask.astype(np.uint8).copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
                areas = [cv2.contourArea(cnt) for cnt in contours]
                rock_size = np.max(areas)
                rock_sizes.append(rock_size)

                i = np.argmax(areas)
                cnt = contours[i]
                ellipse = cv2.fitEllipse(cnt)
                (x, y), (MA, ma), angle = ellipse
                a = ma / 2
                b = MA / 2
                eccentricity = sqrt(pow(a, 2) - pow(b, 2))
                eccentricity = round(eccentricity / a, 2)
                orientations.append([angle, eccentricity])
                major_lengths.append(2*a)
            rock_box["sizes"] = rock_sizes
            rock_box["orientations"] = orientations
            rock_box["major_lengths"] = major_lengths
            rock_statistics.append(rock_box)

        """ histogram"""
        size_bins = np.arange(0, size_max, size_step)
        length_bins = np.arange(0, 200, 10)
        orientation_bins = np.arange(0, 180, 9)
        eccentricity_bins = np.arange(0, 1, 0.05)

        box_nm = len(rock_statistics)
        box = rock_statistics[0]["box"]
        width_v = box[0]-box[3]
        for i, rock_box in enumerate(rock_statistics):
            if mode == "size":
                sizes = rock_box["sizes"]
                if not colormap:
                    if not statistics:
                        hist, _ = np.histogram(sizes, size_bins)
                        hist = (np.array(hist) / np.max(hist) * 255).tolist()  # normalize histogram
                        box = rock_box["box"]
                        self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v)
                    elif statistics:
                        mu = np.array(sizes).mean()
                        sigma = sqrt(np.array(sizes).var())
                        brightness = 255 - sigma/10.0 # you need to print this out and then figure out the formula
                        if brightness<0:
                            brightness = 0
                        #print(mu)
                        hist, _ = np.histogram(mu, size_bins)
                        hist = (np.array(hist) / np.max(hist) * 255).tolist()  # normalize histogram
                        box = rock_box["box"]
                        self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v, brightness=brightness)
                else:
                    mu = np.array(sizes).mean()
                    sigma = sqrt(np.array(sizes).var())
                    hist, _ = np.histogram(mu, size_bins)
                    hist = hist * mu
                    box = rock_box["box"]
                    self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v, colormap_max=size_max)

            elif mode == "major_length":
                major_lengths = rock_box["major_lengths"]
                if not colormap:
                    if not statistics:
                        hist, _ = np.histogram(major_lengths, length_bins)
                        hist = (np.array(hist) / np.max(hist) * 255).tolist()  # normalize histogram
                        box = rock_box["box"]
                        self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v)
                    elif statistics:
                        mu = np.array(major_lengths).mean()
                        sigma = sqrt(np.array(major_lengths).var())
                        brightness = (50-sigma)/50.0*255
                        #print(brightness)
                        print(sigma)
                        if brightness < 0:
                            brightness = 0
                        hist, _ = np.histogram(mu, length_bins)  # hist is one-hot encoded
                        hist = (np.array(hist) / np.max(hist) * 255).tolist()  # normalize histogram
                        box = rock_box["box"]
                        self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v, brightness=brightness)
                else:
                    mu = np.array(major_lengths).mean()
                    sigma = sqrt(np.array(major_lengths).var())
                    hist, _ = np.histogram(mu, length_bins)
                    hist = hist * mu
                    box = rock_box["box"]
                    print(hist)
                    self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v, colormap_max=55)  # the maximum mean of major lengths

            elif mode == "orientation":
                orientations = rock_box["orientations"]

                orientations = np.array(orientations)[:,0]
                hist, _ = np.histogram(orientations, orientation_bins)
                hist = (np.array(hist) / np.max(hist) * 255).tolist()  # normalize histogram
                box = rock_box["box"]
                self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v)

            elif mode == "eccentricity":
                orientations = rock_box["orientations"]

                orientations = np.array(orientations)[:,1]
                hist, _ = np.histogram(orientations, eccentricity_bins)
                hist = (np.array(hist) / np.max(hist) * 255).tolist()  # normalize histogram
                box = rock_box["box"]
                self.__add_hist2colormap(hist, box, side_box, i, box_nm, width_v)
            else:
                print("Undefined mode!")

    def overlap_hist_orth(self, histmap):
        B = self.orth.GetRasterBand(1).ReadAsArray()
        G = self.orth.GetRasterBand(2).ReadAsArray()
        R = self.orth.GetRasterBand(3).ReadAsArray()
        orth = np.zeros((self.h, self.w, 3), dtype=np.uint8)
        orth[:, :, 0], orth[:, :, 1], orth[:, :, 2] = R, G, B
        del R, G, B
        overlap = self.__apply_mask(orth, histmap, alpha=0.8)
        return overlap

def hist_entire_scarp(sa):
    step = 30
    begin = 4146290
    sa.create_heatmap()

    for i in range(16):
        #intersection, norm_v = sa.get_normal_vector(begin - i*step, begin - i*step - step)
        intersection, norm_v = sa.get_normal_vector(begin - (15-i) * step, begin - (15-i) * step - step)

        boxes = sa.generate_boxes(intersection, norm_v, h=200, w=1600, step=200, num=4)
        #sa.rock_hist(boxes, mode="major_length")
        sa.rock_hist(boxes, mode="major_length", side_box=True)
        histmap = sa.histmap.astype(np.uint8).copy()
        cv2.imwrite("../contours/histmap.jpg", histmap)
        print(i)

    return sa.histmap.astype(np.uint8).copy()

if __name__ == "__main__":
    sa = Strike_analysis("../contours/resized_pruned_sk.jpg", "../C3.tif", "../registered_instances_v3.pickle")
    histmap = hist_entire_scarp(sa)
    #intersection, norm_v = sa.get_normal_vector(4146290, 4146260)
    """
    sk = sa.skeleton.copy()
    points = sa.generate_one_box(intersection, norm_v, 100, 30, 100)
    sk = sa.draw_box_in_skeleton(points, sk)
    cv2.imwrite("../box_sk.jpg", sk)
    """

    #boxes = sa.generate_boxes(intersection, norm_v, h=200, w=1600, step=200, num=4)
    """
    sk = sa.skeleton.copy()
    for box in boxes:
        rocks = sa.get_instances_in_box(box)
        sk = sa.draw_box_in_skeleton(box, sk)
        #sk = sa.draw_instances_in_skeleton(rocks, sk)
    #cv2.imwrite("../contours/box_sk.jpg", sk)
    """
    #sa.create_heatmap()
    #sa.rock_hist(boxes)
    #sa.rock_hist(boxes, mode="major_length", side_box=True)
    #histmap = sa.histmap.astype(np.uint8).copy()
    #cv2.imwrite("../contours/histmap.jpg", histmap)

    #del sa
    #sa = Strike_analysis("../contours/resized_pruned_sk.jpg", "../C3.tif", "../registered_instances_v3.pickle")
    #histmap = cv2.imread("../contours/histmap.jpg")
    #overlap = sa.overlap_hist_orth(histmap)
    #cv2.imwrite("../contours/overlap_hist_orth.jpg", overlap)


text:
"""!
@brief Pytorch dataloader for wham dataset for multiple gender combinations.

@author Efthymios Tzinis {etzinis2@illinois.edu}
@copyright University of illinois at Urbana Champaign
"""

import torch
import os
import numpy as np
import pickle
import glob2
import sys

current_dir = os.path.dirname(os.path.abspath('__file__'))
root_dir = os.path.abspath(os.path.join(current_dir, '../../'))
sys.path.append(root_dir)
import approx_ensembles.separation.dataset_loader.abstract_dataset as \
    abstract_dataset
from scipy.io import wavfile
import warnings
from tqdm import tqdm
from time import time

EPS = 1e-8
enh_single = {'mixture': 'mix_single',
              'sources': ['s1', 'noise'],
              'n_sources': 1}
enh_single_white_noise = {
              'mixture': 'source_with_white_noise',
              'sources': ['s1', 'white_noise'],
              'n_sources': 1}
enh_both = {'mixture': 'mix_both',
            'sources': ['mix_clean', 'noise'],
            'n_sources': 1}
sep_clean = {'mixture': 'mix_clean',
             'sources': ['s1', 's2'],
             'n_sources': 2}
sep_noisy = {'mixture': 'mix_both',
             'sources': ['s1', 's2', 'noise'],
             'n_sources': 2}

VALID_GENDER_COMBS = set(['ff', 'mm', 'fm', 'mf'])

WHAM_TASKS = {'enhance_single_white_noise': enh_single_white_noise,
              'enhance_single': enh_single,
              'enhance_both': enh_both,
              'sep_clean': sep_clean,
              'sep_noisy': sep_noisy}
WHAM_TASKS['enh_single'] = WHAM_TASKS['enhance_single']
WHAM_TASKS['enh_both'] = WHAM_TASKS['enhance_both']

def normalize_tensor_wav(wav_tensor, eps=1e-8, std=None):
    mean = wav_tensor.mean(-1, keepdim=True)
    if std is None:
        std = wav_tensor.std(-1, keepdim=True)
    return (wav_tensor - mean) / (std + eps)

class Dataset(torch.utils.data.Dataset, abstract_dataset.Dataset):
    """ Dataset class for WHAM source separation and speech enhancement tasks.

    Example of kwargs:
        root_dirpath='/mnt/data/wham', task='enh_single',
        split='tr', sample_rate=8000, timelength=4.0,
        normalize_audio=False, n_samples=0, zero_pad=False
    """
    def __init__(self, **kwargs):
        super(Dataset, self).__init__()
        warnings.filterwarnings("ignore")
        self.kwargs = kwargs

        self.task = self.get_arg_and_check_validness(
            'task', known_type=str, choices=WHAM_TASKS.keys())

        self.zero_pad = self.get_arg_and_check_validness(
            'zero_pad', known_type=bool)

        self.augment = self.get_arg_and_check_validness(
            'augment', known_type=bool)

        # Gender combination priors for combinations
        # ff, mm, fm/mf
        self.gender_combination_priors = self.get_arg_and_check_validness(
            'gender_combination_priors', known_type=float,
            dict_check={'ff': float, 'mm': float, 'fm': float, 'mf':float},
            extra_lambda_checks=[lambda x: [0 <= y <= 1 for y in x.values()]])

        self.normalize_audio = self.get_arg_and_check_validness(
            'normalize_audio', known_type=bool)

        self.min_or_max = self.get_arg_and_check_validness(
            'min_or_max', known_type=str, choices=['min', 'max'])

        self.split = self.get_arg_and_check_validness(
            'split', known_type=str, choices=['cv', 'tr', 'tt'])

        self.n_samples = self.get_arg_and_check_validness(
            'n_samples', known_type=int, extra_lambda_checks=[lambda x: x >= 0])

        self.sample_rate = self.get_arg_and_check_validness('sample_rate',
                                                            known_type=int)
        self.root_path = self.get_arg_and_check_validness(
            'root_dirpath', known_type=str,
            extra_lambda_checks=[lambda y: os.path.lexists(y)])
        self.dataset_dirpath = self.get_path()

        self.mixtures_info_metadata_path = os.path.join(
            self.dataset_dirpath, 'metadata_v2')

        self.timelength = self.get_arg_and_check_validness(
            'timelength', known_type=float)

        self.time_samples = int(self.sample_rate * self.timelength)

        # Create the indexing for the dataset
        mix_folder_path = os.path.join(self.dataset_dirpath,
                                       WHAM_TASKS[self.task]['mixture'])
        self.file_names = []
        self.available_mixtures = glob2.glob(mix_folder_path + '/*.wav')

        self.mixtures_info = []
        print('Parsing Dataset found at: {}...'.format(self.dataset_dirpath))
        if not os.path.lexists(self.mixtures_info_metadata_path):
            # Parse gender information.
            gender_info_path = os.path.join(
                os.path.dirname(os.path.abspath('__file__')),
                'wham_speaker_info.txt')
            gender_dic = {}
            if os.path.lexists(gender_info_path):
                with open(gender_info_path, 'rb') as filehandle:
                    gender_dic = dict([tuple([x.decode() for x in l.split()])
                                       for l in filehandle.readlines()])

            for file_path in tqdm(self.available_mixtures):
                sample_rate, waveform = wavfile.read(file_path)
                assert sample_rate == self.sample_rate
                numpy_wav = np.array(waveform)

                speaker_info = os.path.basename(file_path).split('.wav')[0]
                speaker_info = [x[:3] for x in speaker_info.split('_')[::2]]

                this_gender_comb = ''
                for speaker in speaker_info:
                    if speaker not in gender_dic:
                        raise ValueError('Speaker with id: {} not '
                                         'found!'.format(speaker))
                    else:
                        this_gender_comb += gender_dic[speaker].lower()

                self.mixtures_info.append([os.path.basename(file_path),
                                           numpy_wav.shape[0],
                                           this_gender_comb])

            print('Dumping metadata in: {}'.format(
                self.mixtures_info_metadata_path))
            with open(self.mixtures_info_metadata_path, 'wb') as filehandle:
                pickle.dump(self.mixtures_info, filehandle)

        if os.path.lexists(self.mixtures_info_metadata_path):
            with open(self.mixtures_info_metadata_path, 'rb') as filehandle:
                self.mixtures_info = pickle.load(filehandle)
                print('Loaded metadata from: {}'.format(
                    self.mixtures_info_metadata_path))

        self.file_names_g_comb = dict([(g, []) for g in VALID_GENDER_COMBS])
        for path, n_samples, gender_comb in self.mixtures_info:
            if (n_samples >= self.time_samples or self.zero_pad):
                self.file_names_g_comb[gender_comb].append((path, n_samples))

        self.file_names = []
        # Apply the priors
        for gender_comb in self.file_names_g_comb:
            percentage = self.gender_combination_priors[gender_comb]
            length = len(self.file_names_g_comb[gender_comb])
            n_requested = int(length * percentage)
            self.file_names += self.file_names_g_comb[gender_comb][:n_requested]
        if self.n_samples > 0:
            self.file_names = self.file_names[:self.n_samples]

        max_time_samples = max([n_s for (_, n_s) in self.file_names])
        self.file_names = [x for (x, _) in self.file_names]

        # for the case that we need the whole audio input
        if self.time_samples <= 0.:
            self.time_samples = max_time_samples

    def get_path(self):
        path = os.path.join(self.root_path,
                            'wav{}k'.format(int(self.sample_rate / 1000)),
                            self.min_or_max, self.split)
        if os.path.lexists(path):
            return path
        else:
            raise IOError('Dataset path: {} not found!'.format(path))

    def safe_pad(self, tensor_wav):
        if self.zero_pad and tensor_wav.shape[0] < self.time_samples:
            appropriate_shape = tensor_wav.shape
            padded_wav = torch.zeros(
                list(appropriate_shape[:-1]) + [self.time_samples],
                dtype=torch.float32)
            padded_wav[:tensor_wav.shape[0]] = tensor_wav
            return padded_wav[:self.time_samples]
        else:
            return tensor_wav[:self.time_samples]

    def __len__(self):
        return len(self.file_names)

    def __getitem__(self, idx):
        if self.augment:
            the_time = int(np.modf(time())[0] * 100000000)
            np.random.seed(the_time)

        filename = self.file_names[idx]

        mixture_path = os.path.join(self.dataset_dirpath,
                                    WHAM_TASKS[self.task]['mixture'],
                                    filename)
        _, waveform = wavfile.read(mixture_path)
        max_len = len(waveform)
        rand_start = 0
        if self.augment and max_len > self.time_samples:
            rand_start = np.random.randint(0, max_len - self.time_samples)
            waveform = waveform[rand_start:rand_start+self.time_samples]
        mixture_wav = np.array(waveform)
        mixture_wav = torch.tensor(mixture_wav, dtype=torch.float32)
        # First normalize the mixture and then pad
        if self.normalize_audio:
            mixture_wav = normalize_tensor_wav(mixture_wav)
        mixture_wav = self.safe_pad(mixture_wav)

        sources_list = []
        for source_name in WHAM_TASKS[self.task]['sources']:
            source_path = os.path.join(self.dataset_dirpath,
                                       source_name, filename)
            try:
                _, waveform = wavfile.read(source_path)
            except Exception as e:
                print(e)
                raise IOError('could not load file from: {}'.format(source_path))
            waveform = waveform[rand_start:rand_start + self.time_samples]
            numpy_wav = np.array(waveform)
            source_wav = torch.tensor(numpy_wav, dtype=torch.float32)
            # First normalize the mixture and then pad
            if self.normalize_audio:
                source_wav = normalize_tensor_wav(source_wav)
            source_wav = self.safe_pad(source_wav)
            sources_list.append(source_wav)

        if self.normalize_audio:
            mix_std = mixture_wav.detach().cpu().numpy().std()
            mixture_wav = normalize_tensor_wav(mixture_wav, std=mix_std)
            sources_list = [normalize_tensor_wav(s, std=mix_std)
                            for s in sources_list]
        sources_wavs = torch.stack(sources_list, dim=0)

        return mixture_wav, sources_wavs

    def get_generator(self, batch_size=4, shuffle=True, num_workers=4):
        generator_params = {'batch_size': batch_size,
                            'shuffle': shuffle,
                            'num_workers': num_workers,
                            'drop_last': True}
        return torch.utils.data.DataLoader(self, **generator_params,
                                           pin_memory=True)

def test_generator():
    wham_root_p = '/mnt/data/wham'
    batch_size = 1
    sample_rate = 8000
    timelength = 4.0
    gender_combination_priors = {
        'ff': 0., 'mm': 0.05, 'fm': 0., 'mf': 0.02
    }
    time_samples = int(sample_rate * timelength)
    data_loader = Dataset(
        root_dirpath=wham_root_p, task='sep_clean',
        gender_combination_priors=gender_combination_priors,
        split='tt', sample_rate=sample_rate, timelength=timelength,
        zero_pad=True, min_or_max='min', augment=True,
        normalize_audio=False, n_samples=10)
    generator = data_loader.get_generator(batch_size=batch_size, num_workers=1)

    for mixture, sources in generator:
        assert mixture.shape == (batch_size, time_samples)
        assert sources.shape == (batch_size, 2, time_samples)

    # test the testing set with batch size 1 only
    data_loader = Dataset(
        root_dirpath=wham_root_p, task='sep_clean',
        gender_combination_priors=gender_combination_priors,
        split='tt', sample_rate=sample_rate, timelength=-1.,
        zero_pad=False, min_or_max='min', augment=False,
        normalize_audio=False, n_samples=10)
    generator = data_loader.get_generator(batch_size=1, num_workers=1)

    for mixture, sources in generator:
        assert mixture.shape[-1] == sources.shape[-1]

if __name__ == "__main__":
    test_generator()

reason: import numpy,from scipy


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

import os
import argparse

import joblib
import matplotlib.pyplot as plt
import numpy as np
from keras_preprocessing.sequence import pad_sequences

from gen_train_captions import EOS_TOKEN, SOS_TOKEN
from train import get_train_captions, max_length
from model import ImgCapModel

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('--dataset_root', type=str, required=True, help='Path to dataset root directory (flickr8k or COCO)')
    parser.add_argument('--img_root', type=str, required=True, help='Path to image root directory')
    parser.add_argument('--model_path', type=str, required=True, help='Path to model.h5')
    return parser.parse_args()

if __name__ == "__main__":
    args = parse_args()

    with open(os.path.join(args.dataset_root, 'test_image_feats.pkl'), 'rb') as f:
        test_img_feats = joblib.load(f)
    test_imgs = list(test_img_feats.keys())
    target_img = np.random.choice(test_imgs, 1)[0]
    target_img_feat = test_img_feats[target_img].reshape(1, 2048)

    with open(os.path.join(args.dataset_root, 'word2idx.pkl'), 'rb') as f:
        word2idx = joblib.load(f)
    with open(os.path.join(args.dataset_root, 'idx2word.pkl'), 'rb') as f:
        idx2word = joblib.load(f)

    # Load trained model
    img_cap_model = ImgCapModel(vocab_size=len(word2idx.keys()) + 1)
    img_cap_model.model.load_weights(args.model_path)
    greedy_text = img_cap_model.greedy_search(target_img_feat, word2idx, idx2word)
    beam_text = img_cap_model.beam_search(target_img_feat, word2idx, idx2word, beam_width=3)

    if args.dataset_root == 'COCO':
        target_img = 'COCO_val2014_' + '%012d.jpg' % (target_img)
    x = plt.imread(os.path.join(args.img_root, target_img))
    plt.imshow(x)
    plt.title('\n'.join(['Greedy: ' + greedy_text, 'Beam:' + beam_text]))
    plt.show()

max_stars_repo_name: satojkovic/DeepCaptioning