DAMASK_EICMD/python/damask/util.py

663 lines
19 KiB
Python
Raw Normal View History

2021-05-07 23:12:23 +05:30
"""Miscellaneous helper functionality."""
2019-05-30 23:32:55 +05:30
import sys
import datetime
2019-05-30 23:32:55 +05:30
import os
import subprocess
import shlex
import re
import fractions
2020-02-21 11:51:45 +05:30
from functools import reduce
2019-05-20 23:24:57 +05:30
2020-02-21 11:51:45 +05:30
import numpy as np
2021-03-20 04:19:41 +05:30
import h5py
2019-05-20 23:24:57 +05:30
from . import version
# limit visibility
__all__=[
'srepr',
'emph','deemph','warn','strikeout',
'execute',
2021-04-03 14:38:22 +05:30
'natural_sort',
'show_progress',
'scale_to_coprime',
'project_stereographic',
2020-09-28 11:10:43 +05:30
'hybrid_IA',
'execution_stamp',
2020-11-15 00:21:15 +05:30
'shapeshifter', 'shapeblender',
2021-03-20 04:19:41 +05:30
'extend_docstring', 'extended_docstring',
'Bravais_to_Miller', 'Miller_to_Bravais',
2021-03-30 23:11:36 +05:30
'DREAM3D_base_group', 'DREAM3D_cell_data_group',
'dict_prune', 'dict_flatten'
]
2021-03-27 12:05:49 +05:30
# https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
# https://stackoverflow.com/questions/287871
_colors = {
'header' : '\033[95m',
'OK_blue': '\033[94m',
'OK_green': '\033[92m',
'warning': '\033[93m',
'fail': '\033[91m',
'end_color': '\033[0m',
'bold': '\033[1m',
'dim': '\033[2m',
'underline': '\033[4m',
'crossout': '\033[9m'
}
####################################################################################################
# Functions
####################################################################################################
def srepr(arg,glue = '\n'):
2020-02-22 04:36:51 +05:30
r"""
2021-04-24 21:30:57 +05:30
Join items with glue string.
Parameters
----------
arg : iterable
Items to join.
glue : str, optional
Glue used for joining operation. Defaults to \n.
2021-04-24 21:30:57 +05:30
Returns
-------
joined : str
String representation of the joined items.
"""
2021-04-03 11:38:48 +05:30
if (not hasattr(arg, 'strip') and
(hasattr(arg, '__getitem__') or
hasattr(arg, '__iter__'))):
return glue.join(str(x) for x in arg)
2021-04-05 20:02:28 +05:30
else:
return arg if isinstance(arg,str) else repr(arg)
def emph(what):
2021-04-24 21:30:57 +05:30
"""
Format with emphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
2021-03-27 12:05:49 +05:30
return _colors['bold']+srepr(what)+_colors['end_color']
def deemph(what):
2021-04-24 21:30:57 +05:30
"""
Format with deemphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
2021-03-27 12:05:49 +05:30
return _colors['dim']+srepr(what)+_colors['end_color']
def warn(what):
2021-04-24 21:30:57 +05:30
"""
Format for warning.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
2021-03-27 12:05:49 +05:30
return _colors['warning']+emph(what)+_colors['end_color']
def strikeout(what):
2021-04-24 21:30:57 +05:30
"""
Format as strikeout.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
2021-03-27 12:05:49 +05:30
return _colors['crossout']+srepr(what)+_colors['end_color']
2021-03-27 12:05:49 +05:30
def execute(cmd,wd='./',env=None):
"""
Execute command.
Parameters
----------
cmd : str
Command to be executed.
wd : str, optional
Working directory of process. Defaults to ./ .
env : dict, optional
Environment for execution.
2021-04-24 21:30:57 +05:30
Returns
-------
stdout, stderr : str
Output of the executed command.
"""
print(f"executing '{cmd}' in '{wd}'")
2021-03-27 12:05:49 +05:30
process = subprocess.run(shlex.split(cmd),
stdout = subprocess.PIPE,
stderr = subprocess.PIPE,
env = os.environ if env is None else env,
cwd = wd,
encoding = 'utf-8')
if process.returncode != 0:
2021-03-27 12:05:49 +05:30
print(process.stdout)
print(process.stderr)
raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}")
2021-03-27 12:05:49 +05:30
return process.stdout, process.stderr
2021-04-03 14:38:22 +05:30
def natural_sort(key):
2021-05-07 23:12:23 +05:30
"""
Natural sort.
For use in python's 'sorted'.
References
----------
https://en.wikipedia.org/wiki/Natural_sort_order
"""
2021-04-03 14:38:22 +05:30
convert = lambda text: int(text) if text.isdigit() else text
return [ convert(c) for c in re.split('([0-9]+)', key) ]
def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
"""
2021-05-11 00:14:58 +05:30
Decorate a loop with a progress bar.
Use similar like enumerate.
Parameters
----------
2021-05-11 00:14:58 +05:30
iterable : iterable or function with yield statement
Iterable (or function with yield statement) to be decorated.
2021-05-07 23:12:23 +05:30
N_iter : int, optional
2021-05-11 00:14:58 +05:30
Total number of iterations. Required unless obtainable as len(iterable).
2021-05-07 23:12:23 +05:30
prefix : str, optional
Prefix string.
bar_length : int, optional
2021-05-11 00:14:58 +05:30
Length of progress bar in characters. Defaults to 50.
"""
if N_iter in [0,1] or (hasattr(iterable,'__len__') and len(iterable) <= 1):
2021-03-31 17:57:36 +05:30
for item in iterable:
yield item
else:
status = _ProgressBar(N_iter if N_iter is not None else len(iterable),prefix,bar_length)
2021-03-31 17:57:36 +05:30
for i,item in enumerate(iterable):
yield item
status.update(i)
def scale_to_coprime(v):
"""
Scale vector to co-prime (relatively prime) integers.
Parameters
----------
v : numpy.ndarray of shape (:)
Vector to scale.
2021-04-24 21:30:57 +05:30
Returns
-------
m : numpy.ndarray of shape (:)
Vector scaled to co-prime numbers.
"""
2020-06-25 11:59:36 +05:30
MAX_DENOMINATOR = 1000000
def get_square_denominator(x):
"""Denominator of the square of a number."""
return fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
def lcm(a,b):
"""Least common multiple."""
try:
return np.lcm(a,b) # numpy > 1.18
except AttributeError:
return a * b // np.gcd(a, b)
2021-02-22 13:05:15 +05:30
m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(int)
m = m//reduce(np.gcd,m)
with np.errstate(invalid='ignore'):
if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]):
raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?')
return m
def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
"""
Apply stereographic projection to vector.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
Vector coordinates to be projected.
direction : str
Projection direction 'x', 'y', or 'z'.
Defaults to 'z'.
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
Maintain three-dimensional output coordinates.
Default two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Returns
-------
coordinates : numpy.ndarray of shape (...,2 | 3)
Projected coordinates.
Examples
--------
>>> project_stereographic(np.ones(3))
[0.3660254, 0.3660254]
>>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
[0, 0.5, 0.5]
>>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False)
[0.41421356, 0]
"""
shift = 'zyx'.index(direction)
v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
shift,axis=-1)
return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def execution_stamp(class_name,function_name=None):
"""Timestamp the execution of a (function within a) class."""
now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
_function_name = '' if function_name is None else f'.{function_name}'
return f'damask.{class_name}{_function_name} v{version} ({now})'
def hybrid_IA(dist,N,rng_seed=None):
"""
Hybrid integer approximation.
Parameters
----------
dist : numpy.ndarray
Distribution to be approximated
N : int
Number of samples to draw.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
"""
N_opt_samples,N_inv_samples = (max(np.count_nonzero(dist),N),0) # random subsampling if too little samples requested
2020-09-28 11:10:43 +05:30
scale_,scale,inc_factor = (0.0,float(N_opt_samples),1.0)
while (not np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
repeats = np.rint(scale*dist).astype(int)
N_inv_samples = np.sum(repeats)
2020-09-28 11:10:43 +05:30
scale_,scale,inc_factor = (scale,scale+inc_factor*0.5*(scale - scale_), inc_factor*2.0) \
if N_inv_samples < N_opt_samples else \
(scale_,0.5*(scale_ + scale), 1.0)
return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]]
2020-09-28 11:10:43 +05:30
def shapeshifter(fro,to,mode='left',keep_ones=False):
"""
Return a tuple that reshapes 'fro' to become broadcastable to 'to'.
Parameters
----------
fro : tuple
Original shape of array.
to : tuple
Target shape of array after broadcasting.
len(to) cannot be less than len(fro).
mode : str, optional
Indicates whether new axes are preferably added to
either 'left' or 'right' of the original shape.
Defaults to 'left'.
keep_ones : bool, optional
Treat '1' in fro as literal value instead of dimensional placeholder.
Defaults to False.
"""
beg = dict(left ='(^.*\\b)',
right='(^.*?\\b)')
sep = dict(left ='(.*\\b)',
right='(.*?\\b)')
end = dict(left ='(.*?$)',
right='(.*$)')
fro = (1,) if not len(fro) else fro
to = (1,) if not len(to) else to
try:
grp = re.match(beg[mode]
+f',{sep[mode]}'.join(map(lambda x: f'{x}'
if x>1 or (keep_ones and len(fro)>1) else
'\\d+',fro))
+f',{end[mode]}',
','.join(map(str,to))+',').groups()
except AttributeError:
raise ValueError(f'Shapes can not be shifted {fro} --> {to}')
fill = ()
for g,d in zip(grp,fro+(None,)):
fill += (1,)*g.count(',')+(d,)
return fill[:-1]
def shapeblender(a,b):
"""
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
Parameters
----------
a : tuple
Shape of first array.
b : tuple
Shape of second array.
Examples
--------
>>> shapeblender((4,4,3),(3,2,1))
(4,4,3,2,1)
>>> shapeblender((1,2),(1,2,3))
(1,2,3)
>>> shapeblender((1,),(2,2,1))
(1,2,2,1)
>>> shapeblender((3,2),(3,2))
(3,2)
"""
i = min(len(a),len(b))
while i > 0 and a[-i:] != b[:i]: i -= 1
return a + b[i:]
def extend_docstring(extra_docstring):
"""
Decorator: Append to function's docstring.
Parameters
----------
extra_docstring : str
Docstring to append.
"""
def _decorator(func):
func.__doc__ += extra_docstring
return func
return _decorator
def extended_docstring(f,extra_docstring):
"""
Decorator: Combine another function's docstring with a given docstring.
2020-11-15 00:21:15 +05:30
Parameters
----------
f : function
Function of which the docstring is taken.
extra_docstring : str
Docstring to append.
2020-11-15 00:21:15 +05:30
"""
def _decorator(func):
func.__doc__ = f.__doc__ + extra_docstring
return func
return _decorator
2020-11-15 00:21:15 +05:30
2021-03-20 04:19:41 +05:30
def DREAM3D_base_group(fname):
2021-03-23 18:58:56 +05:30
"""
Determine the base group of a DREAM.3D file.
The base group is defined as the group (folder) that contains
a 'SPACING' dataset in a '_SIMPL_GEOMETRY' group.
Parameters
----------
2021-04-24 21:30:57 +05:30
fname : str or pathlib.Path
2021-03-23 18:58:56 +05:30
Filename of the DREAM.3D (HDF5) file.
2021-04-24 21:30:57 +05:30
Returns
-------
path : str
Path to the base group.
2021-03-23 18:58:56 +05:30
"""
2021-03-20 04:19:41 +05:30
with h5py.File(fname,'r') as f:
base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None)
2021-03-23 18:58:56 +05:30
2021-03-20 04:19:41 +05:30
if base_group is None:
2021-03-23 19:30:59 +05:30
raise ValueError(f'Could not determine base group in file {fname}.')
2021-03-23 18:58:56 +05:30
2021-03-20 04:19:41 +05:30
return base_group
2021-03-23 18:58:56 +05:30
def DREAM3D_cell_data_group(fname):
"""
Determine the cell data group of a DREAM.3D file.
The cell data group is defined as the group (folder) that contains
a dataset in the base group whose length matches the total number
of points as specified in '_SIMPL_GEOMETRY/DIMENSIONS'.
Parameters
----------
2021-04-24 21:30:57 +05:30
fname : str or pathlib.Path
2021-03-23 18:58:56 +05:30
Filename of the DREAM.3D (HDF5) file.
2021-04-24 21:30:57 +05:30
Returns
-------
path : str
Path to the cell data group.
2021-03-23 18:58:56 +05:30
"""
base_group = DREAM3D_base_group(fname)
with h5py.File(fname,'r') as f:
cells = tuple(f['/'.join([base_group,'_SIMPL_GEOMETRY','DIMENSIONS'])][()][::-1])
2021-03-23 18:58:56 +05:30
cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \
2021-03-23 19:30:59 +05:30
if isinstance(obj,h5py._hl.dataset.Dataset) and np.shape(obj)[:-1] == cells \
2021-03-23 18:58:56 +05:30
else None)
if cell_data_group is None:
2021-03-23 19:30:59 +05:30
raise ValueError(f'Could not determine cell data group in file {fname}/{base_group}.')
2021-03-23 18:58:56 +05:30
return cell_data_group
2021-03-31 14:29:21 +05:30
def Bravais_to_Miller(*,uvtw=None,hkil=None):
"""
Transform 4 MillerBravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
Parameters
----------
uvtw|hkil : numpy.ndarray of shape (...,4)
MillerBravais indices of crystallographic direction [uvtw] or plane normal (hkil).
Returns
-------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of [uvw] direction or (hkl) plane normal.
"""
if (uvtw is not None) ^ (hkil is None):
raise KeyError('Specify either "uvtw" or "hkil"')
axis,basis = (np.array(uvtw),np.array([[1,0,-1,0],
[0,1,-1,0],
[0,0, 0,1]])) \
if hkil is None else \
(np.array(hkil),np.array([[1,0,0,0],
[0,1,0,0],
[0,0,0,1]]))
return np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*,uvw=None,hkl=None):
"""
Transform 3 Miller indices to 4 MillerBravais indices of crystal direction [uvtw] or plane normal (hkil).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of crystallographic direction [uvw] or plane normal (hkl).
Returns
-------
uvtw|hkil : numpy.ndarray of shape (...,4)
MillerBravais indices of [uvtw] direction or (hkil) plane normal.
"""
if (uvw is not None) ^ (hkl is None):
raise KeyError('Specify either "uvw" or "hkl"')
axis,basis = (np.array(uvw),np.array([[ 2,-1, 0],
[-1, 2, 0],
[-1,-1, 0],
[ 0, 0, 3]])/3) \
if hkl is None else \
(np.array(hkl),np.array([[ 1, 0, 0],
[ 0, 1, 0],
[-1,-1, 0],
[ 0, 0, 1]]))
return np.einsum('il,...l',basis,axis)
def dict_prune(d):
2021-03-31 14:29:21 +05:30
"""
Recursively remove empty dictionaries.
2021-03-31 14:29:21 +05:30
Parameters
----------
d : dict
Dictionary to prune.
Returns
-------
pruned : dict
Pruned dictionary.
2021-03-31 14:29:21 +05:30
"""
2021-03-30 23:11:36 +05:30
# https://stackoverflow.com/questions/48151953
2021-03-31 14:29:21 +05:30
new = {}
for k,v in d.items():
2021-03-30 23:11:36 +05:30
if isinstance(v, dict):
v = dict_prune(v)
2021-03-30 23:11:36 +05:30
if not isinstance(v,dict) or v != {}:
2021-03-31 14:29:21 +05:30
new[k] = v
2021-05-07 23:12:23 +05:30
2021-03-31 14:29:21 +05:30
return new
2021-03-30 23:11:36 +05:30
def dict_flatten(d):
2021-03-31 14:29:21 +05:30
"""
Recursively remove keys of single-entry dictionaries.
2021-03-31 14:29:21 +05:30
Parameters
----------
d : dict
Dictionary to flatten.
Returns
-------
flattened : dict
Flattened dictionary.
2021-03-31 14:29:21 +05:30
"""
if isinstance(d,dict) and len(d) == 1:
2021-03-31 17:57:36 +05:30
entry = d[list(d.keys())[0]]
new = dict_flatten(entry.copy()) if isinstance(entry,dict) else entry
else:
new = {k: (dict_flatten(v) if isinstance(v, dict) else v) for k,v in d.items()}
2021-03-31 17:57:36 +05:30
2021-03-31 14:29:21 +05:30
return new
2021-03-30 23:11:36 +05:30
2021-03-27 12:05:49 +05:30
####################################################################################################
# Classes
####################################################################################################
class _ProgressBar:
"""
Report progress of an interation as a status bar.
Works for 0-based loops, ETA is estimated by linear extrapolation.
2019-09-20 01:02:15 +05:30
"""
def __init__(self,total,prefix,bar_length):
"""
Set current time as basis for ETA estimation.
Parameters
----------
total : int
Total # of iterations.
prefix : str
Prefix string.
bar_length : int
Character length of bar.
"""
self.total = total
self.prefix = prefix
self.bar_length = bar_length
2021-07-09 14:59:52 +05:30
self.time_start = self.time_last_update = datetime.datetime.now()
self.fraction_last = 0.0
sys.stderr.write(f"{self.prefix} {''*self.bar_length} 0% ETA n/a")
sys.stderr.flush()
def update(self,iteration):
fraction = (iteration+1) / self.total
filled_length = int(self.bar_length * fraction)
2021-07-09 14:59:52 +05:30
if filled_length > int(self.bar_length * self.fraction_last) or \
datetime.datetime.now() - self.time_last_update > datetime.timedelta(seconds=10):
self.time_last_update = datetime.datetime.now()
bar = '' * filled_length + '' * (self.bar_length - filled_length)
2021-07-09 14:59:52 +05:30
remaining_time = (datetime.datetime.now() - self.time_start) \
* (self.total - (iteration+1)) / (iteration+1)
2021-02-26 11:05:42 +05:30
remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs
2020-06-25 01:04:51 +05:30
sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}')
sys.stderr.flush()
2021-07-09 14:59:52 +05:30
self.fraction_last = fraction
if iteration == self.total - 1:
sys.stderr.write('\n')
sys.stderr.flush()