def normalize_array (solution, prediction):
''' Use min and max of solution as scaling factors to normalize prediction,
then threshold it to [0, 1]. Binarize solution to {0, 1}.
This allows applying classification scores to all cases.
In principle, this should not do anything to properly formatted
classification inputs and outputs.'''
# Binarize solution
sol=np.ravel(solution) # convert to 1-d array
maxi = np.nanmax((filter(lambda x: x != float('inf'), sol))) # Max except NaN and Inf
mini = np.nanmin((filter(lambda x: x != float('-inf'), sol))) # Mini except NaN and Inf
if maxi == mini:
print('Warning, cannot normalize')
return [solution, prediction]
diff = maxi - mini
mid = (maxi + mini)/2.
new_solution = np.copy(solution)
new_solution[solution>=mid] = 1
new_solution[solution<mid] = 0
# Normalize and threshold predictions (takes effect only if solution not in {0, 1})
new_prediction = (np.copy(prediction) - float(mini))/float(diff)
new_prediction[new_prediction>1] = 1 # and if predictions exceed the bounds [0, 1]
new_prediction[new_prediction<0] = 0
# Make probabilities smoother
#new_prediction = np.power(new_prediction, (1./10))
return [new_solution, new_prediction]
python类power()的实例源码
def derivative(self, input=None):
"""The derivative of :meth:`tanh` functions is
.. math:: \\frac{d}{dx} tanh(x) & = \\frac{d}{dx} \\frac{sinh(x)}{cosh(x)} \\\\
& = \\frac{cosh(x) \\frac{d}{dx}sinh(x) - sinh(x) \\frac{d}{dx}cosh(x) }{ cosh^2(x)} \\\\
& = \\frac{ cosh(x) cosh(x) - sinh(x) sinh(x) }{ cosh^2(x)} \\\\
& = 1 - tanh^2(x)
Returns
-------
float32
The derivative of tanh function.
"""
last_forward = self.forward(input) if input else self.last_forward
return 1 - np.power(last_forward, 2)
# tanh-end
# relu-start
def update(self, params, grads):
# init cache and delta
if self.cache is None:
self.cache = [_zero(p.shape) for p in params]
if self.delta is None:
self.delta = [_zero(p.shape) for p in params]
# update parameters
for i, (c, d, p, g) in enumerate(zip(self.cache, self.delta, params, grads)):
c = self.rho * c + (1 - self.rho) * np.power(g, 2)
update = g * np.sqrt(d + self.epsilon) / np.sqrt(c + self.epsilon)
p -= self.lr * update
d = self.rho * d + (1 - self.rho) * np.power(update, 2)
self.cache[i] = c
self.delta[i] = d
def update(self, params, grads):
# init
self.iterations += 1
a_t = self.lr * np.sqrt(1 - np.power(self.beta2, self.iterations)) / \
(1 - np.power(self.beta1, self.iterations))
if self.ms is None:
self.ms = [_zero(p.shape) for p in params]
if self.vs is None:
self.vs = [_zero(p.shape) for p in params]
# update parameters
for i, (m, v, p, g) in enumerate(zip(self.ms, self.vs, params, grads)):
m = self.beta1 * m + (1 - self.beta1) * g
v = self.beta2 * v + (1 - self.beta2) * np.power(g, 2)
p -= a_t * m / (np.sqrt(v) + self.epsilon)
self.ms[i] = m
self.vs[i] = v
def update(self, params, grads):
# init
self.iterations += 1
a_t = self.lr / (1 - np.power(self.beta1, self.iterations))
if self.ms is None:
self.ms = [_zero(p.shape) for p in params]
if self.vs is None:
self.vs = [_zero(p.shape) for p in params]
# update parameters
for i, (m, v, p, g) in enumerate(zip(self.ms, self.vs, params, grads)):
m = self.beta1 * m + (1 - self.beta1) * g
v = np.maximum(self.beta2 * v, np.abs(g))
p -= a_t * m / (v + self.epsilon)
self.ms[i] = m
self.vs[i] = v
def forward(self, outputs, targets):
"""MeanSquaredError forward propagation.
.. math:: L = (p - t)^2
Parameters
----------
outputs, targets : numpy.array
The arrays to compute the squared difference between.
Returns
-------
numpy.array
An expression for the element-wise squared difference.
"""
return 0.5 * np.mean(np.sum(np.power(outputs - targets, 2), axis=1))
def forward(self, outputs, targets):
"""HellingerDistance forward propagation.
Parameters
----------
outputs : numpy 2D array
outputs in (0, 1), such as softmax output of a neural network,
with data points in rows and class probabilities in columns.
targets : numpy 2D array
Either a vector of int giving the correct class index per data point
or a 2D tensor of one-hot encoding of the correct class in the same
layout as predictions (non-binary targets in [0, 1] do not work!)
Returns
-------
numpy 1D array
An expression for the Hellinger Distance
"""
root_difference = np.sqrt(outputs) - np.sqrt(targets)
return np.mean(np.sum(np.power(root_difference, 2), axis=1) / np.sqrt(2))
def evaluation(self, X_test, y_test):
# normalization
X_test = self.normalization(X_test)
# average over the output
pred_y_test = np.zeros([self.M, len(y_test)])
prob = np.zeros([self.M, len(y_test)])
'''
Since we have M particles, we use a Bayesian view to calculate rmse and log-likelihood
'''
for i in range(self.M):
w1, b1, w2, b2, loggamma, loglambda = self.unpack_weights(self.theta[i, :])
pred_y_test[i, :] = self.nn_predict(X_test, w1, b1, w2, b2) * self.std_y_train + self.mean_y_train
prob[i, :] = np.sqrt(np.exp(loggamma)) /np.sqrt(2*np.pi) * np.exp( -1 * (np.power(pred_y_test[i, :] - y_test, 2) / 2) * np.exp(loggamma) )
pred = np.mean(pred_y_test, axis=0)
# evaluation
svgd_rmse = np.sqrt(np.mean((pred - y_test)**2))
svgd_ll = np.mean(np.log(np.mean(prob, axis = 0)))
return (svgd_rmse, svgd_ll)
def plot_beta():
'''plot beta over training
'''
beta = args.beta
scale = args.scale
beta_min = args.beta_min
num_epoch = args.num_epoch
epoch_size = int(float(args.num_examples) / args.batch_size)
x = np.arange(num_epoch*epoch_size)
y = beta * np.power(scale, x)
y = np.maximum(y, beta_min)
epoch_x = np.arange(num_epoch) * epoch_size
epoch_y = beta * np.power(scale, epoch_x)
epoch_y = np.maximum(epoch_y, beta_min)
# plot beta descent curve
plt.semilogy(x, y)
plt.semilogy(epoch_x, epoch_y, 'ro')
plt.title('beta descent')
plt.ylabel('beta')
plt.xlabel('epoch')
plt.show()
def BM_Eval(seq_dict, BMlist, toeholds):
w_exp = np.concatenate([np.zeros((5,)), np.power(2, np.arange(6))])
BM_score = 0
Largest_match = 0
numstrings = len(BMlist);
prog = MyProgress((numstrings**2 - numstrings)/2)
for ctr in range(numstrings):
strand1 = BMlist[ctr];
for strand2 in BMlist[ctr+1:]:
[ismaxmatch, maxmatch, mm_i, mm_j] = \
compare_sequence_notoe(seq_dict[strand1],
seq_dict[strand2],
toeholds)
if maxmatch > Largest_match:
Largest_match = maxmatch
BM_score = BM_score + w_exp[int(min(maxmatch, 10))]
prog.inc()
return [BM_score, Largest_match]
def computeCost(X, y, theta):
inner = np.power(((X * theta.T) - y), 2)
return np.sum(inner) / (2 * len(X))
#def gradientDescent(X, y, theta, alpha, iters):
# temp = np.matrix(np.zeros(theta.shape))
# params = int(theta.ravel().shape[1]) #flattens
# cost = np.zeros(iters)
#
# for i in range(iters):
# err = (X * theta.T) - y
#
# for j in range(params):
# term = np.multiply(err, X[:,j])
# temp[0, j] = theta[0, j] - ((alpha / len(X)) * np.sum(term))
#
# theta = temp
# cost[i] = computeCost(X, y, theta)
#
# return theta, cost
def get_genotype_probability(aln_profile, aln_specificity, sigma=0.12):
# 'aln_specificity' should be a set of unit vectors (at least one of the entry is larger than 1.)
num_haps = len(aln_profile)
aln_vec = unit_vector(aln_profile)
genoprob = []
for i in xrange(num_haps):
v1 = unit_vector(aln_specificity[i])
for j in xrange(i, num_haps):
if j == i:
genoprob.append(sum(np.power(aln_vec - v1, 2))) # homozygotes
else:
v2 = unit_vector(aln_specificity[j])
geno_vec = unit_vector(v1 + v2)
# compute directional similarity
genoprob.append(sum(np.power(aln_vec - geno_vec, 2))) # for heterozygotes
genoprob = np.exp(np.array(genoprob) / (-2 * sigma * sigma))
return np.array(genoprob / sum(genoprob))
def set_params(mo, bparams):
i = 0
for la in mo.layers:
we = bparams[i:i+2]
print len(we)
la.set_weights(we)
i += 2
return mo
#with open("best_model_keras.pkl", 'r') as f:
# b_params = pkl.load(f)
#
#model = set_params(model, b_params)
#out = model.predict(xvl, batch_size=xvl.shape[0], verbose=0)
#error = np.mean(np.mean(np.power(out - yvl, 2), axis=1))
#print "Error vl", error
#sys.exit()
#init_p = get_params(model)
#with open("init_keras_param.pkl", 'w') as f:
# pkl.dump(init_p, f)
def CompLikelihood(X,fx,MCPar,Measurement,Extra):
Sigma=Measurement.Sigma*np.ones((X.shape[0]))
of=np.zeros((fx.shape[0],1))
p=np.zeros((fx.shape[0],1))
log_p=np.zeros((fx.shape[0],1))
for ii in xrange(0,fx.shape[0]):
e=Measurement.MeasData-fx[ii,:]
of[ii,0]=np.sqrt(np.sum(np.power(e,2.0))/e.shape[1])
if MCPar.lik==2: # Compute standard uncorrelated and homoscedastic Gaussian log-likelihood
log_p[ii,0]= - ( Measurement.N / 2.0) * np.log(2.0 * np.pi) - Measurement.N * np.log( Sigma[ii] ) - 0.5 * np.power(Sigma[ii],-2.0) * np.sum( np.power(e,2.0) )
p[ii,0]=(1.0/np.sqrt(2*np.pi* Sigma[ii]**2))**Measurement.N * np.exp(- 0.5 * np.power(Sigma[ii],-2.0) * np.sum( np.power(e,2.0) ))
if MCPar.lik==3: # Box and Tiao (1973) log-likelihood formulation with Sigma integrated out based on prior of the form p(sigma) ~ 1/sigma
log_p[ii,0]= - ( Measurement.N / 2.0) * np.log(np.sum(np.power(e,2.0)))
p[ii,0]=np.exp(log_p[ii,0])
return of, p, log_p
def SLcomputeSNR(X, Xnoisy):
"""
SLcomputeSNR Compute signal to noise ratio (SNR).
Usage:
SNR = SLcomputeSNR(X, Xnoisy)
Input:
X: 2D or 3D signal.
Xnoisy: 2D or 3D noisy signal.
Output:
SNR: The signal to noise ratio (in dB).
"""
if np.linalg.norm(X-Xnoisy) == 0:
return np.Inf
else:
return 10 * np.log10( np.sum(np.power(X,2)) / np.sum(np.power(X-Xnoisy,2)) )
def __init__(self, path, normalize=True, eig=0.0, transpose=False):
if transpose:
ut = np.load(path + '.vt.npy')
self.wi, self.iw = load_vocabulary(path + '.contexts.vocab')
else:
ut = np.load(path + '.ut.npy')
self.wi, self.iw = load_vocabulary(path + '.words.vocab')
s = np.load(path + '.s.npy')
if eig == 0.0:
self.m = ut.T
elif eig == 1.0:
self.m = s * ut.T
else:
self.m = np.power(s, eig) * ut.T
self.dim = self.m.shape[1]
if normalize:
self.normalize()
def run(count_path, out_path, smooth=0, cds=True, normalize=False, neg=1):
counts = create_representation("Explicit", count_path, normalize=False)
old_mat = counts.m
index = counts.wi
smooth = old_mat.sum() * smooth
# getting marginal probs
row_probs = old_mat.sum(1) + smooth
col_probs = old_mat.sum(0) + smooth
if cds:
col_probs = np.power(col_probs, 0.75)
row_probs = row_probs / row_probs.sum()
col_probs = col_probs / col_probs.sum()
# building PPMI matrix
ppmi_mat = make_ppmi_mat(old_mat, row_probs, col_probs, smooth, neg=neg, normalize=normalize)
import pyximport
pyximport.install(setup_args={"include_dirs": np.get_include()})
from representations import sparse_io
sparse_io.export_mat_eff(ppmi_mat.row, ppmi_mat.col, ppmi_mat.data, out_path + ".bin")
util.write_pickle(index, out_path + "-index.pkl")
def __init__(self, path, normalize=True, eig=0.0, **kwargs):
ut = np.load(path + '-u.npy', mmap_mode="c")
s = np.load(path + '-s.npy', mmap_mode="c")
vocabfile = path + '-vocab.pkl'
self.iw = load_pickle(vocabfile)
self.wi = {w:i for i, w in enumerate(self.iw)}
if eig == 0.0:
self.m = ut
elif eig == 1.0:
self.m = s * ut
else:
self.m = np.power(s, eig) * ut
self.dim = self.m.shape[1]
if normalize:
self.normalize()
def dispersion_test(yhat, y, k=100):
""" Implement the regression based dispersion test with k re-sampling.
Args:
yhat (np.array): predicted mutation count
y (np.array): observed mutation count
k (int):
Returns:
float, float: p-value, theta
"""
theta = 0
pval = 0
for i in range(k):
y_sub, yhat_sub = resample(y, yhat, random_state=i)
# (np.power((y - yhat), 2) - y) / yhat for Poisson regression
aux = (np.power((y_sub - yhat_sub), 2) - yhat_sub) / yhat_sub
mod = sm.OLS(aux, yhat_sub)
res = mod.fit()
theta += res.params[0]
pval += res.pvalues[0]
theta = theta/k
pval = pval/k
return pval, theta
def _transform_y(y, lam):
"""Transform a single y, given a single lambda value.
No validation performed.
Parameters
----------
y : array_like, shape (n_samples,)
The vector being transformed
lam : ndarray, shape (n_lambdas,)
The lambda value used for the transformation
"""
# ensure np array
y = np.array(y)
y_prime = np.array([(np.power(x, lam) - 1) / lam if not _eqls(lam, ZERO) else log(x) for x in y])
# rarely -- very rarely -- we can get a NaN. Why?
return y_prime
def _yj_trans_single_x(x, lam):
if x >= 0:
# Case 1: x >= 0 and lambda is not 0
if not _eqls(lam, ZERO):
return (np.power(x + 1, lam) - 1.0) / lam
# Case 2: x >= 0 and lambda is zero
return log(x + 1)
else:
# Case 2: x < 0 and lambda is not two
if not lam == 2.0:
denom = 2.0 - lam
numer = np.power((-x + 1), (2.0 - lam)) - 1.0
return -numer / denom
# Case 4: x < 0 and lambda is two
return -log(-x + 1)
def std(files, batch_size=128):
s = np.zeros(3)
s2 = np.zeros(3)
shape = None
for i in range(0, len(files), batch_size):
print("done with {:>3} / {} images".format(i, len(files)))
images = np.array(data.load_image(files[i : i + batch_size]),
dtype=np.float64)
shape = images.shape
s += images.sum(axis=(0, 2, 3))
s2 += np.power(images, 2).sum(axis=(0, 2, 3))
n = len(files) * shape[2] * shape[3]
var = (s2 - s**2.0 / n) / (n - 1)
print('mean')
print((s / n).astype(np.float32))
print('std')
print(np.sqrt(var))
#return np.sqrt(var)
def test_NotImplemented_not_returned(self):
# See gh-5964 and gh-2091. Some of these functions are not operator
# related and were fixed for other reasons in the past.
binary_funcs = [
np.power, np.add, np.subtract, np.multiply, np.divide,
np.true_divide, np.floor_divide, np.bitwise_and, np.bitwise_or,
np.bitwise_xor, np.left_shift, np.right_shift, np.fmax,
np.fmin, np.fmod, np.hypot, np.logaddexp, np.logaddexp2,
np.logical_and, np.logical_or, np.logical_xor, np.maximum,
np.minimum, np.mod
]
# These functions still return NotImplemented. Will be fixed in
# future.
# bad = [np.greater, np.greater_equal, np.less, np.less_equal, np.not_equal]
a = np.array('1')
b = 1
for f in binary_funcs:
assert_raises(TypeError, f, a, b)
def test_half_coercion(self):
"""Test that half gets coerced properly with the other types"""
a16 = np.array((1,), dtype=float16)
a32 = np.array((1,), dtype=float32)
b16 = float16(1)
b32 = float32(1)
assert_equal(np.power(a16, 2).dtype, float16)
assert_equal(np.power(a16, 2.0).dtype, float16)
assert_equal(np.power(a16, b16).dtype, float16)
assert_equal(np.power(a16, b32).dtype, float16)
assert_equal(np.power(a16, a16).dtype, float16)
assert_equal(np.power(a16, a32).dtype, float32)
assert_equal(np.power(b16, 2).dtype, float64)
assert_equal(np.power(b16, 2.0).dtype, float64)
assert_equal(np.power(b16, b16).dtype, float16)
assert_equal(np.power(b16, b32).dtype, float32)
assert_equal(np.power(b16, a16).dtype, float16)
assert_equal(np.power(b16, a32).dtype, float32)
assert_equal(np.power(a32, a16).dtype, float32)
assert_equal(np.power(a32, b16).dtype, float32)
assert_equal(np.power(b32, a16).dtype, float16)
assert_equal(np.power(b32, b16).dtype, float32)
def __ipow__(self, other):
"""
Raise self to the power other, in place.
"""
other_data = getdata(other)
other_mask = getmask(other)
with np.errstate(divide='ignore', invalid='ignore'):
self._data.__ipow__(np.where(self._mask, self.dtype.type(1),
other_data))
invalid = np.logical_not(np.isfinite(self._data))
if invalid.any():
if self._mask is not nomask:
self._mask |= invalid
else:
self._mask = invalid
np.copyto(self._data, self.fill_value, where=invalid)
new_mask = mask_or(other_mask, invalid)
self._mask = mask_or(self._mask, new_mask)
return self
def normalize_array (solution, prediction):
''' Use min and max of solution as scaling factors to normalize prediction,
then threshold it to [0, 1]. Binarize solution to {0, 1}.
This allows applying classification scores to all cases.
In principle, this should not do anything to properly formatted
classification inputs and outputs.'''
# Binarize solution
sol=np.ravel(solution) # convert to 1-d array
maxi = np.nanmax((filter(lambda x: x != float('inf'), sol))) # Max except NaN and Inf
mini = np.nanmin((filter(lambda x: x != float('-inf'), sol))) # Mini except NaN and Inf
if maxi == mini:
print('Warning, cannot normalize')
return [solution, prediction]
diff = maxi - mini
mid = (maxi + mini)/2.
new_solution = np.copy(solution)
new_solution[solution>=mid] = 1
new_solution[solution<mid] = 0
# Normalize and threshold predictions (takes effect only if solution not in {0, 1})
new_prediction = (np.copy(prediction) - float(mini))/float(diff)
new_prediction[new_prediction>1] = 1 # and if predictions exceed the bounds [0, 1]
new_prediction[new_prediction<0] = 0
# Make probabilities smoother
#new_prediction = np.power(new_prediction, (1./10))
return [new_solution, new_prediction]
def normalize_array (solution, prediction):
''' Use min and max of solution as scaling factors to normalize prediction,
then threshold it to [0, 1]. Binarize solution to {0, 1}.
This allows applying classification scores to all cases.
In principle, this should not do anything to properly formatted
classification inputs and outputs.'''
# Binarize solution
sol=np.ravel(solution) # convert to 1-d array
maxi = np.nanmax((filter(lambda x: x != float('inf'), sol))) # Max except NaN and Inf
mini = np.nanmin((filter(lambda x: x != float('-inf'), sol))) # Mini except NaN and Inf
if maxi == mini:
print('Warning, cannot normalize')
return [solution, prediction]
diff = maxi - mini
mid = (maxi + mini)/2.
new_solution = np.copy(solution)
new_solution[solution>=mid] = 1
new_solution[solution<mid] = 0
# Normalize and threshold predictions (takes effect only if solution not in {0, 1})
new_prediction = (np.copy(prediction) - float(mini))/float(diff)
new_prediction[new_prediction>1] = 1 # and if predictions exceed the bounds [0, 1]
new_prediction[new_prediction<0] = 0
# Make probabilities smoother
#new_prediction = np.power(new_prediction, (1./10))
return [new_solution, new_prediction]
def calc_stoi_from_spec(clean_spec, degraded_spec, analysis_len=30):
freq_bins = np.size(clean_spec, 0)
frames = np.size(clean_spec, 1)
x = np.zeros((freq_bins, frames - analysis_len + 1, analysis_len), dtype=np.float32)
y = np.zeros((freq_bins, frames - analysis_len + 1, analysis_len), dtype=np.float32)
for j in range(0, freq_bins):
for m in range(analysis_len - 1, frames, 1):
x[j, m] = clean_spec[j, m - analysis_len + 1:m + 1]
y[j, m] = degraded_spec[j, m - analysis_len + 1:m + 1]
y[j, m] = np.minimum(np.linalg.norm(x[j,m,:])/np.linalg.norm(y[j,m,:])*y[j,m,:],
(1.+np.power(10., 15./20.))*x[j,m,:]) # y is normalized and clipped
x_mean = np.mean(x, axis=(0, 1))
y_mean = np.mean(y, axis=(0, 1))
score = 0.
for j in range(0, freq_bins):
for m in range(analysis_len - 1, frames, 1):
score += np.dot(x[j, m, :] - x_mean, y[j, m, :] - y_mean) / \
(np.linalg.norm(x[j, m, :] - x_mean) * np.linalg.norm(y[j, m, :] - y_mean))
score /= (freq_bins * analysis_len)
return score
cochleagram_extractor.py 文件源码
项目：speech_feature_extractor
作者: ZhihaoDU
项目源码
文件源码
阅读 24
收藏 0
点赞 0
评论 0
def cochleagram_extractor(xx, sr, win_len, shift_len, channel_number, win_type):
fcoefs, f = make_erb_filters(sr, channel_number, 50)
fcoefs = np.flipud(fcoefs)
xf = erb_frilter_bank(xx, fcoefs)
if win_type == 'hanning':
window = np.hanning(channel_number)
elif win_type == 'hamming':
window = np.hamming(channel_number)
elif win_type == 'triangle':
window = (1 - (np.abs(channel_number - 1 - 2 * np.arange(1, channel_number + 1, 1)) / (channel_number + 1)))
else:
window = np.ones(channel_number)
window = window.reshape((channel_number, 1))
xe = np.power(xf, 2.0)
frames = 1 + ((np.size(xe, 1)-win_len) // shift_len)
cochleagram = np.zeros((channel_number, frames))
for i in range(frames):
one_frame = np.multiply(xe[:, i*shift_len:i*shift_len+win_len], np.repeat(window, win_len, 1))
cochleagram[:, i] = np.sqrt(np.mean(one_frame, 1))
cochleagram = np.where(cochleagram == 0.0, np.finfo(float).eps, cochleagram)
return cochleagram
def log_power_spectrum_extractor(x, win_len, shift_len, win_type, is_log=False):
samples = x.shape[0]
frames = (samples - win_len) // shift_len
stft = np.zeros((win_len, frames), dtype=np.complex64)
spect = np.zeros((win_len // 2 + 1, frames), dtype=np.float64)
if win_type == 'hanning':
window = np.hanning(win_len)
elif win_type == 'hamming':
window = np.hamming(win_len)
elif win_type == 'rectangle':
window = np.ones(win_len)
for i in range(frames):
one_frame = x[i*shift_len: i*shift_len+win_len]
windowed_frame = np.multiply(one_frame, window)
stft[:, i] = np.fft.fft(windowed_frame, win_len)
if is_log:
spect[:, i] = np.log(np.power(np.abs(stft[0: win_len//2+1, i]), 2.))
else:
spect[:, i] = np.power(np.abs(stft[0: win_len//2+1, i]), 2.)
return spect
