diff --git a/lib/damask/util.py b/lib/damask/util.py index b3e35812a..1b239b7b0 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -6,11 +6,13 @@ import numpy as np from optparse import Option class bcolors: - ''' - ASCII Colors (Blender code) - https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py - http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python - ''' + """ + ASCII Colors (Blender code) + + https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py + http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python + """ + HEADER = '\033[95m' OKBLUE = '\033[94m' OKGREEN = '\033[92m' @@ -32,30 +34,28 @@ class bcolors: # ----------------------------- -def srepr(arg, - glue = '\n'): -# ----------------------------- - if (not hasattr(arg, "strip") and - hasattr(arg, "__getitem__") or - hasattr(arg, "__iter__")): - return glue.join(srepr(x) for x in arg) - return arg if isinstance(arg,basestring) else repr(arg) +def srepr(arg,glue = '\n'): + """joins arguments as individual lines""" + if (not hasattr(arg, "strip") and + hasattr(arg, "__getitem__") or + hasattr(arg, "__iter__")): + return glue.join(srepr(x) for x in arg) + return arg if isinstance(arg,basestring) else repr(arg) # ----------------------------- -def croak(what, - newline = True): -# ----------------------------- +def croak(what, newline = True): + """writes formated to stderr""" sys.stderr.write(srepr(what,glue = '\n') + ('\n' if newline else '')) sys.stderr.flush() # ----------------------------- def report(who,what): -# ----------------------------- + """reports script and file name""" croak( (emph(who) if who else '') + (': '+what if what else '') ) # ----------------------------- def emph(what): -# ----------------------------- + """emphasizes string on screen""" return bcolors.BOLD+srepr(what)+bcolors.ENDC # ----------------------------- @@ -68,7 +68,6 @@ for f in ['cos', 'sin', 'tan']: # ----------------------------- def gridLocation(idx,res): -# ----------------------------- return ( idx % res[0], \ ( idx // res[0]) % res[1], \ ( idx // res[0] // res[1]) % res[2] ) @@ -76,7 +75,6 @@ def gridLocation(idx,res): # ----------------------------- def gridIndex(location,res): -# ----------------------------- return ( location[0] % res[0] + \ ( location[1] % res[1]) * res[0] + \ ( location[2] % res[2]) * res[1] * res[0] ) @@ -84,10 +82,12 @@ def gridIndex(location,res): # ----------------------------- class extendableOption(Option): -# ----------------------------- -# used for definition of new option parser action 'extend', which enables to take multiple option arguments -# taken from online tutorial http://docs.python.org/library/optparse.html - + """ + used for definition of new option parser action 'extend', which enables to take multiple option arguments + + taken from online tutorial http://docs.python.org/library/optparse.html + """ + ACTIONS = Option.ACTIONS + ("extend",) STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",) TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",) @@ -102,28 +102,34 @@ class extendableOption(Option): # ----------------------------- class backgroundMessage(threading.Thread): -# ----------------------------- + """reporting with animation to indicate progress""" + choices = {'bounce': ['_','o','O','°','¯','¯','°','O','o',], 'circle': [u'\u25f4',u'\u25f5',u'\u25f6',u'\u25f7'], 'hexagon': [u'\u2b22',u'\u2b23'], 'square': [u'\u2596',u'\u2598',u'\u259d',u'\u2597'], 'triangle': [u'\u140a',u'\u140a',u'\u1403',u'\u1405',u'\u1405',u'\u1403'], - 'amoeba': [u'\u2596',u'\u258f',u'\u2598',u'\u2594',u'\u259d',u'\u2595',u'\u2597',u'\u2582'], - 'beat': [u'\u2581',u'\u2582',u'\u2583',u'\u2585',u'\u2586',u'\u2587',u'\u2587',u'\u2586',u'\u2585',u'\u2583',u'\u2582',], - 'prison': [u'\u168b',u'\u168c',u'\u168d',u'\u168f',u'\u168e',u'\u168d',u'\u168c',u'\u168b',], - 'breath': [u'\u1690',u'\u1691',u'\u1692',u'\u1693',u'\u1694',u'\u1693',u'\u1692',u'\u1691',u'\u1690',], + 'amoeba': [u'\u2596',u'\u258f',u'\u2598',u'\u2594',u'\u259d',u'\u2595', + u'\u2597',u'\u2582'], + 'beat': [u'\u2581',u'\u2582',u'\u2583',u'\u2585',u'\u2586',u'\u2587', + u'\u2587',u'\u2586',u'\u2585',u'\u2583',u'\u2582',], + 'prison': [u'\u168b',u'\u168c',u'\u168d',u'\u168f',u'\u168e',u'\u168d', + u'\u168c',u'\u168b',], + 'breath': [u'\u1690',u'\u1691',u'\u1692',u'\u1693',u'\u1694',u'\u1693', + u'\u1692',u'\u1691',u'\u1690',], 'pulse': [u'·',u'•',u'\u25cf',u'\u25cf',u'•',], - 'ant': [u'\u2801',u'\u2802',u'\u2810',u'\u2820',u'\u2804',u'\u2840',u'\u2880',u'\u2820',u'\u2804',u'\u2802',u'\u2810',u'\u2808'], - 'juggle': [u'\ua708',u'\ua709',u'\ua70a',u'\ua70b',u'\ua70c',u'\ua711',u'\ua710',u'\ua70f',u'\ua70d',], + 'ant': [u'\u2801',u'\u2802',u'\u2810',u'\u2820',u'\u2804',u'\u2840', + u'\u2880',u'\u2820',u'\u2804',u'\u2802',u'\u2810',u'\u2808'], + 'juggle': [u'\ua708',u'\ua709',u'\ua70a',u'\ua70b',u'\ua70c',u'\ua711', + u'\ua710',u'\ua70f',u'\ua70d',], # 'wobbler': [u'\u2581',u'\u25e3',u'\u258f',u'\u25e4',u'\u2594',u'\u25e5',u'\u2595',u'\u25e2',], 'grout': [u'\u2581',u'\u258f',u'\u2594',u'\u2595',], 'partner': [u'\u26ac',u'\u26ad',u'\u26ae',u'\u26af',u'\u26ae',u'\u26ad',], 'classic': ['-', '\\', '|', '/',], } - def __init__(self, - symbol = None, - wait = 0.1): + def __init__(self,symbol = None,wait = 0.1): + """sets animation symbol""" super(backgroundMessage, self).__init__() self._stop = threading.Event() self.message = '' @@ -134,6 +140,7 @@ class backgroundMessage(threading.Thread): self.waittime = wait def __quit__(self): + """cleans output""" length = len(self.symbols[self.counter] + self.gap + self.message) sys.stderr.write(chr(8)*length + ' '*length + chr(8)*length) sys.stderr.write('') @@ -146,8 +153,7 @@ class backgroundMessage(threading.Thread): return self._stop.is_set() def run(self): -# while not threading.enumerate()[0]._Thread__stopped: - while not self.stopped(): + while not threading.enumerate()[0]._Thread__stopped: time.sleep(self.waittime) self.update_message() self.__quit__() @@ -170,24 +176,38 @@ class backgroundMessage(threading.Thread): def animation(self,which = None): return ''.join(self.choices[which]) if which in self.choices else '' -''' -Non-linear least square fitting (Levenberg-Marquardt method) with -bounded parameters. -the codes of transformation between int <-> ext refers to the work of -Jonathan J. Helmus: https://github.com/jjhelmus/leastsqbound-scipy -other codes refers to the source code of minpack.py: -..\Lib\site-packages\scipy\optimize\minpack.py -''' -from numpy import (array, arcsin, asarray, cos, dot, eye, empty_like, - isscalar,finfo, take, triu, transpose, sqrt, sin) -def _check_func(checker, argname, thefunc, x0, args, numinputs, +def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0, + col_deriv=0, ftol=1.49012e-8, xtol=1.49012e-8, + gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None): + from scipy.optimize import _minpack + """ + Non-linear least square fitting (Levenberg-Marquardt method) with + bounded parameters. + the codes of transformation between int <-> ext refers to the work of + Jonathan J. Helmus: https://github.com/jjhelmus/leastsqbound-scipy + other codes refers to the source code of minpack.py: + ..\Lib\site-packages\scipy\optimize\minpack.py + + An internal parameter list is used to enforce contraints on the fitting + parameters. The transfomation is based on that of MINUIT package. + please see: F. James and M. Winkler. MINUIT User's Guide, 2004. + + bounds : list + (min, max) pairs for each parameter, use None for 'min' or 'max' + when there is no bound in that direction. + For example: if there are two parameters needed to be fitting, then + bounds is [(min1,max1), (min2,max2)] + + This function is based on 'leastsq' of minpack.py, the annotation of + other parameters can be found in 'leastsq'. + ..\Lib\site-packages\scipy\optimize\minpack.py + """ + + def _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape=None): - from numpy import atleast_1d, shape, issubdtype, dtype, inexact - ''' - The same as that of minpack.py, - ''' - res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) + """The same as that of minpack.py""" + res = np.atleast_1d(thefunc(*((x0[:numinputs],) + args))) if (output_shape is not None) and (shape(res) != output_shape): if (output_shape[0] != 1): if len(output_shape) > 1: @@ -201,203 +221,171 @@ def _check_func(checker, argname, thefunc, x0, args, numinputs, else: msg += "." raise TypeError(msg) - if issubdtype(res.dtype, inexact): + if np.issubdtype(res.dtype, np.inexact): dt = res.dtype else: dt = dtype(float) return shape(res), dt - -def _int2extGrad(p_int, bounds): - """ - Calculate the gradients of transforming the internal (unconstrained) - to external (constained) parameter. - """ - grad = empty_like(p_int) + + def _int2extGrad(p_int, bounds): + """Calculate the gradients of transforming the internal (unconstrained) to external (constained) parameter.""" + grad = np.empty_like(p_int) for i, (x, bound) in enumerate(zip(p_int, bounds)): lower, upper = bound if lower is None and upper is None: # No constraints grad[i] = 1.0 elif upper is None: # only lower bound - grad[i] = x/sqrt(x*x + 1.0) + grad[i] = x/np.sqrt(x*x + 1.0) elif lower is None: # only upper bound - grad[i] = -x/sqrt(x*x + 1.0) + grad[i] = -x/np.sqrt(x*x + 1.0) else: # lower and upper bounds - grad[i] = (upper - lower)*cos(x)/2.0 + grad[i] = (upper - lower)*np.cos(x)/2.0 return grad - -def _int2extFunc(bounds): - ''' - transform internal parameters into external parameters. - ''' + + def _int2extFunc(bounds): + """transform internal parameters into external parameters.""" local = [_int2extLocal(b) for b in bounds] def _transform_i2e(p_int): - p_ext = empty_like(p_int) + p_ext = np.empty_like(p_int) p_ext[:] = [i(j) for i, j in zip(local, p_int)] return p_ext return _transform_i2e - -def _ext2intFunc(bounds): - ''' - transform external parameters into internal parameters. - ''' + + def _ext2intFunc(bounds): + """transform external parameters into internal parameters.""" local = [_ext2intLocal(b) for b in bounds] def _transform_e2i(p_ext): - p_int = empty_like(p_ext) + p_int = np.empty_like(p_ext) p_int[:] = [i(j) for i, j in zip(local, p_ext)] return p_int return _transform_e2i - -def _int2extLocal(bound): - ''' - transform a single internal parameter to an external parameter. - ''' + + def _int2extLocal(bound): + """transform a single internal parameter to an external parameter.""" lower, upper = bound if lower is None and upper is None: # no constraints return lambda x: x elif upper is None: # only lower bound - return lambda x: lower - 1.0 + sqrt(x*x + 1.0) + return lambda x: lower - 1.0 + np.sqrt(x*x + 1.0) elif lower is None: # only upper bound - return lambda x: upper + 1.0 - sqrt(x*x + 1.0) + return lambda x: upper + 1.0 - np.sqrt(x*x + 1.0) else: - return lambda x: lower + ((upper - lower)/2.0)*(sin(x) + 1.0) - -def _ext2intLocal(bound): - ''' - transform a single external parameter to an internal parameter. - ''' + return lambda x: lower + ((upper - lower)/2.0)*(np.sin(x) + 1.0) + + def _ext2intLocal(bound): + """transform a single external parameter to an internal parameter.""" lower, upper = bound if lower is None and upper is None: # no constraints return lambda x: x elif upper is None: # only lower bound - return lambda x: sqrt((x - lower + 1.0)**2 - 1.0) + return lambda x: np.sqrt((x - lower + 1.0)**2 - 1.0) elif lower is None: # only upper bound - return lambda x: sqrt((x - upper - 1.0)**2 - 1.0) + return lambda x: np.sqrt((x - upper - 1.0)**2 - 1.0) else: - return lambda x: arcsin((2.0*(x - lower)/(upper - lower)) - 1.0) + return lambda x: np.arcsin((2.0*(x - lower)/(upper - lower)) - 1.0) + + i2e = _int2extFunc(bounds) + e2i = _ext2intFunc(bounds) + + x0 = np.asarray(x0).flatten() + n = len(x0) -def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0, - col_deriv=0, ftol=1.49012e-8, xtol=1.49012e-8, - gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None): - from scipy.optimize import _minpack - ''' - An internal parameter list is used to enforce contraints on the fitting - parameters. The transfomation is based on that of MINUIT package. - please see: F. James and M. Winkler. MINUIT User's Guide, 2004. + if len(bounds) != n: + raise ValueError('the length of bounds is inconsistent with the number of parameters ') + + if not isinstance(args, tuple): + args = (args,) + + shape, dtype = _check_func('leastsq', 'func', func, x0, args, n) + m = shape[0] - bounds : list - (min, max) pairs for each parameter, use None for 'min' or 'max' - when there is no bound in that direction. - For example: if there are two parameters needed to be fitting, then - bounds is [(min1,max1), (min2,max2)] - - This function is based on 'leastsq' of minpack.py, the annotation of - other parameters can be found in 'leastsq'. - ..\Lib\site-packages\scipy\optimize\minpack.py - ''' - i2e = _int2extFunc(bounds) - e2i = _ext2intFunc(bounds) - - x0 = asarray(x0).flatten() - n = len(x0) + if n > m: + raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m)) + if epsfcn is None: + epsfcn = np.finfo(dtype).eps - if len(bounds) != n: - raise ValueError('the length of bounds is inconsistent with the number of parameters ') - - if not isinstance(args, tuple): - args = (args,) - - shape, dtype = _check_func('leastsq', 'func', func, x0, args, n) - m = shape[0] + def funcWarp(x, *args): + return func(i2e(x), *args) - if n > m: - raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m)) - if epsfcn is None: - epsfcn = finfo(dtype).eps + xi0 = e2i(x0) + + if Dfun is None: + if maxfev == 0: + maxfev = 200*(n + 1) + retval = _minpack._lmdif(funcWarp, xi0, args, full_output, ftol, xtol, + gtol, maxfev, epsfcn, factor, diag) + else: + if col_deriv: + _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (n, m)) + else: + _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (m, n)) + if maxfev == 0: + maxfev = 100*(n + 1) - # wrapped func - def funcWarp(x, *args): - return func(i2e(x), *args) + def DfunWarp(x, *args): + return Dfun(i2e(x), *args) - xi0 = e2i(x0) - - if Dfun is None: - if maxfev == 0: - maxfev = 200*(n + 1) - retval = _minpack._lmdif(funcWarp, xi0, args, full_output, ftol, xtol, - gtol, maxfev, epsfcn, factor, diag) - else: - if col_deriv: - _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (n, m)) - else: - _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (m, n)) - if maxfev == 0: - maxfev = 100*(n + 1) + retval = _minpack._lmder(funcWarp, DfunWarp, xi0, args, full_output, col_deriv, + ftol, xtol, gtol, maxfev, factor, diag) - # wrapped Dfun - def DfunWarp(x, *args): - return Dfun(i2e(x), *args) + errors = {0: ["Improper input parameters.", TypeError], + 1: ["Both actual and predicted relative reductions " + "in the sum of squares\n are at most %f" % ftol, None], + 2: ["The relative error between two consecutive " + "iterates is at most %f" % xtol, None], + 3: ["Both actual and predicted relative reductions in " + "the sum of squares\n are at most %f and the " + "relative error between two consecutive " + "iterates is at \n most %f" % (ftol, xtol), None], + 4: ["The cosine of the angle between func(x) and any " + "column of the\n Jacobian is at most %f in " + "absolute value" % gtol, None], + 5: ["Number of calls to function has reached " + "maxfev = %d." % maxfev, ValueError], + 6: ["ftol=%f is too small, no further reduction " + "in the sum of squares\n is possible.""" % ftol, + ValueError], + 7: ["xtol=%f is too small, no further improvement in " + "the approximate\n solution is possible." % xtol, + ValueError], + 8: ["gtol=%f is too small, func(x) is orthogonal to the " + "columns of\n the Jacobian to machine " + "precision." % gtol, ValueError], + 'unknown': ["Unknown error.", TypeError]} - retval = _minpack._lmder(funcWarp, DfunWarp, xi0, args, full_output, col_deriv, - ftol, xtol, gtol, maxfev, factor, diag) + info = retval[-1] # The FORTRAN return value + + if info not in [1, 2, 3, 4] and not full_output: + if info in [5, 6, 7, 8]: + np.warnings.warn(errors[info][0], RuntimeWarning) + else: + try: + raise errors[info][1](errors[info][0]) + except KeyError: + raise errors['unknown'][1](errors['unknown'][0]) - errors = {0: ["Improper input parameters.", TypeError], - 1: ["Both actual and predicted relative reductions " - "in the sum of squares\n are at most %f" % ftol, None], - 2: ["The relative error between two consecutive " - "iterates is at most %f" % xtol, None], - 3: ["Both actual and predicted relative reductions in " - "the sum of squares\n are at most %f and the " - "relative error between two consecutive " - "iterates is at \n most %f" % (ftol, xtol), None], - 4: ["The cosine of the angle between func(x) and any " - "column of the\n Jacobian is at most %f in " - "absolute value" % gtol, None], - 5: ["Number of calls to function has reached " - "maxfev = %d." % maxfev, ValueError], - 6: ["ftol=%f is too small, no further reduction " - "in the sum of squares\n is possible.""" % ftol, - ValueError], - 7: ["xtol=%f is too small, no further improvement in " - "the approximate\n solution is possible." % xtol, - ValueError], - 8: ["gtol=%f is too small, func(x) is orthogonal to the " - "columns of\n the Jacobian to machine " - "precision." % gtol, ValueError], - 'unknown': ["Unknown error.", TypeError]} + mesg = errors[info][0] + x = i2e(retval[0]) - info = retval[-1] # The FORTRAN return value - - if info not in [1, 2, 3, 4] and not full_output: - if info in [5, 6, 7, 8]: - warnings.warn(errors[info][0], RuntimeWarning) - else: - try: - raise errors[info][1](errors[info][0]) - except KeyError: - raise errors['unknown'][1](errors['unknown'][0]) - - mesg = errors[info][0] - x = i2e(retval[0]) - - if full_output: - grad = _int2extGrad(retval[0], bounds) - retval[1]['fjac'] = (retval[1]['fjac'].T / take(grad, - retval[1]['ipvt'] - 1)).T - cov_x = None - if info in [1, 2, 3, 4]: - from numpy.dual import inv - from numpy.linalg import LinAlgError - perm = take(eye(n), retval[1]['ipvt'] - 1, 0) - r = triu(transpose(retval[1]['fjac'])[:n, :]) - R = dot(r, perm) - try: - cov_x = inv(dot(transpose(R), R)) - except LinAlgError as inverror: - print inverror - pass - return (x, cov_x) + retval[1:-1] + (mesg, info) - else: - return (x, info) + if full_output: + grad = _int2extGrad(retval[0], bounds) + retval[1]['fjac'] = (retval[1]['fjac'].T / np.take(grad, + retval[1]['ipvt'] - 1)).T + cov_x = None + if info in [1, 2, 3, 4]: + from numpy.dual import inv + from numpy.linalg import LinAlgError + perm = np.take(np.eye(n), retval[1]['ipvt'] - 1, 0) + r = np.triu(np.transpose(retval[1]['fjac'])[:n, :]) + R = np.dot(r, perm) + try: + cov_x = inv(np.dot(np.transpose(R), R)) + except LinAlgError as inverror: + print inverror + pass + return (x, cov_x) + retval[1:-1] + (mesg, info) + else: + return (x, info) def _general_function(params, ydata, xdata, function): return function(xdata, *params) - ydata @@ -405,7 +393,7 @@ def _weighted_general_function(params, ydata, xdata, function, weights): return (function(xdata, *params) - ydata)*weights def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): - ''' Similar as 'curve_fit' in minpack.py''' + """Similar as 'curve_fit' in minpack.py""" if p0 is None: # determine number of parameters by inspecting the function import inspect @@ -418,15 +406,15 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): else: p0 = [1.0] * (len(args)-1) - if isscalar(p0): - p0 = array([p0]) + if np.isscalar(p0): + p0 = np.array([p0]) args = (ydata, xdata, f) if sigma is None: func = _general_function else: func = _weighted_general_function - args += (1.0/asarray(sigma),) + args += (1.0/np.asarray(sigma),) return_full = kw.pop('full_output', False) res = leastsqBound(func, p0, args=args, bounds = bounds, full_output=True, **kw) @@ -440,7 +428,7 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0)) pcov = pcov * s_sq else: - pcov = inf + pcov = np.inf if return_full: return popt, pcov, infodict, errmsg, ier @@ -449,13 +437,11 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): def execute(cmd,streamIn=None,wd='./'): - ''' - executes a command in given directory and returns stdout and stderr for optional stdin - ''' + """executes a command in given directory and returns stdout and stderr for optional stdin""" initialPath=os.getcwd() os.chdir(wd) process = subprocess.Popen(shlex.split(cmd),stdout=subprocess.PIPE,stderr = subprocess.PIPE,stdin=subprocess.PIPE) - if streamIn != None: + if streamIn is not None: out,error = process.communicate(streamIn.read()) else: out,error = process.communicate()