223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

modified background threading because using it in postResults led to an dead lock. 

rest modifications are just style
This commit is contained in:
Martin Diehl 2016-03-04 15:22:01 +01:00
parent 48233d2767
commit 75480bc677
1 changed files with 200 additions and 214 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

414
lib/damask/util.py
View File

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