<commit_before># coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.

from __future__ import division, print_function, unicode_literals
from __future__ import absolute_import

"""
This module provides a base class for tensor-like objects and methods for
basic tensor manipulation.  It also provides a class, SquareTensor,
that provides basic methods for creating and manipulating rank 2 tensors
"""


__author__ = "<NAME>"
__copyright__ = "Copyright 2012, The Materials Project"
__credits__ = ("<NAME>, <NAME>, <NAME>, "
               "<NAME>, <NAME>, <NAME>")
__version__ = "1.0"
__maintainer__ = "<NAME>"
__email__ = "<EMAIL>"
__status__ = "Development"
__date__ = "March 22, 2012"


from scipy.linalg import polar
import numpy as np
import itertools
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.core.operations import SymmOp
from numpy.linalg import norm

class TensorBase(np.ndarray):
    """
    Base class for doing useful general operations on Nth order tensors,
    without restrictions on the type (stress, elastic, strain, piezo, etc.)
    """

    def __new__(cls, input_array):
        """
        Create a TensorBase object.  Note that the constructor uses __new__
        rather than __init__ according to the standard method of
        subclassing numpy ndarrays.

        Args:
            input_array: (3xN array-like): the 3xN array-like representing
                a tensor quantity
        """
        obj = np.asarray(input_array).view(cls)
        obj.rank = len(obj.shape)
        if not all([i == 3 for i in obj.shape]):
            raise ValueError("Pymatgen only supports 3-dimensional tensors")
        
        return obj

    def __array_finalize__(self, obj):
        if obj is None:
            return
        self.rank = getattr(obj, 'rank', None)

    def __array_wrap__(self, obj):
        """
        Overrides __array_wrap__ methods in ndarray superclass to avoid errors
        associated with functions that return scalar values
        """

        if len(obj.shape) == 0:
            return obj[()]
        else:
            return np.ndarray.__array_wrap__(self, obj)
        
    def __hash__(self):
        """
        define a hash function, since numpy arrays
        have their own __eq__ method
        """
        return hash(self.tostring())

    def __repr__(cls):
        return "{}({})".format(cls.__class__.__name__,
                               cls.__str__())

    def zeroed(self, tol = 1e-3):
        """
        returns the matrix with all entries below a certain threshold
        (i.e. tol) set to zero
        """
        new_tensor = self.copy()
        new_tensor[abs(new_tensor) < tol] = 0
        return new_tensor

    def transform(self, symm_op):
        """
        Applies a transformation (via a symmetry operation) to a tensor. 

        Args:
            symm_op (3x3 array-like): a symmetry operation to apply to the tensor
        """
        return self.__class__(symm_op.transform_tensor(self))

    def rotate(self, matrix, tol=1e-5):
        matrix = SquareTensor(matrix)
        if not matrix.is_rotation(tol):
            raise ValueError("Rotation matrix is not valid.")
        sop = SymmOp.from_rotation_and_translation(matrix,
                                                   [0., 0., 0.])
        return self.transform(sop)


    @property
    def symmetrized(self):
        """
        Returns a generally symmetrized tensor, calculated by taking 
        the sum of the tensor and its transpose with respect to all 
        possible permutations of indices
        """
        perms = list(itertools.permutations(range(self.rank)))
        return sum([np.transpose(self, ind) for ind in perms]) / len(perms)

    def is_symmetric(self, tol=1e-5):
        """
        Tests whether a tensor is symmetric or not based on the residual
        with its symmetric part, from self.symmetrized

        Args:
            tol (float): tolerance to test for symmetry
        """
        return (self - self.symmetrized < tol).all()

    def fit_to_structure(self, structure, symprec = 0.1):
        """
        Returns a tensor that is invariant with respect to symmetry
        operations corresponding to a structure

        Args: 
            structure (Structure): structure from which to generate 
                symmetry operations
            symprec (float): symmetry tolerance for the Spacegroup Analyzer
                used to generate the symmetry operations
        """
        sga = SpacegroupAnalyzer(structure, symprec)
        symm_ops = sga.get_symmetry_operations(cartesian=True)
        return sum([self.transform(symm_op)
                    for symm_op in symm_ops]) / len(symm_ops)

    def is_fit_to_structure(self, structure, tol=1e-2):
        """
        Tests whether a tensor is invariant with respect to the
        symmetry operations of a particular structure by testing
        whether the residual of the symmetric portion is below a 
        tolerance
        
        Args:
            tol (float): tolerance for symmetry testing
        """
        return (self - self.fit_to_structure(structure) < tol).all()

    def convert_to_ieee(self, structure, atol = 0.1, ltol = 1e-3):
        """
        Given a structure associated with a tensor, attempts a
        calculation of the tensor in IEEE format according to
        the 1978 IEEE standards.

        Args:
            structure (structure): a structure associated with the
            symprec (float): precision for the symmetry analyzer
        """
        def get_uvec(vec):
            """ Gets a unit vector parallel to input vector"""
            return vec / np.linalg.norm(vec)

        vecs = structure.lattice.matrix
        lengths = np.array(structure.lattice.abc)
        angles = np.array(structure.lattice.angles)

        a = b = c = None
        sga = SpacegroupAnalyzer(structure)
        xtal_sys = sga.get_crystal_system()
        angle_error = "Angle check failure. {} crystal may not be in "\
                      "conventional representation".format(xtal_sys)
        rotation = np.zeros((3,3))

        # IEEE rules: a,b,c || x1,x2,x3
        if xtal_sys == "cubic":
            if (abs(angles - 90.0) > atol).any():
                raise ValueError(angle_error)
            if (abs(lengths - lengths[0]) > ltol).any():
                raise ValueError("Cubic vectors not equal in length")
            rotation = [vecs[i]/lengths[i] for i in range(3)]
        
        # IEEE rules: a=b in length; c,a || x3, x1
        elif xtal_sys == "tetragonal":
            if (abs(angles - 90.0) > atol).any():
                raise ValueError(angle_error)
            # This may be redundant, should be cubic
            if (abs(lengths - lengths[0]) < ltol).all():
                raise ValueError("Tetragonal c vector indistinguishable")

            rotation = [vecs[i]/lengths[i] for i in range(3)]
            if abs(lengths[2] - lengths[0]) < ltol:
                rotation[1], rotation[2] = rotation[2], rotation[1].copy()
            elif abs(lengths[2] - lengths[1]) < ltol:
                rotation[0], rotation[2] = rotation[2], rotation[0].copy()
            else:
                raise ValueError("No latt. vecs equal for tetragonal system")

        # IEEE rules: c<a<b; c,a || x3,x1
        elif xtal_sys == "orthorhombic":
            if (abs(angles - 90.0) > atol).any():
                raise ValueError(angle_error)
            # May be redundant, should be tetragonal
            if len(np.unique(np.floor(lengths/ltol))) < 3:
                raise ValueError("Orthorhombic vectors indistinguishable")
            # TODO: Check this
            # import pdb; pdb.set_trace()
            rotation = [vec/mag for (mag, vec) in sorted(zip(lengths, vecs))]
            rotation = np.roll(rotation, 2, axis = 0)

        # IEEE rules: c,a || x3,x1, c is threefold axis
        elif xtal_sys in ("trigonal", "hexagonal"):
            # Rhombohedral lattice
            if xtal_sys == "trigonal" and (angles - angles[0] < atol).all():
                rotation[0] = get_uvec(vecs[2] - vecs[0])
                rotation[2] = get_uvec(np.sum(vecs, axis=0))
                rotation[1] = get_uvec(np.cross(rotation[2], rotation[0]))
            # Standard hexagonal lattice
            else:
                # find threefold axis:
                tf_mask = np.where(abs(angles-120.0) < atol)
                non_tf_mask = np.logical_not(tf_mask)
                rotation[2] = get_uvec(vecs[tf_mask])
                rotation[0] = get_uvec(vecs[non_tf_mask][0])
                rotation[1] = get_uvec(np.cross(rotation[2], rotation[0]))

        # IEEE rules: b,c || x2,x3; alpha=beta=90, c<a
        elif xtal_sys == "monoclinic":
            # Find unique axis
            umask = np.where(angles - 90.0 > atol)
            n_umask = np.logical_not(umask)
            if sum(umask) != 1:
                raise ValueError(angle_error)
            rotation[1] = get_uvec(vecs[umask])
            # Shorter of remaining lattice vectors for c axis
            rotation[2] = [vec/mag for (mag, vec) in 
                           sorted(zip(lengths[n_umask], vecs[n_umask]))][0]
            rotation[0] = np.cross(rotation[1], rotation[2])

        # IEEE rules: c || x3
        elif xtal_sys == "triclinic":

        else:
            print("{} crystal not implemented".format(xtal_sys))
            return np.zeros(self.shape)
        # import pdb; pdb.set_trace()
        return self.rotate(rotation)


class SquareTensor(TensorBase):
    """
    Base class for doing useful general operations on second rank tensors
    (stress, strain etc.).
    """

    def __new__(cls, input_array):
        """
        Create a SquareTensor object.  Note that the constructor uses __new__
        rather than __init__ according to the standard method of
        subclassing numpy ndarrays.  Error is thrown when the class is
        initialized with non-square matrix.

        Args:
            stress_matrix (3x3 array-like): the 3x3 array-like
                representing the Green-Lagrange strain
        """

        obj = TensorBase(input_array).view(cls)
        if not (len(obj.shape) == 2):
            raise ValueError("SquareTensor only takes 2-D "
                             "tensors as input")
        return obj
        
    @property<commit_msg>added triclinic, last commit before removing checks<commit_after># coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.

from __future__ import division, print_function, unicode_literals
from __future__ import absolute_import

"""
This module provides a base class for tensor-like objects and methods for
basic tensor manipulation.  It also provides a class, SquareTensor,
that provides basic methods for creating and manipulating rank 2 tensors
"""


__author__ = "<NAME>"
__copyright__ = "Copyright 2012, The Materials Project"
__credits__ = ("<NAME>, <NAME>, <NAME>, "
               "<NAME>, <NAME>, <NAME>")
__version__ = "1.0"
__maintainer__ = "<NAME>"
__email__ = "<EMAIL>"
__status__ = "Development"
__date__ = "March 22, 2012"


from scipy.linalg import polar
import numpy as np
import itertools
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.core.operations import SymmOp
from numpy.linalg import norm

class TensorBase(np.ndarray):
    """
    Base class for doing useful general operations on Nth order tensors,
    without restrictions on the type (stress, elastic, strain, piezo, etc.)
    """

    def __new__(cls, input_array):
        """
        Create a TensorBase object.  Note that the constructor uses __new__
        rather than __init__ according to the standard method of
        subclassing numpy ndarrays.

        Args:
            input_array: (3xN array-like): the 3xN array-like representing
                a tensor quantity
        """
        obj = np.asarray(input_array).view(cls)
        obj.rank = len(obj.shape)
        if not all([i == 3 for i in obj.shape]):
            raise ValueError("Pymatgen only supports 3-dimensional tensors")
        
        return obj

    def __array_finalize__(self, obj):
        if obj is None:
            return
        self.rank = getattr(obj, 'rank', None)

    def __array_wrap__(self, obj):
        """
        Overrides __array_wrap__ methods in ndarray superclass to avoid errors
        associated with functions that return scalar values
        """

        if len(obj.shape) == 0:
            return obj[()]
        else:
            return np.ndarray.__array_wrap__(self, obj)
        
    def __hash__(self):
        """
        define a hash function, since numpy arrays
        have their own __eq__ method
        """
        return hash(self.tostring())

    def __repr__(cls):
        return "{}({})".format(cls.__class__.__name__,
                               cls.__str__())

    def zeroed(self, tol = 1e-3):
        """
        returns the matrix with all entries below a certain threshold
        (i.e. tol) set to zero
        """
        new_tensor = self.copy()
        new_tensor[abs(new_tensor) < tol] = 0
        return new_tensor

    def transform(self, symm_op):
        """
        Applies a transformation (via a symmetry operation) to a tensor. 

        Args:
            symm_op (SymmOp): a symmetry operation to apply to the tensor
        """
        return self.__class__(symm_op.transform_tensor(self))

    def rotate(self, matrix, tol=1e-5):
        """
        Applies a rotation directly, and tests input matrix to ensure a valid
        rotation.

        Args:
            matrix (3x3 array-like): rotation matrix to be applied to tensor
            tol (float): tolerance for testing rotation matrix validity
        """
        matrix = SquareTensor(matrix)
        if not matrix.is_rotation(tol):
            raise ValueError("Rotation matrix is not valid.")
        sop = SymmOp.from_rotation_and_translation(matrix,
                                                   [0., 0., 0.])
        return self.transform(sop)


    @property
    def symmetrized(self):
        """
        Returns a generally symmetrized tensor, calculated by taking 
        the sum of the tensor and its transpose with respect to all 
        possible permutations of indices
        """
        perms = list(itertools.permutations(range(self.rank)))
        return sum([np.transpose(self, ind) for ind in perms]) / len(perms)

    def is_symmetric(self, tol=1e-5):
        """
        Tests whether a tensor is symmetric or not based on the residual
        with its symmetric part, from self.symmetrized

        Args:
            tol (float): tolerance to test for symmetry
        """
        return (self - self.symmetrized < tol).all()

    def fit_to_structure(self, structure, symprec = 0.1):
        """
        Returns a tensor that is invariant with respect to symmetry
        operations corresponding to a structure

        Args: 
            structure (Structure): structure from which to generate 
                symmetry operations
            symprec (float): symmetry tolerance for the Spacegroup Analyzer
                used to generate the symmetry operations
        """
        sga = SpacegroupAnalyzer(structure, symprec)
        symm_ops = sga.get_symmetry_operations(cartesian=True)
        return sum([self.transform(symm_op)
                    for symm_op in symm_ops]) / len(symm_ops)

    def is_fit_to_structure(self, structure, tol=1e-2):
        """
        Tests whether a tensor is invariant with respect to the
        symmetry operations of a particular structure by testing
        whether the residual of the symmetric portion is below a 
        tolerance
        
        Args:
            tol (float): tolerance for symmetry testing
        """
        return (self - self.fit_to_structure(structure) < tol).all()

    def convert_to_ieee(self, structure, atol = 0.1, ltol = 1e-3):
        """
        Given a structure associated with a tensor, attempts a
        calculation of the tensor in IEEE format according to
        the 1978 IEEE standards.

        Args:
            structure (structure): a structure associated with the
            symprec (float): precision for the symmetry analyzer
        """
        def get_uvec(vec):
            """ Gets a unit vector parallel to input vector"""
            return vec / np.linalg.norm(vec)

        vecs = structure.lattice.matrix
        lengths = np.array(structure.lattice.abc)
        angles = np.array(structure.lattice.angles)

        a = b = c = None
        sga = SpacegroupAnalyzer(structure)
        xtal_sys = sga.get_crystal_system()
        angle_error = "Angle check failure. {} crystal may not be in "\
                      "conventional representation".format(xtal_sys)
        rotation = np.zeros((3,3))

        # IEEE rules: a,b,c || x1,x2,x3
        if xtal_sys == "cubic":
            if (abs(angles - 90.0) > atol).any():
                raise ValueError(angle_error)
            if (abs(lengths - lengths[0]) > ltol).any():
                raise ValueError("Cubic vectors not equal in length")
            rotation = [vecs[i]/lengths[i] for i in range(3)]
        
        # IEEE rules: a=b in length; c,a || x3, x1
        elif xtal_sys == "tetragonal":
            if (abs(angles - 90.0) > atol).any():
                raise ValueError(angle_error)
            # This may be redundant, should be cubic
            if (abs(lengths - lengths[0]) < ltol).all():
                raise ValueError("Tetragonal c vector indistinguishable")

            rotation = [vecs[i]/lengths[i] for i in range(3)]
            if abs(lengths[2] - lengths[0]) < ltol:
                rotation[1], rotation[2] = rotation[2], rotation[1].copy()
            elif abs(lengths[2] - lengths[1]) < ltol:
                rotation[0], rotation[2] = rotation[2], rotation[0].copy()
            else:
                raise ValueError("No latt. vecs equal for tetragonal system")

        # IEEE rules: c<a<b; c,a || x3,x1
        elif xtal_sys == "orthorhombic":
            if (abs(angles - 90.0) > atol).any():
                raise ValueError(angle_error)
            # May be redundant, should be tetragonal
            if len(np.unique(np.floor(lengths/ltol))) < 3:
                raise ValueError("Orthorhombic vectors indistinguishable")
            # TODO: Check this
            rotation = [vec/mag for (mag, vec) in sorted(zip(lengths, vecs))]
            rotation = np.roll(rotation, 2, axis = 0)

        # IEEE rules: c,a || x3,x1, c is threefold axis
        elif xtal_sys in ("trigonal", "hexagonal"):
            # Rhombohedral lattice
            if xtal_sys == "trigonal" and (angles - angles[0] < atol).all():
                rotation[0] = get_uvec(vecs[2] - vecs[0])
                rotation[2] = get_uvec(np.sum(vecs, axis=0))
                rotation[1] = get_uvec(np.cross(rotation[2], rotation[0]))
            # Standard hexagonal lattice
            else:
                # find threefold axis:
                tf_mask = np.where(abs(angles-120.0) < atol)
                non_tf_mask = np.logical_not(tf_mask)
                rotation[2] = get_uvec(vecs[tf_mask])
                rotation[0] = get_uvec(vecs[non_tf_mask][0])
                rotation[1] = get_uvec(np.cross(rotation[2], rotation[0]))

        # IEEE rules: b,c || x2,x3; alpha=beta=90, c<a
        elif xtal_sys == "monoclinic":
            # Find unique axis
            umask = angles - 90.0 > atol
            n_umask = np.logical_not(umask)
            if np.sum(umask) != 1:
                raise ValueError(angle_error)
            rotation[1] = get_uvec(vecs[umask])
            # Shorter of remaining lattice vectors for c axis
            c = [vec/mag for (mag, vec) in 
                 sorted(zip(lengths[n_umask], vecs[n_umask]))][0]
            rotation[2] = np.array(c)
            rotation[0] = np.cross(rotation[1], rotation[2])
        
        # IEEE rules: c || x3
        elif xtal_sys == "triclinic":
            rotation = [vec/mag for (mag, vec) in sorted(zip(lengths, vecs))]
            rotation = np.roll(rotation, 2, axis = 0)
            rotation[1] = get_uvec(np.cross(rotation[2], rotation[1]))
            rotation[0] = np.cross(rotation[1], rotation[2])

        else:
            print("{} crystal not implemented".format(xtal_sys))
            return np.zeros(self.shape)
        
        return self.rotate(rotation)


class SquareTensor(TensorBase):
    """
    Base class for doing useful general operations on second rank tensors
    (stress, strain etc.).
    """

    def __new__(cls, input_array):
        """
        Create a SquareTensor object.  Note that the constructor uses __new__
        rather than __init__ according to the standard method of
        subclassing numpy ndarrays.  Error is thrown when the class is
        initialized with non-square matrix.

        Args:
            stress_matrix (3x3 array-like): the 3x3 array-like
                representing the Green-Lagrange strain
        """

        obj = TensorBase(input_array).view(cls)
        if not (len(obj.shape) == 2):
            raise ValueError("SquareTensor only takes 2-D "
                             "tensors as input")
        return obj
        
    @property
{'fraction_non_alphanumeric': 0.0966822429906542, 'fraction_numerical': 0.01177570093457944, 'mean_word_length': 2.064728626664757, 'pattern_counts': {'":': 10, '<': 45, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '809036', 'n_tokens_mistral': 6249, 'n_tokens_neox': 5946, 'n_words': 1869}