Here are the examples of the python api numpy.sort taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
152 Examples
3
Example 51
Project: sparkit-learn Source File: test_k_means.py
def test_same_centroids(self):
X, y, X_rdd = self.make_blobs(centers=4, n_samples=200000)
local = KMeans(n_clusters=4, init='k-means++', random_state=42)
dist = SparkKMeans(n_clusters=4, init='k-means++', random_state=42)
local.fit(X)
dist.fit(X_rdd)
local_centers = np.sort(local.cluster_centers_, axis=0)
dist_centers = np.sort(dist.cluster_centers_, axis=0)
assert_array_almost_equal(local_centers, dist_centers, decimal=4)
3
Example 52
def fdr(p, q=.05):
""" Determine FDR threshold given a p value array and desired false
discovery rate q. """
s = np.sort(p)
nvox = p.shape[0]
null = np.array(range(1, nvox + 1), dtype='float') * q / nvox
below = np.where(s <= null)[0]
return s[max(below)] if any(below) else -1
3
Example 53
def _encodeToken(self, token):
"""
Randomly encode an SDR of the input token. We seed the random number
generator such that a given string will return the same SDR each time this
method is called.
@param token (str) String token
@return (list) Numpy arrays, each with a bitmap of the
encoding.
"""
random.seed(token)
return numpy.sort(random.sample(xrange(self.n), self.w))
3
Example 54
Project: OpenRenderManagement Source File: common.py
def lowerQuartile(src):
"""
[15, 15, 14, 13, 14, 14, 15, 15, 14, 13],
13 13 14 14 14 14 15 15 15 15
0 1 2 3 4 5 6 7 8 9
q1 q2
"""
res = []
sortedSrc = np.sort(src)
for i,col in enumerate(sortedSrc):
res.append( col[len(col)/4] )
return res
3
Example 55
Project: OpenRenderManagement Source File: common.py
def higherQuartile(src):
"""
"""
res = []
sortedSrc = np.sort(src)
for i,col in enumerate(sortedSrc):
q1 = col[(len(col)/4)]
q2 = col[3*(len(col)/4)]
print "q1=%d q2 = %d for %r" % (q1, q2, col)
res.append(q2)
return res
3
Example 56
Project: seaborn Source File: test_linearmodels.py
@skipif(_no_statsmodels)
def test_residplot_lowess(self):
ax = lm.residplot("x", "y", self.df, lowess=True)
nt.assert_equal(len(ax.lines), 2)
x, y = ax.lines[1].get_xydata().T
npt.assert_array_equal(x, np.sort(self.df.x))
3
Example 57
Project: scikit-tensor Source File: sptensor.py
def unfold(self, rdims, cdims=None, transp=False):
if isinstance(rdims, type(1)):
rdims = [rdims]
if transp:
cdims = rdims
rdims = setdiff1d(range(self.ndim), cdims)[::-1]
elif cdims is None:
cdims = setdiff1d(range(self.ndim), rdims)[::-1]
if not (arange(self.ndim) == sort(hstack((rdims, cdims)))).all():
raise ValueError(
'Incorrect specification of dimensions (rdims: %s, cdims: %s)'
% (str(rdims), str(cdims))
)
M = prod([self.shape[r] for r in rdims])
N = prod([self.shape[c] for c in cdims])
ridx = _build_idx(self.subs, self.vals, rdims, self.shape)
cidx = _build_idx(self.subs, self.vals, cdims, self.shape)
return unfolded_sptensor((self.vals, (ridx, cidx)), (M, N), rdims, cdims, self.shape)
3
Example 58
Project: nistats Source File: thresholding.py
def fdr_threshold(z_vals, alpha):
""" return the BH fdr for the input z_vals"""
z_vals_ = - np.sort(- z_vals)
p_vals = norm.sf(z_vals_)
n_samples = len(p_vals)
pos = p_vals < alpha * np.linspace(
.5 / n_samples, 1 - .5 / n_samples, n_samples)
if pos.any():
return (z_vals_[pos][-1] - 1.e-8)
else:
return np.infty
3
Example 59
Project: climate Source File: dataset.py
def temporal_boundaries(self):
'''Calculate the temporal range
:returns: The start and end date of the Dataset's temporal range as
a tuple in the form (start_time, end_time).
:rtype: :func:`tuple` of the form (:class:`datetime.datetime`,
:class:`datetime.datetime`)
'''
sorted_time = numpy.sort(self.times)
start_time = sorted_time[0]
end_time = sorted_time[-1]
return (start_time, end_time)
3
Example 60
def scoreatpercentile(data,per,axis=0):
'like the function in scipy.stats but with an axis argument and works on arrays'
a = np.sort(data,axis=axis)
idx = per/100. * (data.shape[axis]-1)
if (idx % 1 == 0):
return a[[slice(None) if ii != axis else idx for ii in range(a.ndim)]]
else:
lowerweight = 1-(idx % 1)
upperweight = (idx % 1)
idx = int(np.floor(idx))
return lowerweight * a[[slice(None) if ii != axis else idx for ii in range(a.ndim)]] \
+ upperweight * a[[slice(None) if ii != axis else idx+1 for ii in range(a.ndim)]]
3
Example 61
Project: numba Source File: test_sort.py
def check_argsort(self, pyfunc, cfunc, val):
orig = copy.copy(val)
expected = pyfunc(val)
got = cfunc(val)
self.assertPreciseEqual(orig[got], np.sort(orig),
msg="the array wasn't argsorted")
# Numba and Numpy results may differ if there are duplicates
# in the array
if not self.has_duplicates(orig):
self.assertPreciseEqual(got, expected)
# The original wasn't mutated
self.assertPreciseEqual(val, orig)
3
Example 62
Project: SciDB-Py Source File: test_afl.py
def test_eval_output(self):
ARR = sdb.random((4))
s = afl.sort(ARR)
out = sdb.new_array()
result = s.eval(out=out)
assert result.name is out.name
expected = np.sort(ARR.toarray())
np.testing.assert_array_almost_equal(s.toarray(), expected)
0
Example 63
def mm(self,X1,X2,p=1):
"""
Given two sampling plans chooses the one with the better space-filling properties
(as per the Morris-Mitchell criterion)
Inputs:
X1,X2-the two sampling plans
p- the distance metric to be used (p=1 rectangular-default, p=2 Euclidean)
Outputs:
Mmplan-if Mmplan=0, identical plans or equally space-
filling, if Mmplan=1, X1 is more space filling, if Mmplan=2,
X2 is more space filling
"""
#thats how two arrays are compared in their sorted form
v = np.sort(X1) == np.sort(X2)
if v.all() == True:#if True, then the designs are the same
# if np.array_equal(X1,X2) == True:
return 0
else:
#calculate the distance and multiplicity arrays
[J1 , d1] = self.jd(X1,p);m1=len(d1)
[J2 , d2] = self.jd(X2,p);m2=len(d2)
#blend the distance and multiplicity arrays together for
#comparison according to definition 1.2B. Note the different
#signs - we are maximising the d's and minimising the J's.
V1 = np.zeros((2*m1))
V1[0:len(V1):2] = d1
V1[1:len(V1):2] = -J1
V2 = np.zeros((2*m2))
V2[0:len(V2):2] = d2
V2[1:len(V2):2] = -J2
#the longer vector can be trimmed down to the length of the shorter one
m = min(m1,m2)
V1 = V1[0:m]
V2 = V2[0:m]
#generate vector c such that c(i)=1 if V1(i)>V2(i), c(i)=2 if V1(i)<V2(i)
#c(i)=0 otherwise
c = np.zeros(m)
for i in xrange(m):
if np.greater(V1[i],V2[i]) == True:
c[i] = 1
elif np.less(V1[i],V2[i]) == True:
c[i] = 2
elif np.equal(V1[i],V2[i]) == True:
c[i] = 0
#If the plans are not identical but have the same space-filling
#properties
if sum(c) == 0:
return 0
else:
#the more space-filling design (mmplan)
#is the first non-zero element of c
i = 0
while c[i] == 0:
i = i+1
return c[i]
0
Example 64
Project: RCN Source File: draw_points_coarse_fine_conv.py
def draw_points_raw(out_path, annotate_kpts=True, high_res=True, max_errors=True, sample_num=20,
plot_colored_kpt=False, indices_given=False, cfNet_path="", mult_probs=False,
merge_sets=False, pickle_path=""):
"""
This method draw points on the test set, when there are two steps for face-detection and
downsampling. So, when the dataset is created, a face-detection is done and the data is downsampled to
orig_downsample size (this is done before pickling the data). Then, later,
at train-time, a second face detection is performed and the data is downsampled
to second_downsample size. In both down_sampling stages, the kpt_norm keeps a value in the range [0,1]. So, the relative
location is the same and just normalized in the downsampled case, compared to the detected face locations.
The process of finding the original keypoint locations is as follows:
for the 2nd face detection, multiply the normalized key-point locations in the bounding_box size
and add to those values to location of the top-left position of the bounding_box. This gives the
positions before the 2nd face detection. Then, normalize the
positions by the downsampling (img) size of the first face-detection to get them normalized. Finally,
multiply the normalized key-point locations in the bounding_box size of the first face detection
and add to those values the location of the top-left position of that bounding_box.
"""
# setting the conv params to the weights
tcdcn, params = create_TCDCN_obejct(pkl_param_file)
tcdcn.load_params(pkl_param_file)
tcdcn_cfNet, params_cfNet = None, None
if cfNet_path != "":
print "loading params of cfNet"
tcdcn_cfNet, params_cfNet = create_TCDCN_obejct(cfNet_path)
tcdcn_cfNet.load_params(cfNet_path)
cfNet_model = (tcdcn_cfNet, params_cfNet)
if params['paral_conv'] in [2, 5, 6] or params['denoise_conv'] in [1, 2]:
params['mask_MTFL'] = 0
params['mask_300W'] = 1
dataSet = '300W'
num_kpt = 68
elif params['paral_conv'] in [1, 3, 4]:
params['mask_300W'] = 0
params['mask_MTFL'] = 1
dataSet = 'MTFL'
num_kpt = 5
##############################################
# getting the 1st bounding box test set data #
##############################################
sets = get_data(**params)
Train, Valid, Test = sets
sets = OrderedDict()
train_set_x, train_set_y = Train[dataSet]
sets['train'] = OrderedDict()
sets['train']['X'] = train_set_x
sets['train']['Y'] = train_set_y
sets['train']['indices'] = []
sets['train']['name'] = '1_train'
valid_set_x, valid_set_y = Valid[dataSet]
sets['valid'] = OrderedDict()
sets['valid']['X'] = valid_set_x
sets['valid']['Y'] = valid_set_y
sets['valid']['indices'] = []
sets['valid']['name'] = '2_valid'
for test_set in np.sort(Test[dataSet].keys()):
test_set_x, test_set_y = Test[dataSet][test_set]
if merge_sets:
if not 'all_sets' in sets.keys():
sets['all_sets'] = OrderedDict()
sets['all_sets']['X'] = []
sets['all_sets']['Y'] = []
sets['all_sets']['name'] = 'all_sets'
sets['all_sets']['X'].extend(test_set_x)
sets['all_sets']['Y'].append(test_set_y)
else:
sets[test_set] = OrderedDict()
sets[test_set]['X'] = test_set_x
sets[test_set]['Y'] = test_set_y
sets[test_set]['name'] = '3_test_%s' %(test_set)
sets[test_set]['indices'] = []
if rotation_file:
with open(rotation_file, 'rb') as fp:
rotation_lists = pickle.load(fp)
rotation_set = None
if merge_sets:
sets['all_sets']['X'] = np.array(sets['all_sets']['X'])
sets['all_sets']['Y'] = append_orderedDict(sets['all_sets']['Y'])
set_names = ['all_sets']
else:
#set_names = Test[dataSet].keys()
set_names = sets.keys()
for sub_set in set_names:
set_x = sets[sub_set]['X']
set_y = sets[sub_set]['Y']
set = sets[sub_set]['name']
set_y_cp = deepcopy(set_y)
set_x_cp = deepcopy(set_x)
if max_errors:
indices, error_kpt_avg_all = eval_test_set(tcdcn, params, set_x, set_y, set,
sample_num, dataSet, cfNet_model)
set_x_indx, set_y_indx = get_indices(set_x_cp, set_y_cp, indices)
elif indices_given:
set_x_indx, set_y_indx = get_indices(set_x_cp, set_y_cp, sets[sub_set]['indices'])
elif pickle_path:
pickle_indices, all_vals = get_pickle_indices(pickle_path)
set_x_indx, set_y_indx = get_indices(set_x_cp, set_y_cp, pickle_indices)
else:
set_x_indx, set_y_indx = get_subsets(set_x_cp, set_y_cp, sample_num)
##############################################
# getting the 2nd bounding box test set data #
##############################################
# preprocessing data
rng_seed = np.random.RandomState(0)
scale_mul = 0.0
translate_mul = 0.0
target_dim = params['target_dim']
td = (target_dim, target_dim)
if rotation_file and 'test' in set:
rotation_list = np.array(rotation_lists[sub_set])
rotation_set = rotation_list[indices]
else:
rotation_set = None
set_y_cp = deepcopy(set_y_indx)
set_x_cp = deepcopy(set_x_indx)
if dataSet == 'MTFL':
dist_ratio = 3.8/4.0
else:
dist_ratio = params['dist_ratio']
set_x2, set_y2 = preprocess_once(set_x_indx, set_y_indx, dist_ratio=dist_ratio, gray_scale=params['gray_scale'])
# using preprocess_iter to prepare set_x2 and set_y2 as the model expect it
set_x2, set_y2 = preprocess_iter(set_x2, set_y2, rng_seed, jitter=False, scale_mul=scale_mul,
translate_mul=translate_mul, target_dim=td, sanity=False, dset=dataSet,
use_lcn=params['use_lcn'], rotation_set=rotation_set)
# using preprocessing with gray_scale=False and use_lcn=True to get set_x_show in RGB format
set_x_show, set_y_show = preprocess_once(set_x_cp, set_y_cp, dist_ratio=dist_ratio, gray_scale=False)
set_x_show, _ = preprocess_iter(set_x_show, set_y_show, rng_seed, jitter=False, scale_mul=scale_mul,
translate_mul=translate_mul, target_dim=td, sanity=False, dset=dataSet,
use_lcn=True, rotation_set=rotation_set)
# getting the keypoint results for the first images in the test set
# face_rect has a vector of four values: rect_start_x, rect_start_y, rect_width (rect_end_x - rect_start_x), rect_height (rect_end_y - rect_start_y)
set_rect_2nd = set_y2['face_rect']
#getting the rectangle size for the detected face
rect_size_2nd = set_rect_2nd[:, 2:]
# getting the start point of the ractangle
rect_start_point_2nd = set_rect_2nd[:, :2]
# broadcasting the rectangle size and the start point by the number of the keypoints
set_rect_size_2nd = np.tile(rect_size_2nd, (1, num_kpt))
set_rect_start_2nd = np.tile(rect_start_point_2nd, (1, num_kpt))
# the keypoint positions in the normalized format
set_kpt_norm = set_y2['kpt_norm']
dim = params['target_dim']
kpt_discret = discretise_y(set_kpt_norm, dim)
set_y2['kpt_norm'] = kpt_discret
# setting batch_size to one to avoid getting different errors for different batch sizes
batch_size = 1
num_batches = int(np.ceil(set_x2.shape[0]/float(batch_size)))
epoch_error_kpt_avg = []
kpt_conv = []
for index in np.arange(num_batches):
x_batch = set_x2[index * batch_size: (index + 1) * batch_size]
if dataSet == 'MTFL':
kpt_conv_batch = tcdcn.get_keypoints_MTFL(x_batch, dropout=0)
else: # dataSet == '300W'
if 'test' in set:
bound_mask = set_y2['bound_mask'][index * batch_size: (index + 1) * batch_size]
pad_ratio = set_y2['pad_ratio'][index * batch_size: (index + 1) * batch_size]
border_pixel = padRatio_to_pixels(pad_ratio, set_x2.shape[1])
if params['denoise_conv'] == 1:
one_hot_Maps = get_one_hot_predictions(tcdcn_cfNet, x_batch, params['target_dim'])
if mult_probs:
ocular_dist = set_y2['ocular_dist'][index * batch_size: (index + 1) * batch_size]
kpt_norm = set_y2['kpt_norm'][index * batch_size: (index + 1) * batch_size]
kpt_conv_batch = get_error_mult_probs_with_iteration(tcdcn=tcdcn, tcdcn_cfNet=tcdcn_cfNet,
x=x_batch, one_hot_Maps=one_hot_Maps,
dim=params['target_dim'], y_kpt_MTFL=kpt_norm,
y_kpt_ocular_dist=ocular_dist, num_kpts=num_kpt,
border_pixel=border_pixel, set_kpts_to_border=True,
return_predictions=True)
else:
kpt_conv_batch = tcdcn.get_keypoints_MTFL(one_hot_Maps, bound_mask, border_pixel, dropout=0)
else: # using coarse_fine_conv models
kpt_conv_batch = tcdcn.get_keypoints_MTFL(x_batch, bound_mask, border_pixel, dropout=0)
else: # train and valid sets
if params['denoise_conv'] == 1:
one_hot_Maps = get_one_hot_predictions(tcdcn_cfNet, x_batch, params['target_dim'])
if mult_probs:
ocular_dist = set_y2['ocular_dist'][index * batch_size: (index + 1) * batch_size]
kpt_norm = set_y2['kpt_norm'][index * batch_size: (index + 1) * batch_size]
kpt_conv_batch = get_error_mult_probs_with_iteration(tcdcn=tcdcn, tcdcn_cfNet=tcdcn_cfNet,
x=x_batch, one_hot_Maps=one_hot_Maps,
dim=params['target_dim'], y_kpt_MTFL=kpt_norm,
y_kpt_ocular_dist=ocular_dist, num_kpts=num_kpt,
border_pixel=None, set_kpts_to_border=False,
return_predictions=True)
else:
kpt_conv_batch = tcdcn.get_keypoints_MTFL_train(one_hot_Maps, dropout=0)
else: # using coarse_fine_conv models
kpt_conv_batch = tcdcn.get_keypoints_MTFL_train(x_batch, dropout=0)
kpt_conv.extend(kpt_conv_batch)
kpt_conv = np.array(kpt_conv)
if plot_colored_kpt:
high_res = True
if high_res:
kpt_conv_shifted = (kpt_conv / float(target_dim)) * set_rect_size_2nd + set_rect_start_2nd
kpt_conv_shifted = kpt_conv_shifted.astype(int)
kpt_true_shifted = set_kpt_norm * set_rect_size_2nd + set_rect_start_2nd
kpt_true_shifted = kpt_true_shifted.astype(int)
else:
kpt_conv_shifted = kpt_conv
kpt_true_shifted = set_kpt_norm * target_dim
n_samples = set_x_indx.shape[0]
for i in xrange(n_samples):
index = i+1
if high_res:
img = set_x_indx[i]
else:
img = set_x_show[i]
if plot_colored_kpt:
img2 = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
img = cv2.cvtColor(img2, cv2.COLOR_GRAY2BGR)
img_kpts = drawpoints(img, kpt_conv_shifted[i], kpt_true_shifted[i], magnify=False, plot_colored_kpt=plot_colored_kpt)
if annotate_kpts:
img_kpts = annotateKpts(img_kpts ,kpt_conv_shifted[i], kpt_true_shifted[i])
out_file_path = "%s/%s_%i_downsampled.png" %(out_path, set, (i+1))
cv2.imwrite(out_file_path, img_kpts)
# tiling the images
print "tiling images"
img = cv2.imread(out_file_path)
row_size, col_size = img.shape[:2]
tile_images(n_samples, set, out_path, row_size, col_size)
0
Example 65
Project: yatsm Source File: prediction.py
def get_prediction(date, result_location, image_ds,
bands='all', prefix='',
after=False, before=False, qa=False,
ndv=-9999, pattern='yatsm_r*', warn_on_empty=False):
""" Output a raster with the predictions from model fit for a given date
Args:
date (int): Ordinal date for prediction image
result_location (str): Location of the results
image_ds (gdal.Dataset): Example dataset
bands (str, list): Bands to predict - 'all' for every band, or specify
a list of bands
prefix (str, optional): Use coef/rmse with refit prefix (default: '')
after (bool, optional): If date intersects a disturbed period, use next
available time segment
before (bool, optional): If date does not intersect a model, use
previous non-disturbed time segment
qa (bool, optional): Add QA flag specifying segment type (intersect,
after, or before)
ndv (int, optional): NoDataValue
pattern (str, optional): filename pattern of saved record results
warn_on_empty (bool, optional): Log warning if result contained no
result records (default: False)
Returns:
np.ndarray: A 3D numpy.ndarray containing the prediction for each band,
for each pixel
"""
# Find results
records = find_results(result_location, pattern)
# Find result attributes to extract
i_bands, _, _, _, design, design_info = find_result_attributes(
records, bands, None, prefix=prefix)
n_bands = len(i_bands)
band_names = ['Band_{0}'.format(b) for b in range(n_bands)]
if qa:
n_bands += 1
band_names.append('SegmentQAQC')
n_i_bands = len(i_bands)
# Create X matrix from date -- ignoring categorical variables
if re.match(r'.*C\(.*\).*', design):
logger.warning('Categorical variable found in design matrix not used'
' in predicted image estimate')
design = re.sub(r'[\+\-][\ ]+C\(.*\)', '', design)
X = patsy.dmatrix(design, {'x': date}).squeeze()
i_coef = []
for k, v in design_info.items():
if not re.match('C\(.*\)', k):
i_coef.append(v)
i_coef = np.sort(i_coef)
logger.debug('Allocating memory')
raster = np.ones((image_ds.RasterYSize, image_ds.RasterXSize, n_bands),
dtype=np.int16) * int(ndv)
logger.debug('Processing results')
for rec in iter_records(records, warn_on_empty=warn_on_empty):
for _qa, index in find_indices(rec, date, after=after, before=before):
if index.shape[0] == 0:
continue
# Calculate prediction
_coef = rec['coef'].take(index, axis=0).\
take(i_coef, axis=1).take(i_bands, axis=2)
raster[rec['py'][index], rec['px'][index], :n_i_bands] = \
np.tensordot(_coef, X, axes=(1, 0))
if qa:
raster[rec['py'][index], rec['px'][index], -1] = _qa
return raster, band_names
0
Example 66
def cdf( d ):
'''
Input: (d) iterable, a data set
Output: (f) NumPy array with (bins) rows and two columns
the first column are values in the range of (d)
the second column are CDF values
alternatively, think of the columns as the
domain and range of the CDF function
(finv) inverse of (f)
---------------------------------------------------------
Calculate the cuemulative distribution function
and its inverse for a given data set
'''
# the number of data points
N = float( len( d ) )
# sorted array of data points
xs = np.sort( d )
# array of unique data points
xu = np.unique( xs )
# number of unique data points
U = len( xu )
# initialize an array of U zeros
cdf = np.zeros((U))
# for each unique data point..
for i in range( U ):
# count the number of points less than
# this point, and then divide by the
# total number of data points
cdf[i] = len( xs[ xs < xu[i] ] ) / N
# f : input value --> output percentage describing
# the number of points less than the input scalar
# in the modeled distribution; if 5 is 20% into
# the distribution, then f[5] = 0.20
f = np.vstack((xu,cdf)).T
# inverse of f
# finv : input percentage --> output value that
# represents the input percentage point of the
# distribution; if 5 is 20% into the distribution,
# then finv[0.20] = 5
finv = np.fliplr(f)
return f, finv
0
Example 67
Project: cmonkey2 Source File: datamatrix.py
def quantile_normalize_scores(matrices, weights=None):
"""quantile normalize scores against each other"""
logging.info("COMPUTING WEIGHTED MEANS...")
start_time = util.current_millis()
# rearranges the scores in the input matrices into a matrix
# with |matrices| columns where the columns contain the values
# of each matrix in sorted order
flat_values = np.transpose(np.asarray([np.sort(matrix.values.flatten())
for matrix in matrices]))
elapsed = util.current_millis() - start_time
logging.info("flattened/sorted score matrices in %f s.", elapsed / 1000.0)
start_time = util.current_millis()
if weights is not None:
# multiply each column of matrix with each component of the
# weight vector: Using matrix multiplication resulted in speedup
# from 125 s. to 0.125 seconds over apply_along_axis() (1000x faster)!
scaled = weights * flat_values
scale = np.sum(np.ma.masked_array(weights, np.isnan(weights)))
tmp_mean = util.row_means(scaled) / scale
else:
tmp_mean = util.row_means(flat_values)
elapsed = util.current_millis() - start_time
logging.info("weighted means in %f s.", elapsed / 1000.0)
start_time = util.current_millis()
result = qm_result_matrices(matrices, tmp_mean)
elapsed = util.current_millis() - start_time
logging.info("result matrices built in %f s.", elapsed / 1000.0)
return result
0
Example 68
Project: pystan Source File: mstats.py
def mquantiles(a, prob=list([.25,.5,.75]), alphap=.4, betap=.4, axis=None,
limit=()):
"""
Computes empirical quantiles for a data array.
Samples quantile are defined by ``Q(p) = (1-gamma)*x[j] + gamma*x[j+1]``,
where ``x[j]`` is the j-th order statistic, and gamma is a function of
``j = floor(n*p + m)``, ``m = alphap + p*(1 - alphap - betap)`` and
``g = n*p + m - j``.
Reinterpreting the above equations to compare to **R** lead to the
equation: ``p(k) = (k - alphap)/(n + 1 - alphap - betap)``
Typical values of (alphap,betap) are:
- (0,1) : ``p(k) = k/n`` : linear interpolation of cdf
(**R** type 4)
- (.5,.5) : ``p(k) = (k - 1/2.)/n`` : piecewise linear function
(**R** type 5)
- (0,0) : ``p(k) = k/(n+1)`` :
(**R** type 6)
- (1,1) : ``p(k) = (k-1)/(n-1)``: p(k) = mode[F(x[k])].
(**R** type 7, **R** default)
- (1/3,1/3): ``p(k) = (k-1/3)/(n+1/3)``: Then p(k) ~ median[F(x[k])].
The resulting quantile estimates are approximately median-unbiased
regardless of the distribution of x.
(**R** type 8)
- (3/8,3/8): ``p(k) = (k-3/8)/(n+1/4)``: Blom.
The resulting quantile estimates are approximately unbiased
if x is normally distributed
(**R** type 9)
- (.4,.4) : approximately quantile unbiased (Cunnane)
- (.35,.35): APL, used with PWM
Parameters
----------
a : array_like
Input data, as a sequence or array of dimension at most 2.
prob : array_like, optional
List of quantiles to compute.
alphap : float, optional
Plotting positions parameter, default is 0.4.
betap : float, optional
Plotting positions parameter, default is 0.4.
axis : int, optional
Axis along which to perform the trimming.
If None (default), the input array is first flattened.
limit : tuple
Tuple of (lower, upper) values.
Values of `a` outside this open interval are ignored.
Returns
-------
mquantiles : MaskedArray
An array containing the calculated quantiles.
Notes
-----
This formulation is very similar to **R** except the calculation of
``m`` from ``alphap`` and ``betap``, where in **R** ``m`` is defined
with each type.
References
----------
.. [1] *R* statistical software at http://www.r-project.org/
Examples
--------
>>> from scipy.stats.mstats import mquantiles
>>> a = np.array([6., 47., 49., 15., 42., 41., 7., 39., 43., 40., 36.])
>>> mquantiles(a)
array([ 19.2, 40. , 42.8])
Using a 2D array, specifying axis and limit.
>>> data = np.array([[ 6., 7., 1.],
[ 47., 15., 2.],
[ 49., 36., 3.],
[ 15., 39., 4.],
[ 42., 40., -999.],
[ 41., 41., -999.],
[ 7., -999., -999.],
[ 39., -999., -999.],
[ 43., -999., -999.],
[ 40., -999., -999.],
[ 36., -999., -999.]])
>>> mquantiles(data, axis=0, limit=(0, 50))
array([[ 19.2 , 14.6 , 1.45],
[ 40. , 37.5 , 2.5 ],
[ 42.8 , 40.05, 3.55]])
>>> data[:, 2] = -999.
>>> mquantiles(data, axis=0, limit=(0, 50))
masked_array(data =
[[19.2 14.6 --]
[40.0 37.5 --]
[42.8 40.05 --]],
mask =
[[False False True]
[False False True]
[False False True]],
fill_value = 1e+20)
"""
def _quantiles1D(data,m,p):
x = np.sort(data.compressed())
n = len(x)
if n == 0:
return ma.array(np.empty(len(p), dtype=float), mask=True)
elif n == 1:
return ma.array(np.resize(x, p.shape), mask=nomask)
aleph = (n*p + m)
k = np.floor(aleph.clip(1, n-1)).astype(int)
gamma = (aleph-k).clip(0,1)
return (1.-gamma)*x[(k-1).tolist()] + gamma*x[k.tolist()]
# Initialization & checks ---------
data = ma.array(a, copy=False)
if data.ndim > 2:
raise TypeError("Array should be 2D at most !")
#
if limit:
condition = (limit[0] < data) & (data < limit[1])
data[~condition.filled(True)] = masked
#
p = np.array(prob, copy=False, ndmin=1)
m = alphap + p*(1.-alphap-betap)
# Computes quantiles along axis (or globally)
if (axis is None):
return _quantiles1D(data, m, p)
return ma.apply_along_axis(_quantiles1D, axis, data, m, p)
0
Example 69
Project: bctpy Source File: visualization.py
def align_matrices(m1, m2, dfun='sqrdiff', verbose=False, H=1e6, Texp=1,
T0=1e-3, Hbrk=10):
'''
This function aligns two matrices relative to one another by reordering
the nodes in M2. The function uses a version of simulated annealing.
Parameters
----------
M1 : NxN np.ndarray
first connection matrix
M2 : NxN np.ndarray
second connection matrix
dfun : str
distance metric to use for matching
'absdiff' : absolute difference
'sqrdiff' : squared difference (default)
'cosang' : cosine of vector angle
verbose : bool
print out cost at each iteration. Default False.
H : int
annealing parameter, default value 1e6
Texp : int
annealing parameter, default value 1. Coefficient of H s.t.
Texp0=1-Texp/H
T0 : float
annealing parameter, default value 1e-3
Hbrk : int
annealing parameter, default value = 10. Coefficient of H s.t.
Hbrk0 = H/Hkbr
Returns
-------
Mreordered : NxN np.ndarray
reordered connection matrix M2
Mindices : Nx1 np.ndarray
reordered indices
cost : float
objective function distance between M1 and Mreordered
Notes
-----
Connection matrices can be weighted or binary, directed or undirected.
They must have the same number of nodes. M1 can be entered in any
node ordering.
Note that in general, the outcome will depend on the initial condition
(the setting of the random number seed). Also, there is no good way to
determine optimal annealing parameters in advance - these parameters
will need to be adjusted "by hand" (particularly H, Texp, T0, and Hbrk).
For large and/or dense matrices, it is highly recommended to perform
exploratory runs varying the settings of 'H' and 'Texp' and then select
the best values.
Based on extensive testing, it appears that T0 and Hbrk can remain
unchanged in most cases. Texp may be varied from 1-1/H to 1-10/H, for
example. H is the most important parameter - set to larger values as
the problem size increases. Good solutions can be obtained for
matrices up to about 100 nodes. It is advisable to run this function
multiple times and select the solution(s) with the lowest 'cost'.
If the two matrices are related it may be very helpful to pre-align them
by reordering along their largest eigenvectors:
[v,~] = eig(M1); v1 = abs(v(:,end)); [a1,b1] = sort(v1);
[v,~] = eig(M2); v2 = abs(v(:,end)); [a2,b2] = sort(v2);
[a,b,c] = overlapMAT2(M1(b1,b1),M2(b2,b2),'dfun',1);
Setting 'Texp' to zero cancels annealing and uses a greedy algorithm
instead.
'''
n = len(m1)
if n < 2:
raise BCTParamError("align_matrix will infinite loop on a singleton "
"or null matrix.")
# define maxcost (greatest possible difference) and lowcost
if dfun in ('absdiff', 'absdff'):
maxcost = np.sum(np.abs(np.sort(m1.flat) - np.sort(m2.flat)[::-1]))
lowcost = np.sum(np.abs(m1 - m2)) / maxcost
elif dfun in ('sqrdiff', 'sqrdff'):
maxcost = np.sum((np.sort(m1.flat) - np.sort(m2.flat)[::-1])**2)
lowcost = np.sum((m1 - m2)**2) / maxcost
elif dfun == 'cosang':
maxcost = np.pi / 2
lowcost = np.arccos(np.dot(m1.flat, m2.flat) /
np.sqrt(np.dot(m1.flat, m1.flat) * np.dot(m2.flat, m2.flat))) / maxcost
else:
raise BCTParamError('dfun must be absdiff or sqrdiff or cosang')
mincost = lowcost
anew = np.arange(n)
amin = np.arange(n)
h = 0
hcnt = 0
# adjust annealing parameters from user provided coefficients
# H determines the maximal number of steps (user-provided)
# Texp determines the steepness of the temperature gradient
Texp = 1 - Texp / H
# T0 sets the initial temperature and scales the energy term (user provided)
# Hbrk sets a break point for the stimulation
Hbrk = H / Hbrk
while h < H:
h += 1
hcnt += 1
# terminate if no new mincost has been found for some time
if hcnt > Hbrk:
break
# current temperature
T = T0 * (Texp**h)
# choose two positions at random and flip them
atmp = anew.copy()
r1, r2 = np.random.randint(n, size=(2,))
while r1 == r2:
r2 = np.random.randint(n)
atmp[r1] = anew[r2]
atmp[r2] = anew[r1]
m2atmp = m2[np.ix_(atmp, atmp)]
if dfun in ('absdiff', 'absdff'):
costnew = np.sum(np.abs(m1 - m2atmp)) / maxcost
elif dfun in ('sqrdiff', 'sqrdff'):
costnew = np.sum((m1 - m2atmp)**2) / maxcost
elif dfun == 'cosang':
costnew = np.arccos(np.dot(m1.flat, m2atmp.flat) / np.sqrt(
np.dot(m1.flat, m1.flat) * np.dot(m2.flat, m2.flat))) / maxcost
if costnew < lowcost or np.random.random() < np.exp(-(costnew - lowcost) / T):
anew = atmp
lowcost = costnew
# is this the absolute best?
if lowcost < mincost:
amin = anew
mincost = lowcost
if verbose:
print('step %i ... current lowest cost = %f' % (h, mincost))
hcnt = 0
# if the cost is 0 we're done
if mincost == 0:
break
if verbose:
print('step %i ... final lowest cost = %f' % (h, mincost))
M_reordered = m2[np.ix_(amin, amin)]
M_indices = amin
cost = mincost
return M_reordered, M_indices, cost
0
Example 70
def kstest(rvs, cdf, args=(), N=20, alternative = 'two_sided', mode='approx',**kwds):
"""
Perform the Kolmogorov-Smirnov test for goodness of fit
This performs a test of the distribution G(x) of an observed
random variable against a given distribution F(x). Under the null
hypothesis the two distributions are identical, G(x)=F(x). The
alternative hypothesis can be either 'two_sided' (default), 'less'
or 'greater'. The KS test is only valid for continuous distributions.
Parameters
----------
rvs : string or array or callable
string: name of a distribution in scipy.stats
array: 1-D observations of random variables
callable: function to generate random variables, requires keyword
argument `size`
cdf : string or callable
string: name of a distribution in scipy.stats, if rvs is a string then
cdf can evaluate to `False` or be the same as rvs
callable: function to evaluate cdf
args : tuple, sequence
distribution parameters, used if rvs or cdf are strings
N : int
sample size if rvs is string or callable
alternative : 'two_sided' (default), 'less' or 'greater'
defines the alternative hypothesis (see explanation)
mode : 'approx' (default) or 'asymp'
defines the distribution used for calculating p-value
'approx' : use approximation to exact distribution of test statistic
'asymp' : use asymptotic distribution of test statistic
Returns
-------
D : float
KS test statistic, either D, D+ or D-
p-value : float
one-tailed or two-tailed p-value
Notes
-----
In the one-sided test, the alternative is that the empirical
cuemulative distribution function of the random variable is "less"
or "greater" than the cuemulative distribution function F(x) of the
hypothesis, G(x)<=F(x), resp. G(x)>=F(x).
Examples
--------
>>> from scipy import stats
>>> import numpy as np
>>> from scipy.stats import kstest
>>> x = np.linspace(-15,15,9)
>>> kstest(x,'norm')
(0.44435602715924361, 0.038850142705171065)
>>> np.random.seed(987654321) # set random seed to get the same result
>>> kstest('norm','',N=100)
(0.058352892479417884, 0.88531190944151261)
is equivalent to this
>>> np.random.seed(987654321)
>>> kstest(stats.norm.rvs(size=100),'norm')
(0.058352892479417884, 0.88531190944151261)
Test against one-sided alternative hypothesis:
>>> np.random.seed(987654321)
Shift distribution to larger values, so that cdf_dgp(x)< norm.cdf(x):
>>> x = stats.norm.rvs(loc=0.2, size=100)
>>> kstest(x,'norm', alternative = 'less')
(0.12464329735846891, 0.040989164077641749)
Reject equal distribution against alternative hypothesis: less
>>> kstest(x,'norm', alternative = 'greater')
(0.0072115233216311081, 0.98531158590396395)
Don't reject equal distribution against alternative hypothesis: greater
>>> kstest(x,'norm', mode='asymp')
(0.12464329735846891, 0.08944488871182088)
Testing t distributed random variables against normal distribution:
With 100 degrees of freedom the t distribution looks close to the normal
distribution, and the kstest does not reject the hypothesis that the sample
came from the normal distribution
>>> np.random.seed(987654321)
>>> stats.kstest(stats.t.rvs(100,size=100),'norm')
(0.072018929165471257, 0.67630062862479168)
With 3 degrees of freedom the t distribution looks sufficiently different
from the normal distribution, that we can reject the hypothesis that the
sample came from the normal distribution at a alpha=10% level
>>> np.random.seed(987654321)
>>> stats.kstest(stats.t.rvs(3,size=100),'norm')
(0.131016895759829, 0.058826222555312224)
"""
if isinstance(rvs, string_types):
#cdf = getattr(stats, rvs).cdf
if (not cdf) or (cdf == rvs):
cdf = getattr(distributions, rvs).cdf
rvs = getattr(distributions, rvs).rvs
else:
raise AttributeError('if rvs is string, cdf has to be the same distribution')
if isinstance(cdf, string_types):
cdf = getattr(distributions, cdf).cdf
if callable(rvs):
kwds = {'size':N}
vals = np.sort(rvs(*args,**kwds))
else:
vals = np.sort(rvs)
N = len(vals)
cdfvals = cdf(vals, *args)
if alternative in ['two_sided', 'greater']:
Dplus = (np.arange(1.0, N+1)/N - cdfvals).max()
if alternative == 'greater':
return Dplus, distributions.ksone.sf(Dplus,N)
if alternative in ['two_sided', 'less']:
Dmin = (cdfvals - np.arange(0.0, N)/N).max()
if alternative == 'less':
return Dmin, distributions.ksone.sf(Dmin,N)
if alternative == 'two_sided':
D = np.max([Dplus,Dmin])
if mode == 'asymp':
return D, distributions.kstwobign.sf(D*np.sqrt(N))
if mode == 'approx':
pval_two = distributions.kstwobign.sf(D*np.sqrt(N))
if N > 2666 or pval_two > 0.80 - N*0.3/1000.0 :
return D, distributions.kstwobign.sf(D*np.sqrt(N))
else:
return D, distributions.ksone.sf(D,N)*2
0
Example 71
Project: scikit-image Source File: test_corner.py
def test_rotated_img():
"""
The harris filter should yield the same results with an image and it's
rotation.
"""
im = img_as_float(data.astronaut().mean(axis=2))
im_rotated = im.T
# Moravec
results = peak_local_max(corner_moravec(im),
min_distance=10, threshold_rel=0)
results_rotated = peak_local_max(corner_moravec(im_rotated),
min_distance=10, threshold_rel=0)
assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all()
assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()
# Harris
results = peak_local_max(corner_harris(im),
min_distance=10, threshold_rel=0)
results_rotated = peak_local_max(corner_harris(im_rotated),
min_distance=10, threshold_rel=0)
assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all()
assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()
# Shi-Tomasi
results = peak_local_max(corner_shi_tomasi(im),
min_distance=10, threshold_rel=0)
results_rotated = peak_local_max(corner_shi_tomasi(im_rotated),
min_distance=10, threshold_rel=0)
assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all()
assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()
0
Example 72
Project: statsmodels Source File: stattools.py
def robust_skewness(y, axis=0):
"""
Calculates the four skewness measures in Kim & White
Parameters
----------
y : array-like
axis : int or None, optional
Axis along which the skewness measures are computed. If `None`, the
entire array is used.
Returns
-------
sk1 : ndarray
The standard skewness estimator.
sk2 : ndarray
Skewness estimator based on quartiles.
sk3 : ndarray
Skewness estimator based on mean-median difference, standardized by
absolute deviation.
sk4 : ndarray
Skewness estimator based on mean-median difference, standardized by
standard deviation.
Notes
-----
The robust skewness measures are defined
.. math ::
SK_{2}=\\frac{\\left(q_{.75}-q_{.5}\\right)
-\\left(q_{.5}-q_{.25}\\right)}{q_{.75}-q_{.25}}
.. math ::
SK_{3}=\\frac{\\mu-\\hat{q}_{0.5}}
{\\hat{E}\\left[\\left|y-\\hat{\\mu}\\right|\\right]}
.. math ::
SK_{4}=\\frac{\\mu-\\hat{q}_{0.5}}{\\hat{\\sigma}}
.. [1] Tae-Hwan Kim and Halbert White, "On more robust estimation of
skewness and kurtosis," Finance Research Letters, vol. 1, pp. 56-73,
March 2004.
"""
if axis is None:
y = y.ravel()
axis = 0
y = np.sort(y, axis)
q1, q2, q3 = np.percentile(y, [25.0, 50.0, 75.0], axis=axis)
mu = y.mean(axis)
shape = (y.size,)
if axis is not None:
shape = list(mu.shape)
shape.insert(axis, 1)
shape = tuple(shape)
mu_b = np.reshape(mu, shape)
q2_b = np.reshape(q2, shape)
sigma = np.mean(((y - mu_b)**2), axis)
sk1 = stats.skew(y, axis=axis)
sk2 = (q1 + q3 - 2.0 * q2) / (q3 - q1)
sk3 = (mu - q2) / np.mean(abs(y - q2_b), axis=axis)
sk4 = (mu - q2) / sigma
return sk1, sk2, sk3, sk4
0
Example 73
Project: statsmodels Source File: _lilliefors.py
def ksstat(x, cdf, alternative='two_sided', args=()):
"""
Calculate statistic for the Kolmogorov-Smirnov test for goodness of fit
This calculates the test statistic for a test of the distribution G(x) of an observed
variable against a given distribution F(x). Under the null
hypothesis the two distributions are identical, G(x)=F(x). The
alternative hypothesis can be either 'two_sided' (default), 'less'
or 'greater'. The KS test is only valid for continuous distributions.
Parameters
----------
x : array_like, 1d
array of observations
cdf : string or callable
string: name of a distribution in scipy.stats
callable: function to evaluate cdf
alternative : 'two_sided' (default), 'less' or 'greater'
defines the alternative hypothesis (see explanation)
args : tuple, sequence
distribution parameters for call to cdf
Returns
-------
D : float
KS test statistic, either D, D+ or D-
See Also
--------
scipy.stats.kstest
Notes
-----
In the one-sided test, the alternative is that the empirical
cuemulative distribution function of the random variable is "less"
or "greater" than the cuemulative distribution function F(x) of the
hypothesis, G(x)<=F(x), resp. G(x)>=F(x).
In contrast to scipy.stats.kstest, this function only calculates the
statistic which can be used either as distance measure or to implement
case specific p-values.
"""
nobs = float(len(x))
if isinstance(cdf, string_types):
cdf = getattr(stats.distributions, cdf).cdf
elif hasattr(cdf, 'cdf'):
cdf = getattr(cdf, 'cdf')
x = np.sort(x)
cdfvals = cdf(x, *args)
if alternative in ['two_sided', 'greater']:
Dplus = (np.arange(1.0, nobs+1)/nobs - cdfvals).max()
if alternative == 'greater':
return Dplus
if alternative in ['two_sided', 'less']:
Dmin = (cdfvals - np.arange(0.0, nobs)/nobs).max()
if alternative == 'less':
return Dmin
D = np.max([Dplus,Dmin])
return D
0
Example 74
def _learnStep(self):
k = len(self.allsamples) % self.windowSize
sample, fit = self._produceNewSample()
self.samples[k], self.fitnesses[k] = sample, fit
self.generation += 1
if len(self.allsamples) < self.windowSize:
return
if self.verbose and len(self.allsamples) % 100 == 0:
print((len(self.allsamples), min(self.fitnesses), max(self.fitnesses)))
# print(len(self.allsamples), min(self.fitnesses), max(self.fitnesses)#, self.alphas)
updateSize = self._computeUpdateSize(self._computeDensities(sample), k)
self.allUpdateSizes.append(deepcopy(updateSize))
if sum(updateSize) > 0:
# update parameters
if self.numberOfCenters > 1:
self._updateAlphas(updateSize)
if not self.mutative:
self._updateMus(updateSize, sample)
self._updateSigmas(updateSize, sample)
else:
self._updateSigmas(updateSize, sample)
self._updateMus(updateSize, sample)
if self.adaptiveShaping:
self._updateShaping()
# storage, e.g. for plotting
self.allalphas.append(deepcopy(self.alphas))
self.allsigmas.append(deepcopy(self.sigmas))
self.allmus.append(deepcopy(self.mus))
if self.oneFifthRule and len(self.allsamples) % 10 == 0 and len(self.allsamples) > 2 * self.windowSize:
lastBatch = self.allfitnesses[-self.windowSize:]
secondLast = self.allfitnesses[-2 * self.windowSize:-self.windowSize]
sortedLast = sort(lastBatch)
sortedSecond = sort(secondLast)
index = int(self.windowSize * 0.8)
if sortedLast[index] >= sortedSecond[index]:
self.sigmas = [1.2 * sigma for sigma in self.sigmas]
#print("+")
else:
self.sigmas = [0.5 * sigma for sigma in self.sigmas]
0
Example 75
def testTiming(self):
reload(test_functions)
ray.init(start_ray_local=True, num_workers=3)
# measure the time required to submit a remote task to the scheduler
elapsed_times = []
for _ in range(1000):
start_time = time.time()
test_functions.empty_function.remote()
end_time = time.time()
elapsed_times.append(end_time - start_time)
elapsed_times = np.sort(elapsed_times)
average_elapsed_time = sum(elapsed_times) / 1000
print "Time required to submit an empty function call:"
print " Average: {}".format(average_elapsed_time)
print " 90th percentile: {}".format(elapsed_times[900])
print " 99th percentile: {}".format(elapsed_times[990])
print " worst: {}".format(elapsed_times[999])
# average_elapsed_time should be about 0.00038
# measure the time required to submit a remote task to the scheduler (where the remote task returns one value)
elapsed_times = []
for _ in range(1000):
start_time = time.time()
test_functions.trivial_function.remote()
end_time = time.time()
elapsed_times.append(end_time - start_time)
elapsed_times = np.sort(elapsed_times)
average_elapsed_time = sum(elapsed_times) / 1000
print "Time required to submit a trivial function call:"
print " Average: {}".format(average_elapsed_time)
print " 90th percentile: {}".format(elapsed_times[900])
print " 99th percentile: {}".format(elapsed_times[990])
print " worst: {}".format(elapsed_times[999])
# average_elapsed_time should be about 0.001
# measure the time required to submit a remote task to the scheduler and get the result
elapsed_times = []
for _ in range(1000):
start_time = time.time()
x = test_functions.trivial_function.remote()
ray.get(x)
end_time = time.time()
elapsed_times.append(end_time - start_time)
elapsed_times = np.sort(elapsed_times)
average_elapsed_time = sum(elapsed_times) / 1000
print "Time required to submit a trivial function call and get the result:"
print " Average: {}".format(average_elapsed_time)
print " 90th percentile: {}".format(elapsed_times[900])
print " 99th percentile: {}".format(elapsed_times[990])
print " worst: {}".format(elapsed_times[999])
# average_elapsed_time should be about 0.0013
# measure the time required to do do a put
elapsed_times = []
for _ in range(1000):
start_time = time.time()
ray.put(1)
end_time = time.time()
elapsed_times.append(end_time - start_time)
elapsed_times = np.sort(elapsed_times)
average_elapsed_time = sum(elapsed_times) / 1000
print "Time required to put an int:"
print " Average: {}".format(average_elapsed_time)
print " 90th percentile: {}".format(elapsed_times[900])
print " 99th percentile: {}".format(elapsed_times[990])
print " worst: {}".format(elapsed_times[999])
# average_elapsed_time should be about 0.00087
ray.worker.cleanup()
0
Example 76
Project: scikit-learn Source File: test_rcv1.py
def test_fetch_rcv1():
try:
data1 = fetch_rcv1(shuffle=False, download_if_missing=False)
except IOError as e:
if e.errno == errno.ENOENT:
raise SkipTest("Download RCV1 dataset to run this test.")
X1, Y1 = data1.data, data1.target
cat_list, s1 = data1.target_names.tolist(), data1.sample_id
# test sparsity
assert_true(sp.issparse(X1))
assert_true(sp.issparse(Y1))
assert_equal(60915113, X1.data.size)
assert_equal(2606875, Y1.data.size)
# test shapes
assert_equal((804414, 47236), X1.shape)
assert_equal((804414, 103), Y1.shape)
assert_equal((804414,), s1.shape)
assert_equal(103, len(cat_list))
# test ordering of categories
first_categories = [u'C11', u'C12', u'C13', u'C14', u'C15', u'C151']
assert_array_equal(first_categories, cat_list[:6])
# test number of sample for some categories
some_categories = ('GMIL', 'E143', 'CCAT')
number_non_zero_in_cat = (5, 1206, 381327)
for num, cat in zip(number_non_zero_in_cat, some_categories):
j = cat_list.index(cat)
assert_equal(num, Y1[:, j].data.size)
# test shuffling and subset
data2 = fetch_rcv1(shuffle=True, subset='train', random_state=77,
download_if_missing=False)
X2, Y2 = data2.data, data2.target
s2 = data2.sample_id
# The first 23149 samples are the training samples
assert_array_equal(np.sort(s1[:23149]), np.sort(s2))
# test some precise values
some_sample_ids = (2286, 3274, 14042)
for sample_id in some_sample_ids:
idx1 = s1.tolist().index(sample_id)
idx2 = s2.tolist().index(sample_id)
feature_values_1 = X1[idx1, :].toarray()
feature_values_2 = X2[idx2, :].toarray()
assert_almost_equal(feature_values_1, feature_values_2)
target_values_1 = Y1[idx1, :].toarray()
target_values_2 = Y2[idx2, :].toarray()
assert_almost_equal(target_values_1, target_values_2)
0
Example 77
Project: pyscf Source File: addons.py
def frac_occ_(mf, tol=1e-3):
assert(isinstance(mf, hf.RHF))
old_get_occ = mf.get_occ
def get_occ(mo_energy, mo_coeff=None):
mol = mf.mol
nocc = mol.nelectron // 2
sort_mo_energy = numpy.sort(mo_energy)
lumo = sort_mo_energy[nocc]
if abs(sort_mo_energy[nocc-1] - lumo) < tol:
mo_occ = numpy.zeros_like(mo_energy)
mo_occ[mo_energy<lumo] = 2
lst = abs(mo_energy-lumo) < tol
degen = int(lst.sum())
frac = 2.*numpy.count_nonzero(lst & (mo_occ == 2))/degen
mo_occ[lst] = frac
logger.warn(mf, 'fraction occ = %6g for orbitals %s',
frac, numpy.where(lst)[0])
logger.info(mf, 'HUMO = %.12g LUMO = %.12g',
sort_mo_energy[nocc-1], sort_mo_energy[nocc])
logger.debug(mf, ' mo_energy = %s', mo_energy)
else:
mo_occ = old_get_occ(mo_energy, mo_coeff)
return mo_occ
mf.get_occ = get_occ
return mf
0
Example 78
def get_data (self, td):
"""
main work is done here: imports the data from datx files into
time_data object td and saves also a '*.h5' file so this import
need not be performed every time the data is needed
"""
if not path.isfile(self.from_file):
# no file there
time_data_import.getdata(self, td)
return
#browse datx information
f0 = open(self.from_file)
config = ConfigParser.ConfigParser()
config.readfp(f0)
sample_rate = float(config.get('keywords', 'sample_rate'))
# reading sample-rate from index-file
channels = []
d_files = {}
# Loop over all channels assigned in index-file
for channel in sort(config.options('channels')):
ch = datx_channel(config, channel)
if ch.label.find('Mic') >= 0:
channels.append(ch)
if not d_files.has_key(ch.d_file):
d_files[ch.d_file] = \
datx_d_file(path.join(path.dirname(self.from_file), \
config.get(ch.d_file, 'fn')), 32)
numchannels = len(channels)
# prepare hdf5
name = td.name
if name == "":
name = path.join(td_dir, \
path.splitext(path.basename(self.from_file))[0]+'.h5')
else:
if td.h5f != None:
td.h5f.close()
# TODO problems with already open h5 files from other instances
f5h = tables.open_file(name, mode = 'w')
ac = f5h.create_earray(f5h.root, 'time_data', \
tables.atom.Float32Atom(), (0, numchannels))
ac.set_attr('sample_freq', sample_rate)
block_data = \
zeros((128*d_files[channels[0].d_file].num_samples_per_block, \
numchannels), 'Float32')
flag = 0
while(not flag):
for i in d_files.values():
flag = i.get_next_blocks()
if flag:
continue
for i in range(numchannels):
data = d_files[channels[i].d_file].data[channels[i].ch_no]
block_data[:data.size, i] = channels[i].scale(data)
ac.append(block_data[:data.size])
f5h.close()
f0.close()
for i in d_files.values():
i.f.close()
td.name = name
td.load_data()
0
Example 79
def median_cihs(data, alpha=0.05, axis=None):
"""
Computes the alpha-level confidence interval for the median of the data.
Uses the Hettmasperger-Sheather method.
Parameters
----------
data : array_like
Input data. Masked values are discarded. The input should be 1D only,
or `axis` should be set to None.
alpha : float, optional
Confidence level of the intervals.
axis : int or None, optional
Axis along which to compute the quantiles. If None, use a flattened
array.
Returns
-------
median_cihs
Alpha level confidence interval.
"""
def _cihs_1D(data, alpha):
data = np.sort(data.compressed())
n = len(data)
alpha = min(alpha, 1-alpha)
k = int(binom._ppf(alpha/2., n, 0.5))
gk = binom.cdf(n-k,n,0.5) - binom.cdf(k-1,n,0.5)
if gk < 1-alpha:
k -= 1
gk = binom.cdf(n-k,n,0.5) - binom.cdf(k-1,n,0.5)
gkk = binom.cdf(n-k-1,n,0.5) - binom.cdf(k,n,0.5)
I = (gk - 1 + alpha)/(gk - gkk)
lambd = (n-k) * I / float(k + (n-2*k)*I)
lims = (lambd*data[k] + (1-lambd)*data[k-1],
lambd*data[n-k-1] + (1-lambd)*data[n-k])
return lims
data = ma.rray(data, copy=False)
# Computes quantiles along axis (or globally)
if (axis is None):
result = _cihs_1D(data.compressed(), alpha)
else:
if data.ndim > 2:
raise ValueError("Array 'data' must be at most two dimensional, "
"but got data.ndim = %d" % data.ndim)
result = ma.apply_along_axis(_cihs_1D, axis, data, alpha)
return result
0
Example 80
Project: TSTools Source File: attach_md.py
@QtCore.pyqtSlot()
def add_metadata(self):
""" """
# Try to match metadata
ts_match_var = (tsm.ts.images['id'] if self.rad_ID.isChecked() is True
else tsm.ts.images['date'])
# Match column
match_col = self.table_metadata.selectedItems()[0].column()
md_match_var = self.md[:, match_col]
# Try to match
if len(ts_match_var) != len(md_match_var):
msg = 'Wrong number of elements to match ({t} vs. {m})'.format(
t=len(ts_match_var), m=len(md_match_var))
qgis_log(msg, logging.WARNING, 5)
return
if not np.all(np.sort(ts_match_var) == np.sort(md_match_var)):
msg = 'Not all elements match'
qgis_log(msg, logging.WARNING, 5)
return
# Perform match
match_ind = []
for i in xrange(len(ts_match_var)):
ind = np.where(md_match_var == ts_match_var[i])[0]
if len(ind) > 1:
msg = 'Multiple index matches for {m}'.format(
m=ts_match_var[i])
qgis_log(msg, logging.WARNING, 5)
return
match_ind.append(ind[0])
match_ind = np.array(match_ind)
# Sort
self.md_sorted = self.md[match_ind, :]
# Add to timeseries
for i, md in enumerate(self.colnames):
# Ignore match column
if i == match_col:
continue
if md not in tsm.ts.metadata and not hasattr(tsm.ts, md):
tsm.ts.metadata.append(md)
tsm.ts.metadata_str.append(md)
setattr(tsm.ts, md, self.md_sorted[:, i])
else:
msg = 'TS already has metadata item {m}'.format(m=md)
qgis_log(msg, logging.WARNING, 5)
# Emit
self.metadata_attached.emit()
0
Example 81
Project: PYPOWER Source File: printpf.py
def printpf(baseMVA, bus=None, gen=None, branch=None, f=None, success=None,
et=None, fd=None, ppopt=None):
"""Prints power flow results.
Prints power flow and optimal power flow results to C{fd} (a file
descriptor which defaults to C{stdout}), with the details of what
gets printed controlled by the optional C{ppopt} argument, which is a
PYPOWER options vector (see L{ppoption} for details).
The data can either be supplied in a single C{results} dict, or
in the individual arguments: C{baseMVA}, C{bus}, C{gen}, C{branch}, C{f},
C{success} and C{et}, where C{f} is the OPF objective function value,
C{success} is C{True} if the solution converged and C{False} otherwise,
and C{et} is the elapsed time for the computation in seconds. If C{f} is
given, it is assumed that the output is from an OPF run, otherwise it is
assumed to be a simple power flow run.
Examples::
ppopt = ppoptions(OUT_GEN=1, OUT_BUS=0, OUT_BRANCH=0)
fd = open(fname, 'w+b')
results = runopf(ppc)
printpf(results)
printpf(results, fd)
printpf(results, fd, ppopt)
printpf(baseMVA, bus, gen, branch, f, success, et)
printpf(baseMVA, bus, gen, branch, f, success, et, fd)
printpf(baseMVA, bus, gen, branch, f, success, et, fd, ppopt)
fd.close()
@author: Ray Zimmerman (PSERC Cornell)
"""
##----- initialization -----
## default arguments
if isinstance(baseMVA, dict):
have_results_struct = 1
results = baseMVA
if gen is None:
ppopt = ppoption() ## use default options
else:
ppopt = gen
if (ppopt['OUT_ALL'] == 0):
return ## nothin' to see here, bail out now
if bus is None:
fd = stdout ## print to stdout by default
else:
fd = bus
baseMVA, bus, gen, branch, success, et = \
results["baseMVA"], results["bus"], results["gen"], \
results["branch"], results["success"], results["et"]
if 'f' in results:
f = results["f"]
else:
f = None
else:
have_results_struct = 0
if ppopt is None:
ppopt = ppoption() ## use default options
if fd is None:
fd = stdout ## print to stdout by default
if ppopt['OUT_ALL'] == 0:
return ## nothin' to see here, bail out now
isOPF = f is not None ## FALSE -> only simple PF data, TRUE -> OPF data
## options
isDC = ppopt['PF_DC'] ## use DC formulation?
OUT_ALL = ppopt['OUT_ALL']
OUT_ANY = OUT_ALL == 1 ## set to true if any pretty output is to be generated
OUT_SYS_SUM = (OUT_ALL == 1) or ((OUT_ALL == -1) and ppopt['OUT_SYS_SUM'])
OUT_AREA_SUM = (OUT_ALL == 1) or ((OUT_ALL == -1) and ppopt['OUT_AREA_SUM'])
OUT_BUS = (OUT_ALL == 1) or ((OUT_ALL == -1) and ppopt['OUT_BUS'])
OUT_BRANCH = (OUT_ALL == 1) or ((OUT_ALL == -1) and ppopt['OUT_BRANCH'])
OUT_GEN = (OUT_ALL == 1) or ((OUT_ALL == -1) and ppopt['OUT_GEN'])
OUT_ANY = OUT_ANY | ((OUT_ALL == -1) and
(OUT_SYS_SUM or OUT_AREA_SUM or OUT_BUS or
OUT_BRANCH or OUT_GEN))
if OUT_ALL == -1:
OUT_ALL_LIM = ppopt['OUT_ALL_LIM']
elif OUT_ALL == 1:
OUT_ALL_LIM = 2
else:
OUT_ALL_LIM = 0
OUT_ANY = OUT_ANY or (OUT_ALL_LIM >= 1)
if OUT_ALL_LIM == -1:
OUT_V_LIM = ppopt['OUT_V_LIM']
OUT_LINE_LIM = ppopt['OUT_LINE_LIM']
OUT_PG_LIM = ppopt['OUT_PG_LIM']
OUT_QG_LIM = ppopt['OUT_QG_LIM']
else:
OUT_V_LIM = OUT_ALL_LIM
OUT_LINE_LIM = OUT_ALL_LIM
OUT_PG_LIM = OUT_ALL_LIM
OUT_QG_LIM = OUT_ALL_LIM
OUT_ANY = OUT_ANY or ((OUT_ALL_LIM == -1) and (OUT_V_LIM or OUT_LINE_LIM or OUT_PG_LIM or OUT_QG_LIM))
ptol = 1e-4 ## tolerance for displaying shadow prices
## create map of external bus numbers to bus indices
i2e = bus[:, BUS_I].astype(int)
e2i = zeros(max(i2e) + 1, int)
e2i[i2e] = arange(bus.shape[0])
## sizes of things
nb = bus.shape[0] ## number of buses
nl = branch.shape[0] ## number of branches
ng = gen.shape[0] ## number of generators
## zero out some data to make printout consistent for DC case
if isDC:
bus[:, r_[QD, BS]] = zeros((nb, 2))
gen[:, r_[QG, QMAX, QMIN]] = zeros((ng, 3))
branch[:, r_[BR_R, BR_B]] = zeros((nl, 2))
## parameters
ties = find(bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] !=
bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA])
## area inter-ties
tap = ones(nl) ## default tap ratio = 1 for lines
xfmr = find(branch[:, TAP]) ## indices of transformers
tap[xfmr] = branch[xfmr, TAP] ## include transformer tap ratios
tap = tap * exp(1j * pi / 180 * branch[:, SHIFT]) ## add phase shifters
nzld = find((bus[:, PD] != 0.0) | (bus[:, QD] != 0.0))
sorted_areas = sort(bus[:, BUS_AREA])
## area numbers
s_areas = sorted_areas[r_[1, find(diff(sorted_areas)) + 1]]
nzsh = find((bus[:, GS] != 0.0) | (bus[:, BS] != 0.0))
allg = find( ~isload(gen) )
ong = find( (gen[:, GEN_STATUS] > 0) & ~isload(gen) )
onld = find( (gen[:, GEN_STATUS] > 0) & isload(gen) )
V = bus[:, VM] * exp(-1j * pi / 180 * bus[:, VA])
out = find(branch[:, BR_STATUS] == 0) ## out-of-service branches
nout = len(out)
if isDC:
loss = zeros(nl)
else:
loss = baseMVA * abs(V[e2i[ branch[:, F_BUS].astype(int) ]] / tap -
V[e2i[ branch[:, T_BUS].astype(int) ]])**2 / \
(branch[:, BR_R] - 1j * branch[:, BR_X])
fchg = abs(V[e2i[ branch[:, F_BUS].astype(int) ]] / tap)**2 * branch[:, BR_B] * baseMVA / 2
tchg = abs(V[e2i[ branch[:, T_BUS].astype(int) ]] )**2 * branch[:, BR_B] * baseMVA / 2
loss[out] = zeros(nout)
fchg[out] = zeros(nout)
tchg[out] = zeros(nout)
##----- print the stuff -----
if OUT_ANY:
## convergence & elapsed time
if success:
fd.write('\nConverged in %.2f seconds' % et)
else:
fd.write('\nDid not converge (%.2f seconds)\n' % et)
## objective function value
if isOPF:
fd.write('\nObjective Function Value = %.2f $/hr' % f)
if OUT_SYS_SUM:
fd.write('\n================================================================================')
fd.write('\n| System Summary |')
fd.write('\n================================================================================')
fd.write('\n\nHow many? How much? P (MW) Q (MVAr)')
fd.write('\n--------------------- ------------------- ------------- -----------------')
fd.write('\nBuses %6d Total Gen Capacity %7.1f %7.1f to %.1f' % (nb, sum(gen[allg, PMAX]), sum(gen[allg, QMIN]), sum(gen[allg, QMAX])))
fd.write('\nGenerators %5d On-line Capacity %7.1f %7.1f to %.1f' % (len(allg), sum(gen[ong, PMAX]), sum(gen[ong, QMIN]), sum(gen[ong, QMAX])))
fd.write('\nCommitted Gens %5d Generation (actual) %7.1f %7.1f' % (len(ong), sum(gen[ong, PG]), sum(gen[ong, QG])))
fd.write('\nLoads %5d Load %7.1f %7.1f' % (len(nzld)+len(onld), sum(bus[nzld, PD])-sum(gen[onld, PG]), sum(bus[nzld, QD])-sum(gen[onld, QG])))
fd.write('\n Fixed %5d Fixed %7.1f %7.1f' % (len(nzld), sum(bus[nzld, PD]), sum(bus[nzld, QD])))
fd.write('\n Dispatchable %5d Dispatchable %7.1f of %-7.1f%7.1f' % (len(onld), -sum(gen[onld, PG]), -sum(gen[onld, PMIN]), -sum(gen[onld, QG])))
fd.write('\nShunts %5d Shunt (inj) %7.1f %7.1f' % (len(nzsh),
-sum(bus[nzsh, VM]**2 * bus[nzsh, GS]), sum(bus[nzsh, VM]**2 * bus[nzsh, BS]) ))
fd.write('\nBranches %5d Losses (I^2 * Z) %8.2f %8.2f' % (nl, sum(loss.real), sum(loss.imag) ))
fd.write('\nTransformers %5d Branch Charging (inj) - %7.1f' % (len(xfmr), sum(fchg) + sum(tchg) ))
fd.write('\nInter-ties %5d Total Inter-tie Flow %7.1f %7.1f' % (len(ties), sum(abs(branch[ties, PF]-branch[ties, PT])) / 2, sum(abs(branch[ties, QF]-branch[ties, QT])) / 2))
fd.write('\nAreas %5d' % len(s_areas))
fd.write('\n')
fd.write('\n Minimum Maximum')
fd.write('\n ------------------------- --------------------------------')
minv = min(bus[:, VM])
mini = argmin(bus[:, VM])
maxv = max(bus[:, VM])
maxi = argmax(bus[:, VM])
fd.write('\nVoltage Magnitude %7.3f p.u. @ bus %-4d %7.3f p.u. @ bus %-4d' % (minv, bus[mini, BUS_I], maxv, bus[maxi, BUS_I]))
minv = min(bus[:, VA])
mini = argmin(bus[:, VA])
maxv = max(bus[:, VA])
maxi = argmax(bus[:, VA])
fd.write('\nVoltage Angle %8.2f deg @ bus %-4d %8.2f deg @ bus %-4d' % (minv, bus[mini, BUS_I], maxv, bus[maxi, BUS_I]))
if not isDC:
maxv = max(loss.real)
maxi = argmax(loss.real)
fd.write('\nP Losses (I^2*R) - %8.2f MW @ line %d-%d' % (maxv, branch[maxi, F_BUS], branch[maxi, T_BUS]))
maxv = max(loss.imag)
maxi = argmax(loss.imag)
fd.write('\nQ Losses (I^2*X) - %8.2f MVAr @ line %d-%d' % (maxv, branch[maxi, F_BUS], branch[maxi, T_BUS]))
if isOPF:
minv = min(bus[:, LAM_P])
mini = argmin(bus[:, LAM_P])
maxv = max(bus[:, LAM_P])
maxi = argmax(bus[:, LAM_P])
fd.write('\nLambda P %8.2f $/MWh @ bus %-4d %8.2f $/MWh @ bus %-4d' % (minv, bus[mini, BUS_I], maxv, bus[maxi, BUS_I]))
minv = min(bus[:, LAM_Q])
mini = argmin(bus[:, LAM_Q])
maxv = max(bus[:, LAM_Q])
maxi = argmax(bus[:, LAM_Q])
fd.write('\nLambda Q %8.2f $/MWh @ bus %-4d %8.2f $/MWh @ bus %-4d' % (minv, bus[mini, BUS_I], maxv, bus[maxi, BUS_I]))
fd.write('\n')
if OUT_AREA_SUM:
fd.write('\n================================================================================')
fd.write('\n| Area Summary |')
fd.write('\n================================================================================')
fd.write('\nArea # of # of Gens # of Loads # of # of # of # of')
fd.write('\n Num Buses Total Online Total Fixed Disp Shunt Brchs Xfmrs Ties')
fd.write('\n---- ----- ----- ------ ----- ----- ----- ----- ----- ----- -----')
for i in range(len(s_areas)):
a = s_areas[i]
ib = find(bus[:, BUS_AREA] == a)
ig = find((bus[e2i[gen[:, GEN_BUS].astype(int)], BUS_AREA] == a) & ~isload(gen))
igon = find((bus[e2i[gen[:, GEN_BUS].astype(int)], BUS_AREA] == a) & (gen[:, GEN_STATUS] > 0) & ~isload(gen))
ildon = find((bus[e2i[gen[:, GEN_BUS].astype(int)], BUS_AREA] == a) & (gen[:, GEN_STATUS] > 0) & isload(gen))
inzld = find((bus[:, BUS_AREA] == a) & logical_or(bus[:, PD], bus[:, QD]))
inzsh = find((bus[:, BUS_AREA] == a) & logical_or(bus[:, GS], bus[:, BS]))
ibrch = find((bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] == a) & (bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA] == a))
in_tie = find((bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] == a) & (bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA] != a))
out_tie = find((bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] != a) & (bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA] == a))
if not any(xfmr + 1):
nxfmr = 0
else:
nxfmr = len(find((bus[e2i[branch[xfmr, F_BUS].astype(int)], BUS_AREA] == a) & (bus[e2i[branch[xfmr, T_BUS].astype(int)], BUS_AREA] == a)))
fd.write('\n%3d %6d %5d %5d %5d %5d %5d %5d %5d %5d %5d' %
(a, len(ib), len(ig), len(igon), \
len(inzld)+len(ildon), len(inzld), len(ildon), \
len(inzsh), len(ibrch), nxfmr, len(in_tie)+len(out_tie)))
fd.write('\n---- ----- ----- ------ ----- ----- ----- ----- ----- ----- -----')
fd.write('\nTot: %6d %5d %5d %5d %5d %5d %5d %5d %5d %5d' %
(nb, len(allg), len(ong), len(nzld)+len(onld),
len(nzld), len(onld), len(nzsh), nl, len(xfmr), len(ties)))
fd.write('\n')
fd.write('\nArea Total Gen Capacity On-line Gen Capacity Generation')
fd.write('\n Num MW MVAr MW MVAr MW MVAr')
fd.write('\n---- ------ ------------------ ------ ------------------ ------ ------')
for i in range(len(s_areas)):
a = s_areas[i]
ig = find((bus[e2i[gen[:, GEN_BUS].astype(int)], BUS_AREA] == a) & ~isload(gen))
igon = find((bus[e2i[gen[:, GEN_BUS].astype(int)], BUS_AREA] == a) & (gen[:, GEN_STATUS] > 0) & ~isload(gen))
fd.write('\n%3d %7.1f %7.1f to %-7.1f %7.1f %7.1f to %-7.1f %7.1f %7.1f' %
(a, sum(gen[ig, PMAX]), sum(gen[ig, QMIN]), sum(gen[ig, QMAX]),
sum(gen[igon, PMAX]), sum(gen[igon, QMIN]), sum(gen[igon, QMAX]),
sum(gen[igon, PG]), sum(gen[igon, QG]) ))
fd.write('\n---- ------ ------------------ ------ ------------------ ------ ------')
fd.write('\nTot: %7.1f %7.1f to %-7.1f %7.1f %7.1f to %-7.1f %7.1f %7.1f' %
(sum(gen[allg, PMAX]), sum(gen[allg, QMIN]), sum(gen[allg, QMAX]),
sum(gen[ong, PMAX]), sum(gen[ong, QMIN]), sum(gen[ong, QMAX]),
sum(gen[ong, PG]), sum(gen[ong, QG]) ))
fd.write('\n')
fd.write('\nArea Disp Load Cap Disp Load Fixed Load Total Load')
fd.write('\n Num MW MVAr MW MVAr MW MVAr MW MVAr')
fd.write('\n---- ------ ------ ------ ------ ------ ------ ------ ------')
Qlim = (gen[:, QMIN] == 0) * gen[:, QMAX] + (gen[:, QMAX] == 0) * gen[:, QMIN]
for i in range(len(s_areas)):
a = s_areas[i]
ildon = find((bus[e2i[gen[:, GEN_BUS].astype(int)], BUS_AREA] == a) & (gen[:, GEN_STATUS] > 0) & isload(gen))
inzld = find((bus[:, BUS_AREA] == a) & logical_or(bus[:, PD], bus[:, QD]))
fd.write('\n%3d %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f' %
(a, -sum(gen[ildon, PMIN]),
-sum(Qlim[ildon]),
-sum(gen[ildon, PG]), -sum(gen[ildon, QG]),
sum(bus[inzld, PD]), sum(bus[inzld, QD]),
-sum(gen[ildon, PG]) + sum(bus[inzld, PD]),
-sum(gen[ildon, QG]) + sum(bus[inzld, QD]) ))
fd.write('\n---- ------ ------ ------ ------ ------ ------ ------ ------')
fd.write('\nTot: %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f' %
(-sum(gen[onld, PMIN]),
-sum(Qlim[onld]),
-sum(gen[onld, PG]), -sum(gen[onld, QG]),
sum(bus[nzld, PD]), sum(bus[nzld, QD]),
-sum(gen[onld, PG]) + sum(bus[nzld, PD]),
-sum(gen[onld, QG]) + sum(bus[nzld, QD])) )
fd.write('\n')
fd.write('\nArea Shunt Inj Branch Series Losses Net Export')
fd.write('\n Num MW MVAr Charging MW MVAr MW MVAr')
fd.write('\n---- ------ ------ -------- ------ ------ ------ ------')
for i in range(len(s_areas)):
a = s_areas[i]
inzsh = find((bus[:, BUS_AREA] == a) & logical_or(bus[:, GS], bus[:, BS]))
ibrch = find((bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] == a) & (bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA] == a) & branch[:, BR_STATUS].astype(bool))
in_tie = find((bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] != a) & (bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA] == a) & branch[:, BR_STATUS].astype(bool))
out_tie = find((bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] == a) & (bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA] != a) & branch[:, BR_STATUS].astype(bool))
fd.write('\n%3d %7.1f %7.1f %7.1f %7.2f %7.2f %7.1f %7.1f' %
(a, -sum(bus[inzsh, VM]**2 * bus[inzsh, GS]),
sum(bus[inzsh, VM]**2 * bus[inzsh, BS]),
sum(fchg[ibrch]) + sum(tchg[ibrch]) + sum(fchg[out_tie]) + sum(tchg[in_tie]),
sum(real(loss[ibrch])) + sum(real(loss[r_[in_tie, out_tie]])) / 2,
sum(imag(loss[ibrch])) + sum(imag(loss[r_[in_tie, out_tie]])) / 2,
sum(branch[in_tie, PT])+sum(branch[out_tie, PF]) - sum(real(loss[r_[in_tie, out_tie]])) / 2,
sum(branch[in_tie, QT])+sum(branch[out_tie, QF]) - sum(imag(loss[r_[in_tie, out_tie]])) / 2 ))
fd.write('\n---- ------ ------ -------- ------ ------ ------ ------')
fd.write('\nTot: %7.1f %7.1f %7.1f %7.2f %7.2f - -' %
(-sum(bus[nzsh, VM]**2 * bus[nzsh, GS]),
sum(bus[nzsh, VM]**2 * bus[nzsh, BS]),
sum(fchg) + sum(tchg), sum(real(loss)), sum(imag(loss)) ))
fd.write('\n')
## generator data
if OUT_GEN:
if isOPF:
genlamP = bus[e2i[gen[:, GEN_BUS].astype(int)], LAM_P]
genlamQ = bus[e2i[gen[:, GEN_BUS].astype(int)], LAM_Q]
fd.write('\n================================================================================')
fd.write('\n| Generator Data |')
fd.write('\n================================================================================')
fd.write('\n Gen Bus Status Pg Qg ')
if isOPF: fd.write(' Lambda ($/MVA-hr)')
fd.write('\n # # (MW) (MVAr) ')
if isOPF: fd.write(' P Q ')
fd.write('\n---- ----- ------ -------- --------')
if isOPF: fd.write(' -------- --------')
for k in range(len(ong)):
i = ong[k]
fd.write('\n%3d %6d %2d ' % (i, gen[i, GEN_BUS], gen[i, GEN_STATUS]))
if (gen[i, GEN_STATUS] > 0) & logical_or(gen[i, PG], gen[i, QG]):
fd.write('%10.2f%10.2f' % (gen[i, PG], gen[i, QG]))
else:
fd.write(' - - ')
if isOPF: fd.write('%10.2f%10.2f' % (genlamP[i], genlamQ[i]))
fd.write('\n -------- --------')
fd.write('\n Total: %9.2f%10.2f' % (sum(gen[ong, PG]), sum(gen[ong, QG])))
fd.write('\n')
if any(onld + 1):
fd.write('\n================================================================================')
fd.write('\n| Dispatchable Load Data |')
fd.write('\n================================================================================')
fd.write('\n Gen Bus Status Pd Qd ')
if isOPF: fd.write(' Lambda ($/MVA-hr)')
fd.write('\n # # (MW) (MVAr) ')
if isOPF: fd.write(' P Q ')
fd.write('\n---- ----- ------ -------- --------')
if isOPF: fd.write(' -------- --------')
for k in range(len(onld)):
i = onld[k]
fd.write('\n%3d %6d %2d ' % (i, gen[i, GEN_BUS], gen[i, GEN_STATUS]))
if (gen[i, GEN_STATUS] > 0) & logical_or(gen[i, PG], gen[i, QG]):
fd.write('%10.2f%10.2f' % (-gen[i, PG], -gen[i, QG]))
else:
fd.write(' - - ')
if isOPF: fd.write('%10.2f%10.2f' % (genlamP[i], genlamQ[i]))
fd.write('\n -------- --------')
fd.write('\n Total: %9.2f%10.2f' % (-sum(gen[onld, PG]), -sum(gen[onld, QG])))
fd.write('\n')
## bus data
if OUT_BUS:
fd.write('\n================================================================================')
fd.write('\n| Bus Data |')
fd.write('\n================================================================================')
fd.write('\n Bus Voltage Generation Load ')
if isOPF: fd.write(' Lambda($/MVA-hr)')
fd.write('\n # Mag(pu) Ang(deg) P (MW) Q (MVAr) P (MW) Q (MVAr)')
if isOPF: fd.write(' P Q ')
fd.write('\n----- ------- -------- -------- -------- -------- --------')
if isOPF: fd.write(' ------- -------')
for i in range(nb):
fd.write('\n%5d%7.3f%9.3f' % tuple(bus[i, [BUS_I, VM, VA]]))
if bus[i, BUS_TYPE] == REF:
fd.write('*')
else:
fd.write(' ')
g = find((gen[:, GEN_STATUS] > 0) & (gen[:, GEN_BUS] == bus[i, BUS_I]) &
~isload(gen))
ld = find((gen[:, GEN_STATUS] > 0) & (gen[:, GEN_BUS] == bus[i, BUS_I]) &
isload(gen))
if any(g + 1):
fd.write('%9.2f%10.2f' % (sum(gen[g, PG]), sum(gen[g, QG])))
else:
fd.write(' - - ')
if logical_or(bus[i, PD], bus[i, QD]) | any(ld + 1):
if any(ld + 1):
fd.write('%10.2f*%9.2f*' % (bus[i, PD] - sum(gen[ld, PG]),
bus[i, QD] - sum(gen[ld, QG])))
else:
fd.write('%10.2f%10.2f ' % tuple(bus[i, [PD, QD]]))
else:
fd.write(' - - ')
if isOPF:
fd.write('%9.3f' % bus[i, LAM_P])
if abs(bus[i, LAM_Q]) > ptol:
fd.write('%8.3f' % bus[i, LAM_Q])
else:
fd.write(' -')
fd.write('\n -------- -------- -------- --------')
fd.write('\n Total: %9.2f %9.2f %9.2f %9.2f' %
(sum(gen[ong, PG]), sum(gen[ong, QG]),
sum(bus[nzld, PD]) - sum(gen[onld, PG]),
sum(bus[nzld, QD]) - sum(gen[onld, QG])))
fd.write('\n')
## branch data
if OUT_BRANCH:
fd.write('\n================================================================================')
fd.write('\n| Branch Data |')
fd.write('\n================================================================================')
fd.write('\nBrnch From To From Bus Injection To Bus Injection Loss (I^2 * Z) ')
fd.write('\n # Bus Bus P (MW) Q (MVAr) P (MW) Q (MVAr) P (MW) Q (MVAr)')
fd.write('\n----- ----- ----- -------- -------- -------- -------- -------- --------')
for i in range(nl):
fd.write('\n%4d%7d%7d%10.2f%10.2f%10.2f%10.2f%10.3f%10.2f' %
(i, branch[i, F_BUS], branch[i, T_BUS],
branch[i, PF], branch[i, QF], branch[i, PT], branch[i, QT],
loss[i].real, loss[i].imag))
fd.write('\n -------- --------')
fd.write('\n Total:%10.3f%10.2f' %
(sum(real(loss)), sum(imag(loss))))
fd.write('\n')
##----- constraint data -----
if isOPF:
ctol = ppopt['OPF_VIOLATION'] ## constraint violation tolerance
## voltage constraints
if (not isDC) & (OUT_V_LIM == 2 | (OUT_V_LIM == 1 &
(any(bus[:, VM] < bus[:, VMIN] + ctol) |
any(bus[:, VM] > bus[:, VMAX] - ctol) |
any(bus[:, MU_VMIN] > ptol) |
any(bus[:, MU_VMAX] > ptol)))):
fd.write('\n================================================================================')
fd.write('\n| Voltage Constraints |')
fd.write('\n================================================================================')
fd.write('\nBus # Vmin mu Vmin |V| Vmax Vmax mu')
fd.write('\n----- -------- ----- ----- ----- --------')
for i in range(nb):
if (OUT_V_LIM == 2) | (OUT_V_LIM == 1 &
((bus[i, VM] < bus[i, VMIN] + ctol) |
(bus[i, VM] > bus[i, VMAX] - ctol) |
(bus[i, MU_VMIN] > ptol) |
(bus[i, MU_VMAX] > ptol))):
fd.write('\n%5d' % bus[i, BUS_I])
if ((bus[i, VM] < bus[i, VMIN] + ctol) |
(bus[i, MU_VMIN] > ptol)):
fd.write('%10.3f' % bus[i, MU_VMIN])
else:
fd.write(' - ')
fd.write('%8.3f%7.3f%7.3f' % tuple(bus[i, [VMIN, VM, VMAX]]))
if (bus[i, VM] > bus[i, VMAX] - ctol) | (bus[i, MU_VMAX] > ptol):
fd.write('%10.3f' % bus[i, MU_VMAX])
else:
fd.write(' - ')
fd.write('\n')
## generator P constraints
if (OUT_PG_LIM == 2) | \
((OUT_PG_LIM == 1) & (any(gen[ong, PG] < gen[ong, PMIN] + ctol) |
any(gen[ong, PG] > gen[ong, PMAX] - ctol) |
any(gen[ong, MU_PMIN] > ptol) |
any(gen[ong, MU_PMAX] > ptol))) | \
((not isDC) & ((OUT_QG_LIM == 2) |
((OUT_QG_LIM == 1) & (any(gen[ong, QG] < gen[ong, QMIN] + ctol) |
any(gen[ong, QG] > gen[ong, QMAX] - ctol) |
any(gen[ong, MU_QMIN] > ptol) |
any(gen[ong, MU_QMAX] > ptol))))):
fd.write('\n================================================================================')
fd.write('\n| Generation Constraints |')
fd.write('\n================================================================================')
if (OUT_PG_LIM == 2) | ((OUT_PG_LIM == 1) &
(any(gen[ong, PG] < gen[ong, PMIN] + ctol) |
any(gen[ong, PG] > gen[ong, PMAX] - ctol) |
any(gen[ong, MU_PMIN] > ptol) |
any(gen[ong, MU_PMAX] > ptol))):
fd.write('\n Gen Bus Active Power Limits')
fd.write('\n # # Pmin mu Pmin Pg Pmax Pmax mu')
fd.write('\n---- ----- ------- -------- -------- -------- -------')
for k in range(len(ong)):
i = ong[k]
if (OUT_PG_LIM == 2) | ((OUT_PG_LIM == 1) &
((gen[i, PG] < gen[i, PMIN] + ctol) |
(gen[i, PG] > gen[i, PMAX] - ctol) |
(gen[i, MU_PMIN] > ptol) | (gen[i, MU_PMAX] > ptol))):
fd.write('\n%4d%6d ' % (i, gen[i, GEN_BUS]))
if (gen[i, PG] < gen[i, PMIN] + ctol) | (gen[i, MU_PMIN] > ptol):
fd.write('%8.3f' % gen[i, MU_PMIN])
else:
fd.write(' - ')
if gen[i, PG]:
fd.write('%10.2f%10.2f%10.2f' % tuple(gen[i, [PMIN, PG, PMAX]]))
else:
fd.write('%10.2f - %10.2f' % tuple(gen[i, [PMIN, PMAX]]))
if (gen[i, PG] > gen[i, PMAX] - ctol) | (gen[i, MU_PMAX] > ptol):
fd.write('%9.3f' % gen[i, MU_PMAX])
else:
fd.write(' - ')
fd.write('\n')
## generator Q constraints
if (not isDC) & ((OUT_QG_LIM == 2) | ((OUT_QG_LIM == 1) &
(any(gen[ong, QG] < gen[ong, QMIN] + ctol) |
any(gen[ong, QG] > gen[ong, QMAX] - ctol) |
any(gen[ong, MU_QMIN] > ptol) |
any(gen[ong, MU_QMAX] > ptol)))):
fd.write('\nGen Bus Reactive Power Limits')
fd.write('\n # # Qmin mu Qmin Qg Qmax Qmax mu')
fd.write('\n--- --- ------- -------- -------- -------- -------')
for k in range(len(ong)):
i = ong[k]
if (OUT_QG_LIM == 2) | ((OUT_QG_LIM == 1) &
((gen[i, QG] < gen[i, QMIN] + ctol) |
(gen[i, QG] > gen[i, QMAX] - ctol) |
(gen[i, MU_QMIN] > ptol) |
(gen[i, MU_QMAX] > ptol))):
fd.write('\n%3d%5d' % (i, gen[i, GEN_BUS]))
if (gen[i, QG] < gen[i, QMIN] + ctol) | (gen[i, MU_QMIN] > ptol):
fd.write('%8.3f' % gen[i, MU_QMIN])
else:
fd.write(' - ')
if gen[i, QG]:
fd.write('%10.2f%10.2f%10.2f' % tuple(gen[i, [QMIN, QG, QMAX]]))
else:
fd.write('%10.2f - %10.2f' % tuple(gen[i, [QMIN, QMAX]]))
if (gen[i, QG] > gen[i, QMAX] - ctol) | (gen[i, MU_QMAX] > ptol):
fd.write('%9.3f' % gen[i, MU_QMAX])
else:
fd.write(' - ')
fd.write('\n')
## dispatchable load P constraints
if (OUT_PG_LIM == 2) | (OUT_QG_LIM == 2) | \
((OUT_PG_LIM == 1) & (any(gen[onld, PG] < gen[onld, PMIN] + ctol) |
any(gen[onld, PG] > gen[onld, PMAX] - ctol) |
any(gen[onld, MU_PMIN] > ptol) |
any(gen[onld, MU_PMAX] > ptol))) | \
((OUT_QG_LIM == 1) & (any(gen[onld, QG] < gen[onld, QMIN] + ctol) |
any(gen[onld, QG] > gen[onld, QMAX] - ctol) |
any(gen[onld, MU_QMIN] > ptol) |
any(gen[onld, MU_QMAX] > ptol))):
fd.write('\n================================================================================')
fd.write('\n| Dispatchable Load Constraints |')
fd.write('\n================================================================================')
if (OUT_PG_LIM == 2) | ((OUT_PG_LIM == 1) &
(any(gen[onld, PG] < gen[onld, PMIN] + ctol) |
any(gen[onld, PG] > gen[onld, PMAX] - ctol) |
any(gen[onld, MU_PMIN] > ptol) |
any(gen[onld, MU_PMAX] > ptol))):
fd.write('\nGen Bus Active Power Limits')
fd.write('\n # # Pmin mu Pmin Pg Pmax Pmax mu')
fd.write('\n--- --- ------- -------- -------- -------- -------')
for k in range(len(onld)):
i = onld[k]
if (OUT_PG_LIM == 2) | ((OUT_PG_LIM == 1) &
((gen[i, PG] < gen[i, PMIN] + ctol) |
(gen[i, PG] > gen[i, PMAX] - ctol) |
(gen[i, MU_PMIN] > ptol) |
(gen[i, MU_PMAX] > ptol))):
fd.write('\n%3d%5d' % (i, gen[i, GEN_BUS]))
if (gen[i, PG] < gen[i, PMIN] + ctol) | (gen[i, MU_PMIN] > ptol):
fd.write('%8.3f' % gen[i, MU_PMIN])
else:
fd.write(' - ')
if gen[i, PG]:
fd.write('%10.2f%10.2f%10.2f' % gen[i, [PMIN, PG, PMAX]])
else:
fd.write('%10.2f - %10.2f' % gen[i, [PMIN, PMAX]])
if (gen[i, PG] > gen[i, PMAX] - ctol) | (gen[i, MU_PMAX] > ptol):
fd.write('%9.3f' % gen[i, MU_PMAX])
else:
fd.write(' - ')
fd.write('\n')
## dispatchable load Q constraints
if (not isDC) & ((OUT_QG_LIM == 2) | ((OUT_QG_LIM == 1) &
(any(gen[onld, QG] < gen[onld, QMIN] + ctol) |
any(gen[onld, QG] > gen[onld, QMAX] - ctol) |
any(gen[onld, MU_QMIN] > ptol) |
any(gen[onld, MU_QMAX] > ptol)))):
fd.write('\nGen Bus Reactive Power Limits')
fd.write('\n # # Qmin mu Qmin Qg Qmax Qmax mu')
fd.write('\n--- --- ------- -------- -------- -------- -------')
for k in range(len(onld)):
i = onld[k]
if (OUT_QG_LIM == 2) | ((OUT_QG_LIM == 1) &
((gen[i, QG] < gen[i, QMIN] + ctol) |
(gen[i, QG] > gen[i, QMAX] - ctol) |
(gen[i, MU_QMIN] > ptol) |
(gen[i, MU_QMAX] > ptol))):
fd.write('\n%3d%5d' % (i, gen(i, GEN_BUS)))
if (gen[i, QG] < gen[i, QMIN] + ctol) | (gen[i, MU_QMIN] > ptol):
fd.write('%8.3f' % gen[i, MU_QMIN])
else:
fd.write(' - ')
if gen[i, QG]:
fd.write('%10.2f%10.2f%10.2f' % gen[i, [QMIN, QG, QMAX]])
else:
fd.write('%10.2f - %10.2f' % gen[i, [QMIN, QMAX]])
if (gen[i, QG] > gen[i, QMAX] - ctol) | (gen[i, MU_QMAX] > ptol):
fd.write('%9.3f' % gen[i, MU_QMAX])
else:
fd.write(' - ')
fd.write('\n')
## line flow constraints
if (ppopt['OPF_FLOW_LIM'] == 1) | isDC: ## P limit
Ff = branch[:, PF]
Ft = branch[:, PT]
strg = '\n # Bus Pf mu Pf |Pmax| Pt Pt mu Bus'
elif ppopt['OPF_FLOW_LIM'] == 2: ## |I| limit
Ff = abs( (branch[:, PF] + 1j * branch[:, QF]) / V[e2i[branch[:, F_BUS].astype(int)]] )
Ft = abs( (branch[:, PT] + 1j * branch[:, QT]) / V[e2i[branch[:, T_BUS].astype(int)]] )
strg = '\n # Bus |If| mu |If| |Imax| |It| |It| mu Bus'
else: ## |S| limit
Ff = abs(branch[:, PF] + 1j * branch[:, QF])
Ft = abs(branch[:, PT] + 1j * branch[:, QT])
strg = '\n # Bus |Sf| mu |Sf| |Smax| |St| |St| mu Bus'
if (OUT_LINE_LIM == 2) | ((OUT_LINE_LIM == 1) &
(any((branch[:, RATE_A] != 0) & (abs(Ff) > branch[:, RATE_A] - ctol)) |
any((branch[:, RATE_A] != 0) & (abs(Ft) > branch[:, RATE_A] - ctol)) |
any(branch[:, MU_SF] > ptol) |
any(branch[:, MU_ST] > ptol))):
fd.write('\n================================================================================')
fd.write('\n| Branch Flow Constraints |')
fd.write('\n================================================================================')
fd.write('\nBrnch From "From" End Limit "To" End To')
fd.write(strg)
fd.write('\n----- ----- ------- -------- -------- -------- ------- -----')
for i in range(nl):
if (OUT_LINE_LIM == 2) | ((OUT_LINE_LIM == 1) &
(((branch[i, RATE_A] != 0) & (abs(Ff[i]) > branch[i, RATE_A] - ctol)) |
((branch[i, RATE_A] != 0) & (abs(Ft[i]) > branch[i, RATE_A] - ctol)) |
(branch[i, MU_SF] > ptol) | (branch[i, MU_ST] > ptol))):
fd.write('\n%4d%7d' % (i, branch[i, F_BUS]))
if (Ff[i] > branch[i, RATE_A] - ctol) | (branch[i, MU_SF] > ptol):
fd.write('%10.3f' % branch[i, MU_SF])
else:
fd.write(' - ')
fd.write('%9.2f%10.2f%10.2f' %
(Ff[i], branch[i, RATE_A], Ft[i]))
if (Ft[i] > branch[i, RATE_A] - ctol) | (branch[i, MU_ST] > ptol):
fd.write('%10.3f' % branch[i, MU_ST])
else:
fd.write(' - ')
fd.write('%6d' % branch[i, T_BUS])
fd.write('\n')
## execute userfcn callbacks for 'printpf' stage
if have_results_struct and 'userfcn' in results:
if not isOPF: ## turn off option for all constraints if it isn't an OPF
ppopt = ppoption(ppopt, 'OUT_ALL_LIM', 0)
run_userfcn(results["userfcn"], 'printpf', results, fd, ppopt)
0
Example 82
Project: brut Source File: bootstrap.py
def bootstrap_negatives(dfs, neg, neg_old):
"""
Resample from a set of negative labels, emphasizng
examples that are currently mis-classified
Parameters
----------
dfs : dict
Decision functions for 'neg', 'cv_pos', and 'neg_old'
neg : list
A set of negative examples
neg_old : list
The negative examples used to train the classifier
Returns
-------
A resampled list of negative values
Notes
-----
This function creates a new version of `neg` that retains
the half of the old examples with the highest decision function.
The other half are drawn randomly from the portion of `neg`
which have decision values >= the 10th percentile of the
decision function of `cv_pos`.
This results in a new set of training data that emphasizes
examples which are currently poorly-classified
"""
df = dfs['neg']
df_pos = dfs['cv_pos']
df_pos = np.sort(df_pos)
df_neg = dfs['neg_old']
#neg dfs > this are problematic
df_thresh = df_pos[.1 * df_pos.shape[0]]
ind = np.argsort(df)[::-1]
nnew = len(neg_old) / 2
nkeep = len(neg_old) - nnew
nresample = max(nnew, (df > df_thresh).sum())
old_ind = np.argsort(df_neg)[:nkeep]
new_ind = np.random.randint(0, nresample - 1, nnew)
new_ind = ind[new_ind]
return [neg_old[i] for i in old_ind] + [neg[i] for i in new_ind]
0
Example 83
Project: RCN Source File: draw_points_denoise_finetune.py
def draw_points_raw(out_path, annotate_kpts=True, high_res=True, max_errors=True, sample_num=20,
plot_colored_kpt=False, indices_given=False, merge_sets=False, pickle_path=""):
"""
This method draw points on the test set, when there are two steps for face-detection and
downsampling. So, when the dataset is created, a face-detection is done and the data is downsampled to
orig_downsample size (this is done before pickling the data). Then, later,
at train-time, a second face detection is performed and the data is downsampled
to second_downsample size. In both down_sampling stages, the kpt_norm keeps a value in the range [0,1]. So, the relative
location is the same and just normalized in the downsampled case, compared to the detected face locations.
The process of finding the original keypoint locations is as follows:
for the 2nd face detection, multiply the normalized key-point locations in the bounding_box size
and add to those values to location of the top-left position of the bounding_box. This gives the
positions before the 2nd face detection. Then, normalize the
positions by the downsampling (img) size of the first face-detection to get them normalized. Finally,
multiply the normalized key-point locations in the bounding_box size of the first face detection
and add to those values the location of the top-left position of that bounding_box.
"""
# setting the conv params to the weights
tcdcn, params = create_TCDCN_obejct(pkl_param_file)
tcdcn.load_params(pkl_param_file)
tcdcn_cfNet, params_cfNet = None, None
cfNet_model = (tcdcn_cfNet, params_cfNet)
if params['paral_conv'] in [2, 5, 6] or params['denoise_conv'] in [1, 2]:
params['mask_MTFL'] = 0
params['mask_300W'] = 1
dataSet = '300W'
num_kpt = 68
elif params['paral_conv'] in [1, 3, 4]:
params['mask_300W'] = 0
params['mask_MTFL'] = 1
dataSet = 'MTFL'
num_kpt = 5
##############################################
# getting the 1st bounding box test set data #
##############################################
sets = get_data(**params)
Train, Valid, Test = sets
sets = OrderedDict()
train_set_x, train_set_y = Train[dataSet]
sets['train'] = OrderedDict()
sets['train']['X'] = train_set_x
sets['train']['Y'] = train_set_y
sets['train']['indices'] = []
sets['train']['name'] = '1_train'
valid_set_x, valid_set_y = Valid[dataSet]
sets['valid'] = OrderedDict()
sets['valid']['X'] = valid_set_x
sets['valid']['Y'] = valid_set_y
sets['valid']['indices'] = []
sets['valid']['name'] = '2_valid'
for test_set in np.sort(Test[dataSet].keys()):
test_set_x, test_set_y = Test[dataSet][test_set]
if merge_sets:
if not 'all_sets' in sets.keys():
sets['all_sets'] = OrderedDict()
sets['all_sets']['X'] = []
sets['all_sets']['Y'] = []
sets['all_sets']['name'] = 'all_sets'
sets['all_sets']['X'].extend(test_set_x)
sets['all_sets']['Y'].append(test_set_y)
else:
sets[test_set] = OrderedDict()
sets[test_set]['X'] = test_set_x
sets[test_set]['Y'] = test_set_y
sets[test_set]['name'] = '3_test_%s' %(test_set)
sets[test_set]['indices'] = []
if merge_sets:
sets['all_sets']['X'] = np.array(sets['all_sets']['X'])
sets['all_sets']['Y'] = append_orderedDict(sets['all_sets']['Y'])
set_names = ['all_sets']
else:
# set_names = Test[dataSet].keys()
set_names = sets.keys()
for sub_set in set_names:
set_x = sets[sub_set]['X']
set_y = sets[sub_set]['Y']
set = sets[sub_set]['name']
set_y_cp = deepcopy(set_y)
set_x_cp = deepcopy(set_x)
if max_errors:
indices, error_kpt_avg_all = eval_test_set(tcdcn, params, set_x, set_y, set,
sample_num, dataSet, cfNet_model)
set_x_indx, set_y_indx = get_indices(set_x_cp, set_y_cp, indices)
elif indices_given:
set_x_indx, set_y_indx = get_indices(set_x_cp, set_y_cp, sets[sub_set]['indices'])
elif pickle_path:
pickle_indices, all_vals = get_pickle_indices(pickle_path)
set_x_indx, set_y_indx = get_indices(set_x_cp, set_y_cp, pickle_indices)
else:
set_x_indx, set_y_indx = get_subsets(set_x_cp, set_y_cp, sample_num)
##############################################
# getting the 2nd bounding box test set data #
##############################################
# preprocessing data
rng_seed = np.random.RandomState(0)
scale_mul = 0.0
translate_mul = 0.0
target_dim = params['target_dim']
td = (target_dim, target_dim)
set_y_orig = deepcopy(set_y_indx)
set_x_orig = deepcopy(set_x_indx)
if dataSet == 'MTFL':
dist_ratio = 3.8/4.0
else:
dist_ratio = params['dist_ratio']
set_x2, set_y2 = preprocess_once(set_x_indx, set_y_indx, dist_ratio=dist_ratio, gray_scale=params['gray_scale'])
# using preprocess_iter to prepare set_x2 and set_y2 as the model expect it
set_x2, set_y2 = preprocess_iter(set_x2, set_y2, rng_seed, jitter=False, scale_mul=scale_mul,
translate_mul=translate_mul, target_dim=td, sanity=False, dset=dataSet,
use_lcn=params['use_lcn'], rotation_set=None)
# getting the keypoint results for the first images in the test set
# face_rect has a vector of four values: rect_start_x, rect_start_y,
# rect_width (rect_end_x - rect_start_x), rect_height (rect_end_y - rect_start_y)
set_rect_2nd = set_y2['face_rect']
#getting the rectangle size for the detected face
rect_size_2nd = set_rect_2nd[:, 2:]
# getting the start point of the ractangle
rect_start_point_2nd = set_rect_2nd[:, :2]
# broadcasting the rectangle size and the start point by the number of the keypoints
set_rect_size_2nd = np.tile(rect_size_2nd, (1, num_kpt))
set_rect_start_2nd = np.tile(rect_start_point_2nd, (1, num_kpt))
# the keypoint positions in the normalized format
set_kpt_norm = set_y2['kpt_norm']
# setting batch_size to one to avoid getting different errors for different batch sizes
batch_size = 1
num_batches = int(np.ceil(set_x2.shape[0]/float(batch_size)))
epoch_error_kpt_avg = []
kpt_conv = []
kpt_conv_cfNet = []
kpt_conv_struc = []
for index in np.arange(num_batches):
x_batch = set_x2[index * batch_size: (index + 1) * batch_size]
assert dataSet == '300W'
assert params['denoise_conv'] == 2:
if 'test' in set:
pad_ratio = set_y2['pad_ratio'][index * batch_size: (index + 1) * batch_size]
border_pixel = padRatio_to_pixels(pad_ratio, set_x2.shape[1])
kpt_conv_batch, kpt_conv_batch_cfNet, kpt_conv_batch_struc = tcdcn.get_keypoints_all(
x_batch, border_pixel, dropout=0)
else: # train and valid sets
kpt_conv_batch, kpt_conv_batch_cfNet, kpt_conv_batch_struc = \
tcdcn.get_keypoints_all_train(x_batch, dropout=0)
kpt_conv.extend(kpt_conv_batch)
kpt_conv_cfNet.extend(kpt_conv_batch_cfNet)
kpt_conv_struc.extend(kpt_conv_batch_struc)
kpt_conv_all = np.array(kpt_conv)
kpt_conv_cfNet = np.array(kpt_conv_cfNet)
kpt_conv_struc = np.array(kpt_conv_struc)
models = [kpt_conv_cfNet, kpt_conv_struc, kpt_conv_all]
model_names = ['cfNet', 'struc', 'all']
for (kpt_conv, model_name) in zip(models, model_names):
set_y_cp = deepcopy(set_y_orig)
set_x_cp = deepcopy(set_x_orig)
# using preprocessing with gray_scale=False and use_lcn=True to get set_x_show in RGB format
set_x_show, set_y_show = preprocess_once(set_x_cp, set_y_cp, dist_ratio=dist_ratio, gray_scale=False)
set_x_show, _ = preprocess_iter(set_x_show, set_y_show, rng_seed, jitter=False, scale_mul=scale_mul,
translate_mul=translate_mul, target_dim=td, sanity=False, dset=dataSet,
use_lcn=True, rotation_set=None)
if plot_colored_kpt:
high_res = True
if high_res:
kpt_conv_shifted = (kpt_conv / float(target_dim)) * set_rect_size_2nd + set_rect_start_2nd
kpt_conv_shifted = kpt_conv_shifted.astype(int)
kpt_true_shifted = set_kpt_norm * set_rect_size_2nd + set_rect_start_2nd
kpt_true_shifted = kpt_true_shifted.astype(int)
else:
kpt_conv_shifted = kpt_conv
kpt_true_shifted = set_kpt_norm * target_dim
n_samples = set_x_cp.shape[0]
for i in xrange(n_samples):
index = i+1
if high_res:
img = set_x_cp[i]
else:
img = set_x_show[i]
if plot_colored_kpt:
img2 = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
img = cv2.cvtColor(img2, cv2.COLOR_GRAY2BGR)
img_kpts = drawpoints(img, kpt_conv_shifted[i], kpt_true_shifted[i], magnify=False, plot_colored_kpt=plot_colored_kpt)
if annotate_kpts:
img_kpts = annotateKpts(img_kpts ,kpt_conv_shifted[i], kpt_true_shifted[i])
out_file_path = "%s/%s_%i_downsampled_%s.png" %(out_path, set, (i+1), model_name)
cv2.imwrite(out_file_path, img_kpts)
# tiling the images
print "tiling images"
img = cv2.imread(out_file_path)
row_size, col_size = img.shape[:2]
tile_images(n_samples, set, out_path, row_size, col_size, model_names)
0
Example 84
Project: qutip Source File: entropy.py
def concurrence(rho):
"""
Calculate the concurrence entanglement measure for a two-qubit state.
Parameters
----------
state : qobj
Ket, bra, or density matrix for a two-qubit state.
Returns
-------
concur : float
Concurrence
References
----------
.. [1] http://en.wikipedia.org/wiki/Concurrence_(quantum_computing)
"""
if rho.isket and rho.dims != [[2, 2], [1, 1]]:
raise Exception("Ket must be tensor product of two qubits.")
elif rho.isbra and rho.dims != [[1, 1], [2, 2]]:
raise Exception("Bra must be tensor product of two qubits.")
elif rho.isoper and rho.dims != [[2, 2], [2, 2]]:
raise Exception("Density matrix must be tensor product of two qubits.")
if rho.isket or rho.isbra:
rho = ket2dm(rho)
sysy = tensor(sigmay(), sigmay())
rho_tilde = (rho * sysy) * (rho.conj() * sysy)
evals = rho_tilde.eigenenergies()
# abs to avoid problems with sqrt for very small negative numbers
evals = abs(sort(real(evals)))
lsum = sqrt(evals[3]) - sqrt(evals[2]) - sqrt(evals[1]) - sqrt(evals[0])
return max(0, lsum)
0
Example 85
Project: spotpy Source File: analyser.py
def sort_like(results):
return np.sort(results,axis=0)
0
Example 86
Project: scikit-learn Source File: test_text.py
def test_vectorizer_unicode():
# tests that the count vectorizer works with cyrillic.
docuement = (
"\xd0\x9c\xd0\xb0\xd1\x88\xd0\xb8\xd0\xbd\xd0\xbd\xd0\xbe\xd0"
"\xb5 \xd0\xbe\xd0\xb1\xd1\x83\xd1\x87\xd0\xb5\xd0\xbd\xd0\xb8\xd0"
"\xb5 \xe2\x80\x94 \xd0\xbe\xd0\xb1\xd1\x88\xd0\xb8\xd1\x80\xd0\xbd"
"\xd1\x8b\xd0\xb9 \xd0\xbf\xd0\xbe\xd0\xb4\xd1\x80\xd0\xb0\xd0\xb7"
"\xd0\xb4\xd0\xb5\xd0\xbb \xd0\xb8\xd1\x81\xd0\xba\xd1\x83\xd1\x81"
"\xd1\x81\xd1\x82\xd0\xb2\xd0\xb5\xd0\xbd\xd0\xbd\xd0\xbe\xd0\xb3"
"\xd0\xbe \xd0\xb8\xd0\xbd\xd1\x82\xd0\xb5\xd0\xbb\xd0\xbb\xd0"
"\xb5\xd0\xba\xd1\x82\xd0\xb0, \xd0\xb8\xd0\xb7\xd1\x83\xd1\x87"
"\xd0\xb0\xd1\x8e\xd1\x89\xd0\xb8\xd0\xb9 \xd0\xbc\xd0\xb5\xd1\x82"
"\xd0\xbe\xd0\xb4\xd1\x8b \xd0\xbf\xd0\xbe\xd1\x81\xd1\x82\xd1\x80"
"\xd0\xbe\xd0\xb5\xd0\xbd\xd0\xb8\xd1\x8f \xd0\xb0\xd0\xbb\xd0\xb3"
"\xd0\xbe\xd1\x80\xd0\xb8\xd1\x82\xd0\xbc\xd0\xbe\xd0\xb2, \xd1\x81"
"\xd0\xbf\xd0\xbe\xd1\x81\xd0\xbe\xd0\xb1\xd0\xbd\xd1\x8b\xd1\x85 "
"\xd0\xbe\xd0\xb1\xd1\x83\xd1\x87\xd0\xb0\xd1\x82\xd1\x8c\xd1\x81\xd1"
"\x8f.")
vect = CountVectorizer()
X_counted = vect.fit_transform([docuement])
assert_equal(X_counted.shape, (1, 15))
vect = HashingVectorizer(norm=None, non_negative=True)
X_hashed = vect.transform([docuement])
assert_equal(X_hashed.shape, (1, 2 ** 20))
# No collisions on such a small dataset
assert_equal(X_counted.nnz, X_hashed.nnz)
# When norm is None and non_negative, the tokens are counted up to
# collisions
assert_array_equal(np.sort(X_counted.data), np.sort(X_hashed.data))
0
Example 87
Project: spotpy Source File: analyser.py
def plot_fast_sensitivity(results,likes=['mean'],like_indices=None,number_of_sensitiv_pars=10):
"""
Example, how to plot the sensitivity for every parameter of your result array, created with the FAST algorithm
:results: Expects an numpy array which should have as first axis an index "like" or "like1".
:type: array
:likes: Optional, header of your objectivefunction
:type: list
:like_indices: Optional, index of objectivefunction to base the sensitivity on, default=None first objectivefunction is taken
:type: int
:number_of_sensitiv_pars: Optional, this number of most sensitive parameters will be shown in the legend
:type: int
:return: Parameter names which are sensitive, Sensitivity indices for every parameter, Parameter names which are not sensitive
:rtype: Three lists
"""
import matplotlib.pyplot as plt
from matplotlib import colors
cnames=list(colors.cnames)
parnames=get_parameternames(results)
fig=plt.figure(figsize=(16,12))
all_names=[]
all_no_names=[]
for i in range(len(likes)):
ax = plt.subplot(len(likes),1,i+1)
if like_indices==None:
Si=get_sensitivity_of_fast(results['like'],parnames)
else:
Si=get_sensitivity_of_fast(results['like'+str(like_indices[i])],parnames)
names=[]
values=[]
no_names=[]
for j in range(len(list(Si.values())[1])):
if list(Si.values())[1][j]>=sorted(np.sort(list(Si.values())[1]),reverse=True)[number_of_sensitiv_pars]:
names.append(parnames[j])
values.append(list(Si.values())[1][j])
else:
no_names.append(parnames[j])
print(names)
ax.plot([sorted(np.sort(list(Si.values())[1]),reverse=True)[number_of_sensitiv_pars]]*len(list(Si.values())[1]),'r--')
#ax.bar(np.arange(0,len(Si.values()[1])),sorted(np.sort(Si.values()[1]),reverse=True),label=str(names))
ax.bar(np.arange(0,len(list(Si.values())[1])),list(Si.values())[1],label=str(names))
ax.set_ylim([0,1])
#ax.set_xlabel(names)
ax.set_ylabel(likes[i])
ax.legend()
all_names.append(names)
all_no_names.append(no_names)
return all_names,values,all_no_names
0
Example 88
Project: pyDOE Source File: build_regression_matrix.py
def build_regression_matrix(H, model, build=None):
"""
Build a regression matrix using a DOE matrix and a list of monomials.
Parameters
----------
H : 2d-array
model : str
build : bool-array
Returns
-------
R : 2d-array
"""
ListOfTokens = model.split(' ')
if H.shape[1]==1:
size_index = len(str(H.shape[0]))
else:
size_index = len(str(H.shape[1]))
if build is None:
build = [True]*len(ListOfTokens)
# Test if the vector has the wrong direction (lines instead of columns)
if H.shape[0]==1:
H = H.T
# Collect the list of monomials
Monom_Index = []
for i in range(len(ListOfTokens)):
if build[i]:
Monom_Index += [grep(ListOfTokens, 'x' + str(0)*(size_index - \
len(str(i))) + str(i))]
Monom_Index = -np.sort(-Monom_Index)
Monom_Index = np.unique(Monom_Index)
if H.shape[1]==1:
nb_var = H.shape[0] # vector "mode": the number of vars is equal to the number of lines of H
VectorMode = True
for i in range(nb_var):
for j in range(ListOfTokens.shape[0]):
ListOfTokens[j] = ListOfTokens[j].replace(
'x' + str(0)*(size_index - len(str(i))) + str(i),
'H(' + str(i) + ')')
else:
nb_var = H.shape[0] # matrix "mode": the number of vars is equal to the number of columns of H
VectorMode = False
for i in range(nb_var):
for j in range(ListOfTokens.shape[0]):
ListOfTokens[j] = ListOfTokens[j].replace(
'x' + str(0)*(size_index - len(str(i))) + str(i),
'H[i,' + str(i) + ')')
# Now build the regression matrix
if VectorMode:
R = np.zeros((len(ListOfTokens), 1))
for j in range(len(ListOfTokens)):
R[j, 0] = eval(ListOfTokens[j])
else:
R = np.zeros((H.shape[0], len(ListOfTokens)))
for i in range(H.shape[0]):
for j in range(len(ListOfTokens)):
R[i, j] = eval(ListOfTokens[j])
return R
0
Example 89
Project: rlpy Source File: GeneralTools.py
def make_colormap(colors):
"""
Define a new color map based on values specified in the dictionary
colors, where colors[z] is the color that value z should be mapped to,
with linear interpolation between the given values of z.
The z values (dictionary keys) are real numbers and the values
colors[z] can be either an RGB list, e.g. [1,0,0] for red, or an
html hex string, e.g. "#ff0000" for red.
"""
from matplotlib.colors import LinearSegmentedColormap, ColorConverter
z = np.sort(colors.keys())
n = len(z)
z1 = min(z)
zn = max(z)
x0 = (z - z1) / ((zn - z1) * 1.)
CC = ColorConverter()
R = []
G = []
B = []
for i in xrange(n):
# i'th color at level z[i]:
Ci = colors[z[i]]
if isinstance(Ci, str):
# a hex string of form '#ff0000' for example (for red)
RGB = CC.to_rgb(Ci)
else:
# assume it's an RGB triple already:
RGB = Ci
R.append(RGB[0])
G.append(RGB[1])
B.append(RGB[2])
cmap_dict = {}
cmap_dict['red'] = [(x0[i], R[i], R[i]) for i in xrange(len(R))]
cmap_dict['green'] = [(x0[i], G[i], G[i]) for i in xrange(len(G))]
cmap_dict['blue'] = [(x0[i], B[i], B[i]) for i in xrange(len(B))]
mymap = LinearSegmentedColormap('mymap', cmap_dict)
return mymap
0
Example 90
Project: scipy Source File: mstats_extras.py
def hdquantiles(data, prob=list([.25,.5,.75]), axis=None, var=False,):
"""
Computes quantile estimates with the Harrell-Davis method.
The quantile estimates are calculated as a weighted linear combination
of order statistics.
Parameters
----------
data : array_like
Data array.
prob : sequence, optional
Sequence of quantiles to compute.
axis : int or None, optional
Axis along which to compute the quantiles. If None, use a flattened
array.
var : bool, optional
Whether to return the variance of the estimate.
Returns
-------
hdquantiles : MaskedArray
A (p,) array of quantiles (if `var` is False), or a (2,p) array of
quantiles and variances (if `var` is True), where ``p`` is the
number of quantiles.
"""
def _hd_1D(data,prob,var):
"Computes the HD quantiles for a 1D array. Returns nan for invalid data."
xsorted = np.squeeze(np.sort(data.compressed().view(ndarray)))
# Don't use length here, in case we have a numpy scalar
n = xsorted.size
hd = np.empty((2,len(prob)), float_)
if n < 2:
hd.flat = np.nan
if var:
return hd
return hd[0]
v = np.arange(n+1) / float(n)
betacdf = beta.cdf
for (i,p) in enumerate(prob):
_w = betacdf(v, (n+1)*p, (n+1)*(1-p))
w = _w[1:] - _w[:-1]
hd_mean = np.dot(w, xsorted)
hd[0,i] = hd_mean
#
hd[1,i] = np.dot(w, (xsorted-hd_mean)**2)
#
hd[0, prob == 0] = xsorted[0]
hd[0, prob == 1] = xsorted[-1]
if var:
hd[1, prob == 0] = hd[1, prob == 1] = np.nan
return hd
return hd[0]
# Initialization & checks
data = ma.array(data, copy=False, dtype=float_)
p = np.array(prob, copy=False, ndmin=1)
# Computes quantiles along axis (or globally)
if (axis is None) or (data.ndim == 1):
result = _hd_1D(data, p, var)
else:
if data.ndim > 2:
raise ValueError("Array 'data' must be at most two dimensional, "
"but got data.ndim = %d" % data.ndim)
result = ma.apply_along_axis(_hd_1D, axis, data, p, var)
return ma.fix_invalid(result, copy=False)
0
Example 91
Project: pyAudioAnalysis Source File: audioSegmentation.py
def silenceRemoval(x, Fs, stWin, stStep, smoothWindow=0.5, Weight=0.5, plot=False):
'''
Event Detection (silence removal)
ARGUMENTS:
- x: the input audio signal
- Fs: sampling freq
- stWin, stStep: window size and step in seconds
- smoothWindow: (optinal) smooth window (in seconds)
- Weight: (optinal) weight factor (0 < Weight < 1) the higher, the more strict
- plot: (optinal) True if results are to be plotted
RETURNS:
- segmentLimits: list of segment limits in seconds (e.g [[0.1, 0.9], [1.4, 3.0]] means that
the resulting segments are (0.1 - 0.9) seconds and (1.4, 3.0) seconds
'''
if Weight >= 1:
Weight = 0.99
if Weight <= 0:
Weight = 0.01
# Step 1: feature extraction
x = audioBasicIO.stereo2mono(x) # convert to mono
ShortTermFeatures = aF.stFeatureExtraction(x, Fs, stWin * Fs, stStep * Fs) # extract short-term features
# Step 2: train binary SVM classifier of low vs high energy frames
EnergySt = ShortTermFeatures[1, :] # keep only the energy short-term sequence (2nd feature)
E = numpy.sort(EnergySt) # sort the energy feature values:
L1 = int(len(E) / 10) # number of 10% of the total short-term windows
T1 = numpy.mean(E[0:L1]) + 0.000000000000001 # compute "lower" 10% energy threshold
T2 = numpy.mean(E[-L1:-1]) + 0.000000000000001 # compute "higher" 10% energy threshold
Class1 = ShortTermFeatures[:, numpy.where(EnergySt <= T1)[0]] # get all features that correspond to low energy
Class2 = ShortTermFeatures[:, numpy.where(EnergySt >= T2)[0]] # get all features that correspond to high energy
featuresSS = [Class1.T, Class2.T] # form the binary classification task and ...
[featuresNormSS, MEANSS, STDSS] = aT.normalizeFeatures(featuresSS) # normalize and ...
SVM = aT.trainSVM(featuresNormSS, 1.0) # train the respective SVM probabilistic model (ONSET vs SILENCE)
# Step 3: compute onset probability based on the trained SVM
ProbOnset = []
for i in range(ShortTermFeatures.shape[1]): # for each frame
curFV = (ShortTermFeatures[:, i] - MEANSS) / STDSS # normalize feature vector
ProbOnset.append(SVM.predict_proba(curFV.reshape(1,-1))[0][1]) # get SVM probability (that it belongs to the ONSET class)
ProbOnset = numpy.array(ProbOnset)
ProbOnset = smoothMovingAvg(ProbOnset, smoothWindow / stStep) # smooth probability
# Step 4A: detect onset frame indices:
ProbOnsetSorted = numpy.sort(ProbOnset) # find probability Threshold as a weighted average of top 10% and lower 10% of the values
Nt = ProbOnsetSorted.shape[0] / 10
T = (numpy.mean((1 - Weight) * ProbOnsetSorted[0:Nt]) + Weight * numpy.mean(ProbOnsetSorted[-Nt::]))
MaxIdx = numpy.where(ProbOnset > T)[0] # get the indices of the frames that satisfy the thresholding
i = 0
timeClusters = []
segmentLimits = []
# Step 4B: group frame indices to onset segments
while i < len(MaxIdx): # for each of the detected onset indices
curCluster = [MaxIdx[i]]
if i == len(MaxIdx)-1:
break
while MaxIdx[i+1] - curCluster[-1] <= 2:
curCluster.append(MaxIdx[i+1])
i += 1
if i == len(MaxIdx)-1:
break
i += 1
timeClusters.append(curCluster)
segmentLimits.append([curCluster[0] * stStep, curCluster[-1] * stStep])
# Step 5: Post process: remove very small segments:
minDuration = 0.2
segmentLimits2 = []
for s in segmentLimits:
if s[1] - s[0] > minDuration:
segmentLimits2.append(s)
segmentLimits = segmentLimits2
if plot:
timeX = numpy.arange(0, x.shape[0] / float(Fs), 1.0 / Fs)
plt.subplot(2, 1, 1)
plt.plot(timeX, x)
for s in segmentLimits:
plt.axvline(x=s[0])
plt.axvline(x=s[1])
plt.subplot(2, 1, 2)
plt.plot(numpy.arange(0, ProbOnset.shape[0] * stStep, stStep), ProbOnset)
plt.title('Signal')
for s in segmentLimits:
plt.axvline(x=s[0])
plt.axvline(x=s[1])
plt.title('SVM Probability')
plt.show()
return segmentLimits
0
Example 92
Project: mystic Source File: symbolic.py
def get_variables(constraints, variables='x'):
"""extract a list of the string variable names from constraints string
Inputs:
constraints -- a string of symbolic constraints, with one constraint
equation per line. Constraints can be equality and/or inequality
constraints. Standard python syntax should be followed (with the
math and numpy modules already imported).
For example:
>>> constraints = '''
... x1 + x2 = x3*4
... x3 = x2*x4'''
>>> get_variables(constraints)
['x1', 'x2', 'x3', 'x4']
Additional Inputs:
variables -- desired variable name. Default is 'x'. A list of variable
name strings is also accepted for when desired variable names
don't have the same base, and can include variables that are not
found in the constraints equation string.
For example:
>>> constraints = '''
... y = min(u,v) - z*sin(x)
... z = x**2 + 1.0
... u = v*z'''
>>> get_variables(constraints, list('pqrstuvwxyz'))
['u', 'v', 'x', 'y', 'z']
"""
if list_or_tuple_or_ndarray(variables):
equations = replace_variables(constraints,variables,'_')
vars = get_variables(equations,'_')
indices = [int(v.strip('_')) for v in vars]
varnamelist = []
from numpy import sort
for i in sort(indices):
varnamelist.append(variables[i])
return varnamelist
import re
target = variables+'[0-9]+'
varnamelist = []
equation_list = constraints.splitlines()
for equation in equation_list:
vars = re.findall(target,equation)
for var in vars:
if var not in varnamelist:
varnamelist.append(var)
return varnamelist
0
Example 93
Project: scipy Source File: mstats_extras.py
def hdquantiles_sd(data, prob=list([.25,.5,.75]), axis=None):
"""
The standard error of the Harrell-Davis quantile estimates by jackknife.
Parameters
----------
data : array_like
Data array.
prob : sequence, optional
Sequence of quantiles to compute.
axis : int, optional
Axis along which to compute the quantiles. If None, use a flattened
array.
Returns
-------
hdquantiles_sd : MaskedArray
Standard error of the Harrell-Davis quantile estimates.
"""
def _hdsd_1D(data,prob):
"Computes the std error for 1D arrays."
xsorted = np.sort(data.compressed())
n = len(xsorted)
#.........
hdsd = np.empty(len(prob), float_)
if n < 2:
hdsd.flat = np.nan
vv = np.arange(n) / float(n-1)
betacdf = beta.cdf
for (i,p) in enumerate(prob):
_w = betacdf(vv, (n+1)*p, (n+1)*(1-p))
w = _w[1:] - _w[:-1]
mx_ = np.fromiter([np.dot(w,xsorted[np.r_[list(range(0,k)),
list(range(k+1,n))].astype(int_)])
for k in range(n)], dtype=float_)
mx_var = np.array(mx_.var(), copy=False, ndmin=1) * n / float(n-1)
hdsd[i] = float(n-1) * np.sqrt(np.diag(mx_var).diagonal() / float(n))
return hdsd
# Initialization & checks
data = ma.array(data, copy=False, dtype=float_)
p = np.array(prob, copy=False, ndmin=1)
# Computes quantiles along axis (or globally)
if (axis is None):
result = _hdsd_1D(data, p)
else:
if data.ndim > 2:
raise ValueError("Array 'data' must be at most two dimensional, "
"but got data.ndim = %d" % data.ndim)
result = ma.apply_along_axis(_hdsd_1D, axis, data, p)
return ma.fix_invalid(result, copy=False).ravel()
0
Example 94
def update_w(self):
# compute distance matrix -> requiresd for the volume
self.init_sivm()
next_sel = list([self.select[0]])
self.select = []
self._v = []
self._t = []
stime = time.time()
for iter in range(self._num_bases-1):
# add new selections to openset
next_sel = list(np.sort(next_sel))
D = pdist(self.data[:, next_sel], self.data[:, next_sel])
V = np.zeros(self.data.shape[1])
d = np.zeros((D.shape[0]+1,D.shape[1]+1))
d[:D.shape[0], :D.shape[1]] = D[:,:]
for i in range(self.data.shape[1]):
# create a temp selection
dtmp = l2_distance(self.data[:,next_sel], self.data[:,i:i+1])
d[:-1,-1] = dtmp
d[-1,:-1] = dtmp
# compute volume for temp selection
V[i] = cmdet(d)
next_index = np.argmax(V)
next_sel.append(next_index)
self._v.append(np.max(V))
self._logger.info('Iter:' + str(iter))
self._logger.info('Current selection:' + str(next_sel))
self._logger.info('Current volume:' + str(self._v[-1]))
self._t.append(time.time() - stime)
# update some values ...
self.select = list(next_sel)
self.W = self.data[:, self.select]
0
Example 95
Project: bhmm Source File: analysis.py
def empirical_confidence_interval(sample, interval=0.95):
"""
Compute specified symmetric confidence interval for empirical sample.
Parameters
----------
sample : numpy.array
The empirical samples.
interval : float, optional, default=0.95
Size of desired symmetric confidence interval (0 < interval < 1)
e.g. 0.68 for 68% confidence interval, 0.95 for 95% confidence interval
Returns
-------
low : float
The lower bound of the symmetric confidence interval.
high : float
The upper bound of the symmetric confidence interval.
Examples
--------
>>> sample = np.random.randn(1000)
>>> [low, high] = empirical_confidence_interval(sample)
>>> [low, high] = empirical_confidence_interval(sample, interval=0.65)
>>> [low, high] = empirical_confidence_interval(sample, interval=0.99)
"""
# Sort sample in increasing order.
sample = np.sort(sample)
# Determine sample size.
N = len(sample)
# Compute low and high indices.
low_index = int(np.round((N-1) * (0.5 - interval/2))) + 1
high_index = int(np.round((N-1) * (0.5 + interval/2))) + 1
# Compute low and high.
low = sample[low_index]
high = sample[high_index]
return [low, high]
0
Example 96
Project: BioSPPy Source File: ecg.py
def _extract_heartbeats(signal=None, rpeaks=None, before=200, after=400):
"""Extract heartbeat templates from an ECG signal, given a list of
R-peak locations.
Parameters
----------
signal : array
Input ECG signal.
rpeaks : array
R-peak location indices.
before : int, optional
Number of samples to include before the R peak.
after : int, optional
Number of samples to include after the R peak.
Returns
-------
templates : array
Extracted heartbeat templates.
rpeaks : array
Corresponding R-peak location indices of the extracted heartbeat
templates.
"""
R = np.sort(rpeaks)
length = len(signal)
templates = []
newR = []
for r in R:
a = r - before
if a < 0:
continue
b = r + after
if b > length:
break
templates.append(signal[a:b])
newR.append(r)
templates = np.array(templates)
newR = np.array(newR, dtype='int')
return templates, newR
0
Example 97
Project: python-control Source File: yottalab.py
def placep(A,B,P):
"""Return the steady state value of the step response os sysmatrix K for
pole placement
Usage
=====
K = placep(A,B,P)
Inputs
------
A : State matrix A
B : INput matrix
P : desired poles
Outputs
-------
K : State gains for pole placement
"""
n = shape(A)[0]
m = shape(B)[1]
tol = 0.0
mode = 1;
wrka = zeros((n,m))
wrk1 = zeros(m)
wrk2 = zeros(m)
iwrk = zeros((m),np.int)
A,B,ncont,indcont,nblk,z = _wrapper.ssxmc(n,m,A,n,B,wrka,wrk1,wrk2,iwrk,tol,mode)
P = sort(P)
wr = real(P)
wi = imag(P)
g = zeros((m,n))
mx = max(2,m)
rm1 = zeros((m,m))
rm2 = zeros((m,mx))
rv1 = zeros(n)
rv2 = zeros(n)
rv3 = zeros(m)
rv4 = zeros(m)
A,B,g,z,ierr,jpvt = _wrapper.polmc(A,B,g,wr,wi,z,indcont,nblk,rm1, rm2, rv1, rv2, rv3, rv4)
return g
0
Example 98
def debug(time, cues, choice, reward):
n = len(cues)
cues = np.sort(cues)
R.append(reward)
if choice == cues[0]:
P.append(1)
else:
P.append(0)
print "Choice: ",
for i in range(n):
if choice == cues[i]:
print "[%d]" % cues[i],
else:
print "%d" % cues[i],
if i < (n-1):
print "/",
if choice == cues[0]:
print " (good)"
else:
print " (bad)"
print "Reward (%3d%%) : %d" % (int(100*CUE["reward"][choice]),reward)
print "Mean performance: %.3f" % np.array(P).mean()
print "Mean reward: %.3f" % np.array(R).mean()
print "Response time: %d ms" % (time)
print "CTX.cog->CTX.ass:", connections["CTX.cog -> CTX.ass"].weights
print
0
Example 99
Project: scipy Source File: mstats_extras.py
def mjci(data, prob=[0.25,0.5,0.75], axis=None):
"""
Returns the Maritz-Jarrett estimators of the standard error of selected
experimental quantiles of the data.
Parameters
----------
data : ndarray
Data array.
prob : sequence, optional
Sequence of quantiles to compute.
axis : int or None, optional
Axis along which to compute the quantiles. If None, use a flattened
array.
"""
def _mjci_1D(data, p):
data = np.sort(data.compressed())
n = data.size
prob = (np.array(p) * n + 0.5).astype(int_)
betacdf = beta.cdf
mj = np.empty(len(prob), float_)
x = np.arange(1,n+1, dtype=float_) / n
y = x - 1./n
for (i,m) in enumerate(prob):
W = betacdf(x,m-1,n-m) - betacdf(y,m-1,n-m)
C1 = np.dot(W,data)
C2 = np.dot(W,data**2)
mj[i] = np.sqrt(C2 - C1**2)
return mj
data = ma.array(data, copy=False)
if data.ndim > 2:
raise ValueError("Array 'data' must be at most two dimensional, "
"but got data.ndim = %d" % data.ndim)
p = np.array(prob, copy=False, ndmin=1)
# Computes quantiles along axis (or globally)
if (axis is None):
return _mjci_1D(data, p)
else:
return ma.apply_along_axis(_mjci_1D, axis, data, p)
0
Example 100
Project: bolt Source File: array.py
def transpose(self, *axes):
"""
Return an array with the axes transposed.
This operation will incur a swap unless the
desiured permutation can be obtained
only by transpoing the keys or the values.
Parameters
----------
axes : None, tuple of ints, or n ints
If None, will reverse axis order.
"""
if len(axes) == 0:
p = arange(self.ndim-1, -1, -1)
else:
p = asarray(argpack(axes))
istransposeable(p, range(self.ndim))
split = self.split
# compute the keys/value axes that need to be swapped
new_keys, new_values = p[:split], p[split:]
swapping_keys = sort(new_values[new_values < split])
swapping_values = sort(new_keys[new_keys >= split])
stationary_keys = sort(new_keys[new_keys < split])
stationary_values = sort(new_values[new_values >= split])
# compute the permutation that the swap causes
p_swap = r_[stationary_keys, swapping_values, swapping_keys, stationary_values]
# compute the extra permutation (p_x) on top of this that
# needs to happen to get the full permutation desired
p_swap_inv = argsort(p_swap)
p_x = p_swap_inv[p]
p_keys, p_values = p_x[:split], p_x[split:]-split
# perform the swap and the the within key/value permutations
arr = self.swap(swapping_keys, swapping_values-split)
arr = arr.keys.transpose(tuple(p_keys.tolist()))
arr = arr.values.transpose(tuple(p_values.tolist()))
return arr