text:
"""
Tests for making unstructured mesh slices

"""

#-----------------------------------------------------------------------------
# Copyright (c) 2013, yt Development Team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
import os
import os.path
import tempfile
import shutil
import numpy as np
import yt
from yt.testing import fake_tetrahedral_ds
from yt.testing import fake_hexahedral_ds, small_fake_hexahedral_ds
from yt.utilities.answer_testing.framework import \
    requires_ds, \
    data_dir_load, \
    GenericImageTest
from yt.utilities.lib.geometry_utils import triangle_plane_intersect
from yt.utilities.lib.mesh_triangulation import triangulate_indices

def setup():
    """Test specific setup."""
    from yt.config import ytcfg
    ytcfg["yt", "__withintesting"] = "True"

def compare(ds, field, test_prefix, decimals=12):
    def slice_image(filename_prefix):
        sl = yt.SlicePlot(ds, 'z', field)
        sl.set_log('all', False)
        image_file = sl.save(filename_prefix)
        return image_file

    slice_image.__name__ = "slice_{}".format(test_prefix)
    test = GenericImageTest(ds, slice_image, decimals)
    test.prefix = test_prefix
    return test

tri2 = "SecondOrderTris/RZ_p_no_parts_do_nothing_bcs_cone_out.e"

@requires_ds(tri2)
def test_tri2():
    ds = data_dir_load(tri2, kwargs={'step':-1})
    for field in ds.field_list:
        yield compare(ds, field, "answers_tri2_%s_%s" % (field[0], field[1]))

quad2 = "SecondOrderQuads/lid_driven_out.e"

@requires_ds(quad2)
def test_quad2():
    ds = data_dir_load(quad2, kwargs={'step':-1})
    field_list = [('all', 'T'), ('all', 'vel_x'), ('all', 'vel_y')]
    for field in field_list:
        yield compare(ds, field, "answers_quad2_%s_%s" % (field[0], field[1]))

multi_region = "MultiRegion/two_region_example_out.e"

@requires_ds(multi_region)
def test_multi_region():
    ds = data_dir_load(multi_region, kwargs={'step':-1})
    for field in ds.field_list:
        yield compare(ds, field, "answers_multi_region_%s_%s" % (field[0], field[1]))

def test_mesh_slices():
    # Perform I/O in safe place instead of yt main dir
    tmpdir = tempfile.mkdtemp()
    curdir = os.getcwd()
    os.chdir(tmpdir)

    np.random.seed(0x4d3d3d3)

    # tetrahedral ds
    ds = fake_tetrahedral_ds()

    for field in ds.field_list:
        for idir in [0, 1, 2]:
            sl = yt.SlicePlot(ds, idir, field)
            sl.annotate_mesh_lines()
            sl.save()

    # hexahedral ds
    ds = fake_hexahedral_ds()

    for field in ds.field_list:
        for idir in [0, 1, 2]:
            sl = yt.SlicePlot(ds, idir, field)
            sl.annotate_mesh_lines()
            sl.save()

    os.chdir(curdir)
    # clean up
    shutil.rmtree(tmpdir)

def test_perfect_element_intersection():
    # This test tests mesh line annotation where a z=0 slice
    # perfectly intersects the top of a hexahedral element with node
    # z-coordinates containing both -0 and +0. Before
    # https://github.com/yt-project/yt/pull/1437 this test falsely
    # yielded three annotation lines, whereas the correct result is four
    # corresponding to the four edges of the top hex face.

    ds = small_fake_hexahedral_ds()
    indices = ds.index.meshes[0].connectivity_indices
    coords = ds.index.meshes[0].connectivity_coords
    tri_indices = triangulate_indices(indices)
    tri_coords = coords[tri_indices]
    lines = triangle_plane_intersect(2, 0, tri_coords)
    non_zero_lines = 0
    for i in range(lines.shape[0]):
        norm = np.linalg.norm(lines[i][0] - lines[i][1])
        if norm > 1e-8:
            non_zero_lines += 1
    np.testing.assert_equal(non_zero_lines, 4)

max_stars_repo_name: kastalpes/yt

max_issues_repo_name: kastalpes/yt

max_forks_repo_name: kastalpes/yt


text:
# Copyright 2018 Amazon.com, Inc. or its affiliates. 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.
#   A copy of the License is located at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   or in the "license" file accompanying this file. This file 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 mxnet as mx
import numpy as np
import mxfusion as mf
import mxnet.gluon.nn as nn
from mxfusion.components.variables.var_trans import PositiveTransformation
from mxnet.gluon import HybridBlock
from ..components.distributions.random_gen import RandomGenerator
from ..components.variables import add_sample_dimension

def prepare_mxnet_array(array, array_has_samples, dtype):
    a_mx = mx.nd.array(array, dtype=dtype)
    if not array_has_samples:
        a_mx = add_sample_dimension(mx.nd, a_mx)
    return a_mx

def numpy_array_reshape(var, isSamples, n_dim):
    """
    Reshape a numpy to a give dimensionality by adding dimensions with size one in front (broadcasting). If *isSamples* is true, keep the first dimension.

    :param var: the variable to be reshaped.
    :type var: numpy.ndarray
    :param isSamples: whether the variable is a sampled variable with the first dimension being the number of samples.
    :type isSamples: boolean
    :param n_dim: the dimensionality that the variable is reshaped to.
    :type n_dim: int
    :returns: the reshaped numpy array
    :rtype: numpy.ndarray
    """
    if len(var.shape) < n_dim:
        if isSamples:
            var_reshape = var.reshape(
                *((var.shape[0],) + (1,) * (n_dim - len(var.shape)) +
                  var.shape[1:]))
        else:
            var_reshape = var.reshape(*((1,) * (n_dim - len(var.shape)) +
                                      var.shape))
    else:
        var_reshape = var
    return var_reshape

class MockMXNetRandomGenerator(RandomGenerator):
    """
    The MXNet pseudo-random number generator.
    """

    def __init__(self, samples):
        self._samples = samples

    def _sample_univariate(self, shape=None, out=None):
        if shape is None:
            shape = (1,)
        res = mx.nd.reshape(self._samples[:np.prod(shape)], shape=shape)
        if out is not None:
            out[:] = res
        return res

    def sample_normal(self, loc=0, scale=1, shape=None, dtype=None, out=None, ctx=None, F=None):
        return self._sample_univariate(shape=shape, out=out)

    def sample_multinomial(self, data, shape=None, get_prob=True, dtype=None, F=None):
        return mx.nd.reshape(self._samples[:np.prod(data.shape[:-1])], shape=data.shape[:-1])

    def sample_bernoulli(self, prob_true=0.5, shape=None, dtype=None, out=None, F=None):
        return mx.nd.reshape(self._samples[:np.prod(shape)], shape=shape)

    def sample_gamma(self, alpha=1, beta=1, shape=None, dtype=None, out=None, ctx=None, F=None):
        res = mx.nd.reshape(self._samples[:np.prod(shape)], shape=alpha.shape)
        if out is not None:
            out[:] = res
        return res

    def sample_uniform(self, low=0., high=1., shape=None, dtype=None, out=None, ctx=None, F=None):
        return self._sample_univariate(shape=shape, out=out)

    def sample_laplace(self, location=0., scale=1., shape=None, dtype=None, out=None, ctx=None, F=None):
        return self._sample_univariate(shape=shape, out=out)

def make_net():
    D=100
    net = nn.HybridSequential(prefix='testnet0_')
    with net.name_scope():
        net.add(nn.Dense(D, activation="tanh"))
        net.add(nn.Dense(D, activation="tanh"))
        net.add(nn.Dense(1, flatten=True))
    net.initialize(mx.init.Xavier(magnitude=3))
    return net

def make_basic_model(finalize=True, verbose=True):
    net = make_net()
    m = mf.models.Model(verbose=verbose)
    m.N = mf.components.Variable()
    m.x = mf.components.Variable(shape=(m.N,))
    m.v = mf.components.Variable(shape=(1,), transformation=PositiveTransformation())
    m.r = mf.components.Variable(shape=(m.N,))
    m.y = mf.components.distributions.Normal.define_variable(mean=m.r, variance=m.v, shape=(m.N,))
    return m

def make_bnn_model(finalize=True, verbose=True):
    net = make_net()
    m = mf.models.Model(verbose=verbose)
    m.N = mf.components.Variable()
    m.f = mf.components.MXFusionGluonFunction(net, num_outputs=1)
    m.x = mf.components.Variable(shape=(m.N,))
    m.v = mf.components.Variable(shape=(1,), transformation=PositiveTransformation())
    m.prior_variance = mf.components.Variable(shape=(1,), transformation=PositiveTransformation())
    m.r = m.f(m.x)
    for _, v in m.r.factor.parameters.items():
        v.set_prior(mf.components.distributions.Normal(mean=mx.nd.array([0]), variance=m.prior_variance))
    m.y = mf.components.distributions.Normal.define_variable(mean=m.r, variance=m.v, shape=(m.N,))
    return m

class DotProduct(HybridBlock):
    def hybrid_forward(self, F, x, *args, **kwargs):
        return F.dot(x, args[0])

class TestBlock(mx.gluon.HybridBlock):
    """
    Block with standard functions for initializing and running an MXNet Gluon block for unit testing.
    """

    def __init__(self):
        """
        Constructor.
        """

        super(TestBlock, self).__init__()
        with self.name_scope():
            self.var1 = self.params.get('var1', shape=(2,))
            self.var2 = self.params.get('var2', shape=(2,))

    def hybrid_forward(self, F, x, var1, var2):
        """
        Simple test function to enable MXFusion tests that use a Gluon block.
        :param F: MXNet computation type <mx.sym, mx.nd>
        :param x: MXNet dummy variable
        :param var1: MXNet dummy variable
        :param var2: MXNet dummy variable
        """

        return x + var1 + var2

def plot_univariate(samples, dist, buffer=0, **kwargs):
    """
    Visual inspection by plotting the distribution: plots a histogram of the samples along with

    :param samples: Samples from the distribution
    :type samples: (mx.nd.NDArray, np.ndarray)
    :param buffer: additional range to plot the distribution over
    :param dist: Distribution that these are samples from (scipy.stats)
    :param kwargs: Keyword arguments for the distribution (e.g. loc, scale)
    """
    if isinstance(samples, mx.nd.NDArray):
        samples = samples.asnumpy().ravel()
    elif isinstance(samples, np.ndarray):
        samples = samples.ravel()
    else:
        raise ValueError("Unexpected type for samples: {}, expected mx.nd.NDArray or np.ndarray".format(type(samples)))

    import matplotlib.pyplot as plt
    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111)
    ax.hist(samples, bins=301, density=True)
    x = np.linspace(samples.min() - buffer, samples.max() + buffer, num=301)
    ax.plot(x, dist.pdf(x, **kwargs).reshape(-1, 1))
    plt.show()

def plot_bivariate(samples, dist, buffer=0, **kwargs):
    """
    Visual inspection by plotting the distribution: plots a scatter plot of samples along with a contour plot

    :param samples: Samples from the distribution
    :type samples: (mx.nd.NDArray, np.ndarray)
    :param buffer: additional range to plot the distribution over
    :param dist: Distribution that these are samples from (scipy.stats)
    :param kwargs: Keyword arguments for the distribution (e.g. loc, scale)
    """
    if isinstance(samples, mx.nd.NDArray):
        samples = samples.asnumpy()
    elif isinstance(samples, np.ndarray):
        pass
    else:
        raise ValueError("Unexpected type for samples: {}, expected mx.nd.NDArray or np.ndarray".format(type(samples)))

    import matplotlib.pyplot as plt
    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111)
    s_min = samples.min(axis=0)
    s_max = samples.max(axis=0)
    x = np.linspace(s_min[0] - buffer, s_max[0] + buffer, num=301)
    y = np.linspace(s_min[1] - buffer, s_max[1] + buffer, num=301)
    x, y = np.meshgrid(x, y)
    z = dist.pdf(x=np.vstack([x.ravel(), y.ravel()]).T, **kwargs).reshape(x.shape)
    ax.contour(x, y, z, levels=10, linewidth=10)
    ax.scatter(samples[:, 0], samples[:, 1], alpha=0.05)
    plt.show()

def make_spd_matrix(dim):
    """
    Generate a random symmetric, positive-definite matrix.

    :param dim: The matrix dimension (matrix is square).
    :type dim: int
    :return X: The random symmetric, positive-definite matrix.
    :rtype: array of shape [n_dim, n_dim]
    """
    A = np.random.rand(dim, dim)
    U, s, V = np.linalg.svd(np.dot(A.T, A))
    X = np.dot(np.dot(U, 1.0 + np.diag(np.random.rand(dim))), V)
    return X

ext: py

max_stars_repo_name: JeremiasKnoblauch/MXFusion

max_issues_repo_name: JeremiasKnoblauch/MXFusion

max_forks_repo_name: JeremiasKnoblauch/MXFusion

reason: import numpy


text:
import copy as cp
import numpy as np

from skmultiflow.core import BaseSKMObject, ClassifierMixin, MetaEstimatorMixin
from skmultiflow.bayes import NaiveBayes

import warnings

def DynamicWeightedMajority(n_estimators=5, base_estimator=NaiveBayes(), period=50, beta=0.5,
                            theta=0.01):     # pragma: no cover
    warnings.warn("'DynamicWeightedMajority' has been renamed to 'DynamicWeightedMajorityClassifier' in v0.5.0.\n"
                  "The old name will be removed in v0.7.0", category=FutureWarning)
    return DynamicWeightedMajorityClassifier(n_estimators=n_estimators,
                                             base_estimator=base_estimator,
                                             period=period,
                                             beta=beta,
                                             theta=theta)

class DynamicWeightedMajorityClassifier(BaseSKMObject, ClassifierMixin, MetaEstimatorMixin):
    """ Dynamic Weighted Majority ensemble classifier.

    Parameters
    ----------
    n_estimators: int (default=5)
        Maximum number of estimators to hold.
    base_estimator: StreamModel or sklearn.BaseEstimator (default=NaiveBayes)
        Each member of the ensemble is an instance of the base estimator.
    period: int (default=50)
        Period between expert removal, creation, and weight update.
    beta: float (default=0.5)
        Factor for which to decrease weights by.
    theta: float (default=0.01)
        Minimum fraction of weight per model.

    Notes
    -----
    The dynamic weighted majority (DWM) [1]_, uses four mechanisms to
    cope with concept drift: It trains online learners of the ensemble,
    it weights those learners based on their performance, it removes them,
    also based on their performance, and it adds new experts based on the
    global performance of the ensemble.

    References
    ----------
    .. [1] Kolter and Maloof. Dynamic weighted majority: An ensemble method
       for drifting concepts. The Journal of Machine Learning Research,
       8:2755-2790, December 2007. ISSN 1532-4435.

    Examples
    --------
    .. code-block:: python

       # Imports
       from skmultiflow.data import SEAGenerator
       from skmultiflow.meta import DynamicWeightedMajorityClassifier
       from skmultiflow.bayes import NaiveBayes

       # Setup a data stream
       stream = SEAGenerator(random_state=1)
       stream.prepare_for_use()

       # Setup Dynamic Weighted Majority Ensemble Classifier
       dynamic_weighted_majority_classifier = DynamicWeightedMajorityClassifier()

       # Setup varibles to control loop and track performance
       n_samples = 0
       correct_cnt = 0
       max_samples = 200

       # Train the classifier with the samples provided by the data stream
       while n_samples < max_samples and stream.has_more_samples():
           X, y = stream.next_sample()
           y_pred = dynamic_weighted_majority_classifier.predict(X)
           if y[0] == y_pred[0]:
               correct_cnt += 1
           dynamic_weighted_majority_classifier = dynamic_weighted_majority_classifier.partial_fit(X, y)
           n_samples += 1

       # Display results
       print('{} samples analyzed.'.format(n_samples))
       print('Dynamic Weighted Majority Classifier: {}'.format(correct_cnt / n_samples))
    """

    class WeightedExpert:
        """
        Wrapper that includes an estimator and its weight.

        Parameters
        ----------
        estimator: StreamModel or sklearn.BaseEstimator
            The estimator to wrap.
        weight: float
            The estimator's weight.
        """
        def __init__(self, estimator, weight):
            self.estimator = estimator
            self.weight = weight

    def __init__(self, n_estimators=5, base_estimator=NaiveBayes(),
                 period=50, beta=0.5, theta=0.01):
        """
        Creates a new instance of DynamicWeightedMajorityClassifier.
        """
        super().__init__()

        self.n_estimators = n_estimators
        self.base_estimator = base_estimator

        self.beta = beta
        self.theta = theta
        self.period = period

        # Following attributes are set later
        self.epochs = None
        self.num_classes = None
        self.experts = None

        self.reset()

    def partial_fit(self, X, y, classes=None, sample_weight=None):
        """ Partially fits the model on the supplied X and y matrices.

        Since it's an ensemble learner, if X and y matrix of more than one
        sample are passed, the algorithm will partial fit the model one sample
        at a time.

        Parameters
        ----------
        X : numpy.ndarray of shape (n_samples, n_features)
            The features to train the model.

        y: numpy.ndarray of shape (n_samples)
            An array-like with the class labels of all samples in X.

        classes: numpy.ndarray, optional (default=None)
            Array with all possible/known class labels. This is an optional parameter, except
            for the first partial_fit call where it is compulsory.

        sample_weight: numpy.ndarray of shape (n_samples), optional (default=None)
            Samples weight. If not provided, uniform weights are assumed. Usage varies depending on the base estimator.

        Returns
        -------
        DynamicWeightedMajorityClassifier
            self
        """
        for i in range(len(X)):
            self.fit_single_sample(
                X[i:i+1, :], y[i:i+1], classes, sample_weight
            )
        return self

    def predict(self, X):
        """ predict

        The predict function will take an average of the predictions of its
        learners, weighted by their respective weights, and return the most
        likely class.

        Parameters
        ----------
        X: numpy.ndarray of shape (n_samples, n_features)
            A matrix of the samples we want to predict.

        Returns
        -------
        numpy.ndarray
            A numpy.ndarray with the label prediction for all the samples in X.
        """
        preds = np.array([np.array(exp.estimator.predict(X)) * exp.weight
                          for exp in self.experts])
        sum_weights = sum(exp.weight for exp in self.experts)
        aggregate = np.sum(preds / sum_weights, axis=0)
        return (aggregate + 0.5).astype(int)    # Round to nearest int

    def predict_proba(self, X):
        raise NotImplementedError

    def fit_single_sample(self, X, y, classes=None, sample_weight=None):
        """
        Fits a single sample of shape `X.shape=(1, n_attributes)` and `y.shape=(1)`

        Aggregates all experts' predictions, diminishes weight of experts whose
        predictions were wrong, and may create or remove experts every _period_
        samples.

        Finally, trains each individual expert on the provided data.

        Train loop as described by Kolter and Maloof in the original paper.

        Parameters
        ----------
        X: numpy.ndarray of shape (n_samples, n_features)
            Features matrix used for partially updating the model.

        y: Array-like
            An array-like of all the class labels for the samples in X.

        classes: list
            List of all existing classes. This is an optional parameter.

        sample_weight: numpy.ndarray of shape (n_samples), optional (default=None)
            Samples weight. If not provided, uniform weights are assumed. Applicability
            depends on the base estimator.

        """
        self.epochs += 1
        self.num_classes = max(
            len(classes) if classes is not None else 0,
            (int(np.max(y)) + 1), self.num_classes)
        predictions = np.zeros((self.num_classes,))
        max_weight = 0
        weakest_expert_weight = 1
        weakest_expert_index = None

        for i, exp in enumerate(self.experts):
            y_hat = exp.estimator.predict(X)
            if np.any(y_hat != y) and (self.epochs % self.period == 0):
                exp.weight *= self.beta

            predictions[y_hat] += exp.weight
            max_weight = max(max_weight, exp.weight)

            if exp.weight < weakest_expert_weight:
                weakest_expert_index = i
                weakest_expert_weight = exp.weight

        y_hat = np.array([np.argmax(predictions)])
        if self.epochs % self.period == 0:
            self._scale_weights(max_weight)
            self._remove_experts()
            if np.any(y_hat != y):
                if len(self.experts) == self.n_estimators:
                    self.experts.pop(weakest_expert_index)
                self.experts.append(self._construct_new_expert())

        # Train individual experts
        for exp in self.experts:
            exp.estimator.partial_fit(X, y, classes, sample_weight)

    def get_expert_predictions(self, X):
        """
        Returns predictions of each class for each expert.
        In shape: (n_experts, n_samples)
        """
        return [exp.estimator.predict(X) for exp in self.experts]

    def reset(self):
        """
        Reset this ensemble learner.
        """
        self.epochs = 0
        self.num_classes = 2    # Minimum of 2 classes
        self.experts = [
            self._construct_new_expert()
        ]

    def _scale_weights(self, max_weight):
        """
        Scales the experts' weights such that the max is 1.
        """
        scale_factor = 1 / max_weight
        for exp in self.experts:
            exp.weight *= scale_factor

    def _remove_experts(self):
        """
        Removes all experts whose weight is lower than self.theta.
        """
        self.experts = [ex for ex in self.experts if ex.weight >= self.theta]

    def _construct_new_expert(self):
        """
        Constructs a new WeightedExpert from the provided base_estimator.
        """
        return self.WeightedExpert(cp.deepcopy(self.base_estimator), 1)

max_stars_repo_name: imran-salim/scikit-multiflow

max_issues_repo_name: imran-salim/scikit-multiflow

max_forks_repo_name: imran-salim/scikit-multiflow


text:
import os
import re

import cv2
import lmdb  # install lmdb by "pip install lmdb"
import numpy as np
from pathlib import Path

home = str(Path.home())

def checkImageIsValid(imageBin):
    if imageBin is None:
        return False
    imageBuf = np.fromstring(imageBin, dtype=np.uint8)
    img = cv2.imdecode(imageBuf, cv2.IMREAD_GRAYSCALE)
    imgH, imgW = img.shape[0], img.shape[1]
    if imgH * imgW == 0:
        return False
    return True

def writeCache(env, cache):
    with env.begin(write=True) as txn:
        for k, v in cache.items():
            if type(v) == str:
                v = v.encode()
            txn.put(k.encode(), v)

def createDataset(outputPath, imagePathList, labelList, lexiconList=None, checkValid=True):
    """
    Create LMDB dataset for CRNN training.

    ARGS:
        outputPath    : LMDB output path
        imagePathList : list of image path
        labelList     : list of corresponding groundtruth texts
        lexiconList   : (optional) list of lexicon lists
        checkValid    : if true, check the validity of every image
    """
    assert (len(imagePathList) == len(labelList))
    nSamples = len(imagePathList)
    env = lmdb.open(outputPath, map_size=1099511627776)
    cache = {}
    cnt = 1
    for i in range(nSamples):
        imagePath = imagePathList[i]
        label = labelList[i]
        if not os.path.exists(imagePath):
            print('%s does not exist' % imagePath)
            continue
        with open(imagePath, 'rb') as f:
            imageBin = f.read()
        if checkValid:
            if not checkImageIsValid(imageBin):
                print('%s is not a valid image' % imagePath)
                continue

        imageKey = 'image-%09d' % cnt
        labelKey = 'label-%09d' % cnt
        cache[imageKey] = imageBin
        cache[labelKey] = label
        if lexiconList:
            lexiconKey = 'lexicon-%09d' % cnt
            cache[lexiconKey] = ' '.join(lexiconList[i])
        if cnt % 1000 == 0:
            writeCache(env, cache)
            cache = {}
            print('Written %d / %d' % (cnt, nSamples))
        cnt += 1
    nSamples = cnt - 1
    cache['num-samples'] = str(nSamples)
    writeCache(env, cache)
    print('Created dataset with %d samples' % nSamples)

def create(data_path, file, path, imagePathList, labelList, vocb=None):
    file = os.path.join(data_path, file)
    with open(file) as f:
        for l in f:
            tokens = l.split(',')
            if len(tokens) < 2:
                continue
            img = tokens[0].strip()
            txt = ''.join(tokens[1:]).strip()
            txt = re.sub('[|]', '', txt)

            path = os.path.join(data_path, path)
            img = os.path.join(path, img) + '.jpg'
            imagePathList.append(img)
            labelList.append(txt)
            if vocb is not None:
                vocb.update(txt)

def COCO(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb):
    data_path = os.path.join(home, 'data/ocr_data/coco')
    create(data_path, 'train_words_gt.txt', 'train_words', trainImagePathList, trainLabelList, vocb)
    create(data_path, 'val_words_gt.txt', 'val_words', valImagePathList, valLabelList)

def Born_Digital(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb):
    data_path = '../dataset/train'
    file = os.path.join(data_path, 'gt.txt')
    imagePathList = []
    labelList = []
    with open(file) as f:
        for l in f:
            tokens = l.split(',')
            img = tokens[0].strip()
            img = os.path.join(data_path, img)
            txt = ''.join(tokens[1:]).strip()
            txt = re.sub('["]', '', txt)

            imagePathList.append(img)
            labelList.append(txt)
            vocb.update(txt)
            # print(img, txt)

    split = int(len(imagePathList) * 0.8)
    trainImagePathList.extend(imagePathList[:split])
    trainLabelList.extend(labelList[:split])
    valImagePathList.extend(imagePathList[split:])
    valLabelList.extend(labelList[split:])

def Gen_Handwritten(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb):
    data_path = os.path.join(home, 'code/handwriting-generation')
    gt_file = os.path.join(data_path, 'gt/ALL')
    img_dir = os.path.join(data_path, 'gen')

    image_path_list = []
    label_list = []
    with open(gt_file) as f:
        for l in f:
            tokens = l.split(',')
            img = tokens[0].strip() + '.png'
            img = os.path.join(img_dir, img)
            txt = ''.join(tokens[1:]).strip()

            if not os.path.isfile(img):
                print('file {} does not exist!'.format(img))
                continue

            image_path_list.append(img)
            label_list.append(txt)
            vocb.update(txt)

            if len(image_path_list) > 50000:
                break

    trainImagePathList.extend(image_path_list)
    trainLabelList.extend(label_list)

    # split = int(len(imagePathList) * 0.95)
    # trainImagePathList.extend(imagePathList[:split])
    # trainLabelList.extend(labelList[:split])
    # valImagePathList.extend(imagePathList[split:])
    # valLabelList.extend(labelList[split:])

def IAM(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb):
    data_path = os.path.join(home, 'data/ocr_data/IAM')
    # img_dir = os.path.join(data_path, 'gen')
    img_dir = os.path.join(data_path, 'word')
    gt_file = os.path.join(data_path, 'words.txt')

    imagePathList = []
    labelList = []
    with open(gt_file) as f:
        for line in f:
            # ignore comment line
            if not line or line[0] == '#':
                continue

            line_split = line.strip().split(' ')
            assert len(line_split) >= 9

            # filename: part1-part2-part3 --> part1/part1-part2/part1-part2-part3.png
            file_name_split = line_split[0].split('-')
            img = '{}/{}-{}/{}.png'.format(file_name_split[0], file_name_split[0], file_name_split[1], line_split[0])
            img = os.path.join(img_dir, img)

            # GT text are columns starting at 9
            txt = ' '.join(line_split[8:])

            # check if image is not empty
            if not os.path.getsize(img):
                print('not found img file: {}'.format(img))
                continue

            imagePathList.append(img)
            labelList.append(txt)
            vocb.update(txt)

    split = int(len(imagePathList) * 0.9)
    trainImagePathList.extend(imagePathList[:split])
    trainLabelList.extend(labelList[:split])
    valImagePathList.extend(imagePathList[split:])
    valLabelList.extend(labelList[split:])

def Gen_Font(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb):
    data_path = os.path.join(home, 'data/ocr_data/font_gen1')
    gt_file = os.path.join(data_path, 'gt/ALL.txt')
    img_dir = os.path.join(data_path, 'img')

    image_path_list = []
    label_list = []
    with open(gt_file) as f:
        for l in f:
            tokens = l.split(',')
            img = tokens[0].strip() + '.png'
            img = os.path.join(img_dir, img)
            txt = ''.join(tokens[1:]).strip()

            if not os.path.isfile(img):
                print('file {} does not exist!'.format(img))
                continue

            image_path_list.append(img)
            label_list.append(txt)
            vocb.update(txt)

    split = int(len(image_path_list) * 0.95)
    trainImagePathList.extend(image_path_list[:split])
    trainLabelList.extend(label_list[:split])
    valImagePathList.extend(image_path_list[split:])
    valLabelList.extend(label_list[split:])

if __name__ == '__main__':
    vocb = set()
    trainImagePathList = []
    trainLabelList = []
    valImagePathList = []
    valLabelList = []

    # COCO(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb)
    # Born_Digital(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb)
    # Gen_Handwritten(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb)
    IAM(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb)
    Gen_Font(trainImagePathList, trainLabelList, valImagePathList, valLabelList, vocb)

    print('train img: {}, label: {}'.format(len(trainImagePathList), len(trainLabelList)))
    print('val img: {}, label: {}'.format(len(valImagePathList), len(valLabelList)))
    vocb = ''.join(vocb)
    vocb = sorted(vocb)
    print('vocb {}: {}'.format(len(vocb), ''.join(vocb)))

    createDataset('train2', trainImagePathList, trainLabelList)
    createDataset('val2', valImagePathList, valLabelList)


text:
"""
.. _extending_pyvista_example:

Extending PyVista
~~~~~~~~~~~~~~~~~

A :class:`pyvista.DataSet`, such as :class:`pyvista.PolyData`, can be extended
by users.  For example, if the user wants to keep track of the location of the
maximum point in the (1, 0, 1) direction on the mesh.

There are two methods by which users can handle subclassing.  One is directly managing
the types objects.  This may require checking types during filter
operations.

The second is automatic managing of types.  Users can control whether user defined
classes are nearly always used for particular types of DataSets.

.. note::
    This is for advanced usage only.  Automatic managing of types
    will not work in all situations, in particular when a builtin dataset is directly
    instantiated.  See examples below.

"""

import numpy as np
import vtk

import pyvista

pyvista.set_plot_theme("document")

###############################################################################
# A user defined subclass of :class:`pyvista.PolyData`, ``FooData`` is defined.
# It includes a property to keep track of the point on the mesh that is
# furthest along in the (1, 0, 1) direction.

class FooData(pyvista.PolyData):
    @property
    def max_point(self):
        """Returns index of point that is furthest along (1, 0, 1) direction."""
        return np.argmax(np.dot(self.points, (1.0, 0.0, 1.0)))

###############################################################################
# Directly Managing Types
# +++++++++++++++++++++++
#
# Now a ``foo_sphere`` object is created of type ``FooData``.
# The index of the point and location of the point of interest can be obtained
# directly. The sphere has a radius of 0.5, so the maximum extent in the
# direction (1, 0, 1) is :math:`0.5\sqrt{0.5}\approx0.354`
#

foo_sphere = FooData(pyvista.Sphere(theta_resolution=100, phi_resolution=100))
print("Original foo sphere:")
print(f"Type: {type(foo_sphere)}")
print(f"Maximum point index: {foo_sphere.max_point}")
print(f"Location of maximum point: {foo_sphere.points[foo_sphere.max_point, :]}")

###############################################################################
# Using an inplace operation like :func:`pyvista.DataSet.rotate_y` does not
# affect the type of the object.

foo_sphere.rotate_y(90, inplace=True)
print("\nRotated foo sphere:")
print(f"Type: {type(foo_sphere)}")
print(f"Maximum point index: {foo_sphere.max_point}")
print(f"Location of maximum point: {foo_sphere.points[foo_sphere.max_point, :]}")

###############################################################################
# However, filter operations can return different ``DataSet`` types including
# ones that differ from the original type.  In this case, the
# :func:`decimate <pyvista.PolyDataFilters.decimate>` method returns a
# :class:`pyvista.PolyData` object.

print("\nDecimated foo sphere:")
decimated_foo_sphere = foo_sphere.decimate(0.5)
print(f"Type: {type(decimated_foo_sphere)}")

###############################################################################
# It is now required to explicitly wrap the object into ``FooData``.

decimated_foo_sphere = FooData(foo_sphere.decimate(0.5))
print(f"Type: {type(decimated_foo_sphere)}")
print(f"Maximum point index: {decimated_foo_sphere.max_point}")
print(f"Location of maximum point: {foo_sphere.points[foo_sphere.max_point, :]}")

###############################################################################
# Automatically Managing Types
# ++++++++++++++++++++++++++++
#
# The default :class:`pyvista.DataSet` type can be set using ``pyvista._wrappers``.
# In general, it is best to use this method when it is expected to primarily
# use the user defined class.
#
# In this example, all objects that would have been created as
# :class:`pyvista.PolyData` would now be created as a ``FooData`` object. Note,
# that the key is the underlying vtk object.

pyvista._wrappers['vtkPolyData'] = FooData

###############################################################################
# It is no longer necessary to specifically wrap :class:`pyvista.PolyData`
# objects to obtain a ``FooData`` object.

foo_sphere = pyvista.Sphere(theta_resolution=100, phi_resolution=100)
print("Original foo sphere:")
print(f"Type: {type(foo_sphere)}")
print(f"Maximum point index: {foo_sphere.max_point}")
print(f"Location of maximum point: {foo_sphere.points[foo_sphere.max_point, :]}")

###############################################################################
# Using an inplace operation like :func:`rotate_y <pyvista.DataSet.rotate_y>` does not
# affect the type of the object.

foo_sphere.rotate_y(90, inplace=True)
print("\nRotated foo sphere:")
print(f"Type: {type(foo_sphere)}")
print(f"Maximum point index: {foo_sphere.max_point}")
print(f"Location of maximum point: {foo_sphere.points[foo_sphere.max_point, :]}")

###############################################################################
# Filter operations that return :class:`pyvista.PolyData` now return
# ``FooData``

print("\nDecimated foo sphere:")
decimated_foo_sphere = foo_sphere.decimate(0.5)
print(f"Type: {type(decimated_foo_sphere)}")
print(f"Maximum point index: {decimated_foo_sphere.max_point}")
print(f"Location of maximum point: {foo_sphere.points[foo_sphere.max_point, :]}")

###############################################################################
# Users can still create a native :class:`pyvista.PolyData` object, but
# using this method may incur unintended consequences.  In this case,
# it is recommended to use the directly managing types method.

poly_object = pyvista.PolyData(vtk.vtkPolyData())
print(f"Type: {type(poly_object)}")
# catch error
try:
    poly_object.rotate_y(90, inplace=True)
except TypeError:
    print("This operation fails")

###############################################################################
# Usage of ``pyvista._wrappers`` may require resetting the default value
# to avoid leaking the setting into cases where it is unused.

pyvista._wrappers['vtkPolyData'] = pyvista.PolyData

###############################################################################
# For instances where a localized usage is preferred, a tear-down method is
# recommended.  One example is a ``try...finally`` block.

try:
    pyvista._wrappers['vtkPolyData'] = FooData
    # some operation that sometimes raises an error
finally:
    pyvista._wrappers['vtkPolyData'] = pyvista.PolyData

max_stars_repo_name: basnijholt/pyvista

max_issues_repo_name: basnijholt/pyvista

max_forks_repo_name: basnijholt/pyvista


text:
"""
------------------------------------------------------------------
FILE:   _cluster_profiler.py
AUTHOR: Edo Altamura
DATE:   12-11-2019
------------------------------------------------------------------
This file is an extension of the cluster.Cluster class. It provides
class methods for performing basic analysis procedures on C-EAGLE
data from the /cosma5 data system.
This file contains a mixin class, affiliated to cluster.Cluster.
Mixins are classes that have no data of their own — only methods —
so although you inherit them, you never have to call super() on them.
They working principle is based on OOP class inheritance.
-------------------------------------------------------------------
"""

import numpy as np
from unyt import hydrogen_mass, boltzmann_constant, gravitational_constant, parsec, solar_mass
from .memory import free_memory
import warnings

# Delete the units from Unyt constants
hydrogen_mass = float(hydrogen_mass.value)
boltzmann_constant = float(boltzmann_constant.value)
G_SI = float(gravitational_constant.value)
G_astro = float(gravitational_constant.in_units('(1e6*pc)*(km/s)**2/(1e10*solar_mass)').value)
parsec = float((1*parsec).in_units('m').value)
solar_mass = float(solar_mass.value)

class Mixin:

    @staticmethod
    def angle_between_vectors(v1, v2):
        """
        :param v1: v1 is your first vector
        :param v2: v2 is your second vector
        :return: Returns the cosine of the angle
        """
        v1 = np.array(v1)
        v2 = np.array(v2)
        # If the vectors are almost equal, then the arccos returns nan because of numerical
        # precision: catch that case and return 0 angle instead.
        if np.abs(np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)) - 1) < 1e-14:
            angle = 1.
        elif np.abs(np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)) + 1) < 1e-14:
            angle = -1.
        else:
            angle = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))
        return angle

    def radial_distance_CoP(self, coords):
        coords = np.asarray(coords)
        coordinates_radial = np.sqrt(
            (coords[:, 0] - self.centre_of_potential[0]) ** 2 +
            (coords[:, 1] - self.centre_of_potential[1]) ** 2 +
            (coords[:, 2] - self.centre_of_potential[2]) ** 2
        )
        return coordinates_radial

    @staticmethod
    def kinetic_energy(mass, vel):
        mass = np.asarray(mass)
        vel = np.asarray(vel)
        ke = 0.5 * mass * np.linalg.norm(vel, axis = 1)**2
        return np.sum(ke)

    @staticmethod
    def thermal_energy(mass, temperature):
        mass = np.asarray(mass)
        temperature = np.asarray(temperature)
        te = 1.5 * boltzmann_constant * temperature * mass / (hydrogen_mass / 1.16)
        return np.sum(te)

    @staticmethod
    def centre_of_mass(mass, coords):
        """
        AIM: reads the FoF group central of mass from the path and file given
        RETURNS: type = np.array of 3 doubles
        ACCESS DATA: e.g. group_CoM[0] for getting the x value
        """
        mass   = np.asarray(mass)
        coords = np.asarray(coords)
        return np.sum(coords*mass[:, None], axis = 0)/np.sum(mass)

    @staticmethod
    def zero_momentum_frame(mass, velocity):
        """
        AIM: reads the FoF group central of mass from the path and file given
        RETURNS: type = np.array of 3 doubles
        """
        mass     = np.asarray(mass)
        velocity = np.asarray(velocity)
        return np.sum(velocity*mass[:, None], axis = 0)/np.sum(mass)

    @staticmethod
    def angular_momentum(mass, coords, velocity):
        """
        Compute angular momentum as r [cross] mv.

        :param mass: the mass array of the particles
        :param coords: the coordinates of the particles, rescaled to the origin of reference
            Usually, coords are rescaled with respect to the centre of potential.
        :param velocity:  the velocity array of the particles, rescaled to the cluster's
            rest frame. I/e/ take out the bulk peculiar velocity to isolate the rotation.
        :return: np.array with the 3D components of the angular momentum vector.
        """
        mass = np.asarray(mass)
        coords = np.asarray(coords)
        velocity = np.asarray(velocity)
        return np.sum(np.cross(coords, velocity*mass[:, None]), axis=0)

    @staticmethod
    def inertia_tensor(mass: np.ndarray, coords: np.ndarray) -> np.ndarray:
        """
		Compute the moment of inertia tensor in matrix form:

		            I_xx    I_xy    I_xz
		            I_yx    I_yy    I_yz
		            I_zx    I_zy    I_zz

		:param mass: the mass array of the particles
		:param coords: the coordinates of the particles, rescaled to the origin of reference
			Usually, coords are rescaled with respect to the centre of potential.

		:return: np.array with the 3x3 component inertia tensor.
		"""
        m = np.asarray(mass)
        coords = np.asarray(coords)
        x = coords[:, 0]
        y = coords[:, 1]
        z = coords[:, 2]
        del coords
        I_xx = np.sum(m*(y**2+z**2))
        I_xy = - np.sum(m*x*y)
        I_xz = - np.sum(m*x*z)
        I_yx = I_xy
        I_yy = np.sum(m*(x**2+z**2))
        I_yz = - np.sum(m*y*z)
        I_zx = I_xz
        I_zy = I_yz
        I_zz = np.sum(m*(x**2+y**2))
        return np.array([
                [I_xx, I_xy, I_xz],
                [I_yx, I_yy, I_yz],
                [I_zx, I_zy, I_zz],
        ])

    @staticmethod
    def principal_axes_ellipsoid(inertia_tensor: np.ndarray, eigenvalues: bool = False) -> np.ndarray:
        """
        Compute the eigenvalues and eigenvectors of the inertia tensor for morphology studies
        :param inertia_tensor:
        :return:
        """
        eigensolution = np.linalg.eig(inertia_tensor)
        return eigensolution if eigenvalues else eigensolution[1]

    def generate_apertures(self):
        physical = [0.1, 0.2, 0.4, 0.8]
        manual = [
                0.1*self.r500,
                self.r2500,
                1.5*self.r2500,
                self.r500
        ]
        auto = np.logspace(np.log10(self.r200), np.log10(5 * self.r200), 15).tolist()
        all_apertures = physical+manual+auto
        return np.asarray(all_apertures)

    def group_thermal_energy(self,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the thermal_energy of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the total thermal_energy.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall thermal_energy, or return the of each particle type separately.
        This toggle is controlled by a boolean.

        The thermal_energy is computed after converting the datasets from Gadget units to
        SI units ([mass] = kg, [velocity] = km/s). The output is given in units of
        (10**46) Joules.

        thermal_energy = 3/2 * k_B * T * (number of H atoms)
        number of H atoms = mass

        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        part_type = '0'
        assert hasattr(self, f'partType{part_type}_coordinates')
        assert hasattr(self, f'partType{part_type}_temperature')
        assert hasattr(self, f'partType{part_type}_mass')
        radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
        aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
        free_memory(['radial_dist'])
        mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
        temperature = getattr(self, f'partType{part_type}_temperature')[aperture_radius_index]
        if mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
        if temperature.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
        mass = self.mass_units(mass, unit_system='SI')
        return self.thermal_energy(mass, temperature)*np.power(10., -46)

    def group_kinetic_energy(self,
                             out_allPartTypes: bool =False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the kinetic_energy of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the total kinetic_energy.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall kinetic_energy, or return the of each particle type separately.
        This toggle is controlled by a boolean.

        The kinetic_energy is computed after converting the datasets from Gadget units to
        SI units ([mass] = kg, [velocity] = km/s). The output is given in units of
        (10**46) Joules.

        kinetic_energy = 1/2 * mass * || velocity ||^2.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            bulk_velocity = self.group_zero_momentum_frame(aperture_radius=aperture_radius)
            kinetic_energy_PartTypes = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_velocity')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass     = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _velocity = getattr(self, f'partType{part_type}_velocity')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _velocity.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                _velocity = np.subtract(_velocity, bulk_velocity)
                _mass     = self.mass_units(_mass, unit_system='SI')
                _velocity = self.velocity_units(_velocity, unit_system='SI')
                kinetic_energy = self.kinetic_energy(_mass, _velocity)
                kinetic_energy_PartTypes = np.append(kinetic_energy_PartTypes, kinetic_energy)

            return kinetic_energy_PartTypes*np.power(10., -46)

        else:
            bulk_velocity = self.group_zero_momentum_frame(aperture_radius=aperture_radius)
            mass = np.zeros(0, dtype=np.float)
            velocity = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_velocity')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass     = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _velocity = getattr(self, f'partType{part_type}_velocity')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _velocity.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                mass = np.concatenate((mass, _mass), axis=0)
                velocity = np.concatenate((velocity, _velocity), axis=0)

            velocity = np.subtract(velocity, bulk_velocity)
            mass     = self.mass_units(mass, unit_system='SI')
            velocity = self.velocity_units(velocity, unit_system='SI')
            return self.kinetic_energy(mass, velocity)*np.power(10., -46)

    def group_mass_aperture(self,
                             out_allPartTypes: bool =False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the mass of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the total particle mass.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall mass, ot return the of each particle type separately.
        This toggle is controlled by a boolean.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            mass_PartTypes = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")

                sum_of_masses = np.sum(_mass)
                mass_PartTypes = np.append(mass_PartTypes, sum_of_masses)

            return mass_PartTypes

        else:

            mass   = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                mass = np.append(mass, _mass)

            return np.sum(mass)

    def group_substructure_mass(self,
                             out_allPartTypes: bool =False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the mass of particles within a specified aperture and belonging to
        subgroups.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the total particle subgroups mass.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall substructure mass, ot return the of each particle type separately.
        This toggle is controlled by a boolean.

        substructure_mass = total mass - fuzz mass

        N.B.: The Fuzz mass is given by subgroupnumber = 0.
        N.B.: The aperture radius could cut subgroups at the outer edge.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            total_mass = self.group_mass_aperture(out_allPartTypes=True, aperture_radius=aperture_radius)
            substructure_mass_PartTypes = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_subgroupnumber')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                subgroupnumber = getattr(self, f'partType{part_type}_subgroupnumber')
                aperture_radius_index = np.where((radial_dist < aperture_radius) & (subgroupnumber == 0))[0]
                free_memory(['radial_dist', 'subgroupnumber'])
                _mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")

                fuzz_mass = np.sum(_mass)
                substructure_mass = total_mass[['0', '1', '4'].index(part_type)] - fuzz_mass
                substructure_mass_PartTypes = np.append(substructure_mass_PartTypes, substructure_mass)

            return substructure_mass_PartTypes

        else:

            total_mass = self.group_mass_aperture(out_allPartTypes=False, aperture_radius=aperture_radius)
            fuzz_mass = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_subgroupnumber')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                subgroupnumber = getattr(self, f'partType{part_type}_subgroupnumber')
                aperture_radius_index = np.where((radial_dist < aperture_radius) & (subgroupnumber == 0))[0]
                free_memory(['radial_dist', 'subgroupnumber'])
                _mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                fuzz_mass = np.append(fuzz_mass, _mass)

            substructure_mass = total_mass - np.sum(fuzz_mass)
            return substructure_mass

    def group_substructure_fraction(self,
                             out_allPartTypes: bool =False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the fraction of mass constituting substructures, taking particles
        within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the total particle substructure fraction.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall substructure mass fraction, ot return the of each particle type separately.
        This toggle is controlled by a boolean.

        substructure_fraction = (total mass - fuzz mass) / total mass

        N.B.: The Fuzz mass is given by subgroupnumber = 0.
        N.B.: The aperture radius could cut subgroups at the outer edge.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            total_mass = self.group_mass_aperture(out_allPartTypes=True, aperture_radius=aperture_radius)
            substructure_frac_PartTypes = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_subgroupnumber')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                subgroupnumber = getattr(self, f'partType{part_type}_subgroupnumber')
                aperture_radius_index = np.where((radial_dist < aperture_radius) & (subgroupnumber == 0))[0]
                free_memory(['radial_dist', 'subgroupnumber'])
                _mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")

                fuzz_mass = np.sum(_mass)
                substructure_mass = 1 - (fuzz_mass/total_mass[['0', '1', '4', '5'].index(part_type)])
                substructure_frac_PartTypes = np.append(substructure_frac_PartTypes, substructure_mass)

            return substructure_frac_PartTypes

        else:

            total_mass = self.group_mass_aperture(out_allPartTypes=False, aperture_radius=aperture_radius)
            fuzz_mass = np.zeros(0, dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_subgroupnumber')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                subgroupnumber = getattr(self, f'partType{part_type}_subgroupnumber')
                aperture_radius_index = np.where((radial_dist < aperture_radius) & (subgroupnumber == 0))[0]
                free_memory(['radial_dist', 'subgroupnumber'])
                _mass = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                fuzz_mass = np.append(fuzz_mass, _mass)

            substructure_fraction = 1 - (np.sum(fuzz_mass)/total_mass)
            return substructure_fraction

    def group_centre_of_mass(self,
                             out_allPartTypes: bool =False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the centre of mass of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the centre of mass using the relative
        static method in this Mixin class.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall centre of mass, ot return the centre of mass of each particle type separately.
        This toggle is controlled by a boolean.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            CoM_PartTypes = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass   = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _coords = getattr(self, f'partType{part_type}_coordinates')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _coords.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")

                centre_of_mass = self.centre_of_mass(_mass, _coords)
                CoM_PartTypes = np.concatenate((CoM_PartTypes, [centre_of_mass]), axis=0)

            return CoM_PartTypes

        else:

            mass   = np.zeros(0, dtype=np.float)
            coords = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass   = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _coords = getattr(self, f'partType{part_type}_coordinates')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _coords.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                mass   = np.concatenate((mass, _mass), axis=0)
                coords = np.concatenate((coords, _coords), axis=0)

            return self.centre_of_mass(mass, coords)

    def group_dynamical_merging_index(self,
                             out_allPartTypes: bool = False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the group_dynamical_merging_index of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the merging index using the relative
        static method in this Mixin class.

        group_dynamical_merging_index = || CoM(aperture) - CoP || / aperture

        It also has the option of combining all the particles of different types and computing
        the overall centre of mass, ot return the centre of mass of each particle type separately.
        This toggle is controlled by a boolean.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        centre_of_mass = self.group_centre_of_mass(out_allPartTypes=out_allPartTypes,
                                                   aperture_radius=aperture_radius)
        if out_allPartTypes:
            dynamical_merging_index = self.radial_distance_CoP(centre_of_mass)/aperture_radius
        else:
            dynamical_merging_index = self.radial_distance_CoP(centre_of_mass[:, None].T)/aperture_radius
            dynamical_merging_index = dynamical_merging_index[0]

        return dynamical_merging_index

    def group_zero_momentum_frame(self,
                             out_allPartTypes: bool = False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the cluster's rest frame of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the rest frame using the relative
        static method in this Mixin class.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall bulk velocity, ot return the bulk velocity of each particle type separately.
        This toggle is controlled by a boolean.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            ZMF_PartTypes = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_velocity')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass     = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _velocity = getattr(self, f'partType{part_type}_velocity')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _velocity.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")

                zmf = self.zero_momentum_frame(_mass, _velocity)
                ZMF_PartTypes = np.concatenate((ZMF_PartTypes, [zmf]), axis=0)

            return ZMF_PartTypes

        else:

            mass     = np.zeros(0, dtype=np.float)
            velocity = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_velocity')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass     = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _velocity = getattr(self, f'partType{part_type}_velocity')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _velocity.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                mass     = np.concatenate((mass, _mass), axis=0)
                velocity = np.concatenate((velocity, _velocity), axis=0)

            return self.zero_momentum_frame(mass, velocity)

    def group_angular_momentum(self,
                             out_allPartTypes: bool = False,
                             aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the cluster's angular momentum from particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the angular momentum using the relative
        static method in this Mixin class.
        The method also checks that the necessary datasets are loaded into the cluster object.
        It also has the option of combining all the particles of different types and computing
        the overall angular momentum, or return the angular momentum of each particle type separately.
        This toggle is controlled by a boolean.

        :param out_allPartTypes: default = False
        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        if out_allPartTypes:

            ANG_PartTypes = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_velocity')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass     = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _coords   = getattr(self, f'partType{part_type}_coordinates')[aperture_radius_index]
                _velocity = getattr(self, f'partType{part_type}_velocity')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _coords.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _velocity.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")

                # Rescale coordinates and velocity
                _coords   = np.subtract(_coords, self.centre_of_potential)
                _velocity = np.subtract(_velocity, self.group_zero_momentum_frame(aperture_radius=aperture_radius))
                ang = self.angular_momentum(_mass, _coords, _velocity)
                ANG_PartTypes = np.concatenate((ANG_PartTypes, [ang]), axis=0)

            return ANG_PartTypes

        else:

            mass     = np.zeros(0, dtype=np.float)
            coords   = np.zeros((0, 3), dtype=np.float)
            velocity = np.zeros((0, 3), dtype=np.float)

            for part_type in ['0', '1', '4']:
                assert hasattr(self, f'partType{part_type}_coordinates')
                assert hasattr(self, f'partType{part_type}_velocity')
                assert hasattr(self, f'partType{part_type}_mass')
                radial_dist = self.radial_distance_CoP(getattr(self, f'partType{part_type}_coordinates'))
                aperture_radius_index = np.where(radial_dist < aperture_radius)[0]
                free_memory(['radial_dist'])
                _mass     = getattr(self, f'partType{part_type}_mass')[aperture_radius_index]
                _coords   = getattr(self, f'partType{part_type}_coordinates')[aperture_radius_index]
                _velocity = getattr(self, f'partType{part_type}_velocity')[aperture_radius_index]
                if _mass.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _coords.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                if _velocity.__len__() == 0: warnings.warn(f"Array PartType{part_type} is empty - check filtering.")
                mass     = np.concatenate((mass, _mass), axis=0)
                coords   = np.concatenate((coords, _coords), axis=0)
                velocity = np.concatenate((velocity, _velocity), axis=0)

            # Rescale coordinates and velocity
            coords   = np.subtract(coords, self.centre_of_potential)
            velocity = np.subtract(velocity, self.group_zero_momentum_frame(aperture_radius=aperture_radius))
            return self.angular_momentum(mass, coords, velocity)

    def group_thermodynamic_merging_index(self,
                                      aperture_radius: float = None) -> np.ndarray:
        """
        Method that computes the group_thermodynamic_merging_index of particles within a specified aperture.
        If the aperture is not specified, it is set by default to the true R500 of the cluster
        radius from the centre of potential and computes the merging index using the relative
        static method in this Mixin class.

        group_thermodynamic_merging_index = kinetic_energy / thermal_energy

        :param aperture_radius: default = None (R500)
        :return: expected a numpy array of dimension 1 if all particletypes are combined, or
            dimension 2 if particle types are returned separately.
        """
        if aperture_radius is None:
            aperture_radius = self.r500
            warnings.warn(f'Aperture radius set to default R_500,true. = {self.r500:.2f} Mpc.')

        kinetic_energy = self.group_kinetic_energy(out_allPartTypes=True, aperture_radius=aperture_radius)[0]
        thermal_energy = self.group_thermal_energy(aperture_radius=aperture_radius)
        return kinetic_energy / thermal_energy

    @staticmethod
    def rotation_matrix_from_vectors(vec1, vec2):
        """
        Find the rotation matrix that aligns vec1 to vec2

        :param vec1: A 3d "source" vector
        :param vec2: A 3d "destination" vector
        :return mat: A transform matrix (3x3) which when applied to vec1, aligns it with vec2.
        """
        a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 / np.linalg.norm(vec2)).reshape(3)
        v = np.cross(a, b)
        c = np.dot(a, b)
        s = np.linalg.norm(v)
        kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
        rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s ** 2))

        # Catch exception where the vectors are perfectly aligned
        if Mixin.angle_between_vectors(vec1, vec2) < 1e-10:
            return np.identity(3)
        elif np.abs(Mixin.angle_between_vectors(vec1, vec2) - 180) < 1e-10:
            return (-1) * np.identity(3)
        else:
            return np.asarray(rotation_matrix)

    @staticmethod
    def rotation_matrix_about_axis(axis, theta):
        """
        Find the rotation matrix that rotates a vector about the origin by delta(theta) in the elevation
        and by delta(phi) about the azimuthal axis.

        :param theta: expect float in degrees
            The displacement angle in elevation.

        :param phi: expect float in degrees
            The displacement angle in azimuth.

        :return: A transform matrix (3x3)
        """
        axis = np.asarray(axis)
        axis = axis / np.sqrt(np.dot(axis, axis))
        a = np.cos(theta / 2.0)
        b, c, d = -axis * np.sin(theta / 2.0)
        aa, bb, cc, dd = a * a, b * b, c * c, d * d
        bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
        return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                         [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                         [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

    def rotation_matrix_about_x(self, angle):
        axis = [1,0,0]
        return self.rotation_matrix_about_axis(axis, angle)

    def rotation_matrix_about_y(self, angle):
        axis = [0,1,0]
        return self.rotation_matrix_about_axis(axis, angle)

    def rotation_matrix_about_z(self, angle):
        axis = [0,0,1]
        return self.rotation_matrix_about_axis(axis, angle)

    @staticmethod
    def _TEST_apply_rotation_matrix(rot_matrix, vectors):
        """Apply this rotation to a set of vectors.
        This function is a self-made test method for comparison with the `scipy`
        method implemented below.

        :param rot_matrix: array_like, shape (3,3) or (N, 3, 3)
            Each matrix represents a rotation in 3D space of the corresponding
            vector.

        :param vectors: array_like, shape (3,) or (N, 3)
            Each `vectors[i]` represents a vector in 3D space. A single vector
            can either be specified with shape `(3, )` or `(1, 3)`. The number
            of rotations and number of vectors given must follow standard numpy
            broadcasting rules: either one of them equals unity or they both
            equal each other.

        :return rotated_vectors : ndarray, shape (3,) or (N, 3)
            Result of applying rotation on input vectors.
            Shape depends on the following cases:
                - If object contains a single rotation (as opposed to a stack
                  with a single rotation) and a single vector is specified with
                  shape ``(3,)``, then `rotated_vectors` has shape ``(3,)``.
                - In all other cases, `rotated_vectors` has shape ``(N, 3)``,
                  where ``N`` is either the number of rotations or vectors.
        """

        vectors = np.asarray(vectors)
        rot_matrix = np.asarray(rot_matrix)

        assert rot_matrix.shape == (3,3)
        assert vectors.__len__() > 0
        return rot_matrix.dot(vectors)

    @staticmethod
    def apply_rotation_matrix(rot_matrix, vectors, inverse=False):
        """Apply this rotation to a set of vectors.
        If the original frame rotates to the final frame by this rotation, then
        its application to a vector can be seen in two ways:
            - As a projection of vector components expressed in the final frame
              to the original frame.
            - As the physical rotation of a vector being glued to the original
              frame as it rotates. In this case the vector components are
              expressed in the original frame before and after the rotation.
        In terms of rotation matricies, this application is the same as
        ``self.as_matrix().dot(vectors)``.

        Parameters
        ----------
        vectors : array_like, shape (3,) or (N, 3)
            Each `vectors[i]` represents a vector in 3D space. A single vector
            can either be specified with shape `(3, )` or `(1, 3)`. The number
            of rotations and number of vectors given must follow standard numpy
            broadcasting rules: either one of them equals unity or they both
            equal each other.
        inverse : boolean, optional
            If True then the inverse of the rotation(s) is applied to the input
            vectors. Default is False.

        Returns
        -------
        rotated_vectors : ndarray, shape (3,) or (N, 3)
            Result of applying rotation on input vectors.
            Shape depends on the following cases:
                - If object contains a single rotation (as opposed to a stack
                  with a single rotation) and a single vector is specified with
                  shape ``(3,)``, then `rotated_vectors` has shape ``(3,)``.
                - In all other cases, `rotated_vectors` has shape ``(N, 3)``,
                  where ``N`` is either the number of rotations or vectors.

        """

        vectors = np.asarray(vectors)
        rot_matrix = np.asarray(rot_matrix)

        if vectors.ndim > 2 or vectors.shape[-1] != 3:
            raise ValueError("Expected input of shape (3,) or (P, 3), "
                             "got {}.".format(vectors.shape))

        single_vector = False
        if vectors.shape == (3,):
            single_vector = True
            vectors = vectors[None, :]

        single_rot = False
        if rot_matrix.shape == (3,3):
            rot_matrix = rot_matrix.reshape((1, 3, 3))
            single_rot = True

        n_vectors = vectors.shape[0]
        n_rotations = len(rot_matrix)

        if n_vectors != 1 and n_rotations != 1 and n_vectors != n_rotations:
            raise ValueError("Expected equal numbers of rotations and vectors "
                             ", or a single rotation, or a single vector, got "
                             "{} rotations and {} vectors.".format(
                             n_rotations, n_vectors))

        if inverse:
            result = np.einsum('ikj,ik->ij', rot_matrix, vectors)
        else:
            result = np.einsum('ijk,ik->ij', rot_matrix, vectors)

        if single_rot and single_vector:
            return result[0]
        else:
            return result

    #####################################################
    #													#
    #				COMOVING UNITS      				#
    # 									 				#
    #####################################################

    def comoving_density(self, density):
        """
        Rescales the density from the comoving coordinates to the physical coordinates
        """
        hubble_par = self.file_hubble_param()
        redshift = self.file_redshift()
        scale_factor = 1 / (redshift + 1)
        return np.multiply(density,  hubble_par ** 2 * scale_factor ** -3)

    def comoving_length(self, coord):
        """
        Rescales the density from the comoving length to the physical length
        """
        hubble_par = self.file_hubble_param()
        redshift = self.file_redshift()
        scale_factor = 1 / (redshift + 1)
        return np.multiply(coord, scale_factor / hubble_par)

    def comoving_velocity(self, vel):
        """
        Rescales the density from the comoving velocity to the physical velocity
        """
        redshift = self.file_redshift()
        scale_factor = 1 / (redshift + 1)
        return np.multiply(vel, np.sqrt(scale_factor))

    def comoving_mass(self, mass):
        """
        Rescales the density from the comoving mass to the physical mass
        """
        hubble_par = self.file_hubble_param()
        return np.divide(mass, hubble_par)

    def comoving_kinetic_energy(self, kinetic_energy):
        """
        Rescales the density from the comoving kinetic_energy to the physical kinetic_energy
        """
        hubble_par = self.file_hubble_param()
        redshift = self.file_redshift()
        scale_factor = 1 / (redshift + 1)
        _kinetic_energy = np.multiply(kinetic_energy, scale_factor)
        _kinetic_energy = np.divide(_kinetic_energy, hubble_par)
        return _kinetic_energy

    def comoving_momentum(self, mom):
        """
        Rescales the momentum from the comoving to the physical
        """
        hubble_par = self.file_hubble_param()
        redshift = self.file_redshift()
        scale_factor = 1 / (redshift + 1)
        return np.multiply(mom, np.sqrt(scale_factor) / hubble_par)

    def comoving_ang_momentum(self, angmom):
        """
        Rescales the angular momentum from the comoving to the physical
        """
        hubble_par = self.file_hubble_param()
        redshift = self.file_redshift()
        scale_factor = 1 / (redshift + 1)
        return np.multiply(angmom, np.sqrt(scale_factor**3) / np.power(hubble_par, 2.))

    #####################################################
    #													#
    #			    UNITS CONEVRSION     				#
    # 									 				#
    #####################################################

    @staticmethod
    def density_units(density, unit_system='SI'):
        """
        CREATED: 12.02.2019
        LAST MODIFIED: 12.02.2019

        INPUTS: density np.array

                metric system used: 'SI' or 'cgs' or astronomical 'astro'
        """
        if unit_system == 'SI':
            # kg*m^-3
            conv_factor = 6.769911178294543 * 10 ** -28
        elif unit_system == 'cgs':
            # g*cm^-3
            conv_factor = 6.769911178294543 * 10 ** -31
        elif unit_system == 'astro':
            # solar masses / (parsec)^3
            conv_factor = 6.769911178294543 * np.power(3.086, 3.) / 1.9891 * 10 ** -10
        elif unit_system == 'nHcgs':
            conv_factor = 6.769911178294543 * 10 ** -31 / (1.674 * 10 ** -24.)
        else:
            raise("[ERROR] Trying to convert SPH density to an unknown metric system.")

        return np.multiply(density, conv_factor)

    @staticmethod
    def velocity_units(velocity, unit_system='SI'):
        """
        CREATED: 14.02.2019
        LAST MODIFIED: 14.02.2019

        INPUTS: velocity np.array

                metric system used: 'SI' or 'cgs' or astronomical 'astro'
        """
        if unit_system == 'SI':
            # m/s
            conv_factor =  10 ** 3
        elif unit_system == 'cgs':
            # cm/s
            conv_factor =  10 ** 5
        elif unit_system == 'astro':
            # km/s
            conv_factor =  1
        else:
            raise("[ERROR] Trying to convert velocity to an unknown metric system.")

        return np.multiply(velocity, conv_factor)

    @staticmethod
    def length_units(len, unit_system='SI'):
        """
		CREATED: 14.02.2019
		LAST MODIFIED: 14.02.2019

		INPUTS: len np.array

				metric system used: 'SI' or 'cgs' or astronomical 'astro'
		"""
        if unit_system == 'SI':
            # m/s
            conv_factor = 1e6 * parsec
        elif unit_system == 'cgs':
            # cm/s
            conv_factor = 1e8 * parsec
        elif unit_system == 'astro':
            # km/s
            conv_factor = 1
        else:
            raise ("[ERROR] Trying to convert length to an unknown metric system.")

        return np.multiply(len, conv_factor)

    @staticmethod
    def mass_units(mass, unit_system='SI'):
        """
        CREATED: 14.02.2019
        LAST MODIFIED: 14.02.2019

        INPUTS: mass np.array

                metric system used: 'SI' or 'cgs' or astronomical 'astro'
        """
        if unit_system == 'SI':
            # m/s
            conv_factor = 1e10 * solar_mass
        elif unit_system == 'cgs':
            # cm/s
            conv_factor = 1e13 * solar_mass
        elif unit_system == 'astro':
            # km/s
            conv_factor = 1e10
        else:
            raise("[ERROR] Trying to convert mass to an unknown metric system.")

        return np.multiply(mass, conv_factor)

    @staticmethod
    def momentum_units(momentum, unit_system='SI'):
        """
        CREATED: 07.03.2019
        LAST MODIFIED: 07.03.2019

        INPUTS: momentum np.array

                metric system used: 'SI' or 'cgs' or astronomical 'astro'
        """
        if unit_system == 'SI':
            # m/s
            conv_factor = 1.9891 * 10 ** 43
        elif unit_system == 'cgs':
            # cm/s
            conv_factor = 1.9891 * 10 ** 48
        elif unit_system == 'astro':
            # km/s
            conv_factor = 10 ** 10
        else:
            raise("[ERROR] Trying to convert mass to an unknown metric system.")

        return np.multiply(momentum, conv_factor)

    @staticmethod
    def energy_units(energy, unit_system='SI'):
        """
		CREATED: 07.03.2019
		LAST MODIFIED: 07.03.2019

		INPUTS: momentum np.array

				metric system used: 'SI' or 'cgs' or astronomical 'astro'
		"""
        if unit_system == 'SI':
            # m/s
            conv_factor = np.power(10., 46)

        else:
            raise ("[ERROR] Trying to convert mass to an unknown metric system.")

        return np.multiply(energy, conv_factor)

max_stars_repo_name: LBJ-Wade/C-Eagle-analysis

max_issues_repo_name: LBJ-Wade/C-Eagle-analysis

max_forks_repo_name: LBJ-Wade/C-Eagle-analysis


text:
import numpy as np
from .taxi import TaxiGrid5x5, TaxiGrid10x10
from ..gridworld.skills import GoToGridPosition

skills5x5 = ['interact'] + list(TaxiGrid5x5.depot_names)
skills10x10 = ['interact'] + list(TaxiGrid10x10.depot_names)

def run_skill(w, name):
    if name in w.depot_locs.keys():
        skill = lambda w, x: GoToDepot(w, x, name)
    else:
        skill = lambda w, x: Interact(w)
    while True:
        can_run, a, term = skill(w, w.agent.position)
        print(can_run, a, term)
        assert (can_run or term)
        if can_run:
            w.step(a)
        if term:
            break

def skill_policy(w, skill_name):
    if skill_name in w.depot_locs.keys():
        skill = lambda w, x: GoToDepot(w, x, skill_name)
    elif skill_name == 'interact':
        skill = lambda w, x: Interact(w)
    else:
        raise ValueError('Invalid skill name' + str(skill_name))
    return skill(w, w.agent.position)

def GoToDepot(gridworld, start, depotname):
    depot = gridworld.depot_locs[depotname]
    return GoToGridPosition(gridworld, start=start, target=depot)[0]

def Interact(gridworld):
    # Check relevant state variables to see if skill can run
    agent_pos = gridworld.agent.position
    at_depot = any(np.all(loc == agent_pos) for _, loc in gridworld.depot_locs.items())
    at_passenger = any(np.all(p.position == agent_pos) for p in gridworld.passengers)
    crowded = (gridworld.passenger is not None and any(
        np.all(p.position == agent_pos) for p in gridworld.passengers if p != gridworld.passenger))

    if at_depot and at_passenger and not crowded:
        can_run = True
        action = 4
    else:
        # Nothing to do
        can_run = False
        action = None

    terminate = True  # skill always terminates in one step
    return (can_run, action, terminate)

max_stars_repo_name: abagaria/visgrid

max_issues_repo_name: abagaria/visgrid

max_forks_repo_name: abagaria/visgrid


text:
from collections import Counter
from itertools import combinations

import clana.io
import clana.visualize_cm
import networkx as nx
import numpy as np
import pandas as pd
import progressbar

# Load the data
df = pd.read_csv("articles.csv")
df["author"] = df["author"].str.split("::")

# Analyze the data
df = df[~df["author"].isna()]
authors = [author for authorset in df["author"].tolist() for author in authorset]
author_count = Counter(authors)

print("* Publications: {}".format(len(df)))
print("* Unique elements: {}".format(len(author_count)))
print("* Most common:")
most_common = sorted(author_count.items(), key=lambda n: n[1], reverse=True)
for name, count in most_common[:10]:
    print(f"    {count:>4}x {name}")

unique_authors = sorted(list(author_count.keys()))

def get_biggest_clusters(edges, n=10):
    G = nx.Graph()
    for authorset in edges.tolist():
        for author in authorset:
            G.add_node(author)

    for authorset in progressbar.progressbar(df["author"].tolist()[:10_000]):
        for author1, author2 in combinations(authorset, 2):
            G.add_edge(author1, author2)

    print("Edges were added")

    components = [c for c in sorted(nx.connected_components(G), key=len, reverse=True)]
    return components[:n]

def create_matrix(nodes, edges):
    n2i = {node: i for i, node in enumerate(sorted(nodes))}
    # node to index
    mat = np.zeros((len(nodes), len(nodes)), dtype=np.int32)
    for edge in edges:
        for a, b in combinations(edge, 2):
            if a not in n2i or b not in n2i:
                continue
            mat[n2i[a]][n2i[b]] += 1
            if a != b:
                mat[n2i[b]][n2i[a]] += 1
    return mat, sorted(nodes)

components = get_biggest_clusters(df["author"])
print("* Biggest clusters: {}".format([len(el) for el in components]))

component_w_publications = [(author, author_count[author]) for author in components[0]]
component_w_publications = sorted(
    component_w_publications, key=lambda n: n[1], reverse=True
)
authors = [author for author, count in component_w_publications[:1_00]]
mat, labels = create_matrix(authors, df["author"].tolist())

clana.visualize_cm.main(
    "coauthors.json",
    perm_file="",
    steps=1_000_000,
    labels_file="labels.json",
    zero_diagonal=False,
    output="cm-ordered.pdf",
)
clana.io.write_cm("coauthors.json", mat)
clana.io.write_labels("labels.json", labels)

max_issues_repo_name: Kerwin-Xie/algorithms

max_forks_repo_name: Kerwin-Xie/algorithms

reason: import numpy,import networkx


text:
# coding: utf-8

# In[1]:

import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
#get_ipython().magic('matplotlib inline')
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 12,10
import random
random.seed('liuyi')

# In[83]:

# In[84]:

#ls

# In[58]:

gene_Data=pd.read_csv("./gene_data.txt")
gene_Data=gene_Data.set_index("Unnamed: 0")
gene_Data=gene_Data.T

# In[59]:

gene_Data

# In[60]:

Symbol=pd.read_csv("./Symbol.txt")
Symbol=Symbol.set_index('Unnamed: 0').T.to_dict('list')

# In[61]:

Symbol.get('1053_at')

# In[62]:

x=gene_Data.values[[0,1,4,5],:]

# In[63]:

x.shape

# # Design Matrix is X

# In[64]:

y=np.array([1,1,0,0])

# In[65]:

y.shape

# In[66]:

from sklearn.svm import LinearSVC
from sklearn.feature_selection import SelectFromModel
from sklearn.linear_model import LogisticRegression

# In[67]:

clf = LogisticRegression(penalty='l1'
                       ,max_iter=2000,C=100).fit(x, y)

# In[68]:

coefs=clf.coef_.reshape(clf.coef_.shape[1])

# In[69]:

nonzero_index=np.nonzero(coefs)

# In[70]:

len(nonzero_index[0])

# In[71]:

np.nonzero(coefs)[0]

# In[38]:

coefs[np.nonzero(coefs)[0]]

# In[72]:

relevant_genes=gene_Data.columns[nonzero_index]
len(relevant_genes)

# In[73]:

len(relevant_genes)

# In[74]:

relevant_genes

# In[75]:

from bioservices.kegg import KEGG
k = KEGG()

# In[76]:

A=list(relevant_genes)

# In[77]:

import sys

# In[78]:

len(A)

# In[79]:

#pwd

# In[81]:

with open('./pathway_baicalin_386_genes.txt','w') as f:
    for i in A:
        try:
            listname=Symbol.get(i)
            n=str(listname[0])
            if (n!='nan'):
                j=k.get_pathway_by_gene(n,"hsa")
                if (j!=None):
                    for keys,values in j.items():
                        print(n,keys,values,sep="\t")
                    #f.write("%s \n" %n)
                    #for keys,values in j.items():
                    #    f.write("%s %s %s\n" %(n,keys,values))
        except AttributeError as err:
            print('error')

# In[134]:

#Print algorithms
#for i in A:
#    try:
#        listname=Symbol.get(i)
#        n=str(listname[0])
#        if (n!='nan'):
#            #print("hello")
#            j=k.get_pathway_by_gene(n,"hsa")
#            if (j!=None):
#                for keys,values in j.items():
#                     print(n,keys,values,sep="\t")
#    except AttributeError as err:
#        print('error')

# # Plot

# In[82]:

plt.rcdefaults()
fig, ax = plt.subplots()

# set x,y axis
genes =A
y_pos = np.arange(len(genes))
performance = coefs[np.nonzero(coefs)[0]]
#error = np.random.rand(len(genes))

ax.barh(y_pos, performance, align='center',
        color='purple', ecolor='black')

ax.set_yticks(y_pos)
ax.set_yticklabels(genes)
ax.invert_yaxis()  # labels read top-to-bottom

ax.set_xlabel('Significant Coefficients')
ax.set_ylabel('Genes')
ax.set_title('Relavant genes')
#plot figues or save figures
#plt.show()
#plt.savefig("gene.png")