text:
import torch
import cv2
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from skimage import transform

torch.manual_seed(0)

def h_point_to_line_intersection(point, line):
    p1, p2 = line
    x1, y1 = p1
    x2, y2 = p2
    x3, y3 = point

    m = (y2 - y1) / (x2 - x1)
    b = y2 - m * x2

    iy = y3
    ix = (y3 - b) / m

    return ix, iy

def line_to_coors(line):
    x_coords = [i[0] for i in line]
    y_coords = [i[1] for i in line]
    return x_coords, y_coords

def eye_mouth_diag_angles(processed_landmarks):

    lm, rm = (processed_landmarks["mouth"][48], processed_landmarks["mouth"][54])

    le, re = (
        list(processed_landmarks["left_eye_center"].values())[0],
        list(processed_landmarks["right_eye_center"].values())[0],
    )

    eye_mouth_diag_angles = np.arccos(
        ((le - rm) * (re - lm))
        / (np.linalg.norm((le - rm)) * np.linalg.norm((re - lm)))
    )
    eye_mouth_diag_angles = np.degrees(eye_mouth_diag_angles)

    return eye_mouth_diag_angles

def nose_mouth_features(processed_landmarks):
    # TODO: Normalize distances by face max x distance and max y distance

    # Get landmarks
    x1, y1 = processed_landmarks["nose"][30]
    lm, rm = (processed_landmarks["mouth"][48], processed_landmarks["mouth"][54])

    x2, y2 = lm
    x3, y3 = rm

    # nose to mouth left corner
    lmnx = x1 - x2
    lmny = y2 - y1
    lmnd = np.sqrt(lmnx ** 2 + lmny ** 2)
    theta = np.degrees(np.arctan(lmny / lmnx))

    # nose to mouth right corner
    rmnx = x3 - x1
    rmny = y3 - y1
    rmnd = np.sqrt(rmnx ** 2 + rmny ** 2)
    phi = np.degrees(np.arctan(rmny / rmnx))
    return lmnx, lmny, lmnd, theta, rmnx, rmny, rmnd, phi

def line_horizontal_deviation(line):

    p1, p2 = line
    x1, y1 = p1
    x2, y2 = p2

    dev = np.degrees(np.arctan(np.abs((y2 - y1) / (x2 - x1))))

    return dev

def h_eyebrow_eyebrow_deviation(processed_landmarks, exclude=[]):

    """
    landmarks indices must be in ascending order from left to right
    and right eyebrow indices larger than left eyebrow
    """

    leb_filtered = {
        k: v for k, v in processed_landmarks["left_eyebrow"].items() if k not in exclude
    }
    reb_filtered = {
        k: v
        for k, v in processed_landmarks["right_eyebrow"].items()
        if k not in exclude
    }

    leb_filtered_sorted_idxs = np.argsort(list(leb_filtered.keys()))
    reb_filtered_sorted_idxs = np.argsort(list(reb_filtered.keys()))[::-1]

    leb_lmks = np.array(list(leb_filtered.values()))[leb_filtered_sorted_idxs]
    reb_lmks = np.array(list(reb_filtered.values()))[reb_filtered_sorted_idxs]

    deviations = []
    for leb, reb in zip(leb_lmks, reb_lmks):

        dev = line_horizontal_deviation([(leb), (reb)])
        deviations.append(dev)

    return deviations

# TODO: include mouth corners horizontal deviation?

class ROIExtractor:
    def __init__(self, frame, processed_landmarks, debug=False):
        self.frame = frame
        self.processed_landmarks = processed_landmarks
        self.debug = debug

    def get_region_of_interest(self, point1, point2, project_to=[]):

        """
        project_to = array of bbox vertices ordered starting at top left
        and continue clockwise
        """

        min_x, max_x = min(point1[0], point2[0]), max(point1[0], point2[0])
        min_y, max_y = min(point1[1], point2[1]), max(point1[1], point2[1])
        roi_bounds = (min_x, max_x, min_y, max_y)

        if project_to != []:

            src = np.array(
                [[min_x, min_y], [max_x, min_y], [max_x, max_y], [min_x, max_y]]
            )

            dst = np.array(project_to)
            dw, dh = dst.max(axis=0) - dst.min(axis=0)

            # dst_2 = dst_2*2 #because image sizes are not the same.
            tform = transform.estimate_transform("projective", src, dst)
            self.frame = transform.warp(self.frame, tform.inverse)

            roi = self.frame[min_y:max_y, min_x:max_x, :]
            roi = cv2.resize(roi, (dh, dw))

        else:
            roi = self.frame[min_y:max_y, min_x:max_x, :]

        if self.debug:
            plt.imshow(roi)
            plt.show()

        return roi, roi_bounds

    def left_mouth(self, **kwargs):
        tr = self.processed_landmarks["nose"][33]
        bl = np.array(
            (
                self.processed_landmarks["head"][2][0],
                self.processed_landmarks["head"][8][1],
            )
        )
        return self.get_region_of_interest(tr, bl, **kwargs)

    def right_mouth(self, **kwargs):
        p1 = self.processed_landmarks["nose"][33]
        p2 = np.array(
            (
                self.processed_landmarks["head"][14][0],
                self.processed_landmarks["head"][8][1],
            )
        )

        return self.get_region_of_interest(p1, p2, **kwargs)

class BarycentricOperations:
    def __init__(self, landmarks_left_to_right_mapping, triangle_indices, constant=100):

        self.landmarks_left_to_right_mapping = dict(landmarks_left_to_right_mapping)
        self.triangle_indices = (
            triangle_indices  # for dlip landmarks: triangle_indices = [17,26,30]
        )

        self.constant = constant

    def order_landmarks(self, landmarks):
        ordered_idxs = list(self.landmarks_left_to_right_mapping.keys()) + list(
            self.landmarks_left_to_right_mapping.values()
        )

        # Place triangle indices to the end
        for i in set(ordered_idxs) & set(self.triangle_indices):
            ordered_idxs.pop(ordered_idxs.index(i))

        ordered_idxs += self.triangle_indices

        return landmarks[ordered_idxs, :]

    def cartesian_to_barycentric(self, coordinates):

        """
        the last 3 coordinates of the array must be the reference triangle coordinates.

        """

        m, n = coordinates.shape

        homogeneous_coordinates = np.hstack(
            [coordinates, np.ones((m, 1)) * self.constant]
        )

        self.triangle = homogeneous_coordinates[-3:, :]

        barycentric_coordinates = homogeneous_coordinates[:-3, :] @ np.linalg.inv(
            self.triangle
        )

        return barycentric_coordinates

    def get_barycentric_landmarks(self, landmarks, order=True):

        # print(landmarks.shape)
        if order:
            ordered_landmarks = self.order_landmarks(landmarks)

        # print(ordered_landmarks.shape)
        return self.cartesian_to_barycentric(ordered_landmarks)

    def landmarks_true_horizontal_flip(self, ordered_landmarks):
        m, n = ordered_landmarks.shape
        c = m // 2
        flip_matrix = np.block(
            [[np.zeros((c, c)), np.eye(c, c)], [np.eye(c, c), np.zeros((c, c))]]
        )

        # print(flip_matrix.shape,ordered_landmarks.shape)
        ordered_landmarks_real_flip = flip_matrix @ ordered_landmarks

        return ordered_landmarks_real_flip

    def landmarks_theoretical_horizontal_flip(self, ordered_landmarks):
        return ordered_landmarks @ np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]])

    def global_asymmetry_index(self, landmarks):

        barycentric_landmarks = self.get_barycentric_landmarks(landmarks)

        true_flipped_landmarks = self.landmarks_true_horizontal_flip(
            barycentric_landmarks
        )
        theoretical_flipped_landmarks = self.landmarks_theoretical_horizontal_flip(
            barycentric_landmarks
        )

        asymmetry_by_landmark = true_flipped_landmarks - theoretical_flipped_landmarks
        return (asymmetry_by_landmark ** 2).sum()

class OpticalFlow:
    def __init__(self, flow_threshold=6):
        self.flow_threshold = flow_threshold
        self.mask = None

    def compute_dense_flow(self, prev_frame, frame):

        prev_gray = cv2.cvtColor(prev_frame, cv2.COLOR_BGR2GRAY)
        # Converts each frame to grayscale - we previously
        # only converted the first frame to grayscale
        gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

        # Calculates dense optical flow by Farneback method
        flow = cv2.calcOpticalFlowFarneback(
            prev_gray, gray, None, 0.5, 3, 15, 3, 5, 1.2, 0
        )

        # Computes the magnitude and angle of the 2D vectors
        magnitude, angle = cv2.cartToPolar(flow[..., 0], flow[..., 1])

        if self.mask is None:
            self.mask = np.zeros_like(prev_frame)

            # Sets image saturation to maximum
            self.mask[..., 1] = 255

        return magnitude, angle

    def filter_flow_magnitude_noise(self, magnitude):

        # Filter noise in magnitude.
        mean_flow_mag = magnitude.mean()
        magnitude[
            magnitude < self.flow_threshold * magnitude.mean()
        ] = 0  # TODO: This should be a hparam

        return magnitude

    def get_rgb_optical_flow(self, magnitude, angle):

        # Sets image hue according to the optical flow
        # direction
        self.mask[..., 0] = angle * 180 / np.pi / 2

        # Sets image value according to the optical flow
        # magnitude (normalized)
        self.mask[..., 2] = cv2.normalize(magnitude, None, 0, 255, cv2.NORM_MINMAX)

        # Converts HSV to RGB (BGR) color representation
        rgb = cv2.cvtColor(self.mask, cv2.COLOR_HSV2BGR)

        self.mask = None

        return rgb

class OpticalFlowSymmetry:
    def __init__(self):
        pass

    def movement_score(self, magnitudes: list):

        magnitudes = np.array(magnitudes)
        return magnitudes.mean()

    def __symmetry_score(
        self, left_movement_score, right_movement_score, lambda_=3.8, ret="abs_diff"
    ):

        abs_dif = abs(left_movement_score - right_movement_score)

        raw_symmetry_score = 1 - lambda_ * abs_dif

        symmetry_score = np.clip(raw_symmetry_score, 0, 1)

        if ret == "abs_diff":
            return abs_dif
        elif ret == "raw_score":
            return raw_symmetry_score
        elif ret == "score":
            return symmetry_score
        else:
            raise ValueError("return type does not exit")

    def symmetry_score(
        self, left_magnitudes, right_magnitudes, lambda_=3.8, ret="abs_diff"
    ):

        left_movement_score = self.movement_score(left_magnitudes)
        right_movement_score = self.movement_score(right_magnitudes)

        return self.__symmetry_score(
            left_movement_score, right_movement_score, lambda_=lambda_, ret=ret
        )

ext: py

reason: import numpy


text:
#Make necessary import
import numpy as np
import pandas as pd
import itertools
from sklearn.model_selection import train_test_split
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.linear_model import PassiveAggressiveClassifier
from sklearn.metrics import accuracy_score, confusion_matrix

#Read the data

#Now, let’s read the data into a DataFrame, and get the shape of the data and the first 5 records.
data_frame=pd.read_csv('data/news.csv')

#Get shape
data_frame.shape

#Get shape and head
print(data_frame.columns)

#after ru in spyder console run data_frame.head()

#Get the labels

#And get the labels from the DataFrame.
labels = data_frame.label
labels.head()

#Split the dataset

#Split the dataset into training and testing sets.
x_train,x_test,y_train,y_test=train_test_split(data_frame['text'], labels, test_size=0.2, random_state=7)

#Initialize a TfidfVectorizer

"""Let’s initialize a TfidfVectorizer with stop words from the English language and a maximum document
frequency of 0.7 (terms with a higher document frequency will be discarded). Stop words are the most common
words in a language that are to be filtered out before processing the natural language data. And a TfidfVectorizer
turns a collection of raw documents into a matrix of TF-IDF features."""
tfidf_vectorizer=TfidfVectorizer(stop_words='english', max_df=0.7)

#Fit and transform train set, transform test set
#Now, fit and transform the vectorizer on the train set, and transform the vectorizer on the test set.
tfidf_train=tfidf_vectorizer.fit_transform(x_train)
tfidf_test=tfidf_vectorizer.transform(x_test)

#Initialize a PassiveAggressiveClassifier
#Next, we’ll initialize a PassiveAggressiveClassifier. This is. We’ll fit this on tfidf_train and y_train.
pac=PassiveAggressiveClassifier(max_iter=50)
pac.fit(tfidf_train,y_train)

#Predict on the test set and calculate accuracy
"""Then, we’ll predict on the test set from the TfidfVectorizer and calculate the accuracy
with accuracy_score() from sklearn.metrics."""
y_pred=pac.predict(tfidf_test)
score=accuracy_score(y_test,y_pred)
print(f'Accuracy: {round(score*100,2)}%')

#Output: Accuracy: 93.13%

#DBuild confusion matrix
"""We got an accuracy of 92.82% with this model. Finally, let’s print out a confusion matrix to gain insight
into the number of false and true negatives and positives."""
confusion_matrix(y_test,y_pred, labels=['FAKE','REAL'])

"""output:
    array([[590,  48],
          [ 39, 590]])"""
#So with this model, we have 590 true positives, 590 true negatives, 39 false positives, and 48 false negatives.


text:
import cv2
import sys
import os
import pickle
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime

try:
    # 拡張モジュールのインポート
    from c_mandel import get_mandel_array
except ModuleNotFoundError as e:
    print(e)
    s = (
        'Run the following command:\n'
        '   \"python setup.py build_ext -i\"'
    )
    print(s)
    sys.exit(1)

DUMP_FILE_NAME = 'param_dict.pickle'    # パラメータを保存するファイル

if __name__ == '__main__':
    W = int(sys.argv[2])    # 分割数
    LIMIT = 4.0             # 計算停止閾値
    INF = 300               # 繰り返し上限

    try:
        f = open(DUMP_FILE_NAME, 'rb')
        param_dict = pickle.load(f)
        X0 = param_dict['X0']   # xの最小値
        X1 = param_dict['X1']   # xの最大値
        Y0 = param_dict['Y0']   # yの最小値
        Y1 = param_dict['Y1']   # yの最大値
        f.close()
    except FileNotFoundError:
        param_dict = {}
        X0 = param_dict['X0'] = -2.25
        X1 = param_dict['X1'] = 0.75
        Y0 = param_dict['Y0'] = -1.50
        Y1 = param_dict['Y1'] = 1.50

    with open(DUMP_FILE_NAME, 'wb') as f:
        param_dict['W'] = W
        pickle.dump(param_dict, f)
        f.close()

    # パラメータの表示
    s = (
        'X0 = {:12.9f}  Y0 = {:12.9f}\n'
        'X1 = {:12.9f}  Y1 = {:12.9f}'
    ).format(X0, Y0, X1, Y1)
    print(s)

    # get_mandel_array -> np.ndarray
    img = get_mandel_array(X0, X1, Y0, Y1, W, INF, LIMIT)

    if sys.argv[1] == 'png':
        dt = datetime.now().strftime('%y%m%d%H%M%S')
        file_dir = './mandel_img'
        filename = file_dir + '/mandelbrot_{}.png'.format(dt)
        if not os.path.exists(file_dir):
            os.mkdir(file_dir)
        cv2.imwrite(filename, img)
        print(filename)
    elif sys.argv[1] == 'show':
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        plt.imshow(img)
        plt.show()

max_stars_repo_name: s-col/cymandel

max_issues_repo_name: s-col/cymandel

max_forks_repo_name: s-col/cymandel


text:
# -*- coding: utf-8 -*-
"""
Clustering Methods
"""
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import normalize

class ClusteringMethods:
    """ Base class of clustering methods """

    @staticmethod
    def _normalize_values(arr, norm=None, axis=None):
        """
        Normalize values in array by column

        Parameters
        ----------
        arr : ndarray
            ndarray of values extracted from meta
            shape (n samples, with m features)
        norm : str
            Normalization method to use (see sklearn.preprocessing.normalize)
            if None range normalize
        axis : int
            Axis to normalize along

        Returns
        ---------
        arr : ndarray
            array with values normalized by column
            shape (n samples, with m features)
        """
        if norm:
            arr = normalize(arr, norm=norm, axis=axis)
        else:
            if np.issubdtype(arr.dtype, np.integer):
                arr = arr.astype(float)

            min_all = arr.min(axis=axis)
            max_all = arr.max(axis=axis)
            range_all = max_all - min_all
            if axis is not None:
                pos = range_all == 0
                range_all[pos] = 1

            arr -= min_all
            arr /= range_all

        return arr

    @staticmethod
    def kmeans(data, **kwargs):
        """ Cluster based on kmeans methodology """

        kmeans = KMeans(random_state=0, **kwargs)
        results = kmeans.fit(data)
        labels = results.labels_

        # Create deterministic cluster labels based on size
        label_n, l_size = np.unique(labels, return_counts=True)
        idx = np.argsort(l_size)
        l_mapping = dict(zip(label_n[idx], label_n))

        sorted_labels = labels.copy()
        for k, v in l_mapping.items():
            sorted_labels[labels == k] = v

        return sorted_labels

max_stars_repo_name: NREL/reVX

max_issues_repo_name: NREL/reVX

max_forks_repo_name: NREL/reVX


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

import itertools

import numpy as np
import pytest
import scipy.linalg

import cirq

I = np.eye(2)
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])
Z = np.array([[1, 0], [0, -1]])
H = np.array([[1, 1], [1, -1]]) * np.sqrt(0.5)
SQRT_X = np.array([[np.sqrt(1j), np.sqrt(-1j)], [np.sqrt(-1j), np.sqrt(1j)]]) * np.sqrt(0.5)
SQRT_Y = np.array([[np.sqrt(1j), -np.sqrt(1j)], [np.sqrt(1j), np.sqrt(1j)]]) * np.sqrt(0.5)
SQRT_Z = np.diag([1, 1j])
E00 = np.diag([1, 0])
E01 = np.array([[0, 1], [0, 0]])
E10 = np.array([[0, 0], [1, 0]])
E11 = np.diag([0, 1])
PAULI_BASIS = cirq.PAULI_BASIS
STANDARD_BASIS = {'a': E00, 'b': E01, 'c': E10, 'd': E11}

def _one_hot_matrix(size: int, i: int, j: int) -> np.ndarray:
    result = np.zeros((size, size))
    result[i, j] = 1
    return result

@pytest.mark.parametrize(
    'basis1, basis2, expected_kron_basis',
    (
        (
            PAULI_BASIS,
            PAULI_BASIS,
            {
                'II': np.eye(4),
                'IX': scipy.linalg.block_diag(X, X),
                'IY': scipy.linalg.block_diag(Y, Y),
                'IZ': np.diag([1, -1, 1, -1]),
                'XI': np.array([[0, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0]]),
                'XX': np.rot90(np.eye(4)),
                'XY': np.rot90(np.diag([1j, -1j, 1j, -1j])),
                'XZ': np.array([[0, 0, 1, 0], [0, 0, 0, -1], [1, 0, 0, 0], [0, -1, 0, 0]]),
                'YI': np.array([[0, 0, -1j, 0], [0, 0, 0, -1j], [1j, 0, 0, 0], [0, 1j, 0, 0]]),
                'YX': np.rot90(np.diag([1j, 1j, -1j, -1j])),
                'YY': np.rot90(np.diag([-1, 1, 1, -1])),
                'YZ': np.array([[0, 0, -1j, 0], [0, 0, 0, 1j], [1j, 0, 0, 0], [0, -1j, 0, 0]]),
                'ZI': np.diag([1, 1, -1, -1]),
                'ZX': scipy.linalg.block_diag(X, -X),
                'ZY': scipy.linalg.block_diag(Y, -Y),
                'ZZ': np.diag([1, -1, -1, 1]),
            },
        ),
        (
            STANDARD_BASIS,
            STANDARD_BASIS,
            {
                'abcd'[2 * row_outer + col_outer]
                + 'abcd'[2 * row_inner + col_inner]: _one_hot_matrix(
                    4, 2 * row_outer + row_inner, 2 * col_outer + col_inner
                )
                for row_outer in range(2)
                for row_inner in range(2)
                for col_outer in range(2)
                for col_inner in range(2)
            },
        ),
    ),
)
def test_kron_bases(basis1, basis2, expected_kron_basis):
    kron_basis = cirq.kron_bases(basis1, basis2)
    assert len(kron_basis) == 16
    assert set(kron_basis.keys()) == set(expected_kron_basis.keys())
    for name in kron_basis.keys():
        assert np.all(kron_basis[name] == expected_kron_basis[name])

@pytest.mark.parametrize(
    'basis1,basis2',
    (
        (PAULI_BASIS, cirq.kron_bases(PAULI_BASIS)),
        (STANDARD_BASIS, cirq.kron_bases(STANDARD_BASIS, repeat=1)),
        (cirq.kron_bases(PAULI_BASIS, PAULI_BASIS), cirq.kron_bases(PAULI_BASIS, repeat=2)),
        (
            cirq.kron_bases(
                cirq.kron_bases(PAULI_BASIS, repeat=2),
                cirq.kron_bases(PAULI_BASIS, repeat=3),
                PAULI_BASIS,
            ),
            cirq.kron_bases(PAULI_BASIS, repeat=6),
        ),
        (
            cirq.kron_bases(
                cirq.kron_bases(PAULI_BASIS, STANDARD_BASIS),
                cirq.kron_bases(PAULI_BASIS, STANDARD_BASIS),
            ),
            cirq.kron_bases(PAULI_BASIS, STANDARD_BASIS, repeat=2),
        ),
    ),
)
def test_kron_bases_consistency(basis1, basis2):
    assert set(basis1.keys()) == set(basis2.keys())
    for name in basis1.keys():
        assert np.all(basis1[name] == basis2[name])

@pytest.mark.parametrize(
    'basis,repeat', itertools.product((PAULI_BASIS, STANDARD_BASIS), range(1, 5))
)
def test_kron_bases_repeat_sanity_checks(basis, repeat):
    product_basis = cirq.kron_bases(basis, repeat=repeat)
    assert len(product_basis) == 4 ** repeat
    for name1, matrix1 in product_basis.items():
        for name2, matrix2 in product_basis.items():
            p = cirq.hilbert_schmidt_inner_product(matrix1, matrix2)
            if name1 != name2:
                assert p == 0
            else:
                assert abs(p) >= 1

@pytest.mark.parametrize(
    'm1,m2,expect_real',
    (
        (X, X, True),
        (X, Y, True),
        (X, H, True),
        (X, SQRT_X, False),
        (I, SQRT_Z, False),
    ),
)
def test_hilbert_schmidt_inner_product_is_conjugate_symmetric(m1, m2, expect_real):
    v1 = cirq.hilbert_schmidt_inner_product(m1, m2)
    v2 = cirq.hilbert_schmidt_inner_product(m2, m1)
    assert v1 == v2.conjugate()

    assert np.isreal(v1) == expect_real
    if not expect_real:
        assert v1 != v2

@pytest.mark.parametrize(
    'a,m1,b,m2',
    (
        (1, X, 1, Z),
        (2, X, 3, Y),
        (2j, X, 3, I),
        (2, X, 3, X),
    ),
)
def test_hilbert_schmidt_inner_product_is_linear(a, m1, b, m2):
    v1 = cirq.hilbert_schmidt_inner_product(H, (a * m1 + b * m2))
    v2 = a * cirq.hilbert_schmidt_inner_product(H, m1) + b * cirq.hilbert_schmidt_inner_product(
        H, m2
    )
    assert v1 == v2

@pytest.mark.parametrize('m', (I, X, Y, Z, H, SQRT_X, SQRT_Y, SQRT_Z))
def test_hilbert_schmidt_inner_product_is_positive_definite(m):
    v = cirq.hilbert_schmidt_inner_product(m, m)
    # Cannot check using np.is_real due to bug in aarch64.
    # See https://github.com/quantumlib/Cirq/issues/4379
    assert np.isclose(np.imag(v), 1e-16)
    assert v.real > 0

@pytest.mark.parametrize(
    'm1,m2,expected_value',
    (
        (X, I, 0),
        (X, X, 2),
        (X, Y, 0),
        (X, Z, 0),
        (H, X, np.sqrt(2)),
        (H, Y, 0),
        (H, Z, np.sqrt(2)),
        (Z, E00, 1),
        (Z, E01, 0),
        (Z, E10, 0),
        (Z, E11, -1),
        (SQRT_X, E00, np.sqrt(-0.5j)),
        (SQRT_X, E01, np.sqrt(0.5j)),
        (SQRT_X, E10, np.sqrt(0.5j)),
        (SQRT_X, E11, np.sqrt(-0.5j)),
    ),
)
def test_hilbert_schmidt_inner_product_values(m1, m2, expected_value):
    v = cirq.hilbert_schmidt_inner_product(m1, m2)
    assert np.isclose(v, expected_value)

@pytest.mark.parametrize(
    'm,basis',
    itertools.product(
        (I, X, Y, Z, H, SQRT_X, SQRT_Y, SQRT_Z),
        (PAULI_BASIS, STANDARD_BASIS),
    ),
)
def test_expand_matrix_in_orthogonal_basis(m, basis):
    expansion = cirq.expand_matrix_in_orthogonal_basis(m, basis)

    reconstructed = np.zeros(m.shape, dtype=complex)
    for name, coefficient in expansion.items():
        reconstructed += coefficient * basis[name]
    assert np.allclose(m, reconstructed)

@pytest.mark.parametrize(
    'expansion',
    (
        {'I': 1},
        {'X': 1},
        {'Y': 1},
        {'Z': 1},
        {'X': 1, 'Z': 1},
        {'I': 0.5, 'X': 0.4, 'Y': 0.3, 'Z': 0.2},
        {'I': 1, 'X': 2, 'Y': 3, 'Z': 4},
    ),
)
def test_matrix_from_basis_coefficients(expansion):
    m = cirq.matrix_from_basis_coefficients(expansion, PAULI_BASIS)

    for name, coefficient in expansion.items():
        element = PAULI_BASIS[name]
        expected_coefficient = cirq.hilbert_schmidt_inner_product(
            m, element
        ) / cirq.hilbert_schmidt_inner_product(element, element)
        assert np.isclose(coefficient, expected_coefficient)

@pytest.mark.parametrize(
    'm1,basis',
    (
        itertools.product(
            (I, X, Y, Z, H, SQRT_X, SQRT_Y, SQRT_Z, E00, E01, E10, E11),
            (PAULI_BASIS, STANDARD_BASIS),
        )
    ),
)
def test_expand_is_inverse_of_reconstruct(m1, basis):
    c1 = cirq.expand_matrix_in_orthogonal_basis(m1, basis)
    m2 = cirq.matrix_from_basis_coefficients(c1, basis)
    c2 = cirq.expand_matrix_in_orthogonal_basis(m2, basis)
    assert np.allclose(m1, m2)
    assert c1 == c2

@pytest.mark.parametrize(
    'coefficients,exponent',
    itertools.product(
        (
            (0, 0, 0, 0),
            (-1, 0, 0, 0),
            (0.5, 0, 0, 0),
            (0.5j, 0, 0, 0),
            (1, 0, 0, 0),
            (2, 0, 0, 0),
            (0, -1, 0, 0),
            (0, 0.5, 0, 0),
            (0, 0.5j, 0, 0),
            (0, 1, 0, 0),
            (0, 2, 0, 0),
            (0, 0, -1, 0),
            (0, 0, 0.5, 0),
            (0, 0, 0.5j, 0),
            (0, 0, 1, 0),
            (0, 0, 2, 0),
            (0, 0, 0, -1),
            (0, 0, 0, 0.5),
            (0, 0, 0, 0.5j),
            (0, 0, 0, 1),
            (0, 0, 0, 2),
            (0, -1, 0, -1),
            (0, 1, 0, 1j),
            (0, 0.5, 0, 0.5),
            (0, 0.5j, 0, 0.5j),
            (0, 0.5, 0, 0.5j),
            (0, 1, 0, 1),
            (0, 2, 0, 2),
            (0, 0.5, 0.5, 0.5),
            (0, 1, 1, 1),
            (0, 1.1j, 0.5 - 0.4j, 0.9),
            (0.7j, 1.1j, 0.5 - 0.4j, 0.9),
            (0.25, 0.25, 0.25, 0.25),
            (0.25j, 0.25j, 0.25j, 0.25j),
            (0.4, 0, 0.5, 0),
            (1, 2, 3, 4),
            (-1, -2, -3, -4),
            (-1, -2, 3, 4),
            (1j, 2j, 3j, 4j),
            (1j, 2j, 3, 4),
        ),
        (0, 1, 2, 3, 4, 5, 100, 101),
    ),
)
def test_pow_pauli_combination(coefficients, exponent):
    i = cirq.PAULI_BASIS['I']
    x = cirq.PAULI_BASIS['X']
    y = cirq.PAULI_BASIS['Y']
    z = cirq.PAULI_BASIS['Z']
    ai, ax, ay, az = coefficients

    matrix = ai * i + ax * x + ay * y + az * z
    expected_result = np.linalg.matrix_power(matrix, exponent)

    bi, bx, by, bz = cirq.pow_pauli_combination(ai, ax, ay, az, exponent)
    result = bi * i + bx * x + by * y + bz * z

    assert np.allclose(result, expected_result)

max_stars_repo_name: dabacon/Cirq

max_issues_repo_name: dabacon/Cirq

max_forks_repo_name: dabacon/Cirq

reason: import numpy,import scipy


text:
from hypergan.gan_component import GANComponent
import numpy as np
import tensorflow as tf

class BaseLoss(GANComponent):
    def __init__(self, gan, config, discriminator=None, generator=None):
        GANComponent.__init__(self, gan, config)
        self.metrics = {}
        self.sample = None
        self.ops = None
        self.discriminator = discriminator
        self.generator = generator

    def create(self, split=2):
        gan = self.gan
        config = self.config
        ops = self.gan.ops

        if self.discriminator is None:
            net = gan.discriminator.sample
        else:
            net = self.discriminator.sample

        if split == 2:
            d_real, d_fake = self.split_batch(net, split)
            d_loss, g_loss = self._create(d_real, d_fake)
        elif split == 3:
            d_real, d_fake, d_fake2 = self.split_batch(net, split)
            d_loss, g_loss = self._create(d_real, d_fake)
            d_loss2, g_loss2 = self._create(d_real, d_fake2)
            g_loss += g_loss2
            d_loss += d_loss2
            #does this double the signal of d_real?

        if d_loss is not None:
            d_loss = ops.squash(d_loss, tf.reduce_mean) #linear doesn't work with this, so we cant pass config.reduce
            self.metrics['d_loss'] = d_loss

            if config.minibatch:
                d_loss += self.minibatch(net)

            if config.gradient_penalty:
                gp = self.gradient_penalty()
                self.metrics['gradient_penalty'] = gp
                print("Gradient penalty applied")
                d_loss += gp

        if g_loss is not None:
            g_loss = ops.squash(g_loss, tf.reduce_mean)
            self.metrics['g_loss'] = g_loss

        self.metrics = self.metrics or sample_metrics

        self.sample = [d_loss, g_loss]
        self.d_loss = d_loss
        self.g_loss = g_loss

        return self.sample

    # This is openai's implementation of minibatch regularization
    def minibatch(self, net):
        discriminator = self.discriminator or self.gan.discriminator
        ops = discriminator.ops
        config = self.config
        batch_size = ops.shape(net)[0]
        single_batch_size = batch_size//2
        n_kernels = config.minibatch_kernels or 300
        dim_per_kernel = config.dim_per_kernel or 50
        print("[discriminator] minibatch from", net, "to", n_kernels*dim_per_kernel)
        x = ops.linear(net, n_kernels * dim_per_kernel)
        activation = tf.reshape(x, (batch_size, n_kernels, dim_per_kernel))

        big = np.zeros((batch_size, batch_size))
        big += np.eye(batch_size)
        big = tf.expand_dims(big, 1)
        big = tf.cast(big,dtype=ops.dtype)

        abs_dif = tf.reduce_sum(tf.abs(tf.expand_dims(activation,3) - tf.expand_dims(tf.transpose(activation, [1, 2, 0]), 0)), 2)
        mask = 1. - big
        masked = tf.exp(-abs_dif) * mask
        def half(tens, second):
            m, n, _ = tens.get_shape()
            m = int(m)
            n = int(n)
            return tf.slice(tens, [0, 0, second * single_batch_size], [m, n, single_batch_size])

        f1 = tf.reduce_sum(half(masked, 0), 2) / tf.reduce_sum(half(mask, 0))
        f2 = tf.reduce_sum(half(masked, 1), 2) / tf.reduce_sum(half(mask, 1))

        return ops.squash(ops.concat([f1, f2]))

    def gradient_penalty(self):
        config = self.config
        gan = self.gan
        gradient_penalty = config.gradient_penalty
        if has_attr(gan.inputs, 'gradient_penalty_label'):
            x = gan.inputs.gradient_penalty_label
        else:
            x = gan.inputs.x
        generator = self.generator or gan.generator
        g = generator.sample
        discriminator = self.discriminator or gan.discriminator
        shape = [1 for t in g.get_shape()]
        shape[0] = gan.batch_size()
        uniform_noise = tf.random_uniform(shape=shape,minval=0.,maxval=1.)
        print("[gradient penalty] applying x:", x, "g:", g, "noise:", uniform_noise)
        interpolates = x + uniform_noise * (g - x)
        reused_d = discriminator.reuse(interpolates)
        gradients = tf.gradients(reused_d, [interpolates])[0]
        penalty = tf.sqrt(tf.reduce_sum(tf.square(gradients), axis=1))
        penalty = tf.reduce_mean(tf.square(penalty - 1.))
        return float(gradient_penalty) * penalty

    def sigmoid_kl_with_logits(self, logits, targets):
       # broadcasts the same target value across the whole batch
       # this is implemented so awkwardly because tensorflow lacks an x log x op
       assert isinstance(targets, float)
       if targets in [0., 1.]:
         entropy = 0.
       else:
         entropy = - targets * np.log(targets) - (1. - targets) * np.log(1. - targets)
         return tf.nn.sigmoid_cross_entropy_with_logits(logits=logits, labels=tf.ones_like(logits) * targets) - entropy

max_stars_repo_name: SlipknotTN/HyperGAN

max_issues_repo_name: SlipknotTN/HyperGAN

max_forks_repo_name: SlipknotTN/HyperGAN


text:
import os
import shutil
import numpy as np
import cv2
from scipy.io import loadmat
import shutil

def flip_points(landmarks):
    result = landmarks.copy()
    result[0:8] = landmarks[9:17][::-1]
    result[9:17] = landmarks[0:8][::-1]
    result[17:22] = landmarks[22:27][::-1]
    result[22:27] = landmarks[17:22][::-1]
    result[36] = landmarks[45]
    result[45] = landmarks[36]
    result[37] = landmarks[44]
    result[44] = landmarks[37]
    result[38] = landmarks[43]
    result[43] = landmarks[38]
    result[39] = landmarks[42]
    result[42] = landmarks[39]
    result[40] = landmarks[47]
    result[47] = landmarks[40]
    result[41] = landmarks[46]
    result[46] = landmarks[41]
    result[31:33] = landmarks[34:36][::-1]
    result[34:36] = landmarks[31:33][::-1]
    result[50] = landmarks[52]
    result[52] = landmarks[50]
    result[49] = landmarks[53]
    result[53] = landmarks[49]
    result[48] = landmarks[54]
    result[54] = landmarks[48]
    result[59] = landmarks[55]
    result[59] = landmarks[55]
    result[50] = landmarks[52]
    result[52] = landmarks[50]
    result[58] = landmarks[56]
    result[56] = landmarks[58]
    result[60] = landmarks[64]
    result[64] = landmarks[60]
    result[61] = landmarks[63]
    result[63] = landmarks[61]
    result[67] = landmarks[65]
    result[65] = landmarks[67]
    return result.astype(np.int64)

subsets = ["AFW", "AFW_Flip", "HELEN", "HELEN_Flip", "IBUG", "IBUG_Flip", "LFPW", "LFPW_Flip"]

base_path = "E:/300W_LP"
output_base = "E:/300W_LP_Out"
for subset in subsets:
    subset_path = os.path.join(base_path, subset)
    img_files = filter(lambda x: ".jpg" in x, os.listdir(subset_path))
    output_subset = os.path.join(output_base, subset)
    if not os.path.exists(output_subset):
        os.makedirs(output_subset)
    for img_file in img_files:
        mat_file = img_file.replace(".jpg", ".mat")
        output_mat_path = os.path.join(output_subset, img_file.replace(".jpg", ".pts"))
        out = open(output_mat_path, "w")
        out.write("version: 1\nn_points: 68\n{\n")
        mat = loadmat(
            os.path.join(base_path, "landmarks", subset.replace("_Flip", ""), img_file.replace(".jpg", "_pts.mat")))[
            'pts_2d']
        if "Flip" in subset:
            mat[:, 0] = 450 - mat[:, 0]
            mat = flip_points(mat)
        for point in mat:
            out.write(str(point[0]))
            out.write(" ")
            out.write(str(point[1]))
            out.write("\n")
        out.write("}\n")
        img_src = os.path.join(subset_path, img_file)
        img_dst = os.path.join(output_subset, img_file)
        shutil.copy(img_src, img_dst)
        # img = cv2.imread(img_dst)
        # for point in mat:
        #     img = cv2.circle(img,tuple(point),2,(255,0,0),-1,1)
        # cv2.imshow("",img)
        # cv2.waitKey()

reason: import numpy,from scipy


text:
from tqdm import tqdm
import numpy as np

def extract_features(window):
    # in: n frames of EOG L, R, H, V

    # window -= np.mean(window, axis=0)  # shift the baseline sometime

    window_mins = np.min(window, axis=0)
    window_maxs = np.max(window, axis=0)

    out = np.mean(window, axis=0)  # mean of each EOG value
    out = np.append( out, window_mins )  # min for each EOG value
    out = np.append( out, window_maxs )  # max for each EOG value
    out = np.append( out, window_maxs - window_mins )  # max - min for each EOG value
    out = np.append( out, window[1] - window[0] )  # diff from start to end for each EOG value
    out = np.append( out, np.percentile(window, q=25, axis=0) )
    out = np.append( out, np.percentile(window, q=50, axis=0) )
    out = np.append( out, np.percentile(window, q=75, axis=0) )

    magnitudes = np.linalg.norm(window, axis=1)  # magnitude of each frame

    out = np.append( out, np.mean(magnitudes) )
    out = np.append( out, np.min(magnitudes) )
    out = np.append( out, np.max(magnitudes) )
    out = np.append( out, out[-1] - out[-2] )  # max - min of magnitudes
    out = np.append( out, magnitudes[-1] - magnitudes[0] )  # diff from start to end magnitude

    # out = np.append( out, np.histogram(window[:,3], normed=True) )

    return out

def slidingWindow(sequence, winSize, stride=1):
    """Returns a generator that will iterate through
    the defined chunks of input sequence.  Input sequence
    must be iterable.
    @author: snoran
    Thanks to https://scipher.wordpress.com/2010/12/02/simple-sliding-window-iterator-in-python/"""

    # Verify the inputs
    try:
        it = iter(sequence)
    except TypeError:
        raise Exception("**ERROR** sequence must be iterable.")
    if not (isinstance(winSize, int) and isinstance(stride, int)):
        raise Exception("**ERROR** type(winSize) and type(stride) must be int.")
    if stride > winSize:
        raise Exception("**ERROR** stride must not be larger than winSize.")
    if winSize > len(sequence):
        raise Exception("**ERROR** winSize must not be larger than sequence length.")

    # Pre-compute number of chunks to emit
    numOfChunks = int((len(sequence)-winSize)/stride)+1

    # Do the work
    for i in range(0,numOfChunks*stride,stride):
        yield i, sequence[i:i+winSize]

# put the data in terms of inputs and outputs
def extract_features_from_dataframe(df, window_size):
    data = df.as_matrix()  # turn it into a numpy array

    inputs = []
    outputs = []

    samples = df.shape[0] - window_size + 1

    with tqdm(total=samples) as pbar:
        for i, window_full in slidingWindow(data, window_size):
            window = window_full[:, 2:-1]  # shave off frame id, time, and blink value
            features = extract_features(window)

            inputs += [ features ]
            # outputs += [ np.mean(window_full[:, -1]) >= 0.6 ]  # turn the blink data into booleans
            outputs += [ window_full[int(window_full.shape[0]/2), -1] ]  # leave the blink regression as-is

            pbar.update(1)

    inputs = np.array(inputs)
    outputs = np.array(outputs)

    return inputs, outputs

max_stars_repo_name: AlexLamson/JINS-blink

max_issues_repo_name: AlexLamson/JINS-blink

max_forks_repo_name: AlexLamson/JINS-blink


text:
import numpy as np
import random

def match_random(a,b):
    state = np.random.get_state()
    np.random.shuffle(a)
    np.random.set_state(state)
    np.random.shuffle(b)

def random_transform(a,alpha):
	return a*random.uniform(1-alpha,1+alpha)

def normlize(npdata,justprice = False):
    _min = np.min(npdata)
    _max = np.max(npdata)
    if justprice:
        _min = 34900.0
        _max = 755000.0
    return (npdata-_min)/(_max-_min)

def convert2price(tensor):
    return tensor*(755000.0-34900.0)+34900

# def normlize(npdata,justprice = False):

#     return np.log1p(npdata)

# def convert2price(tensor):
#     return np.expm1(tensor)


text:
import torch
import random
import numpy as np
import cv2
import torch.nn.functional as F

from .registry import TRANSFORMS

CV2_MODE = {
    'bilinear': cv2.INTER_LINEAR,
    'nearest': cv2.INTER_NEAREST,
    'cubic': cv2.INTER_CUBIC,
    'area': cv2.INTER_AREA,
}

CV2_BORDER_MODE = {
    'constant': cv2.BORDER_CONSTANT,
    'reflect': cv2.BORDER_REFLECT,
    'reflect101': cv2.BORDER_REFLECT101,
    'replicate': cv2.BORDER_REPLICATE,
}

class Compose(object):
    def __init__(self, transforms):
        self.transforms = transforms

    def __call__(self, image, mask):
        for t in self.transforms:
            image, mask = t(image, mask)
        return image, mask

@TRANSFORMS.register_module
class RandomScale(object):
    def __init__(self, min_scale, max_scale, scale_step=0.0, mode='bilinear'):
        self.min_scale = min_scale
        self.max_scale = max_scale
        self.scale_step = scale_step
        self.mode = mode

    def get_scale_factor(self, min_scale, max_scale, scale_step):
        if min_scale == max_scale:
            return min_scale

        if scale_step == 0:
            return random.uniform(min_scale, max_scale)

        num_steps = int((max_scale - min_scale) / scale_step + 1)
        scale_factors = np.linspace(min_scale, max_scale, num_steps)
        scale_factor = np.random.choice(scale_factors).item()

        return scale_factor

    def __call__(self, image, mask):
        scale_factor = self.get_scale_factor(self.min_scale, self.max_scale, self.scale_step)

        if scale_factor == 1.0:
            return image, mask

        h, w, c = image.shape
        new_h = int(h * scale_factor)
        new_w = int(w * scale_factor)

        torch_image = torch.from_numpy(image).permute(2,0,1).unsqueeze(0)
        torch_mask = torch.from_numpy(mask).unsqueeze(0).unsqueeze(0)
        torch_image = F.interpolate(torch_image, size=(new_h, new_w), mode=self.mode, align_corners=True)
        torch_mask = F.interpolate(torch_mask, size=(new_h, new_w), mode='nearest')

        new_image = torch_image.squeeze().permute(1,2,0).numpy()
        new_mask = torch_mask.squeeze().numpy()

        return new_image, new_mask

@TRANSFORMS.register_module
class RandomCrop(object):
    def __init__(self, height, width, image_value, mask_value):
        self.height = height
        self.width = width
        self.image_value = image_value
        self.mask_value = mask_value

    def __call__(self, image, mask):
        h, w, c = image.shape
        target_height = h + max(self.height-h, 0)
        target_width = w + max(self.width-w, 0)

        image_pad_value = np.reshape(np.array(self.image_value, dtype=image.dtype), [1,1,3])
        mask_pad_value = np.reshape(np.array(self.mask_value, dtype=mask.dtype), [1,1])

        new_image = np.tile(image_pad_value, (target_height, target_width, 1))
        new_mask = np.tile(mask_pad_value, (target_height, target_width))

        new_image[:h, :w, :] = image
        new_mask[:h, :w] = mask

        assert np.count_nonzero(mask != self.mask_value) == np.count_nonzero(new_mask != self.mask_value)

        y1 = int(random.uniform(0, target_height-self.height+1))
        y2 = y1 + self.height
        x1 = int(random.uniform(0, target_width-self.width+1))
        x2 = x1 + self.width

        new_image = new_image[y1:y2, x1:x2, :]
        new_mask = new_mask[y1:y2, x1:x2]

        return new_image, new_mask

@TRANSFORMS.register_module
class PadIfNeeded(object):
    def __init__(self, height, width, image_value, mask_value):
        self.height = height
        self.width = width
        self.image_value = image_value
        self.mask_value = mask_value

    def __call__(self, image, mask):
        h, w, c = image.shape

        assert h <= self.height and w <= self.width

        target_height = h + max(self.height-h, 0)
        target_width = w + max(self.width-w, 0)

        image_pad_value = np.reshape(np.array(self.image_value, dtype=image.dtype), [1,1,3])
        mask_pad_value = np.reshape(np.array(self.mask_value, dtype=mask.dtype), [1,1])

        new_image = np.tile(image_pad_value, (target_height, target_width, 1))
        new_mask = np.tile(mask_pad_value, (target_height, target_width))

        new_image[:h, :w, :] = image
        new_mask[:h, :w] = mask

        assert np.count_nonzero(mask!=self.mask_value) == np.count_nonzero(new_mask!=self.mask_value)

        return new_image, new_mask

@TRANSFORMS.register_module
class HorizontalFlip(object):
    def __init__(self, p=0.5):
        self.p = p

    def __call__(self, image, mask):
        if random.random() > self.p:
            image = cv2.flip(image, 1)
            mask = cv2.flip(mask, 1)

        return image, mask

@TRANSFORMS.register_module
class RandomRotate(object):
    def __init__(self, p=0.5, degrees=30, mode='bilinear', border_mode='reflect101', image_value=None, mask_value=None):
        self.p = p
        self.degrees = (-degrees, degrees) if isinstance(degrees, (int, float)) else degrees
        self.mode = CV2_MODE[mode]
        self.border_mode = CV2_BORDER_MODE[border_mode]
        self.image_value = image_value
        self.mask_value = mask_value

    def __call__(self, image, mask):
        if random.random() < self.p:
            h, w, c = image.shape

            angle = random.uniform(*self.degrees)
            matrix = cv2.getRotationMatrix2D((w / 2, h / 2), angle, 1.0)

            image = cv2.warpAffine(image, M=matrix, dsize=(w,h), flags=self.mode, borderMode=self.border_mode, borderValue=self.image_value)
            mask = cv2.warpAffine(mask, M=matrix, dsize=(w,h), flags=cv2.INTER_NEAREST, borderMode=self.border_mode, borderValue=self.mask_value)

        return image, mask

@TRANSFORMS.register_module
class GaussianBlur(object):
    def __init__(self, p=0.5, ksize=7):
        self.p = p
        self.ksize = (ksize, ksize) if isinstance(ksize, int) else ksize

    def __call__(self, image, mask):
        if random.random() < self.p:
            image = cv2.GaussianBlur(image, ksize=self.ksize, sigmaX=0)

        return image, mask

@TRANSFORMS.register_module
class Normalize(object):
    def __init__(self, mean=(123.675, 116.280, 103.530), std=(58.395, 57.120, 57.375)):
        self.mean = mean
        self.std = std

    def __call__(self, image, mask):
        mean = np.reshape(np.array(self.mean, dtype=image.dtype), [1,1,3])
        std = np.reshape(np.array(self.std, dtype=image.dtype), [1,1,3])
        denominator = np.reciprocal(std, dtype=image.dtype)

        new_image = (image - mean) * denominator
        new_mask = mask

        return new_image, new_mask

@TRANSFORMS.register_module
class ToTensor(object):
    def __call__(self, image, mask):
        image = torch.from_numpy(image).permute(2, 0, 1)
        mask = torch.from_numpy(mask)

        return image, mask


text:
"""Functions that help with data processing for human3.6m"""
import numpy as np
import copy
from tqdm.autonotebook import tqdm
import matplotlib as mpl
import xml.etree.ElementTree as ET
from os import path
from lib.utils import add_joints_to_img
from data.base_dataset import BaseDataset
import cv2
from lib.utils import make_joint_img, scale_img
from models.vunets import VunetOrg
import torch
# required for 3d box plot of matplotlib, though not used
# package seems not to be exported
from mpl_toolkits.mplot3d import axes3d, Axes3D

def euler_to_rotation_matrix(angles, deg=True, format="zxy"):
    """

    :param angles: angles around axes (zxy) or xyz
    :return:
    """
    if deg:
        angles = np.radians(angles)

    cx = np.cos(angles[0])
    cy = np.cos(angles[1])
    cz = np.cos(angles[2])
    sx = np.sin(angles[0])
    sy = np.sin(angles[1])
    sz = np.sin(angles[2])

    # init = np.eye(3)
    if format == "zxy":
        R = np.asarray(
            [
                [cy * cz - sx * sy * sz, cy * sz + sx * sy * cz, -sy * cx],
                [-cx * sz, cx * cz, sx],
                [sy * cz + cy * sx * sz, sy * sz - cy * sx * cz, cy * cx],
            ]
        )
    elif format == "xyz":
        R = (
            np.asarray([[cz, sz, 0.0], [-sz, cz, 0.0], [0.0, 0.0, 1.0]])
            @ np.asarray([[cy, 0, -sy], [0.0, 1.0, 0.0], [sy, 0.0, cy]])
            @ np.asarray([[1.0, 0, 0.0], [0.0, cx, sx], [0.0, -sx, cx]])
        )
    else:
        raise NotImplementedError()

    return R

def rotmat2euler(R):
    """
Converts a rotation matrix to Euler angles
Matlab port to python for evaluation purposes
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/mhmublv/Motion/RotMat2Euler.m#L1

Args
  R: a 3x3 rotation matrix
Returns
  eul: a 3x1 Euler angle representation of R
"""
    if R[0, 2] == 1 or R[0, 2] == -1:
        # special case
        E3 = 0  # set arbitrarily
        dlta = np.arctan2(R[0, 1], R[0, 2])

        if R[0, 2] == -1:
            E2 = np.pi / 2
            E1 = E3 + dlta
        else:
            E2 = -np.pi / 2
            E1 = -E3 + dlta

    else:
        E2 = -np.arcsin(R[0, 2])
        E1 = np.arctan2(R[1, 2] / np.cos(E2), R[2, 2] / np.cos(E2))
        E3 = np.arctan2(R[0, 1] / np.cos(E2), R[0, 0] / np.cos(E2))

    eul = np.array([E1, E2, E3])
    return eul

def quat2expmap(q):
    """
Converts a quaternion to an exponential map
Matlab port to python for evaluation purposes
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/mhmublv/Motion/quat2expmap.m#L1

Args
  q: 1x4 quaternion
Returns
  r: 1x3 exponential map
Raises
  ValueError if the l2 norm of the quaternion is not close to 1
"""
    if np.abs(np.linalg.norm(q) - 1) > 1e-3:
        # raise ValueError
        print("corrupteeeeed....")

    sinhalftheta = np.linalg.norm(q[1:])
    coshalftheta = q[0]

    r0 = np.divide(q[1:], (np.linalg.norm(q[1:]) + np.finfo(np.float32).eps))
    theta = 2 * np.arctan2(sinhalftheta, coshalftheta)
    theta = np.mod(theta + 2 * np.pi, 2 * np.pi)

    if theta > np.pi:
        theta = 2 * np.pi - theta
        r0 = -r0

    r = r0 * theta
    return r

def rotmat2quat(R):
    """
Converts a rotation matrix to a quaternion
Matlab port to python for evaluation purposes
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/mhmublv/Motion/rotmat2quat.m#L4

Args
  R: 3x3 rotation matrix
Returns
  q: 1x4 quaternion
"""
    rotdiff = R - R.T

    r = np.zeros(3)
    r[0] = -rotdiff[1, 2]
    r[1] = rotdiff[0, 2]
    r[2] = -rotdiff[0, 1]
    sintheta = np.linalg.norm(r) / 2
    r0 = np.divide(r, np.linalg.norm(r) + np.finfo(np.float32).eps)

    costheta = (np.trace(R) - 1) / 2

    theta = np.arctan2(sintheta, costheta)

    q = np.zeros(4)
    q[0] = np.cos(theta / 2)
    q[1:] = r0 * np.sin(theta / 2)
    return q

def rotmat2expmap(R):
    return quat2expmap(rotmat2quat(R))

def expmap2rotmat(r):
    """
Converts an exponential map angle to a rotation matrix
Matlab port to python for evaluation purposes
I believe this is also called Rodrigues' formula
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/mhmublv/Motion/expmap2rotmat.m

Args
  r: 1x3 exponential map
Returns
  R: 3x3 rotation matrix
"""
    theta = np.linalg.norm(r)
    r0 = np.divide(r, theta + np.finfo(np.float32).eps)
    r0x = np.array([0, -r0[2], r0[1], 0, 0, -r0[0], 0, 0, 0]).reshape(3, 3)
    r0x = r0x - r0x.T
    R = (
        np.eye(3, 3)
        + np.sin(theta) * r0x
        + (1 - np.cos(theta)) * (r0x).dot(r0x)
    )
    return R

def unNormalizeData(normalizedData, data_mean, data_std, dimensions_to_ignore):
    """Borrowed from SRNN code. Reads a csv file and returns a float32 matrix.
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/generateMotionData.py#L12

Args
  normalizedData: nxd matrix with normalized data
  data_mean: vector of mean used to normalize the data
  data_std: vector of standard deviation used to normalize the data
  dimensions_to_ignore: vector with dimensions not used by the model
  actions: list of strings with the encoded actions
  one_hot: whether the data comes with one-hot encoding
Returns
  origData: data originally used to
"""
    T = normalizedData.shape[0]
    D = data_mean.shape[0]

    origData = np.zeros((T, D), dtype=np.float32)
    dimensions_to_use = []
    for i in range(D):
        if i in dimensions_to_ignore:
            continue
        dimensions_to_use.append(i)
    dimensions_to_use = np.array(dimensions_to_use)

    origData[:, dimensions_to_use] = normalizedData

    # potentially ineficient, but only done once per experimentdata_conversions
    stdMat = data_std.reshape((1, D))
    stdMat = np.repeat(stdMat, T, axis=0)
    meanMat = data_mean.reshape((1, D))
    meanMat = np.repeat(meanMat, T, axis=0)
    origData = np.multiply(origData, stdMat) + meanMat
    return origData

def revert_output_format(poses, data_mean, data_std, dim_to_ignore):
    """
Converts the output of the neural network to a format that is more easy to
manipulate for, e.g. conversion to other format or visualization

Args
  poses: Poses, given in an array of shape (n, dim_poses) where n denotes the number of overall samples
Returns
  poses_out: A tensor of size (batch_size, seq_length, dim) output. Each
  batch is an n-by-d sequence of poses.data_conversions
"""
    # seq_len = len(poses)
    # if seq_len == 0:
    #     return []
    #
    # batch_size, dim = poses[0].shape
    #
    # poses_out = np.concatenate(poses)
    # poses_out = np.reshape(poses_out, (seq_len, batch_size, dim))
    # poses_out = np.transpose(poses_out, [1, 0, 2])

    # poses_out_list = []

    poses_unnorm = unNormalizeData(poses, data_mean, data_std, dim_to_ignore)

    return poses_unnorm

def readCSVasFloat(filename):
    """
Borrowed from SRNN code. Reads a csv and returns a float matrix.
https://github.com/asheshjain399/NeuralModels/blob/master/neuralmodels/utils.py#L34

Args
  filename: string. Path to the csv file
Returns
  returnArray: the read data in a float32 matrix
"""
    returnArray = []
    lines = open(filename).readlines()
    for line in lines:
        line = line.strip().split(",")
        if len(line) > 0:
            returnArray.append(np.array([np.float32(x) for x in line]))

    returnArray = np.array(returnArray)
    return returnArray

def load_data(path_to_dataset, subjects, actions, one_hot):
    """
Borrowed from SRNN code. This is how the SRNN code reads the provided .txt files
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/processdata.py#L270

Args
  path_to_dataset: string. directory where the data resides
  subjects: list of numbers. The subjects to load
  actions: list of string. The actions to load
  one_hot: Whether to add a one-hot encoding to the data
Returns
  trainData: dictionary with k:v
    k=(subject, action, subaction, 'even'), v=(nxd) un-normalize.format(
                    , , ,
                )d data
  completeData: nxd matrix with all the data. Used to normlization stats
"""
    nactions = len(actions)

    trainData = {}
    completeData = []
    for subj in tqdm(subjects):
        print(f"processing subject # {subj}")
        for action_idx in tqdm(np.arange(len(actions))):
            print(f"processing action #{action_idx}")
            action = actions[action_idx]

            for subact in [1, 2]:  # subactions

                print(
                    f"Reading subject {subj}, action {action}, subaction {subact}"
                )

                filename = f"{path_to_dataset}/S{subj}/{action}_{subact}.txt"
                action_sequence = readCSVasFloat(filename)

                n, d = action_sequence.shape
                even_list = list(range(0, n, 1))

                if one_hot:
                    # Add a one-hot encoding at the end of the representation
                    the_sequence = np.zeros(
                        (len(even_list), d + nactions), dtype=float
                    )
                    the_sequence[:, 0:d] = action_sequence[even_list, :]
                    the_sequence[:, d + action_idx] = 1
                    trainData[(subj, action, subact, "even")] = the_sequence
                else:
                    trainData[(subj, action, subact, "even")] = action_sequence[
                        even_list, :
                    ]

                if len(completeData) == 0:
                    completeData = copy.deepcopy(action_sequence)
                else:
                    completeData = np.append(
                        completeData, action_sequence, axis=0
                    )

    return trainData, completeData

def normalize_data(data, data_mean, data_std, dim_to_use, actions, one_hot):
    """
Normalize input data by removing unused dimensions, subtracting the mean and
dividing by the standard deviation

Args
  data: nx99 matrix with data to normalize
  data_mean: vector of mean used to normalize the data
  data_std: vector of standard deviation used to normalize the data
  dim_to_use: vector with dimensions used by the model
  actions: list of strings with the encoded actions
  one_hot: whether the data comes with one-hot encoding
Returns
  data_out: the passed data matrix, but normalized
"""
    data_out = {}
    nactions = len(actions)

    if not one_hot:
        # No one-hot encoding... no need to do anything special
        for key in list(data.keys()):
            data_out[key] = np.divide((data[key] - data_mean), data_std)
            data_out[key] = data_out[key][:, dim_to_use]

    else:
        # TODO hard-coding 99 dimensions for un-normalized human poses
        for key in list(data.keys()):
            data_out[key] = np.divide(
                (data[key][:, 0:99] - data_mean), data_std
            )
            data_out[key] = data_out[key][:, dim_to_use]
            data_out[key] = np.hstack((data_out[key], data[key][:, -nactions:]))

    return data_out

def normalization_stats(completeData):
    """"
Also borrowed for SRNN code. Computes mean, stdev and dimensions to ignore.
https://github.com/asheshjain399/RNNexp/blob/srnn/structural_rnn/CRFProblems/H3.6m/processdata.py#L33

Args
  completeData: nx99 matrix with data to normalize
Returns
  data_mean: vector of mean used to normalize the data
  data_std: vector of standard deviation used to normalize the data
  dimensions_to_ignore: vector with dimensions not used by the model
  dimensions_to_use: vector with dimensions used by the model
"""
    data_mean = np.mean(completeData, axis=0)
    data_std = np.std(completeData, axis=0)

    dimensions_to_ignore = []
    dimensions_to_use = []

    dimensions_to_ignore.extend(list(np.where(data_std < 1e-4)[0]))
    dimensions_to_use.extend(list(np.where(data_std >= 1e-4)[0]))

    data_std[dimensions_to_ignore] = 1.0

    return data_mean, data_std, dimensions_to_ignore, dimensions_to_use

def _some_variables(use_posInd=False):
    """
  We define some variables that are useful to run the kinematic tree

  Args
    None
  Returns
    parent: 32-long vector with parent-child relationships in the kinematic tree
    offset: 96-long vector with bone lenghts
    rotInd: 32-long list with indices into angles
    expmapInd: 32-long list with indices into expmap angles
  """

    parent = (
        np.array(
            [
                0,
                1,
                2,
                3,
                4,
                5,
                1,
                7,
                8,
                9,
                10,
                1,
                12,
                13,
                14,
                15,
                13,
                17,
                18,
                19,
                20,
                21,
                20,
                23,
                13,
                25,
                26,
                27,
                28,
                29,
                28,
                31,
            ]
        )
        - 1
    )

    offset = np.array(
        [
            0.000000,
            0.000000,
            0.000000,
            -132.948591,
            0.000000,
            0.000000,
            0.000000,
            -442.894612,
            0.000000,
            0.000000,
            -454.206447,
            0.000000,
            0.000000,
            0.000000,
            162.767078,
            0.000000,
            0.000000,
            74.999437,
            132.948826,
            0.000000,
            0.000000,
            0.000000,
            -442.894413,
            0.000000,
            0.000000,
            -454.206590,
            0.000000,
            0.000000,
            0.000000,
            162.767426,
            0.000000,
            0.000000,
            74.999948,
            0.000000,
            0.100000,
            0.000000,
            0.000000,
            233.383263,
            0.000000,
            0.000000,
            257.077681,
            0.000000,
            0.000000,
            121.134938,
            0.000000,
            0.000000,
            115.002227,
            0.000000,
            0.000000,
            257.077681,
            0.000000,
            0.000000,
            151.034226,
            0.000000,
            0.000000,
            278.882773,
            0.000000,
            0.000000,
            251.733451,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            99.999627,
            0.000000,
            100.000188,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            257.077681,
            0.000000,
            0.000000,
            151.031437,
            0.000000,
            0.000000,
            278.892924,
            0.000000,
            0.000000,
            251.728680,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            99.999888,
            0.000000,
            137.499922,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
        ]
    )

    offset = offset.reshape(-1, 3)

    rotInd = [
        [5, 6, 4],
        [8, 9, 7],
        [11, 12, 10],
        [14, 15, 13],
        [17, 18, 16],
        [],
        [20, 21, 19],
        [23, 24, 22],
        [26, 27, 25],
        [29, 30, 28],
        [],
        [32, 33, 31],
        [35, 36, 34],
        [38, 39, 37],
        [41, 42, 40],
        [],
        [44, 45, 43],
        [47, 48, 46],
        [50, 51, 49],
        [53, 54, 52],
        [56, 57, 55],
        [],
        [59, 60, 58],
        [],
        [62, 63, 61],
        [65, 66, 64],
        [68, 69, 67],
        [71, 72, 70],
        [74, 75, 73],
        [],
        [77, 78, 76],
        [],
    ]

    # definitions are originating from matlab file --> bring them to zero based indexing
    rotInd = [[e - 1 for e in s if len(s) > 0] for s in rotInd]
    posInd = [0, 1, 2] if use_posInd else None

    expmapInd = np.split(np.arange(4, 100) - 1, 32)

    return parent, offset, rotInd, expmapInd, posInd

def apply_affine_transform(x, M):
    """
    Returns affine transformation R*x + t where M = [R,t] \in R^3x4
    :param x:
    :param M:
    :return:
    """
    is1d = len(x.shape) == 1
    if is1d:
        x = np.expand_dims(x, axis=0)

    x_hom = np.concatenate(
        [x, np.ones((x.shape[0], 1), dtype=x.dtype)], axis=-1
    )
    x_out = x_hom @ M.T
    if is1d:
        x_out = np.squeeze(x_out, axis=0)
    return x_out

def fkl(
    angles, parent, offset, rotInd, expmapInd, posInd=None, use_euler=False
):
    """
  Convert joint angles and bone lenghts into the 3d points of a person.

  Args
    angles: vector with 3d position and 3d joint angles in expmap format (99-long) if use_euler=False, or in bvh-motion-capture format (78-long) if use_euler=True
    parent: 32-long vector with parent-child relationships in the kinematic tree
    offset: 96-long vector with bone lenghts
    rotInd: 32-long list with indices into angles
    expmapInd: 32-long list with indices into expmap angles
  Returns
    xyz: 32x3 3d points that represent a person in 3d space
  """
    if use_euler:
        assert len(angles) == 78
    else:
        assert len(angles) == 99

    # Structure that indicates parents for each joint
    njoints = 32
    xyzStruct = [dict() for x in range(njoints)]

    for i in np.arange(njoints):

        if use_euler:
            current_eulers = (
                np.asarray([0.0, 0.0, 0.0])
                if len(rotInd[i]) == 0
                else angles[rotInd[i]]
            )
            thisRotation = euler_to_rotation_matrix(current_eulers, deg=True)
        else:
            r = angles[expmapInd[i]]
            thisRotation = expmap2rotmat(r)

        if parent[i] == -1:  # Root node
            if posInd is not None:
                # note this is currently not in expmap format, since it's a translation and no rotation
                thisPosition = np.array(
                    [angles[posInd[0]], angles[posInd[1]], angles[posInd[2]]]
                )
            else:
                thisPosition = np.asarray([0.0, 0.0, 0.0])
            xyzStruct[i]["rotation"] = thisRotation
            xyzStruct[i]["xyz"] = (
                np.reshape(offset[i, :], (1, 3)) + thisPosition
            )
        else:
            xyzStruct[i]["xyz"] = (offset[i, :]).dot(
                xyzStruct[parent[i]]["rotation"]
            ) + xyzStruct[parent[i]]["xyz"]
            xyzStruct[i]["rotation"] = thisRotation.dot(
                xyzStruct[parent[i]]["rotation"]
            )

    xyz = [xyzStruct[i]["xyz"] for i in range(njoints)]
    xyz = np.array(xyz).squeeze()

    return np.reshape(xyz, [-1])

def revert_coordinate_space(channels, R0, T0):
    """
    Aranges the translation and rotation of the root joint such that the series of poses appears
    to well aranged for visualization in the corrdinate frame
    :param channels: poses, in expmap format; shape (seq_len,99)
    :param R0: initial rotation of root joint
    :param T0: initial translation of root joint
    :param use_euler: whether to use euler angle parametrization or expmap representation
    :return: series of aranged poses
    """
    n, d = channels.shape

    channels_rec = copy.copy(channels)
    R_prev = R0
    T_prev = T0
    rootRotInd = np.arange(3, 6)

    # Loop through the passed posses
    for ii in range(n):

        R_diff = expmap2rotmat(channels[ii, rootRotInd])
        R = R_diff.dot(R_prev)

        channels_rec[ii, rootRotInd] = rotmat2expmap(R)
        T = T_prev + (
            (R_prev.T).dot(np.reshape(channels[ii, :3], [3, 1]))
        ).reshape(-1)
        channels_rec[ii, :3] = T

        T_prev = T
        R_prev = R

    return channels_rec

def kinematic_tree():

    # define mappings
    mappings = {
        "name": str,
        "id": int,
        "offset": float,
        "parent": int,
        "order": str,
        "rotInd": int,
        "children": int,
    }

    p = path.split(path.abspath(__file__))[0]

    tree = ET.parse(path.join(p, "metadata.xml"))
    root = tree.getroot()

    # angles_dict = etree_to_dict(angles)
    skel_tree = root.find("skel_angles")
    kin_tree = {
        "root": "Hips",
        "name": [],
        "id": [],
        "offset": [],
        "parent": [],
        "order": [],
        "rotInd": [],
        "children": [],
    }
    for i, tr in enumerate(skel_tree):
        if not tr.tag == "tree":
            continue
        for i, item in enumerate(tr):
            childs = list(item.iter())

            if i == 0:
                # x,y,z

                posInd = np.asarray(
                    list(map(int, childs[21].text[1:-1].split())), dtype=np.int
                )
                # I HATE MATLAB INDEXING
                kin_tree["posInd"] = {
                    "ids": [posInd[0] - 1, posInd[1] - 1, posInd[2] - 1],
                    "order": "xyz",
                }

            for c in childs:
                if c.tag in mappings.keys():
                    if c.text != "None" and c.text is not None:
                        kin_tree[c.tag].append(
                            list(map(mappings[c.tag], c.text[1:-1].split()))
                            if c.text.startswith("[")
                            else mappings[c.tag](c.text)
                        )
                    else:
                        kin_tree[c.tag].append([])

    kin_tree["order"] = list(
        map(lambda x: None if x == [] else x, kin_tree["order"])
    )
    # same (0-based) indexing is in _some_variables()-function
    kin_tree["parent"] = [e - 1 for e in kin_tree["parent"]]
    kin_tree["parent"] = np.asarray(kin_tree["parent"])
    kin_tree["offset"] = np.asarray(kin_tree["offset"]) * 10
    # kin_tree["offset"] = [[e * 10 for e in r] for r in kin_tree["offset"]]
    kin_tree["expmapInd"] = np.split(np.arange(4, 100) - 1, 32)
    # did I already mention that I HATE MATLAB INDEXING?!!!!!!!!!!!!!!
    kin_tree["rotInd"] = [
        [e - 1 for e in s if len(s) > 0] for s in kin_tree["rotInd"]
    ]
    kin_tree["children"] = [
        [e - 1 for e in s] if isinstance(s, list) else s - 1
        for s in kin_tree["children"]
    ]

    # this offset is more accurate when kps are projected to the image plane
    offset = np.array(
        [
            0.000000,
            0.000000,
            0.000000,
            -132.948591,
            0.000000,
            0.000000,
            0.000000,
            -442.894612,
            0.000000,
            0.000000,
            -454.206447,
            0.000000,
            0.000000,
            0.000000,
            162.767078,
            0.000000,
            0.000000,
            74.999437,
            132.948826,
            0.000000,
            0.000000,
            0.000000,
            -442.894413,
            0.000000,
            0.000000,
            -454.206590,
            0.000000,
            0.000000,
            0.000000,
            162.767426,
            0.000000,
            0.000000,
            74.999948,
            0.000000,
            0.100000,
            0.000000,
            0.000000,
            233.383263,
            0.000000,
            0.000000,
            257.077681,
            0.000000,
            0.000000,
            121.134938,
            0.000000,
            0.000000,
            115.002227,
            0.000000,
            0.000000,
            257.077681,
            0.000000,
            0.000000,
            151.034226,
            0.000000,
            0.000000,
            278.882773,
            0.000000,
            0.000000,
            251.733451,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            99.999627,
            0.000000,
            100.000188,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            257.077681,
            0.000000,
            0.000000,
            151.031437,
            0.000000,
            0.000000,
            278.892924,
            0.000000,
            0.000000,
            251.728680,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
            99.999888,
            0.000000,
            137.499922,
            0.000000,
            0.000000,
            0.000000,
            0.000000,
        ]
    )

    offset = offset.reshape(-1, 3)

    kin_tree["offset"] = offset

    return kin_tree

def camera_projection(poses, camera_parameters):
    """

    :param poses:
    :param camera_parameters: camera parameters as (f_x, x_0, f_y, y_0)
    :return:
    """
    camera_mat = np.asarray(
        [
            [camera_parameters[0], 0.0, camera_parameters[1]],
            [0.0, camera_parameters[2], camera_parameters[3]],
            [0.0, 0.0, 1.0],
        ]
    )
    if len(poses.shape) not in [2, 3]:
        raise Exception("Poses array has to be of dim 2 or 3.")
    # divide by z
    poses3d = poses / np.expand_dims(poses[..., -1], axis=-1)
    poses2d_homogenous = poses3d @ camera_mat.T
    # poses2d_homogenous = np.moveaxis(poses2d_homogenous,[0,1],[1,0])
    return poses2d_homogenous[..., :-1]

class Ax3DPose(object):
    def __init__(
        self,
        ax,
        dataset: BaseDataset = None,
        ticks=True,
        marker_color=None,
        lcolor="#3498db",
        rcolor="#e74c3c",
        limits=None,
    ):
        """
    Create a 3d pose visualizer that can be updated with new poses.

    Args
      ax: 3d axis to plot the 3d pose on
      lcolor: String. Colour for the left part of the body
      rcolor: String. Colour for the right part of the body
    """
        if marker_color is not None:
            self.marker_color = marker_color

        # Start and endpoints of our representation

        # self.I = (
        #     np.array([1, 2, 3, 1, 7, 8, 1, 13, 14, 15, 14, 18, 19, 14, 26, 27])
        #     - 1
        # )
        # self.J = (
        #     np.array([2, 3, 4, 7, 8, 9, 13, 14, 15, 16, 18, 19, 20, 26, 27, 28])
        #     - 1
        # )
        # # Left / right indicator
        # self.LR = np.array(
        #     [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1], dtype=bool
        # )
        self.I = [i[0] for i in dataset.joint_model.total_relative_joints]
        self.J = [i[1] for i in dataset.joint_model.total_relative_joints]
        self.LR = []
        self.ax = ax

        self.n_kps = len(dataset.joint_model.kps_to_use)

        vals = np.zeros((self.n_kps, 3))
        self.limits = limits

        # Make connection matrix
        self.plots = []
        for i in np.arange(len(self.I)):
            x = np.array([vals[self.I[i], 0], vals[self.J[i], 0]])
            y = np.array([vals[self.I[i], 1], vals[self.J[i], 1]])
            z = np.array([vals[self.I[i], 2], vals[self.J[i], 2]])
            if self.marker_color is not None:
                self.plots.append(
                    self.ax.plot(x, y, z, lw=2, c=self.marker_color)
                )
            else:
                self.plots.append(
                    self.ax.plot(
                        x, y, z, lw=2, c=lcolor if self.LR[i] else rcolor
                    )
                )
        self.ticks = ticks
        if self.ticks:
            self.ax.set_xlabel("x")
            self.ax.set_ylabel("y")
            self.ax.set_zlabel("z")
        else:
            self.hide_axes()

    def hide_axes(self):
        # remove grid
        self.ax.grid(False)
        # remove axis
        self.ax.set_axis_off()
        # remomve labels
        self.ax.set_label("")
        # remove surrounding box
        self.ax.set_frame_on(False)

    def update(self, channels, lcolor="#3498db", rcolor="#e74c3c"):
        """
    Update the plotted 3d pose.

    Args
      channels: 96-dim long np array. The pose to plot.
      lcolor: String. Colour for the left part of the body.
      rcolor: String. Colour for the right part of the body.
    Returns
      Nothing. Simply updates the axis with the new pose.
    """

        if hasattr(self, "marker_color"):
            lcolor = rcolor = self.marker_color

        if len(channels.shape) < 3:
            vals = np.reshape(channels, (self.n_kps, -1))
        else:
            vals = channels

        for i in np.arange(len(self.I)):
            x = np.array([vals[self.I[i], 0], vals[self.J[i], 0]])
            y = np.array([vals[self.I[i], 1], vals[self.J[i], 1]])
            z = np.array([vals[self.I[i], 2], vals[self.J[i], 2]])
            self.plots[i][0].set_xdata(x)
            self.plots[i][0].set_ydata(y)
            self.plots[i][0].set_3d_properties(z)
            self.plots[i][0].set_color(lcolor if len(self.LR) == len(self.I) else rcolor)

        r = .75
        xroot, yroot, zroot = vals[0, 0], vals[0, 1], vals[0, 2]
        if self.limits is None:

            self.ax.set_xlim3d([-r + xroot, r + xroot])
            self.ax.set_zlim3d([-r + zroot, r + zroot])
            self.ax.set_ylim3d([-r + yroot, r + yroot])
        else:
            r = 0.2
            self.ax.set_xlim3d(
                [-r + self.limits["x"][0], r + self.limits["x"][1]]
            )
            self.ax.set_ylim3d(
                [-r + self.limits["y"][0], r + self.limits["y"][1]]
            )
            self.ax.set_zlim3d(
                [-r + self.limits["z"][0], r + self.limits["z"][1]]
            )

        if not self.ticks:
            self.hide_axes()

def convert_to_3d(poses_as_angles, kinematic_tree, swap_yz):
    poses_3d = []
    for pose in poses_as_angles:
        poses_3d.append(
            fkl(
                pose,
                kinematic_tree["parent"],
                kinematic_tree["offset"],
                kinematic_tree["rotInd"],
                kinematic_tree["expmapInd"],
                kinematic_tree["posInd"]["ids"],
            )
        )

    poses_3d = np.stack(poses_3d, axis=0)
    poses_3d = poses_3d.reshape((poses_3d.shape[0], 32, -1))
    if swap_yz:
        poses_3d = poses_3d[:, :, [0, 2, 1]]

    return poses_3d

def project_onto_image_plane(
    poses,
    image_sizes,
    color,
    text:str,
    org:tuple,
    extrinsics,
    intrinsics,
    dataset: BaseDataset,
    target_size: tuple = None,
    background_color=None,
    synth_model: VunetOrg = None,  # ,thicken=False
    app_img=None,
    thickness=3,
    font_size=2,
    crop=False,
    target_width=128,
    yield_mask=False,
    overlay=False,
    cond_id=None
):
    """

    :param poses:
    :param image_sizes:
    :param extrinsics:
    :param intrinsics:
    :param color: the color in which to display the generated poses
    :param joint_model:
    :param target_size:
    :param app_img:
    :return:
    """
    if color == None:
        color = [[255, 0, 0],[255, 0, 0],[0, 0, 255],[0, 0, 255],[0, 0, 255],[255, 0, 0],[0, 0, 255],[0, 0, 255],[0, 0, 255],[0, 0, 255],[0, 0, 255],[255, 0, 0],[0, 0, 255],[0, 0, 255],[255, 0, 0],[255, 0, 0]]
        assert len(color) == len(dataset.joint_model.total_relative_joints)

    elif len(color)==2:
        color = [color[0], color[0], color[1], color[1],
                 color[1], color[0], color[1], color[1],
                 color[1], color[1], color[1], color[0],
                 color[1], color[1], color[0], color[0]]
        assert len(color) == len(dataset.joint_model.total_relative_joints)

    if overlay and len(color) != 2:
        color = [[255, 0, 0],[255, 0, 0],[0, 102, 255],[0, 102, 255],[0, 102, 255],[255, 0, 0],[0, 102, 255],[0, 102, 255],[0, 102, 255],[0, 102, 255],[0, 102, 255],[255, 0, 0],[0, 102, 255],[0, 102, 255],[255, 0,    0],[255, 0, 0]]
    yield_mask &= crop
    masks = None

    if background_color is None:
        background_color = 0

    if synth_model is not None:
        assert app_img is not None
        img_size = [synth_model.spatial_size, synth_model.spatial_size]
        size_arr = np.full((1, 2), synth_model.spatial_size,
                           dtype=np.float)
        dev = app_img.get_device() if app_img.get_device() >= 0 else "cpu"

    imgs = []
    imgs_rgb = []
    for id_count,(p, s) in enumerate(zip(poses, image_sizes)):
        pose_c = apply_affine_transform(p, extrinsics)
        pose_i = camera_projection(pose_c, intrinsics)

        img = add_joints_to_img(
            np.full(
                (s[0], s[1], 3), fill_value=background_color, dtype=np.uint8
            ),
            pose_i[dataset.joint_model.kps_to_use] if dataset.small_joint_model else pose_i,
            dataset.joint_model.total_relative_joints,
            color_kps=color,
            color_joints=color,
        )
        if cond_id is not None:
            if id_count < cond_id:
                text = text + "; COND"
            else:
                text = text + "; PRED"

        img = cv2.putText(img,text,org,cv2.FONT_HERSHEY_SIMPLEX,font_size,(0,0,0),thickness=thickness)
        if synth_model is not None:
            # rescale joints
            joint_scaling = size_arr / np.expand_dims(s, axis=0)
            pose_i_rescaled = pose_i * joint_scaling
            stickman = make_joint_img(
                img_size,
                pose_i_rescaled,
                dataset.joint_model,
                line_colors=dataset.line_colors,
                scale_factor=dataset.stickman_scale,
            )
            stickman = (
                dataset.stickman_transforms(stickman).unsqueeze(dim=0).to(dev)
            )
            with torch.no_grad():
                rgb_img = synth_model.transfer(app_img, stickman)
            rgb_img = (
                (scale_img(rgb_img) * 255.0)
                .permute((0, 2, 3, 1))
                .detach()
                .cpu()
                .numpy()
                .astype(np.uint8)
            )
            # rgb_img = cv2.putText(rgb_img, text, org, cv2.FONT_HERSHEY_SIMPLEX, font_size,
            #                   (255, 255, 255), thickness=thickness)
            if crop:
                kps = pose_i[
                    dataset.joint_model.kps_to_use] if dataset.small_joint_model else pose_i
                kps_rescaled = kps * joint_scaling
                center_kps = np.mean(kps_rescaled, axis=0)
                min_x = int(min(max(0, center_kps[0] - target_width // 2),
                                rgb_img.shape[1] - target_width))
                max_x = min_x + target_width
                rgb_img = rgb_img[:, :, min_x:max_x]
            imgs_rgb.append(rgb_img.squeeze(0))

        if target_size is not None:
            scaling = np.asarray(img.shape[:2],dtype=np.float) / np.asarray(target_size,dtype=np.float)
            img = cv2.resize(img, target_size,interpolation=cv2.INTER_NEAREST)
        else:
            scaling = np.ones(2,dtype=np.float)

        if crop:
            kps = pose_i[dataset.joint_model.kps_to_use] if dataset.small_joint_model else pose_i
            kps_rescaled = kps / np.expand_dims(scaling,axis=0)
            center_kps = np.mean(kps_rescaled,axis=0)
            min_x = int(min(max(0, center_kps[0] - target_width// 2),
                        img.shape[1] - target_width))
            max_x = min_x + target_width
            img = img[:,min_x:max_x]
            if masks is None:
                masks = np.expand_dims(np.asarray([min_x,max_x]),axis=0)
            else:
                masks = np.concatenate([masks,np.expand_dims(np.asarray([min_x,max_x]),axis=0)],axis=0)

        imgs.append(img)

    if synth_model is not None and overlay:
        overlays = []
        for stick,r in zip(imgs,imgs_rgb):
            o = np.where(np.expand_dims(np.all(stick==255,axis=2),axis=-1),r,stick)
            overlays.append(o)

    outs = [np.stack(imgs, axis=0)]
    if synth_model is not None:
        outs.append(np.stack(imgs_rgb, axis=0))
    else:
        outs.append(np.zeros_like(outs[0]))

    if yield_mask:
        outs.append(masks)

    if overlay:
        outs.append(np.stack(overlays,axis=0))

    return outs

max_stars_repo_name: CompVis/behavior-driven-video-synthesis

max_issues_repo_name: CompVis/behavior-driven-video-synthesis

max_forks_repo_name: CompVis/behavior-driven-video-synthesis


text:
import time
import random
import numpy as np
from collections import defaultdict
from optparse import OptionParser
import pickle

# Required libraries
import h5py
import keras
from keras.models import Sequential
from keras.layers.core import Dense, Dropout, Activation, Flatten
from keras.layers.convolutional import Convolution2D, MaxPooling2D
from keras.optimizers import SGD
from keras.utils.np_utils import to_categorical
from keras.layers.normalization import BatchNormalization as BN

# Plotting and stuffs (which probably won't work due to X11 issues)
import matplotlib.pylab as plt
import sys

# d-script imports
sys.path.append('..')
import data_iters
from data_iters.hdf5_iterator import Hdf5MiniBatcher
from data_iters.archive.iam_iterator import IAM_MiniBatcher
from data_iters.minibatcher import MiniBatcher

# Denoising stuff
from data_iters.CoffeeStainer import *
from data_iters.NoiseAdder import *

# Which training dataset do we want to train from?
train_dataset='iam-words'
noise_dataset='noise'

# All the images that you require extraction should be in this HDF5 file
if train_dataset=='nmec':
    hdf5authors='/memory/nmec_scaled_author_form.hdf5'
    hdf5authors='/fileserver/nmec-handwriting/nmec_scaled_author_form.hdf5'
    hdf5images='nmecdata/nmec_scaled_flat.hdf5'
elif train_dataset=='iam-words':
    hdf5authors='/fileserver/iam/iam-processed/words/author_words.hdf5'
elif train_dataset=='iam-lines':
    hdf5authors='/fileserver/iam/iam-processed/lines/author_lines.hdf5'
else:
    hdf5authors='/fileserver/iam/iam-processed/forms/author_forms.hdf5'

# Setup HDF5 Files
labels = h5py.File(hdf5authors, 'r')
num_authors=len(labels)
num_forms_per_author=-1
shingle_dim=(120,120)
batch_size=32
load_size=batch_size*1000
iterations = 1000
lr = 0.001
mini_m = Hdf5MiniBatcher(hdf5authors, num_authors, num_forms_per_author,
                            shingle_dim=shingle_dim, default_mode=MiniBatcher.TRAIN,
                            batch_size=batch_size, add_rotation=True)
noise_m = NoiseAdder('noisedir')

(X,Y) = mini_m.get_batch(10)
rng = np.random.RandomState(100)
Xnoise = noise_m.add_noise(X[0], rng=rng, shingle_dim=shingle_dim)


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

# This file was forked from https://github.com/marvis/pytorch-yolo2 ,
# licensed under the MIT License (see LICENSE.external for more details).

from __future__ import division

import dataset
from utils import *
from arg_utils import Yolov2OptionTraining
import os

import nnabla as nn
import nnabla.solvers as S

from train_graph import TrainGraph

from multiprocessing.pool import ThreadPool

import random
random.seed(0)
import numpy as np
np.seterr(**{key: 'raise' for key in ['divide', 'over', 'invalid']})

def get_current_batch(seen, args):
    '''
    Get the current batch number from number of processed examples.
    '''
    return int(seen // (args.batch_size * args.accum_times))

class PrefetchIterator(object):
    '''
    Creator of prefetch iterator from an iterable object.
    '''

    def __init__(self):
        self.pool = ThreadPool(1)

    def create(self, it):
        future_data = self.pool.apply_async(
            raise_info_thread(next), (it, None))
        while True:
            curr_data = future_data
            future_data = self.pool.apply_async(
                raise_info_thread(next), (it, None))
            data = curr_data.get()
            if data is None:
                break
            yield data

class YoloSolver(object):
    def __init__(self, args):

        # Create solvers (only Convolution kernels require weight decay).
        param_convweights = {
            k: v for k, v in nn.get_parameters().items() if k.endswith("conv/W")}
        param_others = {k: v for k, v in nn.get_parameters().items()
                        if not k.endswith("conv/W")}
        convweights = S.Momentum(args.learning_rate, args.momentum)
        others = S.Momentum(args.learning_rate, args.momentum)
        convweights.set_parameters(param_convweights)
        others.set_parameters(param_others)

        # Init parameter gradients.
        convweights.zero_grad()
        others.zero_grad()

        # Set attributes.
        self.convweights = convweights
        self.others = others
        self.args = args
        self.batch_size = args.batch_size * args.accum_times
        self.rate = args.accum_times
        self.count = 0

    def get_current_lr(self, seen):
        '''
        Get the current learning rate by the current number of processed examples.

        Burn-in and step learning rate scheduling are employed.
        '''
        args = self.args
        batch_num = get_current_batch(seen, args)
        lr = float(args.learning_rate)

        # Burn-in (burn_in=1000, burn_in_power=4 by default)
        if batch_num < args.burn_in:
            return lr * ((float(batch_num) / args.burn_in) ** args.burn_in_power)

        # Steps is only implemented
        for step, scale in zip(args.steps, args.scales):
            if step > batch_num:
                return lr
            lr *= scale
        return lr

    def update_at_rate(self, lr, force=False):
        '''

        Returns: int
            Current learning rate
        '''
        if force and self.count > 0:
            self.update(lr)
            return True
        self.count = (self.count + 1) % self.rate
        if self.count == 0:
            self.update(lr)
            return True
        return False

    def update(self, lr):
        # Get current learning rate and set to solvers
        self.convweights.set_learning_rate(lr / self.batch_size)
        self.others.set_learning_rate(lr / self.batch_size)

        # Update
        self.convweights.weight_decay(self.args.decay * self.batch_size)
        self.convweights.update()
        self.others.update()

        # Reset parameter gradients and accumulation count
        self.convweights.zero_grad()
        self.others.zero_grad()
        self.count = 0

def train(args, epoch, max_epochs, train_graph, yolo_solver,
          prefetch_iterator, on_memory_data):

    sample_iter = dataset.listDataset(
        args.train, args,
        shuffle=True,
        train=True,
        seen=train_graph.seen,
        image_sizes=args.size_aug,
        image_size_change_freq=args.batch_size * args.accum_times * 10,
        on_memory_data=on_memory_data,
        use_cv2=not args.disable_cv2)
    batch_iter = dataset.create_batch_iter(
        iter(sample_iter), batch_size=args.batch_size)

    total_loss = []
    epoch_seen = 0
    tic = time.time()
    for batch_idx, ((data_tensor, target_tensor), preprocess_time) \
            in enumerate(prefetch_iterator.create(batch_iter)):
        lr = yolo_solver.get_current_lr(train_graph.seen)
        nB = data_tensor.shape[0]
        print('size={}, '.format(sample_iter.shape[0]), end='')
        stats = train_graph.forward_backward(data_tensor, target_tensor)
        print('s=%d/%d, b=%d/%d, e=%d/%d, lr %g, mIoU %.3f, nGT %d, recall %d,'
              ' proposals %d, loss: %.3f, pre: %.1f [ms], exec: %.1f [ms]' % (
                  stats.seen, len(sample_iter), get_current_batch(
                      stats.seen, args),
                  args.max_batches, epoch, max_epochs, lr, stats.mIoU, stats.nGT,
                  stats.nCorrect, stats.nProposals, stats.loss,
                  preprocess_time, stats.time))
        epoch_seen += nB
        total_loss.append(float(stats.loss))

        # Update parameters for every accum_times iterations.
        updated = yolo_solver.update_at_rate(lr)

    # Update if gradients are computed from previous update.
    yolo_solver.update_at_rate(lr, force=True)
    print()
    if (epoch+1) % args.save_interval == 0:
        logging('save weights to %s/%06d.h5' % (args.output, epoch+1))
        nn.save_parameters('%s/%06d.h5' % (args.output, epoch+1))
    return np.sum(total_loss) / epoch_seen

def main():
    # Training settings
    args = Yolov2OptionTraining().parse_args()

    nsamples = file_lines(args.train)

    set_default_context_by_args(args)

    # Training parameters
    max_epochs = args.max_batches * args.batch_size * args.accum_times / nsamples + 1

    if not os.path.exists(args.output):
        os.mkdir(args.output)

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

    train_graph = TrainGraph(args, (args.size_aug[-1], args.size_aug[-1]))

    # Load parameters
    print("Load", args.weight, "...")
    nn.load_parameters(args.weight)
    yolo_solver = YoloSolver(args)

    if args.on_memory_data:
        on_memory_data = dataset.load_on_memory_data(args.train)
    else:
        on_memory_data = None

    prefetch_iterator = PrefetchIterator()

    from nnabla.monitor import Monitor, MonitorSeries, MonitorTimeElapsed
    monitor = Monitor(args.output)
    monitor_loss = MonitorSeries('Training loss', monitor, interval=1)
    monitor_time = MonitorTimeElapsed('Time per epoch', monitor, interval=1)

    # Epoch loop
    for epoch in range(0, int(max_epochs)):
        loss = train(args, epoch, max_epochs, train_graph,
                     yolo_solver, prefetch_iterator, on_memory_data)
        monitor_loss.add(epoch, loss)
        monitor_time.add(epoch)

    # Save the final parameters
    logging('save weights to %s/%06d.h5' % (args.output, epoch+1))
    nn.save_parameters('%s/%06d.h5' % (args.output, epoch+1))

if __name__ == '__main__':
    main()


text:
import torch
import numpy as np

from mmdet.core.bbox import BaseBBoxCoder
from mmdet.core.bbox.builder import BBOX_CODERS

@BBOX_CODERS.register_module()
class FVNetBBoxCoder(BaseBBoxCoder):
    """Bbox Coder for 3D boxes.

    Args:
        code_size (int): The dimension of boxes to be encoded.
    """
    # TODO: sine difference

    def __init__(self, prior_size=[1.6, 3.9, 1.56], code_size=7):
        super(FVNetBBoxCoder, self).__init__()
        self.prior_size = prior_size
        self.code_size = code_size

    @staticmethod
    def encode(points, boxes, prior_size):
        # normalize
        # size = torch.tensor(prior_size)
        # wa, la, ha = size
        # diagonal = torch.sqrt(wa**2 + la**2)
        # xt = (boxes[:, 0] - points[:, 0]) / diagonal
        # yt = (boxes[:, 1] - points[:, 1]) / diagonal
        # zt = (boxes[:, 2] + boxes[:, 5]/2 - points[:, 2]) / ha
        # wt = torch.log(boxes[:, 3] / prior_size[0])
        # lt = torch.log(boxes[:, 4] / prior_size[1])
        # ht = torch.log(boxes[:, 5] / prior_size[2])
        # rt_cos = torch.cos(boxes[:, 6])
        # rt_sin = torch.sin(boxes[:, 6])

        # not normalize
        xt = boxes[:, 0] - points[:, 0]
        yt = boxes[:, 1] - points[:, 1]
        zt = boxes[:, 2] + boxes[:, 5]/2 - points[:, 2]
        wt = torch.log(boxes[:, 3])
        lt = torch.log(boxes[:, 4])
        ht = torch.log(boxes[:, 5])
        rt_cos = torch.cos(boxes[:, 6])
        rt_sin = torch.sin(boxes[:, 6])

        return torch.stack([xt, yt, zt, wt, lt, ht, rt_cos, rt_sin]).T

    @staticmethod
    def decode(points, deltas, prior_size):
        # normalize
        # size = torch.tensor(prior_size)
        # wa, la, ha = size
        # diagonal = torch.sqrt(wa**2 + la**2)
        # xg = points[:, 0] + deltas[:, 0] * diagonal
        # yg = points[:, 1] + deltas[:, 1] * diagonal
        # zg = points[:, 2] + deltas[:, 2] * ha

        # wg = torch.exp(deltas[:, 3]) * prior_size[0]
        # lg = torch.exp(deltas[:, 4]) * prior_size[1]
        # hg = torch.exp(deltas[:, 5]) * prior_size[2]
        # zg = zg - hg / 2
        # rg = torch.atan2(deltas[:, 7], deltas[:, 6])

        # not normalize
        xg = points[:, 0] + deltas[:, 0]
        yg = points[:, 1] + deltas[:, 1]
        zg = points[:, 2] + deltas[:, 2]

        wg = torch.exp(deltas[:, 3])
        lg = torch.exp(deltas[:, 4])
        hg = torch.exp(deltas[:, 5])
        zg = zg - hg / 2
        rg = torch.atan2(deltas[:, 7], deltas[:, 6])
        return torch.stack([xg, yg, zg, wg, lg, hg, rg]).T


text:
#!python
# Copyright 2019 Vale
# voxelscheme
# v1.0 05/2019 paulo.ernesto

import numpy as np
import csv

# normalized voxel cube around all xyz points of a csv file
class Voxelschema(np.ma.MaskedArray):
  def __new__(cls, size, dimensions):
    if len(size) > 3:
      size = size[0:3]
    self = super().__new__(cls, np.ndarray(dimensions), np.ones(np.prod(dimensions), np.bool_))
    self._bl = np.array(size)
    self._o0 = np.zeros(3);
    return(self)

  # return the ijk of a voxel using our fixed voxel size instead of the voxel own size
  def ijk(self, xyz):
    #(block_centre, block_length) = Voxel.block(self._data.ix[c])
    #return tuple(int((xyz[i] - self._o0[i]) // self._bl[i]) for i in range(3))
    return np.divide(np.subtract(xyz, self._o0), self._bl).astype(np.int_)

  def xyz(self, ijk):
    return np.add(np.multiply(ijk, self._bl), self._o0)

  def block_extent(self, xyzc, xyzl):
    xyzr = np.multiply(xyzl, 0.5)
    return np.subtract(xyzc, xyzr), np.add(xyzc, xyzr)

  def flag_by_block(self, xyzc, xyzl, value, mode = 'first'):
    xyz0, xyz1 = self.block_extent(xyzc, xyzl)
    return self.flag_by_extent(xyz0, xyz1, value, mode)

  # flag grid cells that are touched by a block
  # in the case of multiple blocks touching a cell
  # use the block which intersects with the greater volume
  # input: block corners
  def flag_by_extent(self, xyz0, xyz1, value, mode = 'first'):
    ijk0 = self.ijk(xyz0)
    ijk1 = self.ijk(xyz1)
    # ijkd = np.subtract(ijk1, ijk0, dtype=np.int_, casting='unsafe')
    ijkd = np.subtract(ijk1, ijk0)
    for i in range(ijk0[0], ijk0[0] + ijkd[0]):
      for j in range(ijk0[1], ijk0[1] + ijkd[1]):
        for k in range(ijk0[2], ijk0[2] + ijkd[2]):
          if self[i,j,k] is np.ma.masked:
            pass
          elif mode == 'max' and self[i,j,k] < value:
            pass
          elif mode == 'min' and self[i,j,k] > value:
            pass
          else:
            continue

          self[i,j,k] = value

    return self

  def get_k_minmax(self, i, j, mode = 'max'):
    k = 0
    ks = self[i, j, :]
    if mode == 'max':
      k = ks.argmax()
    else:
      k = ks.argmin()
    return k

  def get_2d_minmax(self, mode = 'max'):
    grid = np.ma.MaskedArray(np.ndarray(self.shape[0:2]), np.ones(np.prod(self.shape[0:2]), np.bool_))
    for i in range(self.shape[0]):
      for j in range(self.shape[1]):
        k = self.get_k_minmax(i, j, mode)
        grid[i,j] = k

    return grid

  def save_2d_grid(self, path, grid):

    with open(path, 'w', newline='') as csvfile:
      csvwriter = csv.writer(csvfile, quoting=csv.QUOTE_MINIMAL)
      csvwriter.writerow(['i', 'j','value'])
      for i in range(grid.shape[0]):
        for j in range(grid.shape[1]):
          csvwriter.writerow([i, j, grid[i,j]])

  def save_grid(self, path):
    with open(path, 'w', newline='') as csvfile:
      csvwriter = csv.writer(csvfile, quoting=csv.QUOTE_MINIMAL)
      csvwriter.writerow(['i', 'j', 'k','value'])
      for i in range(self.shape[0]):
        for j in range(self.shape[1]):
          for k in range(self.shape[2]):
            csvwriter.writerow([i, j, k, self[i,j,k]])