Python numpy 模块，atleast_2d() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.atleast_2d()。
def __init__(self, X, Y, R=None, t=None, s=None, sigma2=None, maxIterations=100, tolerance=0.001, w=0):
if X.shape[1] != Y.shape[1]:
raise 'Both point clouds must have the same number of dimensions!'
self.X = X
self.Y = Y
self.TY = Y
(self.N, self.D) = self.X.shape
(self.M, _) = self.Y.shape
self.R = np.eye(self.D) if R is None else R
self.t = np.atleast_2d(np.zeros((1, self.D))) if t is None else t
self.s = 1 if s is None else s
self.sigma2 = sigma2
self.iteration = 0
self.maxIterations = maxIterations
self.tolerance = tolerance
self.w = w
self.q = 0
self.err = 0
def __init__(self, X, Y, B=None, t=None, sigma2=None, maxIterations=100, tolerance=0.001, w=0):
if X.shape[1] != Y.shape[1]:
raise 'Both point clouds must have the same number of dimensions!'
self.X = X
self.Y = Y
self.TY = Y
(self.N, self.D) = self.X.shape
(self.M, _) = self.Y.shape
self.B = np.eye(self.D) if B is None else B
self.t = np.atleast_2d(np.zeros((1, self.D))) if t is None else t
self.sigma2 = sigma2
self.iteration = 0
self.maxIterations = maxIterations
self.tolerance = tolerance
self.w = w
self.q = 0
self.err = 0
def prepare_2D_traces(data, viz_type=None, fs=None, line_names=None):
data = _np.atleast_2d(data)
N, L = data.shape
x = prepare_2D_x(L, viz_type, fs)
traces = [None] * N
for k in range(0, N):
traces[k] = go.Scatter(
x=x,
y=data[k]
)
try:
traces[k].name = line_names[k]
except TypeError:
pass
return traces
def whiteNoise(fftData, noiseLevel=80):
'''Adds White Gaussian Noise of approx. 16dB crest to a FFT block.
Parameters
----------
fftData : array of complex floats
Input fftData block (e.g. from F/D/T or S/W/G)
noiseLevel : int, optional
Average noise Level in dB [Default: -80dB]
Returns
-------
noisyData : array of complex floats
Output fftData block including white gaussian noise
'''
dimFactor = 10**(noiseLevel / 20)
fftData = _np.atleast_2d(fftData)
channels = fftData.shape[0]
NFFT = fftData.shape[1] * 2 - 2
nNoise = _np.random.rand(channels, NFFT)
nNoise = dimFactor * nNoise / _np.mean(_np.abs(nNoise))
nNoiseSpectrum = _np.fft.rfft(nNoise, axis=1)
return fftData + nNoiseSpectrum
def get_slice(coords, shape, radius):
"""Returns the slice and origin that belong to ``slice_image``"""
# interpret parameters
ndim = len(shape)
radius = validate_tuple(radius, ndim)
coords = np.atleast_2d(np.round(coords).astype(np.int))
# drop features that have no pixels inside the image
in_bounds = np.array([(coords[:, i] >= -r) & (coords[:, i] < sh + r)
for i, sh, r in zip(range(ndim), shape, radius)])
coords = coords[np.all(in_bounds, axis=0)]
# return if no coordinates are left
if len(coords) == 0:
return [slice(None, 0)] * ndim, None
# calculate the box
lower = coords.min(axis=0) - radius
upper = coords.max(axis=0) + radius + 1
# calculate the slices
origin = [None] * ndim
slices = [None] * ndim
for i, sh, low, up in zip(range(ndim), shape, lower, upper):
lower_bound_trunc = max(0, low)
upper_bound_trunc = min(sh, up)
slices[i] = slice(lower_bound_trunc, upper_bound_trunc)
origin[i] = lower_bound_trunc
return slices, origin
def _transform_frame(self, src):
"""Mapping source spectral feature x to target spectral feature y
so that minimize the mean least squared error.
More specifically, it returns the value E(p(y|x)].
Args:
src (array): shape (`order of spectral feature`) source speaker's
spectral feature that will be transformed
Returns:
array: converted spectral feature
"""
D = len(src)
# Eq.(11)
E = np.zeros((self.num_mixtures, D))
for m in range(self.num_mixtures):
xx = np.linalg.solve(self.covarXX[m], src - self.src_means[m])
E[m] = self.tgt_means[m] + self.covarYX[m].dot(xx)
# Eq.(9) p(m|x)
posterior = self.px.predict_proba(np.atleast_2d(src))
# Eq.(13) conditinal mean E[p(y|x)]
return posterior.dot(E).flatten()
def __init__(self,
file_data_source):
self.file_data_source = file_data_source
collected_files = self.file_data_source.collect_files()
# Multiple files
if isinstance(collected_files, tuple):
collected_files = np.asarray(collected_files).T
lengths = np.array([len(files) for files in collected_files])
if not (lengths == lengths[0]).all():
raise RuntimeError(
"""Mismatch of number of collected files {}.
You must collect same number of files when you collect multiple pair of files.""".format(
tuple(lengths)))
else:
collected_files = np.atleast_2d(collected_files).T
if len(collected_files) == 0:
warn("No files are collected. You might have specified wrong data source.")
self.collected_files = collected_files
def _check_freq(f):
"""Check the frequency definition."""
f = np.atleast_2d(np.asarray(f))
#
if len(f.reshape(-1)) == 1:
raise ValueError("The length of f should at least be 2.")
elif 2 in f.shape: # f of shape (N, 2) or (2, N)
if f.shape[1] is not 2:
f = f.T
elif np.squeeze(f).shape == (4,): # (fstart, fend, fwidth, fstep)
f = _pair_vectors(*tuple(np.squeeze(f)))
else: # Sequential
f = f.reshape(-1)
f.sort()
f = np.c_[f[0:-1], f[1::]]
return f
def potential(y,x,v,m,twopiG=1.,omega=None):
"""
NAME:
potential
PURPOSE:
compute the gravitational potential at a set of points
INPUT:
y - positions at which to compute the potential
x - positions of N-body particles [N]
v - velocities of N-body particles [N]
m - masses of N-body particles [N]
twopiG= (1.) value of 2 \pi G
omega= (None) if set, frequency of external harmonic oscillator
OUTPUT:
potential(y)
HISTORY:
2017-05-12 - Written - Bovy (UofT/CCA)
"""
if not omega is None:
out= omega**2.*y**2./2.
else:
out= 0.
return out\
+twopiG\
*numpy.sum(m*numpy.fabs(x-numpy.atleast_2d(y).T),axis=1)
def _control(self, time, trajectory_values=None, feedforward_values=None,
input_values=None, **kwargs):
# input abbreviations
x = input_values
yd = trajectory_values
eq = kwargs.get("eq", None)
if eq is None:
eq = calc_closest_eq_state(self._settings, input_values)
x = x - np.atleast_2d(eq).T
# this is a second version
# x = calc_small_signal_state(self._settings, is_values)
# u corresponds to a force [kg*m/s**2] = [N]
u = - np.dot(self.K, x) + np.dot(self.V, yd[0, 0])
return u
def sparse_to_dense(voxel_data, dims, dtype=np.bool):
if voxel_data.ndim != 2 or voxel_data.shape[0] != 3:
raise ValueError('voxel_data is wrong shape; should be 3xN array.')
if np.isscalar(dims):
dims = [dims] * 3
dims = np.atleast_2d(dims).T
# truncate to integers
xyz = voxel_data.astype(np.int)
# discard voxels that fall outside dims
valid_ix = ~np.any((xyz < 0) | (xyz >= dims), 0)
xyz = xyz[:, valid_ix]
out = np.zeros(dims.flatten(), dtype=dtype)
out[tuple(xyz)] = True
return out
# def get_linear_index(x, y, z, dims):
# """ Assuming xzy order. (y increasing fastest.
# TODO ensure this is right when dims are not all same
# """
# return x*(dims[1]*dims[2]) + z*dims[1] + y
def sort_eigensystem(parameters_dict):
eigenvectors = np.stack(tensor_spherical_to_cartesian(np.squeeze(parameters_dict['theta']),
np.squeeze(parameters_dict['phi']),
np.squeeze(parameters_dict['psi'])), axis=0)
eigenvalues = np.atleast_2d(np.squeeze(np.dstack([parameters_dict['d'],
parameters_dict['dperp0'],
parameters_dict['dperp1']])))
ranking = np.atleast_2d(np.squeeze(np.argsort(eigenvalues, axis=1, kind='mergesort')[:, ::-1]))
voxels_range = np.arange(ranking.shape[0])
sorted_eigenvalues = np.concatenate([eigenvalues[voxels_range, ranking[:, ind], None]
for ind in range(ranking.shape[1])], axis=1)
sorted_eigenvectors = np.stack([eigenvectors[ranking[:, ind], voxels_range, :]
for ind in range(ranking.shape[1])])
return sorted_eigenvalues, sorted_eigenvectors, ranking
def test_RvPCA():
data = get_wine_data()
for d in data:
d.cross_product()
n_datasets = len(data)
expected_output = np.array([[0.10360263],
[0.10363524],
[0.09208477],
[0.10370834],
[0.08063234],
[0.09907428],
[0.09353886],
[0.08881811],
[0.1110871],
[0.12381833]])
output, _, _ = rv_pca(data, n_datasets)
np.testing.assert_array_almost_equal(np.atleast_2d(output).T, expected_output)
def distance_as_discrepancy(dist, *summaries, observed):
"""Evaluate a distance function with signature `dist(summaries, observed)` in ELFI."""
summaries = np.column_stack(summaries)
# Ensure observed are 2d
observed = np.concatenate([np.atleast_2d(o) for o in observed], axis=1)
try:
d = dist(summaries, observed)
except ValueError as e:
raise ValueError('Incompatible data shape for the distance node. Please check '
'summary (XA) and observed (XB) output data dimensions. They '
'have to be at most 2d. Especially ensure that summary nodes '
'outputs 2d data even with batch_size=1. Original error message '
'was: {}'.format(e))
if d.ndim == 2 and d.shape[1] == 1:
d = d.reshape(-1)
return d
def autocov(x, lag=1):
"""Return the autocovariance.
Assumes a (weak) univariate stationary process with mean 0.
Realizations are in rows.
Parameters
----------
x : np.array of size (n, m)
lag : int, optional
Returns
-------
C : np.array of size (n,)
"""
x = np.atleast_2d(x)
# In R this is normalized with x.shape[1]
C = np.mean(x[:, lag:] * x[:, :-lag], axis=1)
return C
def plane2xyz(center, ij, plane):
"""
converts image pixel indices to xyz on the PLANE.
center : 2-tuple
ij : nx2 int array
plane : 4-tuple
return nx3 array.
"""
ij = np.atleast_2d(ij)
n = ij.shape[0]
ij = ij.astype('float')
xy_ray = (ij-center[None,:]) / DepthCamera.f
z = -plane[2]/(xy_ray.dot(plane[:2])+plane[3])
xyz = np.c_[xy_ray, np.ones(n)] * z[:,None]
return xyz
def _normalize_inputs(self, xi):
"""Normalize the inputs."""
xi = np.asarray(xi, dtype=float)
if xi.shape[-1] != self.ndim:
raise ValueError("The requested sample points xi have dimension %d, "
"but this interpolator has dimension %d" % (xi.shape[-1], self.ndim))
xi = np.atleast_2d(xi.copy())
for idx, (offset, scale) in enumerate(self._scale_list):
xi[..., idx] -= offset
xi[..., idx] /= scale
# take extension input account.
xi += self._ext
return xi
def inverse(self, encoded, duration=None):
'''Inverse transformation'''
ann = jams.Annotation(namespace=self.namespace, duration=duration)
for start, end, value in self.decode_intervals(encoded,
duration=duration):
# Map start:end to frames
f_start, f_end = time_to_frames([start, end],
sr=self.sr,
hop_length=self.hop_length)
confidence = np.mean(encoded[f_start:f_end+1, value])
value_dec = self.encoder.inverse_transform(np.atleast_2d(value))[0]
for vd in value_dec:
ann.append(time=start,
duration=end-start,
value=vd,
confidence=confidence)
return ann
def inverse(self, encoded, duration=None):
'''Inverse static tag transformation'''
ann = jams.Annotation(namespace=self.namespace, duration=duration)
if np.isrealobj(encoded):
detected = (encoded >= 0.5)
else:
detected = encoded
for vd in self.encoder.inverse_transform(np.atleast_2d(detected))[0]:
vid = np.flatnonzero(self.encoder.transform(np.atleast_2d(vd)))
ann.append(time=0,
duration=duration,
value=vd,
confidence=encoded[vid])
return ann
def perf_pi_continuous(self, x):
# Use history length 1 (Schreiber k=1), kernel width of 0.5 normalised units
# learnerReward.piCalcC.initialise(40, 1, 0.5);
# learnerReward.piCalcC.initialise(1, 1, 0.5);
# src = np.atleast_2d(x[0:-1]).T # start to end - 1
# dst = np.atleast_2d(x[1:]).T # 1 to end
# learnerReward.piCalcC.setObservations(src, dst)
# print "perf_pi_continuous", x
# learnerReward.piCalcC.initialise(100, 1);
# learnerReward.piCalcC.initialise(50, 1);
learnerReward.piCalcC.initialise(10, 1);
# src = np.atleast_2d(x).T # start to end - 1
# learnerReward.piCalcC.setObservations(src.reshape((src.shape[0],)))
# print "x", x.shape
learnerReward.piCalcC.setObservations(x)
# print type(src), type(dst)
# print src.shape, dst.shape
return learnerReward.piCalcC.computeAverageLocalOfObservations()# * -1
def sparse_to_dense(voxel_data, dims, dtype=np.bool):
if voxel_data.ndim!=2 or voxel_data.shape[0]!=3:
raise ValueError('voxel_data is wrong shape; should be 3xN array.')
if np.isscalar(dims):
dims = [dims]*3
dims = np.atleast_2d(dims).T
# truncate to integers
xyz = voxel_data.astype(np.int)
# discard voxels that fall outside dims
valid_ix = ~np.any((xyz < 0) | (xyz >= dims), 0)
xyz = xyz[:,valid_ix]
out = np.zeros(dims.flatten(), dtype=dtype)
out[tuple(xyz)] = True
return out
#def get_linear_index(x, y, z, dims):
#""" Assuming xzy order. (y increasing fastest.
#TODO ensure this is right when dims are not all same
#"""
#return x*(dims[1]*dims[2]) + z*dims[1] + y
def sparse_to_dense(voxel_data, dims, dtype=np.bool):
if voxel_data.ndim!=2 or voxel_data.shape[0]!=3:
raise ValueError('voxel_data is wrong shape; should be 3xN array.')
if np.isscalar(dims):
dims = [dims]*3
dims = np.atleast_2d(dims).T
# truncate to integers
xyz = voxel_data.astype(np.int)
# discard voxels that fall outside dims
valid_ix = ~np.any((xyz < 0) | (xyz >= dims), 0)
xyz = xyz[:,valid_ix]
out = np.zeros(dims.flatten(), dtype=dtype)
out[tuple(xyz)] = True
return out
#def get_linear_index(x, y, z, dims):
#""" Assuming xzy order. (y increasing fastest.
#TODO ensure this is right when dims are not all same
#"""
#return x*(dims[1]*dims[2]) + z*dims[1] + y
def __init__(self, wls, fls, sigmas, masks=None, orders='all', name=None):
self.wls = np.atleast_2d(wls)
self.fls = np.atleast_2d(fls)
self.sigmas = np.atleast_2d(sigmas)
self.masks = np.atleast_2d(masks) if masks is not None else np.ones_like(self.wls, dtype='b')
self.shape = self.wls.shape
assert self.fls.shape == self.shape, "flux array incompatible shape."
assert self.sigmas.shape == self.shape, "sigma array incompatible shape."
assert self.masks.shape == self.shape, "mask array incompatible shape."
if orders != 'all':
# can either be a numpy array or a list
orders = np.array(orders) #just to make sure
self.wls = self.wls[orders]
self.fls = self.fls[orders]
self.sigmas = self.sigmas[orders]
self.masks = self.masks[orders]
self.shape = self.wls.shape
self.orders = orders
else:
self.orders = np.arange(self.shape[0])
self.name = name
def _run_interface(self, runtime):
out_file = op.abspath(self.inputs.out_file)
self._results['out_file'] = out_file
spikes_list = np.loadtxt(self.inputs.in_spikes, dtype=int).tolist()
# No spikes
if not spikes_list:
with open(out_file, 'w') as f:
f.write('<p>No high-frequency spikes were found in this dataset</p>')
return runtime
spikes_list = [tuple(i) for i in np.atleast_2d(spikes_list).tolist()]
plot_spikes(
self.inputs.in_file, self.inputs.in_fft, spikes_list,
out_file=out_file)
return runtime
def get_knn_score_for(tree, k=5):
tree = tree_copy_with_start(tree)
tree_encoding = encoder.get_encoding([None, tree]) # This makes sure that token-based things fail
tree_str_rep = str(tree)
distances = cdist(np.atleast_2d(tree_encoding), encodings, 'cosine')
knns = np.argsort(distances)[0]
num_non_identical_nns = 0
sum_equiv_nns = 0
current_i = 0
while num_non_identical_nns < k and current_i < len(knns) and eq_class_counts[
tree.symbol] - 1 > num_non_identical_nns:
expression_idx = knns[current_i]
current_i += 1
if eq_class_idx_to_names[expression_data[expression_idx]['eq_class']] == tree.symbol and str(
expression_data[expression_idx]['tree']) == tree_str_rep:
continue # This is an identical expression, move on
num_non_identical_nns += 1
if eq_class_idx_to_names[expression_data[expression_idx]['eq_class']] == tree.symbol:
sum_equiv_nns += 1
return "(%s-nn-stat: %s)" % (k, sum_equiv_nns / k)
def fit(self, X, y, learning_rate=0.2, epochs=10000):
x = np.atleast_2d(X)
# add bias to X
temp = np.ones((x.shape[0], x.shape[1]+1))
temp[:, 0:-1] = x
x = temp
y = np.array(y)
for k in range(epochs):
# random to select one sample
i = np.random.randint(x.shape[0])
a = [x[i]]
for l in range(len(self.weights)):
a.append(self.activation(np.dot(a[l], self.weights[l])))
error = y[i] - a[-1]
deltas = [error*self.activation_derivative(a[-1])]
for l in range(len(a)-2, 0, -1):
deltas.append(deltas[-1].dot(self.weights[l].T)*self.activation_derivative(a[l]))
deltas.reverse()
for i in range(len(self.weights)):
layer = np.atleast_2d(a[i])
delta = np.atleast_2d(deltas[i])
self.weights[i] += learning_rate * layer.T.dot(delta)
def __initialize_simulation(self, x0, time_points, udata, \
integrator_options_user):
self.__x0 = inputchecks.check_states_data(x0, self.__system.nx, 0)
time_points = inputchecks.check_time_points_input(time_points)
number_of_integration_steps = time_points.size - 1
time_steps = time_points[1:] - time_points[:-1]
udata = inputchecks.check_controls_data(udata, self.__system.nu, \
number_of_integration_steps)
self.__simulation_input = ci.vertcat([np.atleast_2d(time_steps), udata])
integrator_options = integrator_options_user.copy()
integrator_options.update({"t0": 0, "tf": 1, "expand": True}) # , "number_of_finite_elements": 1})
# integrator = ci.Integrator("integrator", "rk", \
integrator = ci.Integrator("integrator", "cvodes", \
self.__dae_scaled, integrator_options)
self.__simulation = integrator.mapaccum("simulation", \
number_of_integration_steps)
def __compute_continuity_constraints(self):
integrator = ci.Integrator("integrator", "rk", self.__ode_scaled, {"t0": 0, "tf": 1, "expand": True})
params = ci.vertcat([np.atleast_2d(self.time_points[1:] - self.time_points[:-1]), \
self.optimization_variables["U"], \
ci.repmat(self.optimization_variables["Q"], 1, self.number_of_intervals), \
ci.repmat(self.optimization_variables["P"], 1, self.number_of_intervals), \
self.optimization_variables["EPS_U"]])
shooting = integrator.map("shooting", "openmp", self.number_of_intervals)
X_next = shooting(x0 = self.optimization_variables["X"][:,:-1], \
p = params)["xf"]
self.__continuity_constraints = \
self.optimization_variables["X"][:, 1:] - X_next
def check_controls_data(udata, nu, number_of_controls):
if not nu == 0:
if udata is None:
udata = np.zeros((nu, number_of_controls))
udata = np.atleast_2d(udata)
if udata.shape == (number_of_controls, nu):
udata = udata.T
if not udata.shape == (nu, number_of_controls):
raise ValueError( \
"Time-varying control values provided by user have wrong dimension.")
return udata
else:
return ci.dmatrix(0, number_of_controls)
def check_constant_controls_data(qdata, nq):
if not nq == 0:
if qdata is None:
qdata = np.zeros((nq, 1))
qdata = np.atleast_2d(qdata)
if qdata.shape == (1, nq):
qdata = qdata.T
if not qdata.shape == (nq, 1):
raise ValueError( \
"Time-constant control values provided by user have wrong dimension.")
return qdata
else:
return ci.dmatrix(0, 1)
def check_states_data(xdata, nx, number_of_intervals):
if not nx == 0:
if xdata is None:
xdata = np.zeros((nx, number_of_intervals + 1))
xdata = np.atleast_2d(xdata)
if xdata.shape == (number_of_intervals + 1, nx):
xdata = xdata.T
if not xdata.shape == (nx, number_of_intervals + 1):
raise ValueError( \
"State values provided by user have wrong dimension.")
return xdata
else:
return ci.dmatrix(0,0)
def check_measurement_data(ydata, nphi, number_of_measurements):
if ydata is None:
ydata = np.zeros((nphi, number_of_measurements))
ydata = np.atleast_2d(ydata)
if ydata.shape == (number_of_measurements, nphi):
ydata = ydata.T
if not ydata.shape == (nphi, number_of_measurements):
raise ValueError( \
"Measurement data provided by user has wrong dimension.")
return ydata
def check_measurement_weightings(wv, nphi, number_of_measurements):
if wv is None:
wv = np.ones((nphi, number_of_measurements))
wv = np.atleast_2d(wv)
if wv.shape == (number_of_measurements, nphi):
wv = wv.T
if not wv.shape == (nphi, number_of_measurements):
raise ValueError( \
"Measurement weightings provided by user have wrong dimension.")
return wv
def _reshape(self, array):
"""
checks shapes, eg convert them (2d), raise if not possible
after checks passed, set self._array and return it.
"""
if array.ndim == 1:
array = np.atleast_2d(array).T
elif array.ndim == 2:
pass
else:
shape = array.shape
# hold first dimension, multiply the rest
shape_2d = (shape[0],
functools.reduce(lambda x, y: x * y, shape[1:]))
array = np.reshape(array, shape_2d)
return array
def _add_array_to_storage(self, array):
"""
checks shapes, eg convert them (2d), raise if not possible
after checks passed, add array to self._data
"""
if array.ndim == 1:
array = np.atleast_2d(array).T
elif array.ndim == 2:
pass
else:
shape = array.shape
# hold first dimension, multiply the rest
shape_2d = (shape[0], functools.reduce(lambda x, y: x * y, shape[1:]))
array = np.reshape(array, shape_2d)
self.data.append(array)
def fit(self, x, y, learningRate=0.2, epochs=10000):
x = np.atleast_2d(x)
temp = np.ones([x.shape[0], x.shape[1]+1])
temp[:, 0:-1] = x
x = temp
for k in range(epochs):
i = np.random.randint(x.shape[0])
result = [x[i]]
for l in range(len(self._weights)):
result.append(self._activation(np.dot(result[l], self._weights[l])))
error = y[i] - result[-1]
deltas = [error * self._activationDeriv(result[-1])]
for l in range(len(self._weights)-1, 0, -1):
deltas.append(np.dot(self._weights[l], deltas[-1]) * self._activationDeriv(result[l]))
# deltas.append(deltas[-1].dot(self._weights[l].T) * self._activationDeriv(result[l]))
deltas.reverse()
for i in range(len(self._weights)):
layer = np.atleast_2d(result[i])
delta = np.atleast_2d(deltas[i])
self._weights[i] += learningRate * layer.T.dot(delta)
def register(self, callback):
self.initialize()
while self.iteration < self.maxIterations and self.err > self.tolerance:
self.iterate()
if callback:
callback(iteration=self.iteration, error=self.err, X=self.X, Y=self.TY)
return self.TY, self.R, np.atleast_2d(self.t), self.s
def register(self, callback):
self.initialize()
while self.iteration < self.maxIterations and self.err > self.tolerance:
self.iterate()
if callback:
callback(iteration=self.iteration, error=self.err, X=self.X, Y=self.TY)
return self.TY, self.B, np.atleast_2d(self.t)
def __init__(self, Y, R=None, t=None, maxIterations=100, gamma=0.1, ):
if Y is None:
raise 'Empty list of point clouds!'
dimensions = [cloud.shape[1] for cloud in Y]
if not all(dimension == dimensions[0] for dimension in dimensions):
raise 'All point clouds must have the same number of dimensions!'
self.Y = Y
self.M = [cloud.shape[0] for cloud in self.Y]
self.D = dimensions[0]
if R:
rotations = [rotation.shape for rotation in R]
if not all(rotation[0] == self.D and rotation[1] == self.D for rotation in rotations):
raise 'All rotation matrices need to be %d x %d matrices!' % (self.D, self.D)
self.R = R
else:
self.R = [np.eye(self.D) for cloud in Y]
if t:
translations = [translations.shape for translation in t]
if not all(translations[0] == 1 and translations[1] == self.D for translation in translations):
raise 'All translation vectors need to be 1 x %d matrices!' % (self.D)
self.t = t
else:
self.t = [np.atleast_2d(np.zeros((1, self.D))) for cloud in self.Y]
def plot_prediction(x_test, y_test, prediction, save=False):
import matplotlib
import matplotlib.pyplot as plt
test_size = x_test.shape[0]
fig, ax = plt.subplots(test_size, 3, figsize=(12,12), sharey=True, sharex=True)
x_test = crop_to_shape(x_test, prediction.shape)
y_test = crop_to_shape(y_test, prediction.shape)
ax = np.atleast_2d(ax)
for i in range(test_size):
cax = ax[i, 0].imshow(x_test[i])
plt.colorbar(cax, ax=ax[i,0])
cax = ax[i, 1].imshow(y_test[i, ..., 1])
plt.colorbar(cax, ax=ax[i,1])
pred = prediction[i, ..., 1]
pred -= np.amin(pred)
pred /= np.amax(pred)
cax = ax[i, 2].imshow(pred)
plt.colorbar(cax, ax=ax[i,2])
if i==0:
ax[i, 0].set_title("x")
ax[i, 1].set_title("y")
ax[i, 2].set_title("pred")
fig.tight_layout()
if save:
fig.savefig(save)
else:
fig.show()
plt.show()
def process_contig_chunk( args ):
chunk_id = args[0]
control_pkl = args[1]
cut_CMDs = args[2]
kmers = args[3]
cols_chunk = args[4]
contig_id = args[5]
n_chunks = args[6]
n_contigs = args[7]
opts = args[8]
logging.info(" - Contig %s/%s: chunk %s/%s" % ((contig_id+1), n_contigs, (chunk_id+1), (n_chunks+1)))
control_means = pickle.load(open(control_pkl, "rb"))
contig_motifs = {}
case_motif_Ns = {}
for cut_CMD in cut_CMDs:
sts,stdOutErr = mbin.run_OS_command( cut_CMD )
fns = map(lambda x: x.split("> ")[-1], cut_CMDs)
contig_ipds_sub = np.loadtxt(fns[0], dtype="float")
contig_ipds_N_sub = np.loadtxt(fns[1], dtype="int")
# If there is only one row (read) for this contig, still treat as
# a 2d matrix of many reads
contig_ipds_sub = np.atleast_2d(contig_ipds_sub)
contig_ipds_N_sub = np.atleast_2d(contig_ipds_N_sub)
for j in range(len(cols_chunk)):
motif = kmers[cols_chunk[j]]
case_contig_Ns = contig_ipds_N_sub[:,j]
if control_means.get(motif):
case_contig_means = contig_ipds_sub[:,j]
if np.sum(case_contig_Ns)>0:
case_mean = np.dot(case_contig_means, case_contig_Ns) / np.sum(case_contig_Ns)
else:
case_mean = 0
score = case_mean - control_means[motif]
contig_motifs[motif] = score
case_motif_Ns[motif] = np.sum(case_contig_Ns)
return contig_motifs,case_motif_Ns
def process_contig_chunk( args ):
chunk_id = args[0]
cut_CMDs = args[1]
kmers = args[2]
cols_chunk = args[3]
n_chunks = args[4]
min_motif_count = args[5]
logging.info(" - Control data: chunk %s/%s" % ((chunk_id+1), (n_chunks+1)))
control_means = {}
for cut_CMD in cut_CMDs:
sts,stdOutErr = mbin.run_OS_command( cut_CMD )
fns = map(lambda x: x.split("> ")[-1], cut_CMDs)
control_ipds_sub = np.loadtxt(fns[0], dtype="float")
control_ipds_N_sub = np.loadtxt(fns[1], dtype="int")
# If there is only one row (read) for this contig, still treat as
# a 2d matrix of many reads
control_ipds_sub = np.atleast_2d(control_ipds_sub)
control_ipds_N_sub = np.atleast_2d(control_ipds_N_sub)
not_found = 0
for j in range(len(cols_chunk)):
motif = kmers[cols_chunk[j]]
if np.sum(control_ipds_N_sub[:,j])>=min_motif_count:
if np.sum(control_ipds_N_sub[:,j])>0:
control_mean = np.dot(control_ipds_sub[:,j], control_ipds_N_sub[:,j]) / np.sum(control_ipds_N_sub[:,j])
else:
control_mean = 0
control_means[motif] = control_mean
else:
not_found += 1
return control_means,not_found
def computeSMatrix(self):
for m in range(self.n_tasks):
task_X = self.task_dict[m]['X']
task_Y = self.task_dict[m]['Y']
task_xi = np.array(self.xi[m])
for k in range(self.K):
# Note that transposes are different because we are using different notation than in the paper - specifically we use row vectors where they are using column vectors
# This does all data points (n) at once
inner = np.dot(np.atleast_2d(self.theta[k,:]).T, np.atleast_2d(self.theta[k,:])) + self.gamma[k]
diag_entries = np.einsum('ij,ij->i', np.dot(task_X, inner), task_X)
s_sum = -rhoFunction(task_xi)*diag_entries
s_sum += ((task_Y.T - 0.5)* np.dot(np.atleast_2d(self.theta[k,:]), task_X.T))[0,:]
s_sum += np.log(sigmoid(task_xi))
s_sum += (-0.5)*task_xi
s_sum += rhoFunction(task_xi)*(task_xi**2)
s_sum = np.sum(s_sum)
if k < self.K-1:
s_sum = s_sum + scipy.special.psi(self.small_phi1[k]) \
- scipy.special.psi(self.small_phi1[k] + self.small_phi2[k])
if k > 0:
for i in range(k):
s_sum = s_sum + scipy.special.psi(self.small_phi2[i]) \
- scipy.special.psi(self.small_phi1[i] + self.small_phi2[i])
self.s[m,k] = s_sum
if self.debug: print "s:", self.s
def updatePhi(self):
a = np.array([np.max(self.s, axis=1)]).T #as used in logsumexp trick https://hips.seas.harvard.edu/blog/2013/01/09/computing-log-sum-exp/
self.phi = np.exp(self.s - (a + np.log(np.atleast_2d(np.sum(np.exp(self.s - a),axis=1)).T)))
if self.debug:
print "phi:", self.phi
def updateTheta(self):
for k in range(self.K):
inner_sum = np.zeros((1,self.num_feats))
for m in range(self.n_tasks):
inner_sum = inner_sum + self.phi[m,k] * np.atleast_2d(self.task_vectors[m,:])
self.theta[k,:] = (np.dot(self.gamma[k],(np.dot(la.inv(self.sigma),self.mu.T) + inner_sum.T) )).T
def computeXi(self):
for m in range(self.n_tasks):
task_X = self.task_dict[m]['X']
for n in range(len(task_X)):
inner_sum = 0
for k in range(self.K):
# Note that transposes are different because we are using different notation than in the paper - specifically we use row vectors where they are using column vectors
inner_sum += self.phi[m,k]*np.dot((np.dot(np.atleast_2d(task_X[n,:]),
(np.dot(np.atleast_2d(self.theta[k,:]).T, np.atleast_2d(self.theta[k,:])) + self.gamma[k]))),
np.atleast_2d(task_X[n,:]).T)
assert inner_sum >= 0 # This number can't be negative since we are taking the square root
self.xi[m][n] = np.sqrt(inner_sum[0,0])
if self.xi[m][n]==0:
print m,n
def predictProbability(self, task, X):
prob = 0
for k in range(self.K):
numerator = np.dot(np.atleast_2d(self.theta[k,:]),X.T)
diag_entries = np.einsum('ij,ij->i', np.dot(X, self.gamma[k]), X) ##
denom = np.sqrt(1.0 + np.pi/8 * diag_entries)
prob = prob + self.phi[task,k] * sigmoid(numerator / denom)
return prob
# Code for Predicting for a new task
def dataProb(self,new_task_X,new_task_y,weights):
prod = 1
for i in range(len(new_task_X)):
sig = sigmoid(np.dot(weights,np.atleast_2d(new_task_X[i,:]).T ))
prod = prod*(sig**new_task_y[i]) * (1.0-sig)**(1-new_task_y[i])
return prod
def predictNewTask(self,new_task_X,new_task_y,pred_X,N_sam=1000):
w_dot_array = self.metropolisHastingsAlgorithm(new_task_X,new_task_y,N_sam)
predictions = []
for x_star in pred_X:
predictions.append(sum([sigmoid(np.dot(w,np.atleast_2d(x_star).T))[0,0] for w in w_dot_array])/float(N_sam))
predictions = [1.0 if p>=0.5 else 0.0 for p in predictions]
return predictions
# Helper function
def fit(self, X, y, learning_rate = 0.2, epochs = 10000):
X = np.atleast_2d(X)
# temp.shape=(X.shape[0], X.shape[1] + 1) `+1` is for bais, so X[*][-1] = 1 => numpy.dot(x, weights) + numpy.dot(1 * bais)
temp = np.ones([X.shape[0], X.shape[1] + 1])
temp[:, 0:-1] = X
X = temp
y = np.array(y)
'''
loop operation for epochs times
'''
for k in range(epochs):
# select a random line from X for training
i = np.random.randint(X.shape[0])
x = [X[i]]
# going forward network, for each layer
for l in range(len(self.weights)):
# computer the node value for each layer (O_i) using activation function
x.append(self.activation(np.dot(x[l], self.weights[l])))
# computer the error at the top layer
error = y[i] - x[-1]
deltas = [error * self.activation_deriv(x[-1])] # For output layer, Err calculation (delta is updated error)
# start backprobagation
for l in range(len(x) - 2, 0, -1): # we need to begin at the second to last layer
# compute the updated error (i,e, deltas) for each node going from top layer to input layer
deltas.append(deltas[-1].dot(self.weights[l].T) * self.activation_deriv(x[l]))
deltas.reverse()
for i in range(len(self.weights)):
layer = np.atleast_2d(x[i])
delta = np.atleast_2d(deltas[i])
self.weights[i] += learning_rate * layer.T.dot(delta)