"""
Part of comet/pyhed/misc
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import numpy as np
import scipy
import time
import pygimli as pg
[docs]def fixSingularity(model, drop_tol=1e-2, dipole_z=0.0):
"""
Points in zero get values of drop_tol.
Points on drop_tol get Values of up to 20% the value to the first point out
of the drop_tol.
If all points are in the drop_tol a warning is printed.
"""
from comet.pyhed import log
if drop_tol is None or np.isclose(drop_tol, 0., atol=1e-15):
# no drop_tol
return model
if drop_tol < 0.0:
drop_tol *= -1
index = model[0] <= drop_tol
if np.alltrue(np.logical_not(index)):
# all points outside drop_tol
return model
factor = model[0, index] / drop_tol
index2 = np.where(model[0] > drop_tol)
if len(index2[0]) == 0:
log.warning('No Values outside the drop_tol. Shift points within '
f'{drop_tol} m distance away from dipole.')
Min = drop_tol * 2
else:
Min = np.min(model[0, index2])
twentyPerc = 0.2 * (Min - drop_tol)
model[0, index] = (factor * twentyPerc) + drop_tol
# fixed by drop_tol
return model
[docs]def Ca2Cy(cartesian, dtype=np.float64, drop_tol=1e-2, dipole_z=0.0):
""" Convert cartesian coords to cylindrical coords.
Parameters
----------
cartesian: np.ndarray
Coordinate vector of for N positions of shape (3, N).
dtype: np.dtype or str
Optional choice of output data type.
drop_tol:
Tolerance in m to avoid zeros in horizontal distances (singlularity
removal). All values for the resulting horizontal distance below
the drop_tol are redistributed between drop_tol and 20% of the
distance to the first point outside droptol. Raises Exception if
no point lies outside of drop_tol. drop_tol=None disables the
singularity removal (default).
Returns
-------
cylindrical: np.ndarray
Coordinate Vector in cylindrical cordinates (radius, phi, z) in
given datra type or input data type and with or without
singularities removed.
"""
cylindrical = np.array((np.sqrt(cartesian[0]**2 + cartesian[1]**2), # r
np.arctan2(cartesian[1], cartesian[0]), # phi
cartesian[2]), # z
dtype=cartesian.dtype if dtype is None else dtype)
cylindrical = fixSingularity(cylindrical, drop_tol=drop_tol,
dipole_z=dipole_z)
return cylindrical
[docs]def KtoP(cartesian, dtype=np.float64, drop_tol=1e-2):
'''
# status: depricated, use Ca2Cy instead
'''
return Ca2Cy(cartesian, dtype=dtype, drop_tol=drop_tol)
[docs]def KtoP_all(cartesian, dtype=np.float64, drop_tol=1e-2):
'''
# r = np.sqrt(model[0]**2 + model[1]**2)
# phi = np.arctan2(model[1], model[0])
# z = np.copy(model[2])
'''
polar = np.array((np.sqrt(cartesian[:, 0]**2 + cartesian[:, 1]**2), # r
np.arctan2(cartesian[:, 1], cartesian[:, 0]), # phi
cartesian[:, 2]), # z
dtype=cartesian.dtype if dtype is None else dtype)
polar[0] = fixSingularity(polar[0], drop_tol=drop_tol)
return polar
[docs]def Cy2Ca(cylindrical, dtype=None):
""" Convert cartesian coords to cylindrical coords.
Parameters
----------
cylindrical: np.ndarray
Coordinate vector of for N positions of shape (3, N).
dtype: np.dtype or str
Optional choice of output data type.
Returns
-------
cartesian: np.ndarray
Coordinate Vector in cartesian cordinates (x, y, z) in
given datra type or input data type.
"""
cartesian = np.array((cylindrical[0] * np.cos(cylindrical[1]), # x
cylindrical[0] * np.sin(cylindrical[1]), # y
cylindrical[2]), # z
dtype=cylindrical.dtype if dtype is None else dtype)
return cartesian
[docs]def PtoK(cylindrical, dtype=np.float64):
'''
# status: depricated use Cy2Ca instead
'''
return Cy2Ca(cylindrical, dtype=dtype)
[docs]def Ca2CyField(cartesian, field, dtype=None):
''' Conversion of 3d vector field from cylindrical coords in cartesian.
# cartesian = cartesian coordinate system (x, y, z)
# x = field[x] * cos(phi) - field[y] * sin(phi)
# y = - field[x] * sin(Phi) + field[y] * cos(phi)
# z = field[z]
# output:
# field_cylindrical = data in cylindrical coordinates (r, phi, z)
Parameters
----------
cartesian: np.ndarray
Coordinate vector of for N positions of shape (3, N).
field: np.ndarray
Field vector of for N positions of shape (3, N).
dtype: np.dtype or str
Optional choice of output data type.
Returns
-------
field: np.ndarray
Field vector in cylindrical cordinates (x, y, z) in
given datra type or input data type.
'''
cyl = Ca2Cy(cartesian)
cosphi = np.cos(cyl[1])
sinphi = np.sin(cyl[1])
field_cyl = np.array((field[0] * cosphi + field[1] * sinphi,
-field[0] * sinphi + field[1] * cosphi,
field[2]
), dtype=field.dtype if dtype is None else dtype)
return field_cyl
[docs]def Cy2CaField(cylindrical, field, dtype=None):
''' Conversion of 3d vector field from cartesian coords in cylindrical.
# polar = polar coordinate system
# x = field[r] * cos(phi) - field[phi] * sin(phi)
# y = field[r] * sin(phi) + field[phi] * cos(phi)
# z = field[z]
# output:
# field_cartesian = data in zylindrical coordinates
Parameters
----------
cylindrical: np.ndarray
Coordinate vector of for N positions of shape (3, N).
field: np.ndarray
Field vector of for N positions of shape (3, N).
dtype: np.dtype or str
Optional choice of output data type.
Returns
-------
cartesian field: np.ndarray
Field vector in cartesian cordinates (x, y, z) in
given datra type or input data type.
'''
Xcomp = np.cos(cylindrical[1])
Ycomp = np.sin(cylindrical[1])
field_cartesian = np.array((field[0] * Xcomp - field[1] * Ycomp,
field[0] * Ycomp + field[1] * Xcomp,
field[2]),
dtype=field.dtype if dtype is None else dtype)
return field_cartesian
[docs]def KtoP_field(cartesian, field, dtype=np.float64):
'''
# status: depricated, please use Ca2CyField instead.
'''
return Ca2CyField(cartesian, field, dtype=dtype)
[docs]def PtoK_field(cylindrical, field, dtype=np.float64):
'''
# status: depricated, please use Cy2CaField instead.
'''
return Cy2CaField(cylindrical, field, dtype=dtype)
[docs]def GridtoVector(*args, **kwargs):
'''
# transform the matlab grid to a python vector with the correct shape
# can take vector field data with x, y, z coordinates or simple one
# dimansional vectors
# 2D is not implemented yet
parameters
----------
order: ['F']
order = 'F' -> Fortran style = x varies fastest, instead of z
comp: [3]
number of components
# status: implemented
'''
order = kwargs.pop('order', 'F')
comp = kwargs.pop('comp', 3)
if comp == 3:
if len(args) == 1:
grid_x = args[0][0]
grid_y = args[0][1]
grid_z = args[0][2]
elif len(args) == 3:
grid_x = args[0]
grid_y = args[1]
grid_z = args[2]
vector = np.array((np.ravel(grid_x, order=order),
np.ravel(grid_y, order=order),
np.ravel(grid_z, order=order)))
elif comp == 1:
vector = np.ravel(args[0], order=order)
return vector
[docs]def VectortoGrid(vector, shape, order='F', swap=False):
'''
see 'GridtoVector'
x == VectortoGrid(GridtoVector(x), x.shape) is True
# status: implemented
'''
dx = shape[0]
dy = shape[1]
dz = shape[2]
grid_x = vector[0].reshape(dx, dy, dz, order=order)
grid_y = vector[1].reshape(dx, dy, dz, order=order)
grid_z = vector[2].reshape(dx, dy, dz, order=order)
if swap is True:
grid_x = np.swapaxes(grid_x, 0, 1)
grid_y = np.swapaxes(grid_y, 0, 1)
grid_z = np.swapaxes(grid_z, 0, 1)
return np.array((grid_x, grid_y, grid_z))
[docs]def translate(Vec, x, y, z=0, copy=False):
'''
Translate a given three dimensional vector and returns either a view of it
or a new object.
# status: implemented
'''
if isinstance(Vec, pg.core.PosVector):
Vec = R3VtoNumpy(Vec)
transposed = False
if len(Vec) != 3 and len(Vec[0]) == 3:
transposed = True
Vec = Vec.T
# print('transposed: {}'.format(transposed))
# print(np.shape(Vec))
if copy is True:
totrans = Vec.copy()
if x != 0:
totrans[0] += np.float(x)
if y != 0:
totrans[1] += np.float(y)
if z != 0:
totrans[2] += np.float(z)
if transposed is True:
return totrans.T
return totrans
if x != 0:
Vec[0] += np.float(x)
if y != 0:
Vec[1] += np.float(y)
if z != 0:
Vec[2] += np.float(z)
if transposed is True:
return Vec.T
return Vec
[docs]def translate_all(Vec, coords):
'''
Vec: Dim: num_points, 3
coords: Dim = num_dipoles, 3
output: Dim: num_dipoles, num_points, 3
'''
totrans = Vec.copy()
totrans = totrans.T[np.newaxis, :, :]
totrans = np.repeat(totrans, np.shape(coords)[0], axis=0)
totrans += coords[:, :, np.newaxis]
return np.swapaxes(totrans, 1, 2) # dim: len_coords, num_points, 3
[docs]def translateToDipole(Vec, Dipole, copy=False):
return translate(Vec, -Dipole[0], -Dipole[1], copy=copy)
[docs]def rotate3(Vec, alpha, axis='z', copy=False):
''' Rotates 3 dimensional arrays around a given axis with angle alpha. '''
# case 1/2: pygimli
if isinstance(Vec, pg.core.PosVector):
if axis == 'z':
for i in range(len(Vec)):
Vec[i].rotateZ(alpha)
if axis == 'x':
for i in range(len(Vec)):
Vec[i].rotateX(alpha)
if axis == 'y':
for i in range(len(Vec)):
Vec[i].rotateY(alpha)
return Vec
# case 2/2: numpy
transposed = False
if len(Vec) == 3 and len(Vec[0]) != 3:
transposed = True
Vec = Vec.T
if axis == 'z':
Rotationmatrix = np.array(([np.cos(-alpha), -np.sin(-alpha), 0],
[np.sin(-alpha), np.cos(-alpha), 0],
[0, 0, 1]), dtype=np.float64)
elif axis == 'x':
Rotationmatrix = np.array(([1, 0, 0],
[0, np.cos(-alpha), -np.sin(-alpha)],
[0, np.sin(-alpha), np.cos(-alpha)]),
dtype=np.float64)
elif axis == 'y':
Rotationmatrix = np.array(([np.cos(-alpha), 0, np.sin(-alpha)],
[0, 1, 0],
[-np.sin(-alpha), 0, np.cos(-alpha)]),
dtype=np.float64)
if copy is True:
torotate = Vec.copy()
for i in range(len(torotate)):
torotate[i] = torotate[i].dot(Rotationmatrix)
if transposed is True:
return torotate.T
return torotate
for i in range(len(Vec)):
Vec[i] = Vec[i].dot(Rotationmatrix)
if transposed is True:
return Vec.T
return Vec
[docs]def rotate3_all(Vec, alpha, axis='z', copy=False):
'''
'ijk,k...i->k...j'
... = number of points, broadcasted dimension
i = 3 (3 coords per point)
j = 3 (3 coords per point)
k = number of different dipoles, number of alphas
'''
RotMat = np.zeros((3, 3, len(alpha)), np.float64)
if axis == 'z':
a, b = 0, 1
c = 2
elif axis == 'x':
a, b, = 1, 2
c = 0
elif axis == 'y':
a, b = 2, 0
c = 1
RotMat[[a, b], [a, b], :] = np.cos(alpha)
RotMat[b, a, :] = np.sin(-alpha)
RotMat[a, b, :] = -RotMat[b, a, :]
RotMat[c, c, :] += 1.
# matrix multiplication for the rotation matrix and the coordinate points
if copy is True:
torotate = Vec.copy()
torotate = np.einsum('ijk,k...i->k...j', RotMat, torotate)
return torotate
Vec = np.einsum('ijk,k...i->k...j', RotMat, Vec)
return Vec
[docs]def R3VtoNumpy(R3Vector, **kwargs):
''' Creates a numpy vector from a pygimli R3Vector. '''
Array = np.array((np.array(pg.x(R3Vector)),
np.array(pg.y(R3Vector)),
np.array(pg.z(R3Vector))), **kwargs)
return Array.T # shape = (3, len(R3Vector)), order='c'
[docs]def regular_slice(dim1, dim2, direction, value):
'''
out: regular pygimli 2D-mesh object with given discretisation and
orientation with respect to a 3D coordinate system.
for now there are only x-, y- and z-orientated slices possible
# status: implemented
'''
vec1 = np.linspace(dim1[0], dim1[1], dim1[2])
vec2 = np.linspace(dim2[0], dim2[1], dim2[2])
Slice = pg.createMesh2D(pg.Vector(vec1),
pg.Vector(vec2))
Slice.translate(pg.RVector3(0, 0, value))
if direction == 'x':
Slice.swapCoordinates(0, 1)
Slice.swapCoordinates(0, 2)
# output (const value, x, y)
elif direction == 'y':
Slice.swapCoordinates(1, 2)
# output (x, const value, y)
return Slice
[docs]def regular_sliceFrom3DMesh(mesh, discretisation1, discretisation2,
direction, value):
'''
Cut a slice with regular grid discretisation from an arbitrary shaped
irregular mesh. Used for plotting purposes with matplotlib.
input: mesh
x discretisation
y discreatisation
normal direction of the slice
value for position of the slice on the normal axis
# status: implemented
'''
Pos = mesh.positions()
if direction == 'x':
dim1 = (np.min(pg.y(Pos)), np.max(pg.y(Pos)), discretisation1)
dim2 = (np.min(pg.z(Pos)), np.max(pg.z(Pos)), discretisation2)
if direction == 'y':
dim1 = (np.min(pg.x(Pos)), np.max(pg.x(Pos)), discretisation1)
dim2 = (np.min(pg.z(Pos)), np.max(pg.z(Pos)), discretisation2)
if direction == 'z':
dim1 = (np.min(pg.x(Pos)), np.max(pg.x(Pos)), discretisation1)
dim2 = (np.min(pg.y(Pos)), np.max(pg.y(Pos)), discretisation2)
return regular_slice(dim1, dim2, direction, value)
[docs]def interpolateField_Matrix(Field, InterpolationMatrix, verbose=False):
'''
# status: implemented
'''
from comet.pyhed import log
tick = time.time()
NewField = np.zeros((Field.shape[0], InterpolationMatrix.rows()),
dtype=np.complex128)
for j in range(np.shape(Field)[0]): # 0...3: x, y, z
data_real = InterpolationMatrix * pg.Vector(Field[j].real)
data_imag = InterpolationMatrix * pg.Vector(Field[j].imag)
NewField[j] = np.array(data_real) + 1j * np.array(data_imag)
log.debug(
'{:.2f} seconds for interpolation ({} * {} to {} * {} complex values)'
.format(time.time() - tick, Field.shape[0], Field.shape[1],
NewField.shape[0], NewField.shape[1]))
return NewField
[docs]def interpolateField_rotatedMatrix(Field, base_mat=None, sin_mat=None,
cos_mat=None):
'''
# status: in testing
'''
from comet.pyhed import log
tickiF_rM = time.time()
if len(Field) != 3:
if len(Field[0]) != 3:
raise Exception('Abort. Wrong field shape ({})'
.format(Field.shape))
else:
Field = Field.T
NewField = np.zeros((len(Field), base_mat.rows()),
dtype=np.complex128)
fXreal = Field[0].real
fYreal = Field[1].real
fXimag = Field[0].imag
fYimag = Field[1].imag
xreal = cos_mat * fXreal + sin_mat * fYreal
ximag = cos_mat * fXimag + sin_mat * fYimag
NewField[0] = np.array(xreal) + 1j * np.array(ximag)
yreal = sin_mat * (-fXreal) + cos_mat * fYreal
yimag = sin_mat * (-fXimag) + cos_mat * fYimag
NewField[1] = np.array(yreal) + 1j * np.array(yimag)
zreal = base_mat * Field[2].real
zimag = base_mat * Field[2].imag
NewField[2] = np.array(zreal) + 1j * np.array(zimag)
if log.getEffectiveLevel() <= 10: # debug
print('%.03f seconds for interpolation \
(%d * %d to %d * %d complex values)' %
(time.time() - tickiF_rM, Field.shape[0], Field.shape[1],
NewField.shape[0], NewField.shape[1]))
return NewField
[docs]def interpolateField(Mesh, positions, Field, interpolationMatrix=None,
verbose=False):
'''
simple case: (meshInput, meshOutput, fieldFromInputMesh) ready
# status: implemented
'''
if interpolationMatrix is None:
# six seperate calcualtions !
if isinstance(positions, pg.Mesh):
positions = positions.positions()
if not isinstance(positions, pg.core.PosVector):
positions = pg.core.PosVector(positions)
# type positions = pg.core.PosVector, regardless of input type
interpolationMatrix = Mesh.interpolationMatrix(positions)
NewField = interpolateField_Matrix(Field, interpolationMatrix,
verbose=verbose)
return NewField
[docs]def interpolateVector(Mesh, Slice, Vector, verbose=False):
'''
Interpolates a given vectorfield(Vector) based on the given Mesh to a
second mesh or slice. The field can either be real or complex.
# status: implemented
'''
tick = time.time()
NewVector = np.zeros((1, len(Slice.positions())), dtype=type(Vector[0]))
data_real = pg.Vector()
data_imag = pg.Vector()
add = ''
if type(Vector[0]) == np.complex128:
add = 'complex '
pg.interpolate(srcMesh=Mesh, inVec=Vector.real,
destPos=Slice.positions(), outVec=data_real)
pg.interpolate(srcMesh=Mesh, inVec=Vector.imag,
destPos=Slice.positions(), outVec=data_imag)
NewVector = np.asarray(data_real, dtype=np.complex128) + \
1j * np.asarray(data_imag, dtype=np.complex128)
else:
pg.interpolate(srcMesh=Mesh, inVec=Vector,
destPos=Slice.positions(), outVec=data_real)
NewVector = np.asarray(data_real)
if verbose is True:
print('%.03f seconds for interpolation (%d to %d %svalues)' %
(time.time() - tick, len(Vector), len(Slice.positions()), add))
return NewVector
[docs]def pointDataToCellData_np(mesh, field, mixed=False, weight=True):
"""
Interpolates vector- or skalarfield data defined on the nodes of the given
mesh to its cell midpoints. For now it has to be either a uniform mesh with
field (3d, complex or real) or scalar (1d, complex or real) datasets or a
mixed mesh with a simple real scalar data set (takes more time).
Parameters
----------
mesh: pg.Mesh
For now the algorithm takes only pygimli meshes.
field: array of shape (n) or (n, 3) or pg.Vector
Data set with n = number of nodes in the mesh.
mixed: bool [False]
Flag to determine if the mesh is either of mixed (True) shape (cells
can consist of variable number of nodes) or uniform (False).
weight: bool [True]
The cell data can be calculated as simple average of the surrounding
node values (False) or additionally weighted by their inverse distance
from the nodes (True).
Output
------
newfield: array of shape (c) or (3, c)
Data set with c = number of cells in the mesh.
"""
def pointDataToCellData_vector(mesh, field, weight=True):
TYPE = type(field[0, 0])
tick = time.time()
result = np.array(mesh.cellCount())
np_pos = R3VtoNumpy(mesh.positions())
np_cc = R3VtoNumpy(mesh.cellCenters())
nodespercell = len(mesh.cell(0).nodes())
idMat = np.zeros((nodespercell, len(np_cc)), dtype=int)
DisMat = np.zeros((nodespercell, mesh.cellCount()))
FieldMat = np.zeros((3, nodespercell, mesh.cellCount()), dtype=TYPE)
for ic, c in enumerate(mesh.cells()):
for iN, n in enumerate(c.nodes()):
idMat[iN, ic] = n.id()
for i in range(0, nodespercell):
DisMat[i, :] = np.sqrt(np.sum((np_pos[idMat[i]] - np_cc)**2, 1))
FieldMat[:, i, :] = field[:, idMat[i]]
weight = DisMat / np.sum(DisMat, 0)
result = np.sum(FieldMat * weight, 1)
print('pointDataToCellData %0.2f sec' % (time.time() - tick))
return result
def pointDataToCellData_scalar(mesh, field, weight=True):
TYPE = type(field[0])
tick = time.time()
result = np.array(mesh.cellCount())
np_pos = R3VtoNumpy(mesh.positions())
np_cc = R3VtoNumpy(mesh.cellCenters())
nodespercell = len(mesh.cell(0).nodes())
idMat = np.zeros((nodespercell, len(np_cc)), dtype=int)
DisMat = np.zeros((nodespercell, mesh.cellCount()))
FieldMat = np.zeros((nodespercell, mesh.cellCount()), dtype=TYPE)
for ic, c in enumerate(mesh.cells()):
for iN, n in enumerate(c.nodes()):
idMat[iN, ic] = n.id()
for i in range(0, nodespercell):
DisMat[i, :] = np.sqrt(np.sum((np_pos[idMat[i]] - np_cc)**2, 1))
FieldMat[i, :] = field[idMat[i]]
weight = DisMat / np.sum(DisMat, 0)
result = np.sum(FieldMat * weight, 0)
print('pointDataToCellData %0.2f sec' % (time.time() - tick))
return result
def pointDataToCellData_mixed(mesh, field, weight=True):
tick = time.time()
result = pg.Vector(mesh.cellCount())
for ic, c in enumerate(mesh.cells()):
nodes = [n.id() for n in c.nodes()]
w = pg.Vector(len(nodes), 1.0)
if weight is True: # simple weighting (inverse distance)
cc = c.center()
for i in range(len(nodes)):
w[i] = 1 / mesh.node(nodes[i]).pos().distance(cc)
result[ic] = pg.sum(field[nodes] * w) / pg.sum(w)
print('pointDataToCellData %0.2f sec' % (time.time() - tick))
return np.array(result)
if mixed is True:
return pointDataToCellData_mixed(mesh, field, weight=weight)
DIM = np.shape(field)[0]
if DIM > 3:
if isinstance(field, pg.Vector):
field = np.array(field)
return pointDataToCellData_scalar(mesh, field)
elif DIM == 3:
return pointDataToCellData_vector(mesh, field)
[docs]def getRSparseValues(sparseMapMatrix, indices=True, getInCRS=False):
'''
Get CRS Arrays (Row Index, Column Start_End, Values)
from SparseMatrix (CRS format).
'''
if isinstance(sparseMapMatrix, pg.core.SparseMapMatrix):
mat = pg.core.SparseMatrix(sparseMapMatrix)
else:
mat = sparseMapMatrix
vals = np.array(mat.vecVals())
if indices is True:
rows = mat.vecRowIdx()
cols = mat.vecColPtr()
if getInCRS:
return list(rows), list(cols), vals
else:
rr, cc = convertCRStoMap(list(rows), list(cols))
return rr, cc, vals
else:
return vals
[docs]def getIndicesFromConstraintMatrix(mat):
import tempfile as temp
with temp.NamedTemporaryFile(delete=True) as temp:
mat.save(temp.name)
indices = np.array(np.genfromtxt(temp.name, usecols=1), dtype=int)
id1 = indices[::2]
id2 = indices[1::2]
return (id1, id2)
[docs]def getConstraints(inv):
cmap = inv.fop().constraints()
node0, node1 = getIndicesFromConstraintMatrix(cmap)
vals = inv.cWeight().array()
return (node0, node1, vals)
[docs]def convertCRStoMap(rowIdx, colPtr):
''' Converts CRS indices to map indices. '''
ii = []
jj = []
for i in range(len(colPtr) - 1):
for j in range(colPtr[i], colPtr[i + 1]):
ii.append(i)
jj.append(rowIdx[j])
return ii, jj
[docs]def fillCRS(crsMat, rowIdx, colPtr, vals):
''' Fill CRS format SparseMatrix with values. Very Slow. '''
if not isinstance(vals, pg.Vector):
vals = pg.Vector(vals)
for i in range(crsMat.rows()):
for j in range(colPtr[i], colPtr[i + 1]):
crsMat.addVal(i, rowIdx[j], vals[j])
return crsMat
[docs]def uniqueAndSum(indices, to_sum, return_index=False, verbose=False):
"""Summs double values found by indices in a various number of arrays.
Returns the sorted unique elements of a column_stacked array of indices.
Another column_stacked array is returned with values at the unique
indices, while values at double indices are properly summed.
Parameters
----------
ar : array_like
Input array. This will be flattened if it is not already 1-D.
to_sum : array_like
Input array to be summed over axis 0. Other exsisting axes will be
broadcasted remain untouched.
return_index : bool, optional
If True, also return the indices of `ar` (along the specified axis,
if provided, or in the flattened array) that result in the unique
array.
Returns
-------
unique : ndarray
The sorted unique values.
summed_array : ndarray
The summed array, whereas all values for a specific index is the sum
over all corresponding nonunique values.
unique_indices : ndarray, optional
The indices of the first occurrences of the unique values in the
original array. Only provided if `return_index` is True.
Example
-------
>>> import numpy as np
>>> from comet.pyhed.misc.vec import uniqueAndSum
>>> idx1 = np.array([0, 0, 1, 1, 2, 2])
>>> idx2 = np.array([0, 0, 1, 2, 3, 3])
>>> # indices at positions 0 and 1 and at positions 5 and 6 are not unique
>>> to_sort = np.column_stack((idx1, idx2))
>>> # its possible to stack more than two array
>>> # you need for example 3 array to find unique node positions in a mesh
>>> values = np.arange(0.1, 0.7, 0.1)
>>> print(values)
[ 0.1 0.2 0.3 0.4 0.5 0.6]
>>> # some values to be summed together (for example attributes of nodes)
>>> unique_idx, summed_vals = uniqueAndSum(to_sort, values)
>>> print(unique_idx)
[[0 0]
[1 1]
[1 2]
[2 3]]
>>> print(summed_vals)
[ 0.3 0.3 0.4 1.1]
>>> # [0.1 + 0.2, 03., 0.4, 0.5 + 0.6]
"""
flag_mult = len(indices) != indices.size
if verbose:
print('Get {} indices for sorting'.format(np.shape(indices)))
if flag_mult:
ar = indices.ravel().view(
np.dtype((np.void,
indices.dtype.itemsize * indices.shape[1]))).flatten()
else:
ar = np.asanyarray(indices).flatten()
to_sum = np.asanyarray(to_sum)
if ar.size == 0:
ret = (ar,)
ret += (to_sum)
if return_index:
ret += (np.empty(0, np.bool),)
return ret
if verbose:
print('Performing argsort...')
perm = ar.argsort(kind='mergesort')
aux = ar[perm]
flag = np.concatenate(([True], aux[1:] != aux[:-1]))
if flag_mult:
ret = (indices[perm[flag]],)
else:
ret = (aux[flag],) # unique indices
if verbose:
print('Identified {} unique indices'.format(np.shape(ret)))
if verbose:
print('Performing reduceat...')
summed = np.add.reduceat(to_sum[perm], np.nonzero(flag)[0])
ret += (summed,) # summed values
if return_index:
ret += (perm[flag],) # optional: indices
return ret
[docs]def perimeterFromPolyPoints(points, circle_radius=None, closed=True):
""" Get perimeter of a polygon. """
if circle_radius is not None:
perimeter = np.pi * circle_radius * 2.0
else:
p2 = np.row_stack((points[1:], points[0]))
perimeter = np.sum(np.sqrt((p2[:, 0] - points[:, 0])**2 +
(p2[:, 1] - points[:, 1])**2))
if not closed:
perimeter -= np.sqrt((points[-1, 0] - points[0, 0])**2 +
(points[-1, 1] - points[0, 1])**2)
return perimeter
[docs]def areaFromPolyPoints(points):
""" Get perimeter of a polygon. """
points = np.asarray(points)
if len(points) == 1:
# assuming dipole
return 0.0
if len(points) == 2:
# assuming unrotated rectangle
return (np.diff(points[:, 0]) * np.diff(points[:, 1]))[0]
# n polygon area formula, not for figure of eight
p2 = np.row_stack((points[1:], points[0]))
raw = np.sum(points[:, 0] * p2[:, 1]) - np.sum(points[:, 1] * p2[:, 0])
area = np.abs(raw/2.0)
return area
[docs]def sumBetweenIndices(array, indices, use_thickness=False, axis=0):
"""
Split array at given axis and indices and sum up the parts.
If use_thickness is True, the difference of the absolute values are used,
usfull for array representing a depth for example.
"""
if axis != 0: # swap, so that axis = 0 is to be reduced
array = np.swapaxes(array, 0, axis)
new_array = []
if use_thickness:
new_array.append(array[0])
array = np.diff(array, axis=0)
for i in range(0, len(indices) - 1):
new_array.append(np.sum(array[indices[i]:indices[i + 1]], axis=0))
if use_thickness:
return np.cumsum(np.array(new_array), axis=0)
new_array = np.array(new_array)
if axis != 0: # swap back
new_array = np.swapaxes(new_array, 0, axis)
return new_array
[docs]def cumsumDepth(a, min_thk=0.5):
"""
Summs part of a array, until all layers have a given minimum thickness.
only use on array with increasing thickness.
"""
thk = np.diff(np.abs(a))
indices = [0]
new_val = 0.0
for i in range(len(thk)):
new_val += thk[i]
if i == len(thk) - 1:
indices[-1] = i + 1 # last layer to small, append to second last
elif new_val > min_thk:
new_val = 0.0
indices.append(i + 1)
if i <= len(a) - 1:
if thk[i + 1] >= min_thk:
indices.extend(np.arange(i + 2, len(a), 1))
break
else:
continue
return indices, sumBetweenIndices(a, indices, use_thickness=True)
[docs]def sinhspace(start, stop, num_step):
""" Like linspace but using a hyperbolic sine function. """
off = start
start2 = 0
end2 = stop - start
vec = np.sinh(np.linspace(
np.arcsinh(start2), np.arcsinh(end2), num_step))
final = vec + off
return final
[docs]def sinhZVolumeFunction(z, z_range=[0, -100], area_range=[0.1, 100]):
""" Maps values from z per z_range to area_range using a sinh
function instead of linear interpolation.
Values outside z_range are assigned the limits of area_range.
Example
-------
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> area = []
>>> zrange = np.linspace(30, -130, 160)
>>> for z in zrange:
>>> area.append(sinhZVolumeFunction([z]))
>>> fig, ax = plt.subplots(1, 1)
>>> ax.plot(zrange, area)
"""
if z_range[0] < z_range[1]:
z_range = np.array(z_range)[::-1]
area_range = np.array(area_range)[::-1]
inverse = False
if area_range[0] > area_range[1]:
inverse = True
if z >= z_range[0]:
area = area_range[0]
elif z <= z_range[1]:
area = area_range[1]
else:
z_width = z_range[1] - z_range[0]
z_perc = (z - z_range[0]) / z_width
trans_max = np.arcsinh(np.abs(area_range[1] - area_range[0]))
if inverse:
area = np.sinh((1 - z_perc) * trans_max) + area_range[1]
else:
area = np.sinh(z_perc * trans_max) + area_range[0]
return area
[docs]def linspace2D(Point1, Point2, num):
""" Internal function. Like linspace but for twodimensional points. """
return np.array((np.linspace(Point1[0], Point2[0], num),
np.linspace(Point1[1], Point2[1], num)),
np.float64).T
[docs]def linspace3D(Point1, Point2, num):
""" Internal function. Like linspace but for threedimensional points. """
return np.array((np.linspace(Point1[0], Point2[0], num),
np.linspace(Point1[1], Point2[1], num),
np.linspace(Point1[2], Point2[2], num)),
np.float64).T
[docs]def rotationMatrix(axis, theta):
"""
Return the rotation matrix associated with counterclockwise rotation about
the given axis (3, n) by theta radians (n,). Supports broadcasting along
second axis of input **axis**.
Outputs:
--------
Array with rotation matrices of shape (3, 3, n) or (3, 3) if n==1.
"""
# norm axis, otherwise the rotation warps the vector
axis = np.array(axis)
assert axis.shape[0] == 3, '{} != {}'.format(axis.shape, '(3, 3, ...)')
if axis.shape == (3,):
axis = axis[:, np.newaxis]
theta = np.atleast_1d(theta)
normed = scipy.linalg.norm(axis, axis=0)
coaxial = np.isclose(normed, 0)
noco = np.logical_not(coaxial)
# mats for coaxial edge case
rotmats = np.eye(3)
if len(axis.shape) == 2:
rotmats = np.tile(rotmats[:, :, np.newaxis], (1, 1, axis.shape[1]))
if not np.all(coaxial):
# calc real rotation matrices where needed
axis_normed = axis[:, noco] / normed[noco]
a = np.cos(theta[noco] / 2.0)
b, c, d = -axis_normed * np.sin(theta[noco] / 2.0)
aa = a * a
bb = b * b
cc = c * c
dd = d * d
bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
rm_noco = np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
[2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
[2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
# exchange dummy rot matrices with real ones where needed
rotmats[:, :, noco] = rm_noco
return np.squeeze(rotmats)
[docs]def angle(ax1, ax2):
"""
Returns angle between two arbitrary vectors of shape (3, ...).
Allows broadcasting.
"""
if len(np.shape(ax1)) == 2:
a1n = ax1 / scipy.linalg.norm(ax1, axis=-1)[:, np.newaxis]
a2n = ax2 / scipy.linalg.norm(ax2, axis=-1)[:, np.newaxis]
else:
a1n = ax1 / scipy.linalg.norm(ax1)
a2n = ax2 / scipy.linalg.norm(ax2)
dot = np.einsum('ij, ij->i',
np.atleast_2d(a1n),
np.atleast_2d(a2n))
rot_angle = np.arccos(np.clip(dot, -1.0, 1.0))
return rot_angle
# u = A(:)/norm(A);
# % a and b must be column vectors
# v = B(:)/norm(B); % of equal length
# N = length(u);
# S = reflection( eye(N), v+u ); % S*u = -v, S*v = -u
# R = reflection( S, v ); % v = R*u
# function v = reflection(u, n ) % Reflection of u on hyperplane n.
# % u can be a matrix. u and v must have the same number of rows.
# v = u - 2 * n * (n'*u) / (n'*n);
# end
[docs]def rotFromAtoB(vec, ax1, ax2):
"""
Rotates input vector **vec** from one direction **ax1** (x, y, z) to
another direction **ax2** (x, y, z).
"""
if len(np.shape(ax1)) == 2:
a1n = ax1 / scipy.linalg.norm(ax1, axis=-1)[:, np.newaxis]
a2n = ax2 / scipy.linalg.norm(ax2, axis=-1)[:, np.newaxis]
else:
a1n = ax1 / scipy.linalg.norm(ax1)
a2n = ax2 / scipy.linalg.norm(ax2)
assert ax1.shape == ax2.shape, '{} != {}'.format(ax1.shape, ax2.shape)
rax = np.cross(a1n, a2n)
rang = angle(a1n, a2n)
rotmats = rotationMatrix(rax.T, -rang)
if rotmats.shape == (3, 3):
rotated_vec = np.dot(vec, rotmats)
else:
rotated_vec = np.squeeze(np.einsum('ki,ijk->kj',
vec,
rotmats))
return rotated_vec
# THE END