The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def _select_edge_sur(self, edges, k):
"""
Select the five cell indices surrounding each edge cell.
"""
i, j = edges[k]['k']
if k == 'n':
return ([i + 0, i + 1, i + 1, i + 1, i + 0],
[j + 1, j + 1, j + 0, j - 1, j - 1])
elif k == 'e':
return ([i - 1, i + 1, i + 1, i + 0, i - 1],
[j + 0, j + 0, j - 1, j - 1, j - 1])
elif k == 's':
return ([i - 1, i - 1, i + 0, i + 0, i - 1],
[j + 0, j + 1, j + 1, j - 1, j - 1])
elif k == 'w':
return ([i - 1, i - 1, i + 0, i + 1, i + 1],
[j + 0, j + 1, j + 1, j + 1, j + 0]) Example 2
def take(self, indices, axis=None, out=None, mode='raise'):
"""
"""
(_data, _mask) = (self._data, self._mask)
cls = type(self)
# Make sure the indices are not masked
maskindices = getattr(indices, '_mask', nomask)
if maskindices is not nomask:
indices = indices.filled(0)
# Get the data
if out is None:
out = _data.take(indices, axis=axis, mode=mode).view(cls)
else:
np.take(_data, indices, axis=axis, mode=mode, out=out)
# Get the mask
if isinstance(out, MaskedArray):
if _mask is nomask:
outmask = maskindices
else:
outmask = _mask.take(indices, axis=axis, mode=mode)
outmask |= maskindices
out.__setmask__(outmask)
return out
# Array methods Example 3
def put(a, indices, values, mode='raise'):
"""
Set storage-indexed locations to corresponding values.
This function is equivalent to `MaskedArray.put`, see that method
for details.
See Also
--------
MaskedArray.put
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode) Example 4
def fixOffset(self, offset, img):
size = img.shape
finalImg = np.ndarray(size)
indices = np.indices((self.videoSize[0],self.videoSize[1])).swapaxes(0,2).swapaxes(0,1)
indices = np.around(indices, decimals=1)
indices.shape = (self.videoSize[1] * self.videoSize[0], 2)
phi = 2 * np.arctan(np.exp(indices[:, 1] / self.videoSize[1])) - 1/2 * np.pi - offset[0]
lamb = indices[:, 0] - offset[1]
x = lamb
y = np.log(np.tan(np.pi / 4 + 1/2 * phi)) * self.videoSize[1]
finalIdx = np.ndarray((self.videoSize[1] * self.videoSize[0], 2))
finalIdx = np.around(finalIdx, decimals=1).astype(int)
finalIdx[:, 1] = y % self.videoSize[1]
finalIdx[:, 0] = x % self.videoSize[0]
finalImg[indices[:,1], indices[:,0]] = img[finalIdx[:,1], finalIdx[:,0]]
return finalImg Example 5
def test_encode_data_roundtrip():
minrand, maxrand = np.sort(np.random.randint(-427, 8848, 2))
testdata = np.round((np.sum(
np.dstack(
np.indices((512, 512),
dtype=np.float64)),
axis=2) / (511. + 511.)) * maxrand, 2) + minrand
baseval = -1000
interval = 0.1
rtripped = _decode(data_to_rgb(testdata.copy(), baseval, interval), baseval, interval)
assert testdata.min() == rtripped.min()
assert testdata.max() == rtripped.max() Example 6
def _parse_mask(mask):
r"""
Interprets a string mask to return the number of coefficients to be kept
and the indices of the first and last ones in the zig-zagged flattened DCT matrix
Example: '1-44' returns 44, first=1, last=44
Parameters
----------
mask
Returns
-------
"""
tmp = mask.split('-')
first = int(tmp[0])
last = int(tmp[1])
ncoeffs = last-first+1
return ncoeffs, first, last Example 7
def diff_approx(self, fields, pars, eps=1E-8):
nvar, N = len(fields.dependent_variables), fields.size
fpars = {key: pars[key] for key in self.pars}
fpars['dx'] = (fields['x'][-1] - fields['x'][0]) / fields['x'].size
J = np.zeros((N * nvar, N * nvar))
indices = np.indices(fields.uarray.shape)
for i, (var_index, node_index) in enumerate(zip(*map(np.ravel,
indices))):
fields_plus = fields.copy()
fields_plus.uarray[var_index, node_index] += eps
fields_moins = fields.copy()
fields_moins.uarray[var_index, node_index] -= eps
Fplus = self(fields_plus, pars)
Fmoins = self(fields_moins, pars)
J[i] = (Fplus - Fmoins) / (2 * eps)
return J.T Example 8
def continuous_loss(self, y, y_hat):
if isinstance(y_hat, DocLabel):
raise ValueError("continuous loss on discrete input")
if isinstance(y_hat[0], tuple):
y_hat = y_hat[0]
prop_marg, link_marg = y_hat
y_nodes = self.prop_encoder_.transform(y.nodes)
y_links = self.link_encoder_.transform(y.links)
prop_ix = np.indices(y.nodes.shape)
link_ix = np.indices(y.links.shape)
# relies on prop_marg and link_marg summing to 1 row-wise
prop_loss = np.sum(self.prop_cw_[y_nodes] *
(1 - prop_marg[prop_ix, y_nodes]))
link_loss = np.sum(self.link_cw_[y_links] *
(1 - link_marg[link_ix, y_links]))
loss = prop_loss + link_loss
return loss Example 9
def find_beam_position_blur(z, sigma=30):
"""Estimate direct beam position by blurring the image with a large
Gaussian kernel and finding the maximum.
Parameters
----------
sigma : float
Sigma value for Gaussian blurring kernel.
Returns
-------
center : np.array
np.array containing indices of estimated direct beam positon.
"""
blurred = ndi.gaussian_filter(z, sigma)
center = np.unravel_index(blurred.argmax(), blurred.shape)
return np.array(center) Example 10
def make_gisaxs_grid( qr_w= 10, qz_w = 12, dim_r =100,dim_z=120):
''' Dec 16, 2015, [email protected]
'''
y, x = np.indices( [dim_z,dim_r] )
Nr = int(dim_r/qp_w)
Nz = int(dim_z/qz_w)
noqs = Nr*Nz
ind = 1
for i in range(0,Nr):
for j in range(0,Nz):
y[ qr_w*i: qr_w*(i+1), qz_w*j:qz_w*(j+1)]= ind
ind += 1
return y
###########################################
#for Q-map, convert pixel to Q
########################################### Example 11
def get_reflected_angles(inc_x0, inc_y0, refl_x0, refl_y0, thetai=0.0,
pixelsize=[75,75], Lsd=5.0,dimx = 2070.,dimy=2167.):
''' Dec 16, 2015, [email protected]
giving: incident beam center: bcenx,bceny
reflected beam on detector: rcenx, rceny
sample to detector distance: Lsd, in mm
pixelsize: 75 um for Eiger4M detector
detector image size: dimx = 2070,dimy=2167 for Eiger4M detector
get reflected angle alphaf (outplane)
reflected angle thetaf (inplane )
'''
#if Lsd>=1000:#it should be something wrong and the unit should be meter
#convert Lsd from mm to m
if Lsd>=1000:
Lsd = Lsd/1000.
alphai, phi = get_incident_angles( inc_x0, inc_y0, refl_x0, refl_y0, pixelsize, Lsd)
print ('The incident_angle (alphai) is: %s'%(alphai* 180/np.pi))
px,py = pixelsize
y, x = np.indices( [int(dimy),int(dimx)] )
#alphaf = np.arctan2( (y-inc_y0)*py*10**(-6), Lsd )/2 - alphai
alphaf = np.arctan2( (y-inc_y0)*py*10**(-6), Lsd ) - alphai
thetaf = np.arctan2( (x-inc_x0)*px*10**(-6), Lsd )/2 - thetai
return alphaf,thetaf, alphai, phi Example 12
def take(self, indices, axis=None, out=None, mode='raise'):
"""
"""
(_data, _mask) = (self._data, self._mask)
cls = type(self)
# Make sure the indices are not masked
maskindices = getattr(indices, '_mask', nomask)
if maskindices is not nomask:
indices = indices.filled(0)
# Get the data
if out is None:
out = _data.take(indices, axis=axis, mode=mode).view(cls)
else:
np.take(_data, indices, axis=axis, mode=mode, out=out)
# Get the mask
if isinstance(out, MaskedArray):
if _mask is nomask:
outmask = maskindices
else:
outmask = _mask.take(indices, axis=axis, mode=mode)
outmask |= maskindices
out.__setmask__(outmask)
return out
# Array methods Example 13
def put(a, indices, values, mode='raise'):
"""
Set storage-indexed locations to corresponding values.
This function is equivalent to `MaskedArray.put`, see that method
for details.
See Also
--------
MaskedArray.put
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode) Example 14
def _parse_output(self):
unique_ids = np.unique(self.tags)
counts = np.bincount(self.tags + 1)
sort_indices = np.argsort(self.tags)
grab_indices = np.indices(self.tags.shape).ravel()[sort_indices]
dens = self.densities[sort_indices]
cp = 0
for i in unique_ids:
cp_c = cp + counts[i + 1]
if i == -1:
cp += counts[i + 1]
continue
group_indices = grab_indices[cp:cp_c]
self._groups.append(self._halo_class(self, i, group_indices,
ptype=self.ptype))
md_i = np.argmax(dens[cp:cp_c])
px, py, pz = \
[self.particle_fields['particle_position_%s' % ax][group_indices]
for ax in 'xyz']
self._max_dens[i] = (dens[cp:cp_c][md_i], px[md_i],
py[md_i], pz[md_i])
cp += counts[i + 1] Example 15
def _parse_halolist(self, threshold_adjustment):
groups = []
max_dens = {}
hi = 0
LE, RE = self.bounds
for halo in self._groups:
this_max_dens = halo.maximum_density_location()
# if the most dense particle is in the box, keep it
if np.all((this_max_dens >= LE) & (this_max_dens <= RE)):
# Now we add the halo information to OURSELVES, taken from the
# self.hop_list
# We need to mock up the HOPHaloList thingie, so we need to
# set self._max_dens
max_dens_temp = list(self._max_dens[halo.id])[0] / \
threshold_adjustment
max_dens[hi] = [max_dens_temp] + \
list(self._max_dens[halo.id])[1:4]
groups.append(self._halo_class(self, hi, ptype=self.ptype))
groups[-1].indices = halo.indices
self.comm.claim_object(groups[-1])
hi += 1
del self._groups, self._max_dens # explicit >> implicit
self._groups = groups
self._max_dens = max_dens Example 16
def test_fill_region():
for level in range(2):
rf = 2**level
output_fields = [np.zeros((NDIM*rf,NDIM*rf,NDIM*rf), "float64")
for i in range(3)]
input_fields = [np.empty(NDIM**3, "float64")
for i in range(3)]
v = np.mgrid[0.0:1.0:NDIM*1j, 0.0:1.0:NDIM*1j, 0.0:1.0:NDIM*1j]
input_fields[0][:] = v[0].ravel()
input_fields[1][:] = v[1].ravel()
input_fields[2][:] = v[2].ravel()
left_index = np.zeros(3, "int64")
ipos = np.empty((NDIM**3, 3), dtype="int64")
ind = np.indices((NDIM,NDIM,NDIM))
ipos[:,0] = ind[0].ravel()
ipos[:,1] = ind[1].ravel()
ipos[:,2] = ind[2].ravel()
ires = np.zeros(NDIM*NDIM*NDIM, "int64")
ddims = np.array([NDIM, NDIM, NDIM], dtype="int64") * rf
fill_region(input_fields, output_fields, level,
left_index, ipos, ires, ddims,
np.array([2, 2, 2], dtype="i8"))
for r in range(level + 1):
for o, i in zip(output_fields, v):
assert_equal( o[r::rf,r::rf,r::rf], i) Example 17
def weighted_distances( dx=10, dy=10, c=(5,5)):
'''
Map with weighted distances to a point
args: Dimension maps and point
'''
a = np.zeros((dx,dy))
a[c]=1
indr = np.indices(a.shape)[0,:]
indc = np.indices(a.shape)[1,:]
difr = indr-c[0]
difc = indc-c[1]
map_diff = np.sqrt((difr**2)+(difc**2))
map_diff = 1.0 - (map_diff/ map_diff.flatten().max())
# Return inverse distance map
return map_diff Example 18
def moments(data):
total = data.sum()
X, Y = np.indices(data.shape)
x = (X*data).sum()/total
y = (Y*data).sum()/total
col = data[:, int(y)]
width_x = np.sqrt(abs((np.arange(col.size)-y)**2*col).sum()/col.sum())
row = data[int(x), :]
width_y = np.sqrt(abs((np.arange(row.size)-x)**2*row).sum()/row.sum())
height = data.max()
# Return the parameters
return height, x, y, width_x, width_y
# ----------------------------------------------------------------- Example 19
def moments(data):
total = data.sum()
X, Y = np.indices(data.shape)
x = (X*data).sum()/total
y = (Y*data).sum()/total
col = data[:, int(y)]
width_x = np.sqrt(abs((np.arange(col.size)-y)**2*col).sum()/col.sum())
row = data[int(x), :]
width_y = np.sqrt(abs((np.arange(row.size)-x)**2*row).sum()/row.sum())
height = data.max()
# Return the parameters
return height, x, y, width_x, width_y
# ----------------------------------------------------------------- Example 20
def compute_dt_stats(self):
self.datestack = True
print("Computing date stats")
allmask = np.ma.getmaskarray(self.ma_stack).all(axis=0)
minidx = np.argmin(np.ma.getmaskarray(self.ma_stack), axis=0)
maxidx = np.argmin(np.ma.getmaskarray(self.ma_stack[::-1]), axis=0)
dt_stack_min = np.zeros(minidx.shape, dtype=self.dtype)
dt_stack_max = np.zeros(maxidx.shape, dtype=self.dtype)
for n, dt_o in enumerate(self.date_list_o):
dt_stack_min[minidx == n] = dt_o
dt_stack_max[maxidx == (len(self.date_list_o)-1 - n)] = dt_o
self.dt_stack_min = np.ma.array(dt_stack_min, mask=allmask)
self.dt_stack_max = np.ma.array(dt_stack_max, mask=allmask)
self.dt_stack_ptp = np.ma.masked_equal((self.dt_stack_max - self.dt_stack_min), 0)
self.dt_stack_center = self.dt_stack_min + self.dt_stack_ptp.filled(0)/2.0
#Should pull out unmasked indices at each pixel along axis 0
#Take min index along axis 0
#Then create grids by pulling out corresponding value from date_list_o Example 21
def take(self, indices, axis=None, out=None, mode='raise'):
"""
"""
(_data, _mask) = (self._data, self._mask)
cls = type(self)
# Make sure the indices are not masked
maskindices = getattr(indices, '_mask', nomask)
if maskindices is not nomask:
indices = indices.filled(0)
# Get the data
if out is None:
out = _data.take(indices, axis=axis, mode=mode).view(cls)
else:
np.take(_data, indices, axis=axis, mode=mode, out=out)
# Get the mask
if isinstance(out, MaskedArray):
if _mask is nomask:
outmask = maskindices
else:
outmask = _mask.take(indices, axis=axis, mode=mode)
outmask |= maskindices
out.__setmask__(outmask)
return out
# Array methods Example 22
def put(a, indices, values, mode='raise'):
"""
Set storage-indexed locations to corresponding values.
This function is equivalent to `MaskedArray.put`, see that method
for details.
See Also
--------
MaskedArray.put
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode) Example 23
def take(self, indices, axis=None, out=None, mode='raise'):
"""
"""
(_data, _mask) = (self._data, self._mask)
cls = type(self)
# Make sure the indices are not masked
maskindices = getattr(indices, '_mask', nomask)
if maskindices is not nomask:
indices = indices.filled(0)
# Get the data
if out is None:
out = _data.take(indices, axis=axis, mode=mode).view(cls)
else:
np.take(_data, indices, axis=axis, mode=mode, out=out)
# Get the mask
if isinstance(out, MaskedArray):
if _mask is nomask:
outmask = maskindices
else:
outmask = _mask.take(indices, axis=axis, mode=mode)
outmask |= maskindices
out.__setmask__(outmask)
return out
# Array methods Example 24
def put(a, indices, values, mode='raise'):
"""
Set storage-indexed locations to corresponding values.
This function is equivalent to `MaskedArray.put`, see that method
for details.
See Also
--------
MaskedArray.put
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode) Example 25
def _parse_halolist(self, threshold_adjustment):
groups = []
max_dens = {}
hi = 0
LE, RE = self.bounds
for halo in self._groups:
this_max_dens = halo.maximum_density_location()
# if the most dense particle is in the box, keep it
if np.all((this_max_dens >= LE) & (this_max_dens <= RE)):
# Now we add the halo information to OURSELVES, taken from the
# self.hop_list
# We need to mock up the HOPHaloList thingie, so we need to
# set self._max_dens
max_dens_temp = list(self._max_dens[halo.id])[0] / \
threshold_adjustment
max_dens[hi] = [max_dens_temp] + \
list(self._max_dens[halo.id])[1:4]
groups.append(self._halo_class(self, hi, ptype=self.ptype))
groups[-1].indices = halo.indices
self.comm.claim_object(groups[-1])
hi += 1
del self._groups, self._max_dens # explicit >> implicit
self._groups = groups
self._max_dens = max_dens Example 26
def put(a, indices, values, mode='raise'):
"""
Set storage-indexed locations to corresponding values.
This function is equivalent to `MaskedArray.put`, see that method
for details.
See Also
--------
MaskedArray.put
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode) Example 27
def put(a, indices, values, mode='raise'):
"""
Set storage-indexed locations to corresponding values.
This function is equivalent to `MaskedArray.put`, see that method
for details.
See Also
--------
MaskedArray.put
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode) Example 28
def transform_to_2d(data, max_axis):
"""
Projects 3d data cube along one axis using maximum intensity with
preservation of the signs. Adapted from nilearn.
"""
import numpy as np
# get the shape of the array we are projecting to
new_shape = list(data.shape)
del new_shape[max_axis]
# generate a 3D indexing array that points to max abs value in the
# current projection
a1, a2 = np.indices(new_shape)
inds = [a1, a2]
inds.insert(max_axis, np.abs(data).argmax(axis=max_axis))
# take the values where the absolute value of the projection
# is the highest
maximum_intensity_data = data[inds]
return np.rot90(maximum_intensity_data) Example 29
def gaussian_image(label):
label = tf.reshape(label, [-1, 2])
indices = np.indices([368, 368])[:, ::8, ::8].astype(np.float32)
coords = tf.constant(indices)
stretch = tf.reshape(tf.to_float(label), [-1, 2, 1, 1])
stretch = tf.tile(stretch, [1, 1, 46, 46])
# pdf = 1.0/(np.sqrt(2*(sigma**2)*np.pi)) * tf.exp(-tf.pow(coords-stretch,2)/(2*sigma**2))
pdf = tf.pow(coords - stretch, 2) / (2 * sigma ** 2)
pdf = tf.reduce_sum(pdf, [1])
# pdf = tf.reduce_prod(pdf,[1])
# print debug
pdf = tf.expand_dims(pdf, 3)
debug = tf.exp(-pdf) # 1.0 / (np.sqrt(2 * (sigma ** 2) * np.pi)) *
pdf_debug_img('super', debug, sigma)
return debug Example 30
def endmembers_by_query(rast, query, gt, wkt, dd=False):
'''
Returns a list of endmember locations based on a provided query, e.g.:
> query = rast[1,...] < -25 # Band 2 should be less than -25
> endmembers_by_query(rast, query, gt, wkt)
Arguments:
rast The raster array to find endmembers within
query A NumPy boolean array representing a query in the feature space
gt The GDAL GeoTransform
wkt The GDAL WKT projection
dd True for coordinates in decimal degrees
'''
assert isinstance(rast, np.ndarray), 'Requires a NumPy array'
shp = rast.shape
idx = np.indices((shp[-2], shp[-1]))
# Execute query on the indices (pixel locations), then return the coordinates
return list(pixel_to_xy([
(x, y) for y, x in idx[:,query].T
], gt, wkt, dd=dd)) Example 31
def mae(reference, predictions, idx=None, n=1):
'''
Mean absolute error (MAE) for (p x n) raster arrays, where p is the number
of bands and n is the number of pixels. Arguments:
reference Raster array of reference ("truth" or measured) data
predictions Raster array of predictions
idx Optional array of indices at which to sample the arrays
n A normalizing constant for residuals; e.g., the number
of endmembers when calculating RMSE for modeled reflectance
'''
if idx is None:
r = reference.shape[1]
residuals = reference - predictions
else:
r = len(idx)
residuals = reference[:, idx] - predictions[:, idx]
# Divide the MSE by the number of bands before taking the root
return np.apply_along_axis(lambda x: np.divide(np.abs(x).sum(), n), 0,
residuals) Example 32
def measure_background(image, Fibers, width=30, niter=3, order=3):
t = []
a,b = image.shape
ygrid,xgrid = np.indices(image.shape)
ygrid = 1. * ygrid.ravel() / a
xgrid = 1. * xgrid.ravel() / b
image = image.ravel()
s = np.arange(a*b)
for fiber in Fibers:
t.append(fiber.D*fiber.yind + fiber.xind)
t = np.hstack(t)
t = np.array(t, dtype=int)
ind = np.setdiff1d(s,t)
mask = np.zeros((a*b))
mask[ind] = 1.
mask[ind] = 1.-is_outlier(image[ind])
sel = np.where(mask==1.)[0]
for i in xrange(niter):
V = polyvander2d(xgrid[sel],ygrid[sel],[order,order])
sol = np.linalg.lstsq(V, image[sel])[0]
vals = np.dot(V,sol) - image[sel]
sel = sel[~is_outlier(vals)]
V = polyvander2d(xgrid,ygrid,[order,order])
back = np.dot(V, sol).reshape(a,b)
return back Example 33
def test_xor():
# Check on a XOR problem
y = np.zeros((10, 10))
y[:5, :5] = 1
y[5:, 5:] = 1
gridx, gridy = np.indices(y.shape)
X = np.vstack([gridx.ravel(), gridy.ravel()]).T
y = y.ravel()
for name, Tree in CLF_TREES.items():
clf = Tree(random_state=0)
clf.fit(X, y)
assert_equal(clf.score(X, y), 1.0,
"Failed with {0}".format(name))
clf = Tree(random_state=0, max_features=1)
clf.fit(X, y)
assert_equal(clf.score(X, y), 1.0,
"Failed with {0}".format(name)) Example 34
def apply_SL2C_elt_to_image(M_SL2C, src_image, out_size=None):
s_im = np.atleast_3d(src_image)
in_size = s_im.shape[:-1]
if out_size is None:
out_size = in_size
#We are going to find the location in the source image that each pixel in the output image comes from
#least squares matrix inversion (find X such that M @ X = I ==> X = inv(M) @ I = inv(M))
Minv = np.linalg.lstsq(M_SL2C, np.eye(2))[0]
#all of the x,y pairs in o_im:
pts_out = np.indices(out_size).reshape((2,-1)) #results in a 2 x (num pixels) array of indices
pts_out_a = angles_from_pixel_coords(pts_out, out_size)
pts_out_s = sphere_from_angles(pts_out_a)
pts_out_c = CP1_from_sphere(pts_out_s)
pts_in_c = np.dot(Minv, pts_out_c) # (2x2) @ (2xn) => (2xn)
pts_in_s = sphere_from_CP1(pts_in_c)
pts_in_a = angles_from_sphere(pts_in_s)
pts_in = pixel_coords_from_angles(pts_in_a, in_size)
#reshape pts into 2 x image_shape for the interpolation
o_im = get_interpolated_pixel_color(pts_in.reshape((2,)+out_size), s_im, in_size)
return o_im Example 35
def improve_ipopt(x0, prob, *args, **kwargs):
try:
import pyipopt
except ImportError:
raise Exception("PyIpopt package is not installed.")
lb = pyipopt.NLP_LOWER_BOUND_INF
ub = pyipopt.NLP_UPPER_BOUND_INF
g_L = np.zeros(prob.m)
for i in range(prob.m):
if prob.fs[i].relop == '<=':
g_L[i] = lb
g_U = np.zeros(prob.m)
def eval_grad_f(x, user_data = None):
return 2*prob.f0.P.dot(x) + prob.f0.qarray
def eval_g(x, user_data = None):
return np.array([f.eval(x) for f in prob.fs])
jac_grid = np.indices((prob.m, prob.n))
jac_r = jac_grid[0].ravel()
jac_c = jac_grid[1].ravel()
def eval_jac_g(x, flag, user_data = None):
if flag:
return (jac_r, jac_c)
else:
return np.vstack([2*f.P.dot(x)+f.qarray for f in prob.fs])
nlp = pyipopt.create(
prob.n, lb*np.ones(prob.n), ub*np.ones(prob.n),
prob.m, g_L, g_U, prob.m*prob.n, 0,
prob.f0.eval, eval_grad_f,
eval_g, eval_jac_g
)
try:
x, zl, zu, constraint_multipliers, obj, status = nlp.solve(x0)
except:
pass
return x Example 36
def check(p1, p2, base_array):
''' Checks if the values in the base array fall inside of the triangle
enclosed in the points (p1, p2, (0,0)).
Args:
p1 (`iterable`): iterable containing (x,y) coordinates of a point.
p2 (`iterable`): iterable containing (x,y) coordinates of a point.
base_array (`numpy.ndarray`): a logical array.
Returns:
`numpy.ndarray`: array with True value inside and False value outside bounds
'''
# Create 3D array of indices
idxs = np.indices(base_array.shape)
# ensure points are floats
p1 = p1.astype(float)
p2 = p2.astype(float)
# Calculate max column idx for each row idx based on interpolated line between two points
max_col_idx = (idxs[0] - p1[0]) / (p2[0] - p1[0]) * (p2[1] - p1[1]) + p1[1]
sign = np.sign(p2[0] - p1[0])
return idxs[1] * sign <= max_col_idx * sign Example 37
def _select_surround(self, i, j):
"""
Select the eight indices surrounding a given index.
"""
return ([i - 1, i - 1, i + 0, i + 1, i + 1, i + 1, i + 0, i - 1],
[j + 0, j + 1, j + 1, j + 1, j + 0, j - 1, j - 1, j - 1]) Example 38
def _select_surround_ravel(self, i, shape):
"""
Select the eight indices surrounding a flattened index.
"""
offset = shape[1]
return np.array([i + 0 - offset,
i + 1 - offset,
i + 1 + 0,
i + 1 + offset,
i + 0 + offset,
i - 1 + offset,
i - 1 + 0,
i - 1 - offset]).T Example 39
def rebin(a, newshape):
"""Rebin an array to a new shape."""
assert len(a.shape) == len(newshape)
slices = [slice(0, old, float(old) / new)
for old, new in zip(a.shape, newshape)]
coordinates = np.mgrid[slices]
indices = coordinates.astype('i')
return a[tuple(indices)] Example 40
def rotatedCrystal(V, size=(2, 2, 1), a=1.3968418, cType='gr'):
"""
Generates a triangular crystal lattice of the given size and rotates it so that the new unit vectors
align with the columns of V. The positions are set so that the center atom is at the
origin. Size is expected to be even in all directions.
'a' is the atomic distance between the atoms of the hexagonal lattice daul to this crystal.
In other words, a*sqrt(3) is the lattice constant of the triangular lattice.
The returned object is of ase.Atoms type
"""
if cType == 'gr':
cr = GB.grapheneCrystal(1, 1, 'armChair').aseCrystal(ccBond=a)
elif cType == 'tr':
numbers = [6.0]
cell = numpy.array([[a * (3.0 ** 0.5), 0, 0], [0.5 * a * (3.0 ** 0.5), 1.5 * a, 0], [0, 0, 10 * a]])
positions = numpy.array([[0, 0, 0]])
cr = ase.Atoms(numbers=numbers, positions=positions, cell=cell, pbc=[True, True, True])
elif cType == 'tr-or':
numbers = [6.0, 6.0]
cell = numpy.array([[a * (3.0 ** 0.5), 0, 0], [0, 3.0 * a, 0], [0, 0, 10 * a]])
positions = numpy.array([[0, 0, 0], [0.5 * a * (3.0 ** 0.5), 1.5 * a, 0]])
cr = ase.Atoms(numbers=numbers, positions=positions, cell=cell, pbc=[True, True, True]) # Repeating
ix = numpy.indices(size, dtype=int).reshape(3, -1)
tvecs = numpy.einsum('ki,kj', ix, cr.cell)
rPos = numpy.ndarray((len(cr) * len(tvecs), 3))
for i in range(len(cr)):
rPos[i * len(tvecs):(i + 1) * len(tvecs)] = tvecs + cr.positions[i]
# New cell size
for i in range(3):
cr.cell[i] *= size[i]
cr = Atoms(symbols=['C'] * len(rPos), positions=rPos, cell=cr.cell, pbc=[True, True, True])
center = numpy.sum(cr.cell, axis=0) * 0.5
cr.positions = cr.positions - center
cr.cell = numpy.einsum('ik,jk', cr.cell, V)
cr.positions = numpy.einsum('ik,jk', cr.positions, V)
return cr Example 41
def rotatedCrystal(V,size=(2,2,1),a=1.3968418):
"""
Generates a triangular crystal lattice of the given size and rotates it so that the new unit vectors
align with the columns of V. The positions are set so that the center atom is at the
origin. Size is expected to be even in all directions.
'a' is the atomic distance between the atoms of the hexagonal lattice daul to this crystal.
In other words, a*sqrt(3) is the lattice constant of the triangular lattice.
The returned object is of ase.Atoms type
"""
numbers = [6.0]
cell = numpy.array([[a*(3.0**0.5),0,0],[0.5*a*(3.0**0.5),1.5*a,0],[0,0,10*a]])
positions = numpy.array([[0,0,0]])
cr = ase.Atoms(numbers=numbers,positions=positions,cell=cell,pbc=[True,True,True])
# Repeating
ix = numpy.indices(size, dtype=int).reshape(3,-1)
tvecs = numpy.einsum('ki,kj',ix,cr.cell)
rPos = numpy.ndarray((len(cr)*len(tvecs),3))
for i in range(len(cr)):
rPos[i*len(tvecs):(i+1)*len(tvecs)] = tvecs + cr.positions[i]
# New cell size
for i in range(3):
cr.cell[i]*=size[i]
cr = Atoms(symbols=['C']*len(rPos), positions=rPos, cell = cr.cell, pbc=[True,True,True])
center = numpy.sum(cr.cell,axis=0)*0.5
cr.positions = cr.positions - center
cr.cell = numpy.einsum('ik,jk',cr.cell,V)
cr.positions = numpy.einsum('ik,jk',cr.positions,V)
return cr Example 42
def test_copy_detection_zero_dim(self, level=rlevel):
# Ticket #658
np.indices((0, 3, 4)).T.reshape(-1, 3) Example 43
def test_copy_detection_corner_case(self, level=rlevel):
# Ticket #658
np.indices((0, 3, 4)).T.reshape(-1, 3)
# Cannot test if NPY_RELAXED_STRIDES_CHECKING changes the strides.
# With NPY_RELAXED_STRIDES_CHECKING the test becomes superfluous,
# 0-sized reshape itself is tested elsewhere. Example 44
def test_copy_detection_corner_case2(self, level=rlevel):
# Ticket #771: strides are not set correctly when reshaping 0-sized
# arrays
b = np.indices((0, 3, 4)).T.reshape(-1, 3)
assert_equal(b.strides, (3 * b.itemsize, b.itemsize)) Example 45
def test_take(self):
tgt = [2, 3, 5]
indices = [1, 2, 4]
a = [1, 2, 3, 4, 5]
out = np.take(a, indices)
assert_equal(out, tgt) Example 46
def test_results(self):
a = np.arange(1*2*3*4).reshape(1, 2, 3, 4).copy()
aind = np.indices(a.shape)
assert_(a.flags['OWNDATA'])
for (i, j) in self.tgtshape:
# positive axis, positive start
res = np.rollaxis(a, axis=i, start=j)
i0, i1, i2, i3 = aind[np.array(res.shape) - 1]
assert_(np.all(res[i0, i1, i2, i3] == a))
assert_(res.shape == self.tgtshape[(i, j)], str((i,j)))
assert_(not res.flags['OWNDATA'])
# negative axis, positive start
ip = i + 1
res = np.rollaxis(a, axis=-ip, start=j)
i0, i1, i2, i3 = aind[np.array(res.shape) - 1]
assert_(np.all(res[i0, i1, i2, i3] == a))
assert_(res.shape == self.tgtshape[(4 - ip, j)])
assert_(not res.flags['OWNDATA'])
# positive axis, negative start
jp = j + 1 if j < 4 else j
res = np.rollaxis(a, axis=i, start=-jp)
i0, i1, i2, i3 = aind[np.array(res.shape) - 1]
assert_(np.all(res[i0, i1, i2, i3] == a))
assert_(res.shape == self.tgtshape[(i, 4 - jp)])
assert_(not res.flags['OWNDATA'])
# negative axis, negative start
ip = i + 1
jp = j + 1 if j < 4 else j
res = np.rollaxis(a, axis=-ip, start=-jp)
i0, i1, i2, i3 = aind[np.array(res.shape) - 1]
assert_(np.all(res[i0, i1, i2, i3] == a))
assert_(res.shape == self.tgtshape[(4 - ip, 4 - jp)])
assert_(not res.flags['OWNDATA']) Example 47
def test_swapaxes(self):
a = np.arange(1*2*3*4).reshape(1, 2, 3, 4).copy()
idx = np.indices(a.shape)
assert_(a.flags['OWNDATA'])
b = a.copy()
# check exceptions
assert_raises(ValueError, a.swapaxes, -5, 0)
assert_raises(ValueError, a.swapaxes, 4, 0)
assert_raises(ValueError, a.swapaxes, 0, -5)
assert_raises(ValueError, a.swapaxes, 0, 4)
for i in range(-4, 4):
for j in range(-4, 4):
for k, src in enumerate((a, b)):
c = src.swapaxes(i, j)
# check shape
shape = list(src.shape)
shape[i] = src.shape[j]
shape[j] = src.shape[i]
assert_equal(c.shape, shape, str((i, j, k)))
# check array contents
i0, i1, i2, i3 = [dim-1 for dim in c.shape]
j0, j1, j2, j3 = [dim-1 for dim in src.shape]
assert_equal(src[idx[j0], idx[j1], idx[j2], idx[j3]],
c[idx[i0], idx[i1], idx[i2], idx[i3]],
str((i, j, k)))
# check a view is always returned, gh-5260
assert_(not c.flags['OWNDATA'], str((i, j, k)))
# check on non-contiguous input array
if k == 1:
b = c Example 48
def __getslice__(self, i, j):
"""
x.__getslice__(i, j) <==> x[i:j]
Return the slice described by (i, j). The use of negative indices
is not supported.
"""
return self.__getitem__(slice(i, j)) Example 49
def argmax(self, axis=None, fill_value=None, out=None):
"""
Returns array of indices of the maximum values along the given axis.
Masked values are treated as if they had the value fill_value.
Parameters
----------
axis : {None, integer}
If None, the index is into the flattened array, otherwise along
the specified axis
fill_value : {var}, optional
Value used to fill in the masked values. If None, the output of
maximum_fill_value(self._data) is used instead.
out : {None, array}, optional
Array into which the result can be placed. Its type is preserved
and it must be of the right shape to hold the output.
Returns
-------
index_array : {integer_array}
Examples
--------
>>> a = np.arange(6).reshape(2,3)
>>> a.argmax()
5
>>> a.argmax(0)
array([1, 1, 1])
>>> a.argmax(1)
array([2, 2])
"""
if fill_value is None:
fill_value = maximum_fill_value(self._data)
d = self.filled(fill_value).view(ndarray)
return d.argmax(axis, out=out) Example 50
def take(a, indices, axis=None, out=None, mode='raise'):
"""
"""
a = masked_array(a)
return a.take(indices, axis=axis, out=out, mode=mode) 