numpy.sort

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 7

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)

Example 52

Project: neurosynth Source File: stats.py
Function: fdr
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

Example 53

Project: nupic.research Source File: classify_keywords.py
Function: encode_token
  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))

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

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

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))

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)

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

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)

Example 60

Project: pybasicbayes Source File: general.py
Function: scoreatpercentile
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)]]

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)

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)

Example 63

Project: pyKriging Source File: samplingplan.py
Function: mm
    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]

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)

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

Example 66

Project: geostatsmodels Source File: zscoretrans.py
Function: cdf
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

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

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)

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

Example 70

Project: statsmodels Source File: gof_new.py
Function: ks_test
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

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()

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

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

Example 74

Project: pybrain Source File: fem.py
Function: learnstep
    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]

Example 75

Project: ray Source File: microbenchmarks.py
Function: test_timing
  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()

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)

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

Example 78

Project: acoular Source File: fileimport.py
Function: get_data
    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()

Example 79

Project: scipy Source File: mstats_extras.py
Function: median_cihs
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

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()

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)

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]

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)

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)

Example 85

Project: spotpy Source File: analyser.py
def sort_like(results):
    return np.sort(results,axis=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))

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

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

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

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)

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

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

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()

Example 94

Project: msaf Source File: sivm_sgreedy.py
Function: update_w
    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]

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]

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

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

Example 98

Project: Neurosciences Source File: 250-trials.py
Function: debug
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

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)

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
See More Examples - Go to Next Page
Page 1 Page 2 Selected Page 3 Page 4