diff --git a/groupy/garray/C4h_array.py b/groupy/garray/C4h_array.py new file mode 100644 index 0000000..62a2130 --- /dev/null +++ b/groupy/garray/C4h_array.py @@ -0,0 +1,148 @@ +import numpy as np +from groupy.garray.finitegroup import FiniteGroup +from groupy.garray.matrix_garray import MatrixGArray +from groupy.garray.C4ht_array import C4htArray +from groupy.garray.Z3_array import Z3Array + +""" +Implementation of the finite group C4h, the cyclic group of rectangular cuboid symmetry, consisting of 8 elements in total: + - Identity element + - one 4-fold axis over opposing square faces (90 degree rotations) + - two 2-fold axes over opposing rectangular faces (180 degree rotations) + - two 2-fold axes over opposing long edges (180 degree rotations) +As C4h is abelian, each transformation can be described in terms of number of rotations around the various fold axes. +However, all elements from C4h can also be generated by two elements: 180 degree rotation over the y-axis and 90 +degree rotation over the z-axis. + +The int parameterization is in the form of (y, z) where y represents the number of 180 degree rotations +over y (in {0, 1}) and z the number of 90 degree rotations over z (in {0, 1, 2, 3}) +""" + + +class C4hArray(MatrixGArray): + parameterizations = ['int', 'mat', 'hmat'] + _g_shapes = {'int': (2,), 'mat': (3, 3), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'C4h' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + + # classes C4hArray can be multiplied with + self._left_actions[C4hArray] = self.__class__.left_action_hmat + self._left_actions[C4htArray] = self.__class__.left_action_hmat + self._left_actions[Z3Array] = self.__class__.left_action_vec + + super(C4hArray, self).__init__(data, p) + self.elements = self.get_elements() + + def mat2int(self, mat_data): + ''' + Transforms 3x3 matrix representation to int representation. + To handle any size and shape of mat_data, the original mat_data + is reshaped to a long list of 3x3 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + + mat-2-int is achieved by taking the matrix, looking up the index in the + element list, and converting that index to two numbers: y and z. The index + is the result of (y * 4) + z. + ''' + + input = mat_data.reshape((-1, 3, 3)) + data = np.zeros((input.shape[0], 2), dtype=np.int) + for i in xrange(input.shape[0]): + index = self.elements.index(input[i].tolist()) + z = int(index % 4) + y = int((index - z) / 4) + data[i, 0] = y + data[i, 1] = z + data = data.reshape(mat_data.shape[:-2] + (2,)) + return data + + def int2mat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + + The element is located in the list at index (y * 4) + z. + ''' + y = int_data[..., 0].flatten() + z = int_data[..., 1].flatten() + data = np.zeros((len(y),) + (3, 3), dtype=np.int) + + for j in xrange(len(y)): + index = (y[j] * 4) + z[j] + mat = self.elements[index] + data[j, 0:3, 0:3] = mat + + data = data.reshape(int_data.shape[:-1] + (3, 3)) + return data + + def _multiply(self, element, generator, times): + element = np.array(element) + for i in range(times): + element = np.dot(element, np.array(generator)) + return element + + def get_elements(self): + ''' + Elements are stored as lists rather than numpy arrays to enable + lookup through self.elements.index(x) and sorting. + + All elements are found by multiplying the identity matrix with all + possible combinations of the generators, i.e. 0 or 1 rotations over y + and 0, 1, 2, or 3 rotations over z. + ''' + # specify generators + mode = 'zyx' + if mode == 'xyz': + g1 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # 180 degrees over y + g2 = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) # 90 degrees over z + elif mode == 'zyx': + g1 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # 180 degrees over y + g2 = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) # 90 degrees over z + + + element_list = [] + element = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + for i in range(0, 2): + element = self._multiply(element, g1, i) + for j in range(0, 4): + element = self._multiply(element, g2, j) + element_list.append(element.tolist()) + return element_list + + +class C4hGroup(FiniteGroup, C4hArray): + def __init__(self): + C4hArray.__init__( + self, + data=np.array([[i, j] for i in xrange(2) for j in xrange(4)]), + p='int' + ) + FiniteGroup.__init__(self, C4hArray) + + def factory(self, *args, **kwargs): + return C4hArray(*args, **kwargs) + +C4h = C4hGroup() + +def rand(size=()): + ''' + Returns an C4hArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (2,), dtype=np.int) + data[..., 0] = np.random.randint(0, 2, size) # rotations over y + data[..., 1] = np.random.randint(0, 4, size) # rotations over z + return C4hArray(data=data, p='int') + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] + e = C4hArray(data=np.array(li, dtype=np.int), p='mat') + return e.reparameterize(p) diff --git a/groupy/garray/C4ht_array.py b/groupy/garray/C4ht_array.py new file mode 100644 index 0000000..d7b2e30 --- /dev/null +++ b/groupy/garray/C4ht_array.py @@ -0,0 +1,151 @@ +import numpy as np +from groupy.garray.matrix_garray import MatrixGArray + +""" +Implementation of the non-orientation perserving variant of group C4h -- C4h with translations. +The int parameterisation is similar to that of C4h, but with the added 3D translation (u, v, w) to indicate +translation in Z3. + +4x4 homogeneous matrices (hmat) are used to represent the transformation in matrix format. +""" + +class C4htArray(MatrixGArray): + parameterizations = ['int', 'hmat'] + _g_shapes = {'int': (5,), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'C4ht' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + self._left_actions[C4htArray] = self.__class__.left_action_hmat + super(C4htArray, self).__init__(data, p) + self.elements = self.get_elements() + + def hmat2int(self, hmat_data): + ''' + Transforms 4x4 matrix representation to int representation. + To handle any size and shape of hmat_data, the original hmat_data + is reshaped to a long list of 4x4 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + + hmat-2-int is achieved by taking the matrix, looking up the index in the + element list, and converting that index to two numbers: y and z. The index + is the result of (y * 4) + z. u, v, w are retrieved by looking at the last + column of the hmat. + ''' + + input = hmat_data.reshape((-1, 4, 4)) + data = np.zeros((input.shape[0], 5), dtype=np.int) + for i in xrange(input.shape[0]): + hmat = input[i] + mat = [elem[0:3] for elem in hmat.tolist()][0:3] + index = self.elements.index(mat) + z = int(index % 4) + y = int((index - z) / 4) + u, v, w, _ = hmat[:, 3] + data[i, 0] = y + data[i, 1] = z + data[i, 2] = u + data[i, 3] = v + data[i, 4] = w + data = data.reshape(hmat_data.shape[:-2] + (5,)) + return data + + def int2hmat(self, int_data): + ''' + Transforms integer representation to 4x4 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + # rotations over y and z + y = int_data[..., 0].flatten() + z = int_data[..., 1].flatten() + + # translations + u = int_data[..., 2].flatten() + v = int_data[..., 3].flatten() + w = int_data[..., 4].flatten() + data = np.zeros((len(y),) + (4, 4), dtype=np.int) + + for j in xrange(len(y)): + index = (y[j] * 4) + z[j] + mat = self.elements[index] + + data[j, 0:3, 0:3] = mat + data[j, 0, 3] = u[j] + data[j, 1, 3] = v[j] + data[j, 2, 3] = w[j] + data[j, 3, 3] = 1 + + data = data.reshape(int_data.shape[:-1] + (4, 4)) + return data + + def _multiply(self, element, generator, times): + ''' + Helper function to multiply an _element_ with a _generator_ + _times_ number of times. Used in self.get_elements() + ''' + element = np.array(element) + for i in range(times): + element = np.dot(element, np.array(generator)) + return element + + def get_elements(self): + ''' + Function to generate a list containing elements of group C4ht, + similar to get_elements() of BArray. + + These are the base elements in 3x3 matrix notation without translations. + ''' + # specify generators + mode = 'zyx' + if mode == 'xyz': + g1 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # 180 degrees over y + g2 = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) # 90 degrees over z + elif mode == 'zyx': + g1 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # 180 degrees over y + g2 = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) # 90 degrees over z + + element_list = [] + element = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + for i in range(0, 2): + element = self._multiply(element, g1, i) + for j in range(0, 4): + element = self._multiply(element, g2, j) + element_list.append(element.tolist()) + return element_list + + +def rand(minu=0, maxu=5, minv=0, maxv=5, minw=0, maxw=5, size=()): + ''' + Returns an C4htArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (5,), dtype=np.int64) + data[..., 0] = np.random.randint(0, 2, size) # rotations over y + data[..., 1] = np.random.randint(0, 4, size) # rotations over x + data[..., 2] = np.random.randint(minu, maxu, size) # translation on x + data[..., 3] = np.random.randint(minv, maxv, size) # translation on y + data[..., 4] = np.random.randint(minw, maxw, size) # translation on z + return C4htArray(data=data, p='int') + + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] + e = C4htArray(data=np.array(li, dtype=np.int), p='hmat') + return e.reparameterize(p) + + +def meshgrid(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2): + ''' + Creates a meshgrid of all elements of the group, within the given + translation parameters. + ''' + li = [[i, m, u, v, w] for i in xrange(2) for m in xrange(4) for u in xrange(minu, maxu) for v in xrange(minv, maxv) + for + w in xrange(minw, maxw)] + return C4htArray(li, p='int') diff --git a/groupy/garray/D4_array.py b/groupy/garray/D4_array.py index 3fe6fb2..1410c78 100644 --- a/groupy/garray/D4_array.py +++ b/groupy/garray/D4_array.py @@ -6,7 +6,6 @@ from groupy.garray.matrix_garray import MatrixGArray - class D4Array(MatrixGArray): parameterizations = ['int', 'mat', 'hmat'] diff --git a/groupy/garray/D4h_array.py b/groupy/garray/D4h_array.py new file mode 100644 index 0000000..1c0b4b9 --- /dev/null +++ b/groupy/garray/D4h_array.py @@ -0,0 +1,148 @@ +import numpy as np +from groupy.garray.finitegroup import FiniteGroup +from groupy.garray.matrix_garray import MatrixGArray +from groupy.garray.D4ht_array import D4htArray +from groupy.garray.Z3_array import Z3Array + +""" +Implementation of dihedral finite group D4h, consisting of 16 elements in total. +These are the elements of C4h, with added reflection. +Int parameterisation contains an extra parameter, m (in {0, 1}) to represent this reflection. +""" + + +class D4hArray(MatrixGArray): + parameterizations = ['int', 'mat', 'hmat'] + _g_shapes = {'int': (3,), 'mat': (3, 3), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'D4h' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + + # classes OArray can be multiplied with + self._left_actions[D4hArray] = self.__class__.left_action_hmat + self._left_actions[D4htArray] = self.__class__.left_action_hmat + self._left_actions[Z3Array] = self.__class__.left_action_vec + + super(D4hArray, self).__init__(data, p) + self.elements = self.get_elements() + + def mat2int(self, mat_data): + ''' + Transforms 3x3 matrix representation to int representation. + To handle any size and shape of mat_data, the original mat_data + is reshaped to a long list of 3x3 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + + mat-2-int is achieved by taking the matrix, and looking up whether it + exists in the element list. If not, the matrix should be multiplied with -1 + to retrieve the reflection. The resulting matrix can be looked up in the + element list, and that index can be converted to y and z. + ''' + + input = mat_data.reshape((-1, 3, 3)) + data = np.zeros((input.shape[0], 3), dtype=np.int) + for i in xrange(input.shape[0]): + mat = input[i] + # check for reflection + if mat.tolist() not in self.elements: + mat = np.array(mat) * -1 + data[i, 2] = 1 + + # determine z and y + index = self.elements.index(mat.tolist()) + z = int(index % 4) + y = int((index - z) / 4) + data[i, 0] = y + data[i, 1] = z + data = data.reshape(mat_data.shape[:-2] + (3,)) + return data + + def int2mat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + # rotations over y, z and reflection + y = int_data[..., 0].flatten() + z = int_data[..., 1].flatten() + m = int_data[..., 2].flatten() + data = np.zeros((len(y),) + (3, 3), dtype=np.int) + + for j in xrange(len(y)): + index = (y[j] * 4) + z[j] + mat = self.elements[index] + mat = np.array(mat) * ((-1) ** m[j]) # mirror if reflection + data[j, 0:3, 0:3] = mat.tolist() + + data = data.reshape(int_data.shape[:-1] + (3, 3)) + return data + + def _multiply(self, element, generator, times): + ''' + Helper function to multiply an _element_ with a _generator_ + _times_ number of times. + ''' + element = np.array(element) + for i in range(times): + element = np.dot(element, np.array(generator)) + return element + + def get_elements(self): + ''' + Function to generate a list containing elements of group D4h, + similar to get_elements() of BArray. + + Elements are stored as lists rather than numpy arrays to enable + lookup through self.elements.index(x). + ''' + # specify generators + g1 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # 180 degrees over y + g2 = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) # 90 degrees over z + + element_list = [] + element = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) # starting point = identity matrix + for i in range(0, 2): + element = self._multiply(element, g1, i) + for j in range(0, 4): + element = self._multiply(element, g2, j) + element_list.append(element.tolist()) + return element_list + + +class D4hGroup(FiniteGroup, D4hArray): + def __init__(self): + D4hArray.__init__( + self, + data=np.array([[i, j, m] for i in xrange(2) for j in xrange(4) for m in xrange(2)]), + p='int' + ) + FiniteGroup.__init__(self, D4hArray) + + def factory(self, *args, **kwargs): + return D4hArray(*args, **kwargs) + + +D4h = D4hGroup() + +def rand(size=()): + ''' + Returns an D4hArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (3,), dtype=np.int) + data[..., 0] = np.random.randint(0, 2, size) + data[..., 1] = np.random.randint(0, 4, size) + data[..., 2] = np.random.randint(0, 2, size) + return D4hArray(data=data, p='int') + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] + e = D4hArray(data=np.array(li, dtype=np.int), p='mat') + return e.reparameterize(p) diff --git a/groupy/garray/D4ht_array.py b/groupy/garray/D4ht_array.py new file mode 100644 index 0000000..32ed725 --- /dev/null +++ b/groupy/garray/D4ht_array.py @@ -0,0 +1,175 @@ +import numpy as np +from groupy.garray.matrix_garray import MatrixGArray + +''' +Implementation of space group D4h that allows translations. It has no official name, and is therefor referred to as D4ht. +Implementation is similar to that of group D4h. However, to represent the translations in 3D space, +the int parameterization is now in the form of (y, z, m, u, v, w) + +Implementation of the space group Oh that allows translations. +It has no official name, and is therefore now referred to as Oht. + +Implementation is similar to that of group Oh. However, to represent +the translations in a 3D space, the int parameterization is now +in the form of (i, m, u, v, w) representing the index of the unmirrored +element in the element list, the mirror (1 or -1) and the translation +in the x, y and z direction respectively. + +To accurately represent the translation, we use 4x4 homogeneous matrices +(hmat) instead of the 3x3 matrix. + +Note: self.base_elements are 3x3 matrices. +''' + + +class D4htArray(MatrixGArray): + parameterizations = ['int', 'hmat'] + _g_shapes = {'int': (6,), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'D4ht' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + self._left_actions[D4htArray] = self.__class__.left_action_hmat + super(D4htArray, self).__init__(data, p) + self.elements = self.get_elements() + + def hmat2int(self, hmat_data): + ''' + Transforms 4x4 matrix representation to int representation. + To handle any size and shape of hmat_data, the original hmat_data + is reshaped to a long list of 4x4 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + + + hmat-2-int is achieved by taking the matrix, and looking up whether it + exists in the element list. If not, the matrix should be multiplied with -1 + to retrieve the reflection. The resulting matrix can be looked up in the + element list, and that index can be converted to y and z. u, v, and w + are retrieved by looking at the last column in the matrix. + ''' + + input = hmat_data.reshape((-1, 4, 4)) + data = np.zeros((input.shape[0], 6), dtype=np.int) + for i in xrange(input.shape[0]): + hmat = input[i] + mat = [elem[0:3] for elem in hmat.tolist()][0:3] + # check for reflection + if mat not in self.elements: + mat = (np.array(mat) * -1).tolist() + data[i, 2] = 1 # reflection + + # retrieve values + index = self.elements.index(mat) + z = int(index % 4) + y = int((index - z) / 4) + u, v, w, _ = hmat[:, 3] + + # rotations over y and z + data[i, 0] = y + data[i, 1] = z + + # translation + data[i, 3] = u + data[i, 4] = v + data[i, 5] = w + + data = data.reshape(hmat_data.shape[:-2] + (6,)) + return data + + def int2hmat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + # rotatiins over y, z and reflection + y = int_data[..., 0].flatten() + z = int_data[..., 1].flatten() + m = int_data[..., 2].flatten() + + # translations + u = int_data[..., 3].flatten() + v = int_data[..., 4].flatten() + w = int_data[..., 5].flatten() + + data = np.zeros((len(y),) + (4, 4), dtype=np.int) + + for j in xrange(len(y)): + index = (y[j] * 4) + z[j] + mat = self.elements[index] + mat = np.array(mat) * ((-1) ** m[j]) # mirror if reflection + + data[j, 0:3, 0:3] = mat + data[j, 0, 3] = u[j] + data[j, 1, 3] = v[j] + data[j, 2, 3] = w[j] + data[j, 3, 3] = 1 + + data = data.reshape(int_data.shape[:-1] + (4, 4)) + return data + + def _multiply(self, element, generator, times): + ''' + Helper function to multiply an _element_ with a _generator_ + _times_ number of times. + ''' + element = np.array(element) + for i in range(times): + element = np.dot(element, np.array(generator)) + return element + + def get_elements(self): + ''' + Function to generate a list containing elements of group D4ht, + similar to get_elements() of BArray. + + These are the base elements in 3x3 matrix notation without translations. + ''' + # specify generators + g1 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # 180 degrees over y + g2 = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) # 90 degrees over z + + element_list = [] + element = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + for i in range(0, 2): + element = self._multiply(element, g1, i) + for j in range(0, 4): + element = self._multiply(element, g2, j) + element_list.append(element.tolist()) + return element_list + + +def rand(minu=0, maxu=5, minv=0, maxv=5, minw=0, maxw=5, size=()): + ''' + Returns an D4htArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (6,), dtype=np.int64) + data[..., 0] = np.random.randint(0, 2, size) # 180 degree rotation over y + data[..., 1] = np.random.randint(0, 4, size) # 90 degree rotation over z + data[..., 2] = np.random.randint(0, 2, size) # reflection or not + data[..., 3] = np.random.randint(minu, maxu, size) # translation on x + data[..., 4] = np.random.randint(minv, maxv, size) # translation on y + data[..., 5] = np.random.randint(minw, maxw, size) # translation on z + return D4htArray(data=data, p='int') + + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] + e = D4htArray(data=np.array(li, dtype=np.int), p='hmat') + return e.reparameterize(p) + + +def meshgrid(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2): + ''' + Creates a meshgrid of all elements of the group, within the given + translation parameters. + ''' + li = [[y, z, m, u, v, w] for y in xrange(2) for z in xrange(4) for m in xrange(2) for u in xrange(minu, maxu) for v + in xrange(minv, maxv) for w in xrange(minw, maxw)] + return D4htArray(li, p='int') diff --git a/groupy/garray/O_array.py b/groupy/garray/O_array.py new file mode 100644 index 0000000..32ac17c --- /dev/null +++ b/groupy/garray/O_array.py @@ -0,0 +1,135 @@ +import random +import numpy as np +from groupy.garray.Ot_array import OtArray +from groupy.garray.Z3_array import Z3Array +from groupy.garray.finitegroup import FiniteGroup +from groupy.garray.matrix_garray import MatrixGArray + +""" +Implementation of finite group O, the group of octahedral symmetry. +Group O represents the 24 symmetries of a octahedron and can be constructed +by combining its generators (90 degree rotation over the x- and y-axis). + +All elements belonging to the group are generated by self.get_elements(), +and saved to a list (self.elements). Because O is not commutative, +it cannot simply be parameterised in terms of number of times the generators +are combined. Therefore, the int representation is the index of the transformation +in the list of elements generated at initialisation of this class. +""" + + +class OArray(MatrixGArray): + parameterizations = ['int', 'mat', 'hmat'] + _g_shapes = {'int': (1,), 'mat': (3, 3), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'O' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + + # classes OArray can be multiplied with + self._left_actions[OArray] = self.__class__.left_action_hmat + self._left_actions[OtArray] = self.__class__.left_action_hmat + self._left_actions[Z3Array] = self.__class__.left_action_vec + + super(OArray, self).__init__(data, p) + self.elements = self.get_elements() + + def mat2int(self, mat_data): + ''' + Transforms 3x3 matrix representation to int representation. + To handle any size and shape of mat_data, the original mat_data + is reshaped to a long list of 3x3 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + ''' + + input = mat_data.reshape((-1, 3, 3)) + data = np.zeros((input.shape[0], 1), dtype=np.int) + for i in xrange(input.shape[0]): + mat = input[i] + index = self.elements.index(mat.tolist()) + data[i, 0] = index + data = data.reshape(mat_data.shape[:-2] + (1,)) + return data + + def int2mat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + i = int_data[..., 0].flatten() + data = np.zeros((len(i),) + (3, 3), dtype=np.int) + + for j in xrange(len(i)): + mat = self.elements[i[j]] + data[j, 0:3, 0:3] = mat + + data = data.reshape(int_data.shape[:-1] + (3, 3)) + return data + + def get_elements(self): + ''' + Function to generate a list containing all 24 elements of + group O. This list is necessary for the integer representation: + the index of an element in this list represents the group element. + + Elements are stored as lists rather than numpy arrays to enable + lookup through self.elements.index(x) and sorting. + + Basic principle to find all elements in the list: + specify the generators (rotations over x and y axis) + while not all elements have been found, repeat: + choose random element from current list of elements as multiplier + multiply this element with the last found element + if this element is new, add it to the list of elements. + ''' + # specify generators + g1 = [[1, 0, 0], [0, 0, -1], [0, 1, 0]] # 90o degree rotation over x + g2 = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]] # 90o degree rotation over y + # initialise element list + element_list = [g1, g2] + current = g1 + while len(element_list) < 24: + multiplier = random.choice(element_list) + current = np.dot(np.array(current), np.array(multiplier)).tolist() + if current not in element_list: + element_list.append(current) + element_list.sort() + return element_list + + +class OGroup(FiniteGroup, OArray): + def __init__(self): + OArray.__init__( + self, + data=np.arange(24)[:, None], # 24 elements + p='int' + ) + FiniteGroup.__init__(self, OArray) + + def factory(self, *args, **kwargs): + return OArray(*args, **kwargs) + + +O = OGroup() + + +def rand(size=()): + ''' + Returns an OArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (1,), dtype=np.int64) + data[..., 0] = np.random.randint(0, 24, size) + return OArray(data=data, p='int') + + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] + e = OArray(data=np.array(li, dtype=np.int), p='mat') + return e.reparameterize(p) diff --git a/groupy/garray/Oh_array.py b/groupy/garray/Oh_array.py new file mode 100644 index 0000000..b6c09cd --- /dev/null +++ b/groupy/garray/Oh_array.py @@ -0,0 +1,159 @@ +import copy +import random + +import numpy as np +from groupy.garray.Oht_array import OhtArray +from groupy.garray.Z3_array import Z3Array +from groupy.garray.finitegroup import FiniteGroup +from groupy.garray.matrix_garray import MatrixGArray + +""" +Implementation of finite group Oh, the group of octahedral symmetry. +Group O represents the 24 symmetries of a octahedron and the mirrors thereof. +This group can be constructed by combining its generators (90 degree rotation +over the x- and y-axis) with all mirros of these elements. + +All elements belonging to the group are generated by self.get_elements(), +and saved to a list (self.elements). The int representation is (i, m) where +i represents the index of the unmirrored element (similar to O) in this list +and m (in {0, 1}) whether or not the element is a mirror. +""" + + +class OhArray(MatrixGArray): + parameterizations = ['int', 'mat', 'hmat'] + _g_shapes = {'int': (2,), 'mat': (3, 3), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'Oh' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + self._left_actions[OhArray] = self.__class__.left_action_hmat + self._left_actions[OhtArray] = self.__class__.left_action_hmat + self._left_actions[Z3Array] = self.__class__.left_action_vec + super(OhArray, self).__init__(data, p) + self.elements = self.get_elements() + + def mat2int(self, mat_data): + ''' + Transforms 3x3 matrix representation to int representation. + To handle any size and shape of mat_data, the original mat_data + is reshaped to a long list of 3x3 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + ''' + + input = mat_data.reshape((-1, 3, 3)) + data = np.zeros((input.shape[0], 2), dtype=np.int) + for i in xrange(input.shape[0]): + mat = input[i] + index, mirror = self.get_int(mat) + data[i, 0] = index + data[i, 1] = mirror + data = data.reshape(mat_data.shape[:-2] + (2,)) + return data + + def get_int(self, mat_data): + ''' + Return int (index, mirror) representation of given mat + by mirroring if necessary to find the original mat and + looking up the index in the list of base elements + ''' + assert mat_data.shape == (3, 3) + orig_data = copy.deepcopy(mat_data) + m = 0 if orig_data.tolist() in self.elements else 1 + orig_data[0:3] = orig_data[0:3] * ((-1) ** m) + i = self.elements.index(orig_data.tolist()) + return i, m + + def int2mat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + index = int_data[..., 0].flatten() + m = int_data[..., 1].flatten() + data = np.zeros((len(index),) + (3, 3), dtype=np.int) + + for j in xrange(len(index)): + hmat = self.get_mat(index[j], m[j]) + data[j, 0:3, 0:3] = hmat + + data = data.reshape(int_data.shape[:-1] + (3, 3)) + return data + + def get_mat(self, index, mirror): + ''' + Return matrix representation of a given int parameterization (index, mirror) + by determining looking up the mat by index and mirroring if necessary + (note: deepcopy to avoid alterations to original self.base_elements) + ''' + element = copy.deepcopy(self.elements[index]) + element = np.array(element, dtype=np.int) + element[0:3] = element[0:3] * ((-1) ** mirror) + element = element.astype(dtype=np.int) + return element + + def get_elements(self): + ''' + Function to generate a list containing elements of group Oh, + similar to get_elements() of OArray. However, group Oh also + includes the mirrors of these elements. + ''' + g1 = [[1, 0, 0], [0, 0, -1], [0, 1, 0]] # 90o degree rotation over x + g2 = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]] # 90o degree rotation over y + + element_list = [g1, g2] + current = g1 + while len(element_list) < 24: + multiplier = random.choice(element_list) + current = np.dot(np.array(current), np.array(multiplier)).tolist() + if current not in element_list: + element_list.append(current) + element_list.sort() + return element_list + + +class OhGroup(FiniteGroup, OhArray): + def __init__(self): + OhArray.__init__( + self, + data=np.array([[i, j] for i in xrange(24) for j in xrange(2)]), + p='int' + ) + FiniteGroup.__init__(self, OhArray) + + def factory(self, *args, **kwargs): + return OhArray(*args, **kwargs) + + +Oh = OhGroup() + + +def rand(size=()): + ''' + Returns an OhArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (2,), dtype=np.int) + data[..., 0] = np.random.randint(0, 24, size) + data[..., 1] = np.random.randint(0, 2, size) + return OhArray(data=data, p='int') + + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] + e = OhArray(data=np.array(li, dtype=np.int), p='mat') + return e.reparameterize(p) + + +def meshgrid(i=24, m=2): + ''' + Creates a meshgrid of all elements of the group, within the given + translation parameters. + ''' + return OhArray([[[k, l] for l in xrange(m)] for k in xrange(i)]) diff --git a/groupy/garray/Oht_array.py b/groupy/garray/Oht_array.py new file mode 100644 index 0000000..9fd52de --- /dev/null +++ b/groupy/garray/Oht_array.py @@ -0,0 +1,152 @@ +import copy +import random + +import numpy as np +from groupy.garray.matrix_garray import MatrixGArray + +""" +Implementation of the non-orientation perserving variant of group Oh -- O hwith translations. +The int parameterisation is similar to that of Oh, but with the added 3D translation (u, v, w) to indicate +translation in Z3 (i.e. i, m, u, v, w). + +4x4 homogeneous matrices (hmat) are used to represent the transformation in matrix format. +""" + + +class OhtArray(MatrixGArray): + parameterizations = ['int', 'hmat'] + _g_shapes = {'int': (5,), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'Oht' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + self._left_actions[OhtArray] = self.__class__.left_action_hmat + super(OhtArray, self).__init__(data, p) + self.base_elements = self.get_base_elements() + + def hmat2int(self, hmat_data): + ''' + Transforms 4x4 matrix representation to int representation. + To handle any size and shape of hmat_data, the original hmat_data + is reshaped to a long list of 4x4 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + ''' + + input = hmat_data.reshape((-1, 4, 4)) + data = np.zeros((input.shape[0], 5), dtype=np.int) + for i in xrange(input.shape[0]): + hmat = input[i] + mat = [elem[0:3] for elem in hmat.tolist()][0:3] + index, mirror = self.get_int(mat) + u, v, w, _ = hmat[:, 3] + data[i, 0] = index + data[i, 1] = mirror + data[i, 2] = u + data[i, 3] = v + data[i, 4] = w + data = data.reshape(hmat_data.shape[:-2] + (5,)) + return data + + def int2hmat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + i = int_data[..., 0].flatten() + m = int_data[..., 1].flatten() + u = int_data[..., 2].flatten() + v = int_data[..., 3].flatten() + w = int_data[..., 4].flatten() + data = np.zeros((len(i),) + (4, 4), dtype=np.int) + + for j in xrange(len(i)): + mat = self.get_mat(i[j], m[j]) + data[j, 0:3, 0:3] = mat + data[j, 0, 3] = u[j] + data[j, 1, 3] = v[j] + data[j, 2, 3] = w[j] + data[j, 3, 3] = 1 + + data = data.reshape(int_data.shape[:-1] + (4, 4)) + return data + + def get_mat(self, index, mirror): + ''' + Return matrix representation of a given int parameterization (index, mirror) + by determining looking up the mat by index and mirroring if necessary + (note: deepcopy to avoid alterations to original self.base_elements) + ''' + element = copy.deepcopy(self.base_elements[index]) + element = np.array(element, dtype=np.int) + element = element * ((-1) ** mirror) + return element + + def get_int(self, hmat_data): + ''' + Return int (index, mirror) representation of given mat + by mirroring if necessary to find the original mat and + looking up the index in the list of base elements + ''' + orig_data = copy.deepcopy(hmat_data) + m = 0 if orig_data in self.base_elements else 1 + orig_data = np.array(orig_data) * ((-1) ** m) + i = self.base_elements.index(orig_data.tolist()) + return i, m + + def get_base_elements(self): + ''' + Function to generate a list containing elements of group Oh, + similar to get_elements() of OArray. However, group Oh also + includes the mirrors of these elements. + + These are the base elements in 3x3 matrix notation without translations. + ''' + g1 = [[1, 0, 0], [0, 0, -1], [0, 1, 0]] # 90o degree rotation over x + g2 = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]] # 90o degree rotation over y + + element_list = [g1, g2] + current = g1 + while len(element_list) < 24: + multiplier = random.choice(element_list) + current = np.dot(np.array(current), np.array(multiplier)).tolist() + if current not in element_list: + element_list.append(current) + element_list.sort() + return element_list + + +def rand(minu, maxu, minv, maxv, minw, maxw, size=()): + ''' + Returns an OhtArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (5,), dtype=np.int64) + data[..., 0] = np.random.randint(0, 24, size) + data[..., 1] = np.random.randint(0, 2, size) + data[..., 2] = np.random.randint(minu, maxu, size) + data[..., 3] = np.random.randint(minv, maxv, size) + data[..., 4] = np.random.randint(minw, maxw, size) + return OhtArray(data=data, p='int') + + +def identity(p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] + e = OhtArray(data=np.array(li, dtype=np.int), p='hmat') + return e.reparameterize(p) + + +def meshgrid(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2): + ''' + Creates a meshgrid of all elements of the group, within the given + translation parameters. + ''' + li = [[i, m, u, v, w] for i in xrange(24) for m in xrange(2) for u in xrange(minu, maxu) for v in xrange(minv, maxv) + for + w in xrange(minw, maxw)] + return OhtArray(li, p='int') diff --git a/groupy/garray/Ot_array.py b/groupy/garray/Ot_array.py new file mode 100644 index 0000000..73bed83 --- /dev/null +++ b/groupy/garray/Ot_array.py @@ -0,0 +1,131 @@ +import random + +import numpy as np +from groupy.garray.matrix_garray import MatrixGArray + +""" +Implementation of the non-orientation perserving variant of group O -- O with translations. +The int parameterisation is similar to that of O, but with the added 3D translation (u, v, w) to indicate +translation in Z3 (i.e. i, u, v, w). + +4x4 homogeneous matrices (hmat) are used to represent the transformation in matrix format. +""" + +class OtArray(MatrixGArray): + ''' + Implementation of space group Ot. + ''' + parameterizations = ['int', 'hmat'] + _g_shapes = {'int': (4,), 'hmat': (4, 4)} + _left_actions = {} + _reparameterizations = {} + _group_name = 'Ot' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + self._left_actions[OtArray] = self.__class__.left_action_hmat + super(OtArray, self).__init__(data, p) + self.base_elements = self.get_base_elements() + + def hmat2int(self, hmat_data): + ''' + Transforms 4x4 matrix representation to int representation. + To handle any size and shape of hmat_data, the original hmat_data + is reshaped to a long list of 4x4 matrices, converted to a list of + int representations, and reshaped back to the original mat_data shape. + ''' + + input = hmat_data.reshape((-1, 4, 4)) + data = np.zeros((input.shape[0], 4), dtype=np.int) + for i in xrange(input.shape[0]): + hmat = input[i] + mat = [elem[0:3] for elem in hmat.tolist()][0:3] + index = self.base_elements.index(mat) + u, v, w, _ = hmat[:, 3] + data[i, 0] = index + data[i, 1] = u + data[i, 2] = v + data[i, 3] = w + data = data.reshape(hmat_data.shape[:-2] + (4,)) + return data + + def int2hmat(self, int_data): + ''' + Transforms integer representation to 3x3 matrix representation. + Original int_data is flattened and later reshaped back to its original + shape to handle any size and shape of input. + ''' + i = int_data[..., 0].flatten() + u = int_data[..., 1].flatten() + v = int_data[..., 2].flatten() + w = int_data[..., 3].flatten() + data = np.zeros((len(i),) + (4, 4), dtype=np.int) + + for j in xrange(len(i)): + mat = self.base_elements[i[j]] + data[j, 0:3, 0:3] = mat + data[j, 0, 3] = u[j] + data[j, 1, 3] = v[j] + data[j, 2, 3] = w[j] + data[j, 3, 3] = 1 + + data = data.reshape(int_data.shape[:-1] + (4, 4)) + return data + + def get_base_elements(self): + ''' + Function to generate a list containing all 24 elements of + group O. This list is necessary for the integer representation: + the index of an element in this list represents the group element. + + Elements are stored as lists rather than numpy arrays to enable + lookup through self.elements.index(x) and sorting. + + Basic principle to find all elements in the list: + specify the generators (rotations over x and y axis) + while not all elements have been found, repeat: + choose random element from current list of elements as multiplier + multiply this element with the last found element + if this element is new, add it to the list of elements. + + These are the base elements in 3x3 matrix notation without translations. + ''' + g1 = [[1, 0, 0], [0, 0, -1], [0, 1, 0]] # 90o degree rotation over x + g2 = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]] # 90o degree rotation over y + element_list = [g1, g2] + current = g1 + while len(element_list) < 24: + multiplier = random.choice(element_list) + current = np.dot(np.array(current), np.array(multiplier)).tolist() + if current not in element_list: + element_list.append(current) + element_list.sort() + return element_list + + +def identity(shape=(), p='int'): + ''' + Returns the identity element: a matrix with 1's on the diagonal. + ''' + li = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] + e = OtArray(data=np.array(li, dtype=np.int), p='hmat') + return e.reparameterize(p) + + +def rand(minu, maxu, minv, maxv, minw, maxw, size=()): + ''' + Returns an OtArray of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (4,), dtype=np.int64) + data[..., 0] = np.random.randint(0, 24, size) + data[..., 1] = np.random.randint(minu, maxu, size) + data[..., 2] = np.random.randint(minv, maxv, size) + data[..., 3] = np.random.randint(minw, maxw, size) + return OtArray(data=data, p='int') + + +def meshgrid(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2): + li = [[i, u, v, w] for i in xrange(24) for u in xrange(minu, maxu) for v in xrange(minv, maxv) for + w in xrange(minw, maxw)] + return OtArray(li, p='int') diff --git a/groupy/garray/Z3_array.py b/groupy/garray/Z3_array.py new file mode 100644 index 0000000..bccca0b --- /dev/null +++ b/groupy/garray/Z3_array.py @@ -0,0 +1,59 @@ +import numpy as np + +from groupy.garray.garray import GArray + + +class Z3Array(GArray): + parameterizations = ['int'] + _left_actions = {} + _reparameterizations = {} + _g_shapes = {'int': (3,)} + _group_name = 'Z3' + + def __init__(self, data, p='int'): + data = np.asarray(data) + assert data.dtype == np.int + self._left_actions[Z3Array] = self.__class__.z3_composition + super(Z3Array, self).__init__(data, p) + + def z3_composition(self, other): + return Z3Array(self.data + other.data) + + def inv(self): + return Z3Array(-self.data) + + def __repr__(self): + return 'Z3\n' + self.data.__repr__() + + def reparameterize(self, p): + assert p == 'int' + return self + + +def identity(shape=()): + ''' + Returns the identity element: an array of 3 zeros. + ''' + e = Z3Array(np.zeros(shape + (3,), dtype=np.int), 'int') + return e + + +def rand(minu, maxu, minv, maxv, minw, maxw, size=()): + ''' + Returns an Z3Array of shape size, with randomly chosen elements in int parameterization. + ''' + data = np.zeros(size + (3,), dtype=np.int64) + data[..., 0] = np.random.randint(minu, maxu, size) + data[..., 1] = np.random.randint(minv, maxv, size) + data[..., 2] = np.random.randint(minw, maxw, size) + return Z3Array(data=data, p='int') + + +def meshgrid(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2): + ''' + Creates a meshgrid of all elements of the group, within the given + translation parameters. + ''' + li = [[u, v, w] for u in xrange(minu, maxu) for v in xrange(minv, maxv) for + w in xrange(minw, maxw)] + return Z3Array(li, p='int') diff --git a/groupy/garray/__init__.py b/groupy/garray/__init__.py index 725516f..e9435db 100644 --- a/groupy/garray/__init__.py +++ b/groupy/garray/__init__.py @@ -1,6 +1,16 @@ - from groupy.garray.Z2_array import Z2Array +from groupy.garray.Z3_array import Z3Array from groupy.garray.p4_array import P4Array from groupy.garray.p4m_array import P4MArray from groupy.garray.C4_array import C4Array, C4Group from groupy.garray.D4_array import D4Array, D4Group +from groupy.garray.C4h_array import C4hArray +from groupy.garray.C4ht_array import C4htArray +from groupy.garray.D4h_array import D4hArray +from groupy.garray.D4ht_array import D4htArray +from groupy.garray.O_array import OArray +from groupy.garray.Ot_array import OtArray +from groupy.garray.Oh_array import OhArray +from groupy.garray.Oht_array import OhtArray + + diff --git a/groupy/garray/test_garray.py b/groupy/garray/test_garray.py index 820e558..be742ad 100644 --- a/groupy/garray/test_garray.py +++ b/groupy/garray/test_garray.py @@ -1,4 +1,3 @@ - # TODO: reshaping / flattening tests, check updating of shape, g_shape, ndim, g_ndim # TODO: test all left_actions, not just composition in group @@ -18,6 +17,11 @@ def test_z2_array(): check_wallpaper_group(Z2_array, Z2_array.Z2Array) +def test_z3_array(): + from groupy.garray import Z3_array + check_space_group(Z3_array, Z3_array.Z3Array) + + def test_c4_array(): from groupy.garray import C4_array check_finite_group(C4_array, C4_array.C4Array, C4_array.C4) @@ -28,8 +32,62 @@ def test_d4_array(): check_finite_group(D4_array, D4_array.D4Array, D4_array.D4) -def check_wallpaper_group(garray_module, garray_class): +def test_o_array(): + from groupy.garray import O_array + check_finite_group(O_array, O_array.OArray, O_array.O) + + +def test_c4h_array(): + from groupy.garray import C4h_array + check_finite_group(C4h_array, C4h_array.C4hArray, C4h_array.C4h) + + +def test_c4ht_array(): + from groupy.garray import C4ht_array + check_space_group(C4ht_array, C4ht_array.C4htArray) +# + +def test_d4h_array(): + from groupy.garray import D4h_array + check_finite_group(D4h_array, D4h_array.D4hArray, D4h_array.D4h) + + +def test_d4ht_array(): + from groupy.garray import D4ht_array + check_space_group(D4ht_array, D4ht_array.D4htArray) + + +def test_oh_array(): + from groupy.garray import Oh_array + check_finite_group(Oh_array, Oh_array.OhArray, Oh_array.Oh) + + +def test_ot_array(): + from groupy.garray import Ot_array + check_space_group(Ot_array, Ot_array.OtArray) + +def test_oht_array(): + from groupy.garray import Oht_array + check_space_group(Oht_array, Oht_array.OhtArray) + + +def check_space_group(garray_module, garray_class): + a = garray_module.rand(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2, size=(2, 3)) + b = garray_module.rand(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2, size=(2, 3)) + c = garray_module.rand(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2, size=(2, 3)) + + check_associative(a, b, c) + check_identity(garray_module, a) + check_inverse(garray_module, a) + + check_reparameterize_invertible(garray_class, a) + + m = garray_module.meshgrid(minu=-1, maxu=2, minv=-1, maxv=2, minw=-1, maxw=2) + check_closed_inverse(m) + + +def check_wallpaper_group(garray_module, garray_class): a = garray_module.rand(minu=-1, maxu=2, minv=-1, maxv=2, size=(2, 3)) b = garray_module.rand(minu=-1, maxu=2, minv=-1, maxv=2, size=(2, 3)) c = garray_module.rand(minu=-1, maxu=2, minv=-1, maxv=2, size=(2, 3)) @@ -48,7 +106,6 @@ def check_wallpaper_group(garray_module, garray_class): def check_finite_group(garray_module, garray_class, G): - a = garray_module.rand() b = garray_module.rand() c = garray_module.rand() diff --git a/groupy/gconv/make_gconv_indices.py b/groupy/gconv/make_gconv_indices.py index 570ebac..67b7a6d 100644 --- a/groupy/gconv/make_gconv_indices.py +++ b/groupy/gconv/make_gconv_indices.py @@ -1,18 +1,23 @@ - # Code for generating indices used in G-convolutions for various groups G. # The indices created by these functions are used to rotate and flip filters on the plane or on a group. # These indices depend only on the filter size, so they are created only once at the beginning of training. import numpy as np - from groupy.garray.C4_array import C4 from groupy.garray.D4_array import D4 +from groupy.garray.O_array import O +from groupy.garray.Oh_array import Oh +from groupy.garray.C4h_array import C4h +from groupy.garray.D4h_array import D4h from groupy.garray.p4_array import C4_halfshift -from groupy.gfunc.z2func_array import Z2FuncArray +from groupy.gfunc.otfunc_array import OtFuncArray +from groupy.gfunc.ohtfunc_array import OhtFuncArray +from groupy.gfunc.c4htfunc_array import C4htFuncArray +from groupy.gfunc.d4htfunc_array import D4htFuncArray from groupy.gfunc.p4func_array import P4FuncArray from groupy.gfunc.p4mfunc_array import P4MFuncArray - - +from groupy.gfunc.z2func_array import Z2FuncArray +from groupy.gfunc.z3func_array import Z3FuncArray def make_c4_z2_indices(ksize): x = np.random.randn(1, ksize, ksize) f = Z2FuncArray(v=x) @@ -20,7 +25,8 @@ def make_c4_z2_indices(ksize): if ksize % 2 == 0: uv = f.left_translation_indices(C4_halfshift[:, None, None, None]) else: - uv = f.left_translation_indices(C4[:, None, None, None]) + a = C4[:, None, None, None] + uv = f.left_translation_indices(a) r = np.zeros(uv.shape[:-1] + (1,)) ruv = np.c_[r, uv] return ruv.astype('int32') @@ -37,6 +43,82 @@ def make_c4_p4_indices(ksize): return li.astype('int32') +def make_d4h_z3_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(1, ksize, ksize, ksize) + f = Z3FuncArray(v=x) + a = D4h[:, None, None, None, None] + uvw = f.left_translation_indices(a) + i = np.zeros(uvw.shape[:-1] + (1,)) + iuvw = np.c_[i, uvw] + return iuvw.astype('int32') + + +def make_d4h_d4ht_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(16, ksize, ksize, ksize) + f = D4htFuncArray(v=x) + li = f.left_translation_indices(D4h[:, None, None, None, None]) + return li.astype('int32') + + +def make_c4h_z3_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(1, ksize, ksize, ksize) + f = Z3FuncArray(v=x) + a = C4h[:, None, None, None, None] + uvw = f.left_translation_indices(a) + i = np.zeros(uvw.shape[:-1] + (1,)) + iuvw = np.c_[i, uvw] + return iuvw.astype('int32') + + +def make_c4h_c4ht_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(8, ksize, ksize, ksize) + f = C4htFuncArray(v=x) + li = f.left_translation_indices(C4h[:, None, None, None, None]) + return li.astype('int32') + + +def make_o_z3_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(1, ksize, ksize, ksize) + f = Z3FuncArray(v=x) + a = O[:, None, None, None, None] + uvw = f.left_translation_indices(a) + i = np.zeros(uvw.shape[:-1] + (1,)) + iuvw = np.c_[i, uvw] + return iuvw.astype('int32') + + +def make_o_ot_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(24, ksize, ksize, ksize) + f = OtFuncArray(v=x) + li = f.left_translation_indices(O[:, None, None, None, None]) + return li.astype('int32') + + +def make_oh_z3_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(1, ksize, ksize, ksize) + f = Z3FuncArray(v=x) + a = Oh[:, None, None, None, None] + uvw = f.left_translation_indices(a) + i = np.zeros(uvw.shape[:-1] + (1,)) + iuvw = np.c_[i, uvw] + return iuvw.astype('int32') + + +def make_oh_oht_indices(ksize): + assert ksize % 2 == 1 + x = np.random.randn(48, ksize, ksize, ksize) + f = OhtFuncArray(v=x) + li = f.left_translation_indices(Oh[:, None, None, None, None]) + return li.astype('int32') + + def make_d4_z2_indices(ksize): assert ksize % 2 == 1 # TODO x = np.random.randn(1, ksize, ksize) @@ -60,10 +142,8 @@ def flatten_indices(inds): The Chainer implementation of G-Conv uses indices into a 5D filter tensor (with an additional axis for the transformations H. For the tensorflow implementation it was more convenient to flatten the filter tensor into a 3D tensor with shape (output channels, input channels, transformations * width * height). - This function takes indices in the format required for Chainer and turns them into indices into the flat array used by tensorflow. - :param inds: np.ndarray of shape (output transformations, input transformations, n, n, 3), as output by the functions like make_d4_p4m_indices(n). :return: np.ndarray of shape (output transformations, input transformations, n, n) @@ -75,4 +155,26 @@ def flatten_indices(inds): V = inds[..., 2] # shape (nto, nti, n, n) # inds_flat = T * n * n + U * n + V inds_flat = U * n * nti + V * nti + T - return inds_flat \ No newline at end of file + + return inds_flat + + +def flatten_indices_3d(inds): + """ + The Chainer implementation of G-Conv uses indices into a 5D filter tensor (with an additional axis for the + transformations H. For the tensorflow implementation it was more convenient to flatten the filter tensor into + a 3D tensor with shape (output channels, input channels, transformations * width * height). + This function takes indices in the format required for Chainer and turns them into indices into the flat array + used by tensorflow. + :param inds: np.ndarray of shape (output transformations, input transformations, n, n, 3), as output by + the functions like make_d4_p4m_indices(n). + :return: np.ndarray of shape (output transformations, input transformations, n, n) + """ + n = inds.shape[-2] + nti = inds.shape[1] + T = inds[..., 0] # shape (nto, nti, n, n, n) + U = inds[..., 1] # shape (nto, nti, n, n, n) + V = inds[..., 2] # shape (nto, nti, n, n, n) + W = inds[..., 3] # shape (nto, nti, n, n, n) + inds_flat = U * n * n * nti + V * n * nti + W * nti + T + return inds_flat diff --git a/groupy/gconv/tensorflow_gconv/check_gconv3d.py b/groupy/gconv/tensorflow_gconv/check_gconv3d.py new file mode 100644 index 0000000..d07accb --- /dev/null +++ b/groupy/gconv/tensorflow_gconv/check_gconv3d.py @@ -0,0 +1,104 @@ +import groupy.garray.O_array as O +import groupy.garray.Oh_array as Oh +import groupy.garray.C4h_array as C4h +import groupy.garray.D4h_array as D4h +import numpy as np +import tensorflow as tf +from groupy.gconv.tensorflow_gconv.splitgconv3d import gconv3d_util, gconv3d +from groupy.gfunc.otfunc_array import OtFuncArray +from groupy.gfunc.c4htfunc_array import C4htFuncArray +from groupy.gfunc.d4htfunc_array import D4htFuncArray +from groupy.gfunc.ohtfunc_array import OhtFuncArray +from groupy.gfunc.z3func_array import Z3FuncArray + + +def check_o_z3_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 1) + x, y = make_graph('Z3', 'O', ksize) + check_equivariance(im, x, y, Z3FuncArray, OtFuncArray, O) + + +def check_o_o_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 24) + x, y = make_graph('O', 'O', ksize) + check_equivariance(im, x, y, OtFuncArray, OtFuncArray, O) + + +def check_oh_z3_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 1) + x, y = make_graph('Z3', 'OH', ksize) + check_equivariance(im, x, y, Z3FuncArray, OhtFuncArray, Oh) + + +def check_oh_oh_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 48) + x, y = make_graph('OH', 'OH', ksize) + check_equivariance(im, x, y, OhtFuncArray, OhtFuncArray, Oh) + + +def check_c4h_z3_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 1) + x, y = make_graph('Z3', 'C4H', ksize) + check_equivariance(im, x, y, Z3FuncArray, C4htFuncArray, C4h) + + +def check_c4h_c4h_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 8) + x, y = make_graph('C4H', 'C4H', ksize) + check_equivariance(im, x, y, C4htFuncArray, C4htFuncArray, C4h) + + +def check_d4h_z3_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 1) + x, y = make_graph('Z3', 'D4H', ksize) + check_equivariance(im, x, y, Z3FuncArray, D4htFuncArray, D4h) + + +def check_d4h_d4h_conv_equivariance(): + ksize = 3 + im = np.random.randn(2, ksize, ksize, ksize, 16) + x, y = make_graph('D4H', 'D4H', ksize) + check_equivariance(im, x, y, Z3FuncArray, D4htFuncArray, D4h) + + +def make_graph(h_input, h_output, ksize): + gconv_indices, gconv_shape_info, w_shape = gconv3d_util( + h_input=h_input, h_output=h_output, in_channels=1, out_channels=1, ksize=ksize) + nti = gconv_shape_info[-2] + + x = tf.placeholder(tf.float32, [None, ksize, ksize, ksize, 1 * nti]) + w = tf.Variable(tf.truncated_normal(w_shape, stddev=1.)) + + y = gconv3d(input=x, filter=w, strides=[1, 1, 1, 1, 1], padding='SAME', + gconv_indices=gconv_indices, gconv_shape_info=gconv_shape_info) + return x, y + + +def check_equivariance(im, input, output, input_array, output_array, point_group): + # Transform the image + f = input_array(im.transpose((0, 4, 1, 2, 3))) + g = point_group.rand() + gf = g * f + im1 = gf.v.transpose((0, 2, 3, 4, 1)) + + # Compute + init = tf.global_variables_initializer() + sess = tf.Session() + sess.run(init) + yx = sess.run(output, feed_dict={input: im}) + yrx = sess.run(output, feed_dict={input: im1}) + sess.close() + + # Transform the computed feature maps + fmap1_garray = output_array(yrx.transpose((0, 4, 1, 2, 3))) + r_fmap1_data = (g.inv() * fmap1_garray).v.transpose((0, 2, 3, 4, 1)) + + assert np.allclose(yx, r_fmap1_data, rtol=1e-5, atol=1e-3) + diff --git a/groupy/gconv/tensorflow_gconv/splitgconv3d.py b/groupy/gconv/tensorflow_gconv/splitgconv3d.py new file mode 100644 index 0000000..f78dc5f --- /dev/null +++ b/groupy/gconv/tensorflow_gconv/splitgconv3d.py @@ -0,0 +1,112 @@ +import tensorflow as tf + +from groupy.gconv.make_gconv_indices import make_o_z3_indices, make_o_ot_indices, make_c4h_z3_indices, \ + make_c4h_c4ht_indices, make_d4h_z3_indices, make_d4h_d4ht_indices, make_oh_z3_indices, make_oh_oht_indices, \ + flatten_indices_3d +from groupy.gconv.tensorflow_gconv.transform_filter import transform_filter_3d_nhwc + + +def gconv3d(input, filter, strides, padding, gconv_indices, gconv_shape_info, + use_cudnn_on_gpu=None, data_format='NHWC', name=None): + """ + Implements the g-convolution. Similar interface as the standard 3D convolution in tensorflow, with gconv_indices + and gconv_shape_info as additional parameters. These can be obtained using gconv3d_util. + Args: + input: tensor with (b, z, y, x, c) axes + filter: tensor with (ksize, ksize, ksize, in channels * transformations, out_channels) axes + strides: A list of ints. 1-D of length 5. The stride of the sliding window for each dimension of input. + padding: A string from: "SAME", "VALID". The type of padding algorithm to use. + gconv_indices: indices used in the filter transformation step of the G-Conv. + gconv_shape_info: a tuple containing information for the gconv: in/out channels, in/out transformations, ksize + data_format: only nhwc supported + name: a name for the operation + + Returns: + conv: tensor with (batch, z, y, x, c) axes + + """ + if data_format != 'NHWC': + raise NotImplemented('Currently only NHWC data_format is supported. Got:' + str(data_format)) + + # Transform the filters + transformed_filter = transform_filter_3d_nhwc(w=filter, flat_indices=gconv_indices, shape_info=gconv_shape_info) + + # Convolve input with transformed filters + conv = tf.nn.conv3d(input=input, filter=transformed_filter, strides=strides, padding=padding, name=name) + + return conv + + +def gconv3d_util(h_input, h_output, in_channels, out_channels, ksize): + """ + Convenience function for setting up static data required for the G-Conv. The number of 3D channels will be + 1, 8, 16, 24 or 48 times larger depending on the value of h_input and h_output. + + Args: + h_input: Z3, C4H, D4H, O, OH -- use one. Z3 for first layer. + h_output: Z3, C4H, D4H, O, OH -- use one. + in_channels: number of input channels of the 3D channels on the group. + out_channels: number of output channels of the 3D channels on the group. + ksize: the spatial size of filter kernels, typicall 3, 5 or 7. Only uneven ksize is supported. + + Returns: + gconv_indices: an array of indices used in the filter transformation step of gconv3d + w_shape: the shape of the filter tensor to be allocated and passed to gconv3d + gconv_shape_info: shape information required by gconv3d + -- (nr. out channels, nr. out transformations, nr. in channels, nr. in tranformations, ksize) + """ + + # uppercase for consistency + h_input = h_input.upper() + h_output = h_output.upper() + + # get number of transformations in and out + mapping = {'Z3': 1, 'C4': 4, 'D4': 8, 'O': 24, 'C4H': 8, 'D4H': 16, 'OH': 48} + nti = mapping[h_input] + nto = mapping[h_output] + + # get gconv_indices + if h_input == 'Z3' and h_output == 'O': + gconv_indices = make_o_z3_indices(ksize=ksize) + elif h_input == 'O' and h_output == 'O': + gconv_indices = make_o_ot_indices(ksize=ksize) + elif h_input == 'Z3' and h_output == 'C4H': + gconv_indices = make_c4h_z3_indices(ksize=ksize) + elif h_input == 'C4H' and h_output == 'C4H': + gconv_indices = make_c4h_c4ht_indices(ksize=ksize) + elif h_input == 'Z3' and h_output == 'D4H': + gconv_indices = make_d4h_z3_indices(ksize=ksize) + elif h_input == 'D4H' and h_output == 'D4H': + gconv_indices = make_d4h_d4ht_indices(ksize=ksize) + elif h_input == 'Z3' and h_output == 'OH': + gconv_indices = make_oh_z3_indices(ksize=ksize) + elif h_input == 'OH' and h_output == 'OH': + gconv_indices = make_oh_oht_indices(ksize=ksize) + else: + raise ValueError('Unknown (h_input, h_output) pair:' + str((h_input, h_output))) + + # flatten and get shape information and filter tensor shape + gconv_indices = flatten_indices_3d(gconv_indices) + w_shape = (ksize, ksize, ksize, in_channels * nti, out_channels) + gconv_shape_info = (out_channels, nto, in_channels, nti, ksize) + + return gconv_indices, gconv_shape_info, w_shape + + +def gconv2d_addbias(input, bias, nti=8): + """ + In a G-CNN, the feature maps are interpreted as functions on a group G instead of functions on the plane Z^2. + Just like how we use a single scalar bias per 2D feature map, in a G-CNN we should use a single scalar bias per + G-feature map. Failing to do this breaks the equivariance and typically hurts performance. + A G-feature map usually consists of a number (e.g. 4 or 8) adjacent channels. + This function will add a single bias vector to a stack of feature maps that has e.g. 4 or 8 times more 2D channels + than G-channels, by replicating the bias across adjacent groups of 2D channels. + + :param input: tensor of shape (n, h, w, ni * nti), where n is the batch dimension, (h, w) are the height and width, + ni is the number of input G-channels, and nti is the number of transformations in H. + :param bias: tensor of shape (ni,) + :param nti: number of transformations, e.g. 4 for C4/p4 or 8 for D4/p4m. + :return: input with bias added + """ + # input = tf.reshape(input, ()) + pass # TODO diff --git a/groupy/gconv/tensorflow_gconv/transform_filter.py b/groupy/gconv/tensorflow_gconv/transform_filter.py index 9f94697..a09b5b0 100644 --- a/groupy/gconv/tensorflow_gconv/transform_filter.py +++ b/groupy/gconv/tensorflow_gconv/transform_filter.py @@ -1,7 +1,22 @@ - import tensorflow as tf +def transform_filter_3d_nhwc(w, flat_indices, shape_info, validate_indices=True): + no, nto, ni, nti, n = shape_info + w_flat = tf.reshape(w, [n * n * n * nti, ni, no]) # shape (n * n * n * nti, ni, no) + + # Do the transformation / indexing operation. + transformed_w = tf.gather(w_flat, flat_indices, + validate_indices=validate_indices) # shape (nto, nti, n, n, n, ni, no) + + # Put the axes in the right order, and collapse them to get a standard shape filter bank + transformed_w = tf.transpose(transformed_w, [2, 3, 4, 5, 1, 6, 0]) # shape (n, n, n, ni, nti, no, nto) + + transformed_w = tf.reshape(transformed_w, [n, n, n, ni * nti, no * nto]) # shape (n, n, n, ni * nti, no * nto) + + return transformed_w + + def transform_filter_2d_nhwc(w, flat_indices, shape_info, validate_indices=True): """ Transform a set of filters defined on a split plane group G. @@ -27,15 +42,15 @@ def transform_filter_2d_nhwc(w, flat_indices, shape_info, validate_indices=True) # The indexing is done using tf.gather. This function can only do integer indexing along the first axis. # We want to index the spatial and transformation axes of our filter, so we must flatten them into one axis. no, nto, ni, nti, n = shape_info - w_flat = tf.reshape(w, [n * n * nti, ni, no]) # shape (n * n * nti, ni, no) + w_flat = tf.reshape(w, [n * n * nti, ni, no]) # shape (n * n * nti, ni, no) # Do the transformation / indexing operation. transformed_w = tf.gather(w_flat, flat_indices, - validate_indices=validate_indices) # shape (nto, nti, n, n, ni, no) + validate_indices=validate_indices) # shape (nto, nti, n, n, ni, no) # Put the axes in the right order, and collapse them to get a standard shape filter bank - transformed_w = tf.transpose(transformed_w, [2, 3, 4, 1, 5, 0]) # shape (n, n, ni, nti, no, nto) - transformed_w = tf.reshape(transformed_w, [n, n, ni * nti, no * nto]) # shape (n, n, ni * nti, no * nto) + transformed_w = tf.transpose(transformed_w, [2, 3, 4, 1, 5, 0]) # shape (n, n, ni, nti, no, nto) + transformed_w = tf.reshape(transformed_w, [n, n, ni * nti, no * nto]) # shape (n, n, ni * nti, no * nto) return transformed_w @@ -66,14 +81,14 @@ def transform_filter_2d_nchw(w, flat_indices, shape_info, validate_indices=True) # We want to index the spatial and transformation axes of our filter, so we must flatten them into one axis, # and bring them to the first axis no, nto, ni, nti, n = shape_info - w_flat = tf.transpose(tf.reshape(w, [no, ni, nti * n * n]), [2, 0, 1]) # shape (nti * n * n, no, ni) + w_flat = tf.transpose(tf.reshape(w, [no, ni, nti * n * n]), [2, 0, 1]) # shape (nti * n * n, no, ni) # Do the transformation / indexing operation. transformed_w = tf.gather(w_flat, flat_indices, - validate_indices=validate_indices) # shape (nto, nti, n, n, no, ni) + validate_indices=validate_indices) # shape (nto, nti, n, n, no, ni) # Put the axes in the right order, and collapse them to get a standard-shape filter bank - transformed_w = tf.transpose(transformed_w, [4, 0, 5, 1, 2, 3]) # shape (no, nto, ni, nti, n, n) - transformed_w = tf.reshape(transformed_w, (no * nto, ni * nti, n, n)) # shape (no * nto, ni * nti, n, n) + transformed_w = tf.transpose(transformed_w, [4, 0, 5, 1, 2, 3]) # shape (no, nto, ni, nti, n, n) + transformed_w = tf.reshape(transformed_w, (no * nto, ni * nti, n, n)) # shape (no * nto, ni * nti, n, n) return transformed_w diff --git a/groupy/gfunc/__init__.py b/groupy/gfunc/__init__.py index b552da0..460374e 100644 --- a/groupy/gfunc/__init__.py +++ b/groupy/gfunc/__init__.py @@ -1,4 +1,8 @@ - -from groupy.gfunc.p4func_array import P4FuncArray -from groupy.gfunc.p4mfunc_array import P4MFuncArray -from groupy.gfunc.z2func_array import Z2FuncArray +from ohtfunc_array import OhtFuncArray +from otfunc_array import OtFuncArray +from c4htfunc_array import C4htFuncArray +from d4htfunc_array import D4htFuncArray +from p4func_array import P4FuncArray +from p4mfunc_array import P4MFuncArray +from z2func_array import Z2FuncArray +from z3func_array import Z3FuncArray diff --git a/groupy/gfunc/c4htfunc_array.py b/groupy/gfunc/c4htfunc_array.py new file mode 100644 index 0000000..9fa3b88 --- /dev/null +++ b/groupy/gfunc/c4htfunc_array.py @@ -0,0 +1,57 @@ +import groupy.garray.C4ht_array as c4ht +from groupy.gfunc.gfuncarray import GFuncArray + + +class C4htFuncArray(GFuncArray): + def __init__(self, v, umin=None, umax=None, vmin=None, vmax=None, wmin=None, wmax=None): + + # TODO: error message + if umin is None or umax is None or vmin is None or vmax is None: + if not (umin is None and umax is None and vmin is None and vmax is None): + raise ValueError('Either all or none of umin, umax, vmin, vmax must equal None') + + # If (u, v, w) ranges are not given, determine them from the shape of v, + # assuming the grid is centered. + nu, nv, nw = v.shape[-3:] + + hnu = nu // 2 + hnv = nv // 2 + hnw = nw // 2 + + umin = -hnu + umax = hnu - (nu % 2 == 0) + vmin = -hnv + vmax = hnv - (nv % 2 == 0) + wmin = -hnw + wmax = hnw - (nw % 2 == 0) + + self.umin = umin + self.umax = umax + 1 + self.vmin = vmin + self.vmax = vmax + 1 + self.wmin = wmin + self.wmax = wmax + 1 + + i2g = c4ht.meshgrid( + minu=self.umin, + maxu=self.umax, + minv=self.vmin, + maxv=self.vmax, + minw=self.wmin, + maxw=self.wmax + ) + i2g = i2g.reshape(v.shape[-4:]) + + super(C4htFuncArray, self).__init__(v=v, i2g=i2g) + + def g2i(self, g): + gint = g.reparameterize('int').data.copy() + gint[..., 2] -= self.umin + gint[..., 3] -= self.vmin + gint[..., 4] -= self.wmin + + # flat stabilizer: instead of (4, 2, ...) use (8, ...) + gint[..., 1] += gint[..., 0] * 4 + gint = gint[..., 1:] + + return gint diff --git a/groupy/gfunc/d4htfunc_array.py b/groupy/gfunc/d4htfunc_array.py new file mode 100644 index 0000000..8211a74 --- /dev/null +++ b/groupy/gfunc/d4htfunc_array.py @@ -0,0 +1,59 @@ +import groupy.garray.D4ht_array as d4ht +from groupy.gfunc.gfuncarray import GFuncArray + + +class D4htFuncArray(GFuncArray): + def __init__(self, v, umin=None, umax=None, vmin=None, vmax=None, wmin=None, wmax=None): + + # TODO: error message + if umin is None or umax is None or vmin is None or vmax is None: + if not (umin is None and umax is None and vmin is None and vmax is None): + raise ValueError('Either all or none of umin, umax, vmin, vmax must equal None') + + # If (u, v, w) ranges are not given, determine them from the shape of v, + # assuming the grid is centered. + nu, nv, nw = v.shape[-3:] + + hnu = nu // 2 + hnv = nv // 2 + hnw = nw // 2 + + umin = -hnu + umax = hnu - (nu % 2 == 0) + vmin = -hnv + vmax = hnv - (nv % 2 == 0) + wmin = -hnw + wmax = hnw - (nw % 2 == 0) + + self.umin = umin + self.umax = umax + 1 + self.vmin = vmin + self.vmax = vmax + 1 + self.wmin = wmin + self.wmax = wmax + 1 + + i2g = d4ht.meshgrid( + minu=self.umin, + maxu=self.umax, + minv=self.vmin, + maxv=self.vmax, + minw=self.wmin, + maxw=self.wmax + ) + i2g = i2g.reshape(v.shape[-4:]) + + super(D4htFuncArray, self).__init__(v=v, i2g=i2g) + + def g2i(self, g): + gint = g.reparameterize('int').data.copy() + gint[..., 3] -= self.umin + gint[..., 4] -= self.vmin + gint[..., 5] -= self.wmin + + # flat stabilizer: instead of (2, 4, 2, ...) use (16, ...) + gint[..., 1] += gint[..., 0] * 4 + gint[..., 2] += gint[..., 1] * 2 + + gint = gint[..., 2:] + + return gint diff --git a/groupy/gfunc/gfuncarray.py b/groupy/gfunc/gfuncarray.py index b9c9489..232feb7 100644 --- a/groupy/gfunc/gfuncarray.py +++ b/groupy/gfunc/gfuncarray.py @@ -9,8 +9,8 @@ class GFuncArray(object): def __init__(self, v, i2g): """ A GFunc is a discretely sampled function on a group or homogeneous space G. - The GFuncArray stores an array of GFuncs, - together with a map from G to an index set I (the set of sampling points) and the inverse of this map. + The GFuncArray stores an array of GFuncs,together with a map from G to an + index set I (the set of sampling points) and the inverse of this map. The ndarray v can be thought of as a map v : J x I -> R @@ -39,8 +39,8 @@ def __init__(self, v, i2g): So v implicitly defines a function v' on G: v'(g) = v(g2i(g)) - If we have a map T: G - > G (e.g. left multiplication by g^-1), that we want to precompose with v', - w'(g) = v'(T(g)) + If we have a map T: G - > G (e.g. left multiplication by g^-1), that we want to + precompose with v', w'(g) = v'(T(g)) we can get the corresponding map v by composing maps like this: I ---> G ---> G ---> I ---> R diff --git a/groupy/gfunc/ohtfunc_array.py b/groupy/gfunc/ohtfunc_array.py new file mode 100644 index 0000000..644ba8a --- /dev/null +++ b/groupy/gfunc/ohtfunc_array.py @@ -0,0 +1,57 @@ +import groupy.garray.Oht_array as oht +from groupy.gfunc.gfuncarray import GFuncArray + + +class OhtFuncArray(GFuncArray): + def __init__(self, v, umin=None, umax=None, vmin=None, vmax=None, wmin=None, wmax=None): + + # TODO: error message + if umin is None or umax is None or vmin is None or vmax is None: + if not (umin is None and umax is None and vmin is None and vmax is None): + raise ValueError('Either all or none of umin, umax, vmin, vmax must equal None') + + # If (u, v, w) ranges are not given, determine them from the shape of v, + # assuming the grid is centered. + nu, nv, nw = v.shape[-3:] + + hnu = nu // 2 + hnv = nv // 2 + hnw = nw // 2 + + umin = -hnu + umax = hnu - (nu % 2 == 0) + vmin = -hnv + vmax = hnv - (nv % 2 == 0) + wmin = -hnw + wmax = hnw - (nw % 2 == 0) + + self.umin = umin + self.umax = umax + 1 + self.vmin = vmin + self.vmax = vmax + 1 + self.wmin = wmin + self.wmax = wmax + 1 + + i2g = oht.meshgrid( + minu=self.umin, + maxu=self.umax, + minv=self.vmin, + maxv=self.vmax, + minw=self.wmin, + maxw=self.wmax + ) + i2g = i2g.reshape(v.shape[-4:]) + + super(OhtFuncArray, self).__init__(v=v, i2g=i2g) + + def g2i(self, g): + gint = g.reparameterize('int').data.copy() + + gint[..., 2] -= self.umin + gint[..., 3] -= self.vmin + gint[..., 4] -= self.wmin + + # flat stabilizer: instead of (24, 2, ...) use (48, ...) + gint[..., 1] += gint[..., 0] * 2 + gint = gint[..., 1:] + return gint diff --git a/groupy/gfunc/otfunc_array.py b/groupy/gfunc/otfunc_array.py new file mode 100644 index 0000000..16bb37c --- /dev/null +++ b/groupy/gfunc/otfunc_array.py @@ -0,0 +1,53 @@ +import groupy.garray.Ot_array as ot +from groupy.gfunc.gfuncarray import GFuncArray + + +class OtFuncArray(GFuncArray): + def __init__(self, v, umin=None, umax=None, vmin=None, vmax=None, wmin=None, wmax=None): + + # TODO: error message + if umin is None or umax is None or vmin is None or vmax is None: + if not (umin is None and umax is None and vmin is None and vmax is None): + raise ValueError('Either all or none of umin, umax, vmin, vmax must equal None') + + # If (u, v, w) ranges are not given, determine them from the shape of v, + # assuming the grid is centered. + nu, nv, nw = v.shape[-3:] + + hnu = nu // 2 + hnv = nv // 2 + hnw = nw // 2 + + umin = -hnu + umax = hnu - (nu % 2 == 0) + vmin = -hnv + vmax = hnv - (nv % 2 == 0) + wmin = -hnw + wmax = hnw - (nw % 2 == 0) + + self.umin = umin + self.umax = umax + 1 + self.vmin = vmin + self.vmax = vmax + 1 + self.wmin = wmin + self.wmax = wmax + 1 + + i2g = ot.meshgrid( + minu=self.umin, + maxu=self.umax, + minv=self.vmin, + maxv=self.vmax, + minw=self.wmin, + maxw=self.wmax + ) + + i2g = i2g.reshape(v.shape[-4:]) + + super(OtFuncArray, self).__init__(v=v, i2g=i2g) + + def g2i(self, g): + gint = g.reparameterize('int').data.copy() + gint[..., 1] -= self.umin + gint[..., 2] -= self.vmin + gint[..., 3] -= self.wmin + return gint diff --git a/groupy/gfunc/p4func_array.py b/groupy/gfunc/p4func_array.py index da02753..d70a3b1 100644 --- a/groupy/gfunc/p4func_array.py +++ b/groupy/gfunc/p4func_array.py @@ -1,4 +1,3 @@ - import numpy as np import groupy.garray.p4_array as p4a from groupy.gfunc.gfuncarray import GFuncArray @@ -49,7 +48,7 @@ def g2i(self, g): def tst(): - from groupy.garray.p4_array import P4Array, meshgrid, u_range, v_range, rotation, translation + from groupy.garray.p4_array import meshgrid, u_range, v_range, rotation x = np.random.randn(4, 3, 3) c = meshgrid(u=u_range(-1, 2), v=v_range(-1, 2)) @@ -65,4 +64,5 @@ def tst(): gf = g * f gfi = f.v[li[..., 0], li[..., 1], li[..., 2]] - return x, c, f, li, gf, gfp, gfi \ No newline at end of file + return x, c, f, li, gf, gfp, gfi + diff --git a/groupy/gfunc/test_gfuncarray.py b/groupy/gfunc/test_gfuncarray.py index 62687d5..f09af49 100644 --- a/groupy/gfunc/test_gfuncarray.py +++ b/groupy/gfunc/test_gfuncarray.py @@ -1,6 +1,68 @@ import numpy as np +def test_ot_func(): + from groupy.gfunc.otfunc_array import OtFuncArray + import groupy.garray.O_array as o + + v = np.random.randn(2, 6, 24, 5, 5, 5) + f = OtFuncArray(v=v) + + g = o.rand(size=(1,)) + h = o.rand(size=(1,)) + + check_associative(g, h, f) + check_identity(o, f) + check_invertible(g, f) + check_i2g_g2i_invertible(f) + + +def test_oht_func(): + from groupy.gfunc.ohtfunc_array import OhtFuncArray + import groupy.garray.Oh_array as oh + + v = np.random.randn(2, 6, 48, 1, 1, 1) + f = OhtFuncArray(v=v) + + g = oh.rand(size=(1,)) + h = oh.rand(size=(1,)) + + check_associative(g, h, f) + check_identity(oh, f) + check_invertible(g, f) + check_i2g_g2i_invertible(f) + +def test_c4ht_func(): + from groupy.gfunc.c4htfunc_array import C4htFuncArray + import groupy.garray.C4h_array as c4h + + v = np.random.randn(2, 6, 8, 3, 3, 3) + f = C4htFuncArray(v=v) + + g = c4h.rand(size=(1,)) + h = c4h.rand(size=(1,)) + + check_associative(g, h, f) + check_identity(c4h, f) + check_invertible(g, f) + check_i2g_g2i_invertible(f) + +def test_d4ht_func(): + from groupy.gfunc.d4htfunc_array import D4htFuncArray + import groupy.garray.D4h_array as d4h + + v = np.random.randn(2, 6, 16, 3, 3, 3) + f = D4htFuncArray(v=v) + + g = d4h.rand(size=(1,)) + h = d4h.rand(size=(1,)) + + check_associative(g, h, f) + check_identity(d4h, f) + check_invertible(g, f) + check_i2g_g2i_invertible(f) + + def test_p4_func(): from groupy.gfunc.p4func_array import P4FuncArray import groupy.garray.C4_array as c4a @@ -51,7 +113,33 @@ def test_z2_func(): g = d4a.rand(size=(1,)) h = d4a.rand(size=(1,)) check_associative(g, h, f) - check_identity(c4a, f) + check_identity(d4a, f) + check_invertible(g, f) + check_i2g_g2i_invertible(f) + + +def test_z3_func(): + from groupy.gfunc.z3func_array import Z3FuncArray + import groupy.garray.O_array as o + import groupy.garray.Oh_array as oh + + v = np.random.randn(2, 6, 5, 5, 5) + f = Z3FuncArray(v=v) + + g = o.rand(size=(1,)) + h = o.rand(size=(1,)) + check_associative(g, h, f) + check_identity(o, f) + check_invertible(g, f) + check_i2g_g2i_invertible(f) + + v = np.random.randn(2, 6, 5, 5, 5) + f = Z3FuncArray(v=v) + + g = oh.rand(size=(1,)) + h = oh.rand(size=(1,)) + check_associative(g, h, f) + check_identity(oh, f) check_invertible(g, f) check_i2g_g2i_invertible(f) @@ -78,8 +166,3 @@ def check_i2g_g2i_invertible(f): i = f.g2i(i2g) inds = [i[..., j] for j in range(i.shape[-1])] assert (i2g[inds] == i2g).all() - - - - - diff --git a/groupy/gfunc/z3func_array.py b/groupy/gfunc/z3func_array.py new file mode 100644 index 0000000..a091ef9 --- /dev/null +++ b/groupy/gfunc/z3func_array.py @@ -0,0 +1,54 @@ +import groupy.garray.Z3_array as z3a +from groupy.gfunc.gfuncarray import GFuncArray + + +class Z3FuncArray(GFuncArray): + def __init__(self, v, umin=None, umax=None, vmin=None, vmax=None): + + if umin is None or umax is None or vmin is None or vmax is None: + if not (umin is None and umax is None and vmin is None and vmax is None): + raise ValueError('Either all or none of umin, umax, vmin, vmax must equal None') + + # If (u, v) ranges are not given, determine them from the shape of v, + # assuming the grid is centered. + nu, nv, nw = v.shape[-3:] + + hnu = nu // 2 + hnv = nv // 2 + hnw = nw // 2 + + umin = -hnu + umax = hnu - (nu % 2 == 0) + vmin = -hnv + vmax = hnv - (nv % 2 == 0) + wmin = -hnw + wmax = hnw - (nw % 2 == 0) + + self.umin = umin + self.umax = umax + 1 + self.vmin = vmin + self.vmax = vmax + 1 + self.wmin = wmin + self.wmax = wmax + 1 + + i2g = z3a.meshgrid( + minu=self.umin, + maxu=self.umax, + minv=self.vmin, + maxv=self.vmax, + minw=self.wmin, + maxw=self.wmax + ) + + i2g = i2g.reshape(v.shape[-3:]) + super(Z3FuncArray, self).__init__(v=v, i2g=i2g) + + def g2i(self, g): + # TODO: check validity of indices and wrap / clamp if necessary + # (or do this in a separate function, so that this function can be more easily tested?) + + gint = g.reparameterize('int').data.copy() + gint[..., 0] -= self.umin + gint[..., 1] -= self.vmin + gint[..., 2] -= self.wmin + return gint