def check_finite_real(M):
'''
Check that all entries in array M are finite and real-valued
'''
if np.any(~np.isreal(M)):
raise ValueError("Complex value encountered for real vector")
if np.any(~np.isfinite(M)):
raise ValueError("Non-finite number encountered")
# need a faster covariance matrix checker
python类isreal()的实例源码
def check_covmat_fast(C,N=None,eps=1e-6):
'''
Verify that matrix M is a size NxN precision or covariance matirx
'''
if not type(C)==np.ndarray:
raise ValueError("Covariance matrix should be a 2D numpy array")
if not len(C.shape)==2:
raise ValueError("Covariance matrix should be a 2D numpy array")
if N is None:
N = C.shape[0]
if not C.shape==(N,N):
raise ValueError("Expected size %d x %d matrix"%(N,N))
if np.any(~np.isreal(C)):
raise ValueError("Covariance matrices should not contain complex numbers")
C = np.real(C)
if np.any(~np.isfinite(C)):
raise ValueError("Covariance matrix contains NaN or ±inf!")
if not np.all(np.abs(C-C.T)<eps):
raise ValueError("Covariance matrix is not symmetric up to precision %0.1e"%eps)
try:
ch = chol(C)
except numpy.linalg.linalg.LinAlgError:
# Check smallest eigenvalue if cholesky fails
mine = np.real(scipy.linalg.decomp.eigh(C,eigvals=(0,0))[0][0])
if np.any(mine<-eps):
raise ValueError('Covariance matrix contains eigenvalue %0.3e<%0.3e'%(mine,-eps))
if (mine<eps):
C = C + np.eye(N)*(eps-mine)
C = 0.5*(C+C.T)
return C
def _strekklikevekt(self, L, G_sno, T, H_0=None):
"""Finner kabelstrekk under gitte forhold for fastavspent ledning.
Følgende likevektsligning ligger til grunn for beregningene:
:math:`H_x^2 [H_x - H_0 + \\frac{EA(G_0 L)^2}{24H_0^2} + EA\\alpha \\Delta_T]
= \\frac{EA(G_x L)^2}{24}`
Ligningen for kabellikevekten er hentet fra «KL-bibelen» («Contact Lines for
Electric Railways» av Kiessling, Puschmann etc.) ligning (5.57) side 282.
Løsningen finnes ved å finne den reelle, positive egenverdien
tilhørende "companion matrix" for residualfunksjonens koeffisienter.
:param float H_0: Initiell spennkraft i kabel :math:`[N]`
:param float E: Kabelens E-modul :math:`[\\frac{N}{mm^2}]`
:param float A: Kabelens tverrsnittsareal :math:`[mm^2]`
:param float G_0: Kabelens egenvekt :math:`[\\frac{N}{m}]`
:param float G_sno: Egenvekt snølast :math:`[\\frac{N}{m}]`
:param float L: Masteavstand :math:`[m]`
:param float alpha: Lengdeutvidelseskoeffisient :math:`[\\frac{1}{^{\\circ}C}]`
:param float T: Lufttemperatur :math:`[^{\\circ}C]`
:return: Endelig kabelstrekk ``H_x`` :math:`[N]`
:rtype: :class:`float`
"""
# Inngangsparametre
if H_0 is not None:
H_0 = H_0
else:
H_0 = self.temperaturdata["5C"]["s"]
E = self.E
A = self.A
G_0 = self.G_0
alpha = self.alpha
G_x = G_0 + G_sno
delta_T = T - 5
# Konstanter
a = E * A * (G_x * L) ** 2 / 24
b = - H_0 + E * A * (G_0 * L) ** 2 / (24 * H_0 ** 2) + E * A * alpha * delta_T
roots = numpy.roots([-1, -b, 0, a])
H_x = 0
for r in roots:
if numpy.isreal(r) and r > 0:
H_x = numpy.real(r)
break
return H_x
def _make_widget(self):
self.widget = pg.GraphicsWindow(title="Curve")
self.plot_item = self.widget.addPlot(title="Curve")
self.plot_item_phase = self.widget.addPlot(row=1, col=0, title="Phase (deg)")
self.plot_item_phase.setXLink(self.plot_item)
self.plot_item.showGrid(y=True, alpha=1.)
self.plot_item_phase.showGrid(y=True, alpha=1.)
self.curve = self.plot_item.plot(pen='g')
self.curve_phase = self.plot_item_phase.plot(pen='g')
self._is_real = True
self._set_real(True)
#def _set_widget_value(self, new_value):
# data = new_value
# if data is None:
# return
# shape = np.shape(new_value)
# if len(shape)>2:
# raise ValueError("Shape of data should be (1) or (2, 1)")
# if len(shape)==1:
# x = np.linspace(0, len(data), len(data))
# y = [data]
# if len(shape)==2:
# if shape[0] == 1:
# x = np.linspace(0, len(data), len(data[0]))
# y = [data[0]]
# if shape[0] >= 2:
# x = data[0]
# y = data[1:]
# self._set_real(np.isreal(y).all())
# for i, values in enumerate(y):
# self._display_curve_index(x, values, i)
# while (i + 1 < len(self.curves)): # delete remaining curves
# i += 1
# self.curves[i].hide()
#def _display_curve_index(self, x, values, i):
# y_mag = values if self._is_real else self._magnitude(values)
# y_phase = np.zeros(len(values)) if self._is_real else \
# self._phase(values)
# if len(self.curves)<=i:
# color = self._defaultcolors[i%len(self._defaultcolors)]
# self.curves.append(self.plot_item.plot(pen=color))
# self.curves_phase.append(self.plot_item_phase.plot(pen=color))
# self.curves[i].setData(x, y_mag)
# self.curves_phase[i].setData(x, y_phase)
def __init__(self, value, shape=None, dtype=None):
"""
Create a new lazy array.
`value` : may be an int, long, float, bool, NumPy array, iterator,
generator or a function, `f(i)` or `f(i,j)`, depending on the
dimensions of the array.
`f(i,j)` should return a single number when `i` and `j` are integers,
and a 1D array when either `i` or `j` or both is a NumPy array (in the
latter case the two arrays must have equal lengths).
"""
self.dtype = dtype
self.operations = []
if isinstance(value, basestring):
raise TypeError("An larray cannot be created from a string")
elif isinstance(value, larray):
if shape is not None and value.shape is not None:
assert shape == value.shape
self._shape = shape or value.shape
self.base_value = value.base_value
self.dtype = dtype or value.dtype
self.operations = value.operations # should deepcopy?
elif isinstance(value, collections.Sized): # False for numbers, generators, functions, iterators
if have_scipy and sparse.issparse(value): # For sparse matrices
self.dtype = dtype or value.dtype
elif not isinstance(value, numpy.ndarray):
value = numpy.array(value, dtype=dtype)
elif dtype is not None:
assert numpy.can_cast(value.dtype, dtype, casting='safe') # or could convert value to the provided dtype
if shape and value.shape != shape:
raise ValueError("Array has shape %s, value has shape %s" % (shape, value.shape))
self._shape = value.shape
self.base_value = value
else:
assert numpy.isreal(value) # also True for callables, generators, iterators
self._shape = shape
if dtype is None or isinstance(value, dtype):
self.base_value = value
else:
try:
self.base_value = dtype(value)
except TypeError:
self.base_value = value
def get_PGM_from_numpy_arr(arr, window_center, window_width,
lut_min=0, lut_max=255):
'''real-valued numpy input -> PGM-image formatted byte string
arr: real-valued numpy array to display as grayscale image
window_center, window_width: to define max/min values to be mapped to the
lookup-table range. WC/WW scaling is done
according to DICOM-3 specifications.
lut_min, lut_max: min/max values of (PGM-) grayscale table: do not change
'''
if np.isreal(arr).sum() != arr.size:
raise ValueError
# currently only support 8-bit colors
if lut_max != 255:
raise ValueError
if arr.dtype != np.float64:
arr = arr.astype(np.float64)
# LUT-specific array scaling
# width >= 1 (DICOM standard)
window_width = max(1, window_width)
wc, ww = np.float64(window_center), np.float64(window_width)
lut_range = np.float64(lut_max) - lut_min
minval = wc - 0.5 - (ww - 1.0) / 2.0
maxval = wc - 0.5 + (ww - 1.0) / 2.0
min_mask = (minval >= arr)
to_scale = (arr > minval) & (arr < maxval)
max_mask = (arr >= maxval)
if min_mask.any():
arr[min_mask] = lut_min
if to_scale.any():
arr[to_scale] = ((arr[to_scale] - (wc - 0.5)) /
(ww - 1.0) + 0.5) * lut_range + lut_min
if max_mask.any():
arr[max_mask] = lut_max
# round to next integer values and convert to unsigned int
arr = np.rint(arr).astype(np.uint8)
# return PGM byte-data string
return get_PGM_bytedata_string(arr)
def get_PGM_from_numpy_arr(arr, window_center, window_width,
lut_min=0, lut_max=255):
'''real-valued numpy input -> PGM-image formatted byte string
arr: real-valued numpy array to display as grayscale image
window_center, window_width: to define max/min values to be mapped to the
lookup-table range. WC/WW scaling is done
according to DICOM-3 specifications.
lut_min, lut_max: min/max values of (PGM-) grayscale table: do not change
'''
if np.isreal(arr).sum() != arr.size:
raise ValueError
# currently only support 8-bit colors
if lut_max != 255:
raise ValueError
if arr.dtype != np.float64:
arr = arr.astype(np.float64)
# LUT-specific array scaling
# width >= 1 (DICOM standard)
window_width = max(1, window_width)
wc, ww = np.float64(window_center), np.float64(window_width)
lut_range = np.float64(lut_max) - lut_min
minval = wc - 0.5 - (ww - 1.0) / 2.0
maxval = wc - 0.5 + (ww - 1.0) / 2.0
min_mask = (minval >= arr)
to_scale = (arr > minval) & (arr < maxval)
max_mask = (arr >= maxval)
if min_mask.any():
arr[min_mask] = lut_min
if to_scale.any():
arr[to_scale] = ((arr[to_scale] - (wc - 0.5)) /
(ww - 1.0) + 0.5) * lut_range + lut_min
if max_mask.any():
arr[max_mask] = lut_max
# round to next integer values and convert to unsigned int
arr = np.rint(arr).astype(np.uint8)
# return PGM byte-data string
return get_PGM_bytedata_string(arr)
def check_covmat(C,N=None,eps=1e-6):
'''
Verify that matrix M is a size NxN precision or covariance matirx
'''
if not type(C)==np.ndarray:
raise ValueError("Covariance matrix should be a 2D numpy array")
if not len(C.shape)==2:
raise ValueError("Covariance matrix should be a 2D numpy array")
if N is None:
N = C.shape[0]
if not C.shape==(N,N):
raise ValueError("Expected size %d x %d matrix"%(N,N))
if np.any(~np.isreal(C)):
raise ValueError("Covariance matrices should not contain complex numbers")
C = np.real(C)
if np.any(~np.isfinite(C)):
raise ValueError("Covariance matrix contains NaN or ±inf!")
if not np.all(np.abs(C-C.T)<eps):
raise ValueError("Covariance matrix is not symmetric up to precision %0.1e"%eps)
# Get just highest eigenvalue
maxe = np.real(scipy.linalg.decomp.eigh(C,eigvals=(N-1,N-1))[0][0])
# Get just lowest eigenvalue
mine = np.real(scipy.linalg.decomp.eigh(C,eigvals=(0,0))[0][0])
#if np.any(w<-eps):
# raise ValueError('Covariance matrix contains eigenvalue %0.3e<%0.3e'%(np.min(w),-eps))
if mine<0:
raise ValueError('Covariance matrix contains negative eigenvalue %0.3e'%mine)
if (mine<eps):
C = C + eye(N)*(eps-mine)
# trucate spectrum at some small value
# w[w<eps]=eps
# Very large eigenvalues can also cause numeric problems
# w[w>1./eps]=1./eps;
# maxe = np.max(np.abs(w))
# if maxe>10./eps:
# raise ValueError('Covariance matrix eigenvalues %0.2e is larger than %0.2e'%(maxe,10./eps))
# Rebuild matrix
# C = v.dot(np.diag(w)).dot(v.T)
# Ensure symmetry (only occurs as a numerical error for very large matrices?)
C = 0.5*(C+C.T)
return C
# need a faster covariance matrix checker
