Set the centre of the Gaussian basis 
        functions be spaced evenly throughout run time
187321
def gen_centers(self):
        """"""

        c = np.linspace(0, 2*np.pi, self.bfs+1)
        c = c[0:-1]
        self.c = c
1125515
def fit(self, X):
        
        self.centers_, self.labels_, self.sse_arr_, self.n_iter_ = \
              _kmedoids(X, self.n_clusters, self.distance, self.max_iter, self.n_trials, self.tol, self.rng)
280255
def _update_centers(X, membs, n_clusters):
    
    centers = np.empty(shape=(n_clusters, X.shape[1]), dtype=float)
    sse = np.empty(shape=n_clusters, dtype=float)
    for clust_id in range(n_clusters):
        memb_ids = np.where(membs == clust_id)[0]

        if memb_ids.shape[0] == 0:
            memb_ids = np.random.choice(X.shape[0], size=1)
            #print("Empty cluster replaced with ", memb_ids)
        centers[clust_id,:] = np.mean(X[memb_ids,:], axis=0)
        
        sse[clust_id] = _cal_dist2center(X[memb_ids,:], centers[clust_id,:]) 
    return(centers, sse)
1049172
def overlap(alpha1,lmn1,A,alpha2,lmn2,B):
    
    l1,m1,n1 = lmn1
    l2,m2,n2 = lmn2
    rab2 = norm2(A-B)
    gamma = alpha1+alpha2
    P = gaussian_product_center(alpha1,A,alpha2,B)

    pre = pow(pi/gamma,1.5)*exp(-alpha1*alpha2*rab2/gamma)

    wx = overlap1d(l1,l2,P[0]-A[0],P[0]-B[0],gamma)
    wy = overlap1d(m1,m2,P[1]-A[1],P[1]-B[1],gamma)
    wz = overlap1d(n1,n2,P[2]-A[2],P[2]-B[2],gamma)
    return pre*wx*wy*wz
1332731
def centering_centroid(data, xi, yi, box, nloop=10, toldist=1e-3,
                       maxdist=10.0):
    '''
        
    '''

    # Store original center
    cxy = (xi, yi)
    origin = (xi, yi)
    # initial background
    back = 0.0

    if nloop == 0:
        return xi, yi, 0.0, 0, 'not recentering'

    for i in range(nloop):
        nxy, back = _centering_centroid_loop_xy(data, cxy, box)
        # _logger.debug('new center is %s', ncenter)
        # if we are to far away from the initial point, break
        dst = distance.euclidean(origin, nxy)
        if dst > maxdist:
            msg = 'maximum distance (%5.2f) from origin reached' % maxdist
            return cxy[0], cxy[1], back, 2, msg

        # check convergence
        dst = distance.euclidean(nxy, cxy)
        if dst < toldist:
            return nxy[0], nxy[1], back, 1, 'converged in iteration %i' % i
        else:
            cxy = nxy

    return nxy[0], nxy[1], back, 3, 'not converged in %i iterations' % nloop
992713
def gaussian(x,  a, b, c, d=0):
    '''
    
    '''
    return a * np.exp(  -(((x-b)**2 )/ (2*(c**2)))  ) + d
417841
private static void balbak(DenseMatrix V, double[] scale) {
        int n = V.nrows();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                V.mul(i, j, scale[i]);
            }
        }
    }
635257
def init_prior(self, R):
        
        centers, widths = self.init_centers_widths(R)
        # update prior
        prior = np.zeros(self.K * (self.n_dim + 1))
        self.set_centers(prior, centers)
        self.set_widths(prior, widths)
        self.set_prior(prior)
        return self
695018
def linear_positions(n_channels):
    """"""
    return np.c_[np.zeros(n_channels),
                 np.linspace(0., 1., n_channels)]
1049851
def set(self, opts, popsize=None, ccovfac=1, verbose=True):
        

        alpha_cc = 1.0  # cc-correction for mueff, was zero before

        def conedf(df, mu, N):
            """used for computing separable learning rate"""
            return 1. / (df + 2.*sqrt(df) + float(mu) / N)

        def cmudf(df, mu, alphamu):
            """used for computing separable learning rate"""
            return (alphamu + mu - 2. + 1. / mu) / (df + 4.*sqrt(df) + mu / 2.)

        sp = self
        N = sp.N
        if popsize:
            opts.evalall({'N':N, 'popsize':popsize})
        else:
            popsize = opts.evalall({'N':N})['popsize']  # the default popsize is computed in CMAOptions()
        ## meta_parameters.lambda_exponent == 0.0
        popsize = int(popsize + N** 0.0 - 1)
        sp.popsize = popsize
        if opts['CMA_mirrors'] < 0.5:
            sp.lam_mirr = int(0.5 + opts['CMA_mirrors'] * popsize)
        elif opts['CMA_mirrors'] > 1:
            sp.lam_mirr = int(0.5 + opts['CMA_mirrors'])
        else:
            sp.lam_mirr = int(0.5 + 0.16 * min((popsize, 2 * N + 2)) + 0.29)  # 0.158650... * popsize is optimal
            # lam = arange(2,22)
            # mirr = 0.16 + 0.29/lam
            # print(lam); print([int(0.5 + l) for l in mirr*lam])
            # [ 2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21]
            # [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4]
        ## meta_parameters.parent_fraction == 0.5
        sp.mu_f = 0.5 * sp.popsize  # float value of mu
        if opts['CMA_mu'] is not None:
            sp.mu_f = opts['CMA_mu']
        sp.mu = int(sp.mu_f + 0.499999)  # round down for x.5
        sp.mu = max((sp.mu, 1))
        # in principle we have mu_opt = popsize/2 + lam_mirr/2,
        # which means in particular weights should only be negative for q > 0.5+mirr_frac/2
        if sp.mu > sp.popsize - 2 * sp.lam_mirr + 1:
            _print_warning("pairwise selection is not implemented, therefore " +
                  " mu = %d > %d = %d - 2*%d + 1 = popsize - 2*mirr + 1 can produce a bias" % (
                    sp.mu, sp.popsize - 2 * sp.lam_mirr + 1, sp.popsize, sp.lam_mirr))
        if sp.lam_mirr > sp.popsize // 2:
            raise _Error("fraction of mirrors in the population as read from option CMA_mirrors cannot be larger 0.5, " +
                         "theoretically optimal is 0.159")
        sp.weights = log(max([sp.mu, sp.popsize / 2.0]) + 0.5) - log(1 + np.arange(sp.mu))
        sp.weights /= sum(sp.weights)
        sp.mueff = 1 / sum(sp.weights**2)
        # TODO: this will disappear, as it is done in class CMAAdaptSigmaCSA
        ## meta_parameters.cs_exponent == 1.0
        b = 1.0
        ## meta_parameters.cs_multiplier == 1.0
        sp.cs = 1.0 * (sp.mueff + 2)**b / (N + (sp.mueff + 3)**b)  # TODO: this doesn't change dependency of dimension
        # sp.cs = (sp.mueff + 2) / (N + 1.5*sp.mueff + 1)
        ## meta_parameters.cc_exponent == 1.0
        b = 1.0
        ## meta_parameters.cc_multiplier == 1.0
        sp.cc = 1.0 * \
                (4 + alpha_cc * sp.mueff / N)**b / \
                (N**b + (4 + alpha_cc * 2 * sp.mueff / N)**b)
        sp.cc_sep = (1 + 1 / N + alpha_cc * sp.mueff / N) / (N**0.5 + 1 / N + alpha_cc * 2 * sp.mueff / N)  # \not\gg\cc
        if hasattr(opts['vv'], '__getitem__') and opts['vv'][0] == 'sweep_ccov1':
            ## meta_parameters.cc_multiplier == 1.0
            sp.cc = 1.0 * (4 + sp.mueff / N)**0.5 / ((N + 4)**0.5 + (2 * sp.mueff / N)**0.5)
        sp.rankmualpha = opts['CMA_rankmualpha']
        # sp.rankmualpha = _evalOption(opts['CMA_rankmualpha'], 0.3)
        ## meta_parameters.c1_multiplier == 1.0
        sp.c1 = ( 1.0 * ccovfac * min(1, sp.popsize / 6) *
                 ## meta_parameters.c1_exponent == 2.0
                 2 / ((N + 1.3)** 2.0 + sp.mueff))
        # 1/0
        sp.c1_sep = ccovfac * conedf(N, sp.mueff, N)
        if opts['CMA_rankmu'] != 0:  # also empty
            ## meta_parameters.cmu_multiplier == 2.0
            alphacov, mu = 2.0 , sp.mueff
            sp.cmu = min(1 - sp.c1, ccovfac * alphacov *
                         ## meta_parameters.cmu_exponent == 2.0
                         (sp.rankmualpha + mu - 2 + 1 / mu) / ((N + 2)** 2.0 + alphacov * mu / 2))
            if hasattr(opts['vv'], '__getitem__') and opts['vv'][0] == 'sweep_ccov':
                sp.cmu = opts['vv'][1]
            sp.cmu_sep = min(1 - sp.c1_sep, ccovfac * cmudf(N, sp.mueff, sp.rankmualpha))
        else:
            sp.cmu = sp.cmu_sep = 0
        if hasattr(opts['vv'], '__getitem__') and opts['vv'][0] == 'sweep_ccov1':
            sp.c1 = opts['vv'][1]

        sp.neg = _BlancClass()
        if opts['CMA_active'] and opts['CMA_on']:
            # in principle we have mu_opt = popsize/2 + lam_mirr/2,
            # which means in particular weights should only be negative for q > 0.5+mirr_frac/2
            if 1 < 3: # seems most natural: continuation of log(lambda/2) - log(k) qqqqqqqqqqqqqqqqqqqqqqqqqq
                sp.neg.mu_f = popsize // 2  # not sure anymore what this is good for
                sp.neg.weights = array([log(k) - log(popsize/2 + 1/2) for k in np.arange(np.ceil(popsize/2 + 1.1/2), popsize + .1)])
            sp.neg.mu = len(sp.neg.weights)
            sp.neg.weights /= sum(sp.neg.weights)
            sp.neg.mueff = 1 / sum(sp.neg.weights**2)
            ## meta_parameters.cact_exponent == 1.5
            sp.neg.cmuexp = opts['CMA_active'] * 0.3 * sp.neg.mueff / ((N + 2)** 1.5 + 1.0 * sp.neg.mueff)
            if hasattr(opts['vv'], '__getitem__') and opts['vv'][0] == 'sweep_ccov_neg':
                sp.neg.cmuexp = opts['vv'][1]
            # reasoning on learning rate cmuexp: with sum |w| == 1 and
            #   length-normalized vectors in the update, the residual
            #   variance in any direction exceeds exp(-N*cmuexp)
            assert sp.neg.mu >= sp.lam_mirr  # not really necessary
            # sp.neg.minresidualvariance = 0.66  # not it use, keep at least 0.66 in all directions, small popsize is most critical
        else:
            sp.neg.cmuexp = 0

        sp.CMA_on = sp.c1 + sp.cmu > 0
        # print(sp.c1_sep / sp.cc_sep)

        if not opts['CMA_on'] and opts['CMA_on'] not in (None, [], (), ''):
            sp.CMA_on = False
            # sp.c1 = sp.cmu = sp.c1_sep = sp.cmu_sep = 0
        mueff_exponent = 0.5
        if 1 < 3:
            mueff_exponent = opts['CSA_damp_mueff_exponent']
        # TODO: this will disappear, as it is done in class CMAAdaptSigmaCSA
        sp.damps = opts['CSA_dampfac'] * (0.5 +
                                          0.5 * min([1, (sp.lam_mirr / (0.159 * sp.popsize) - 1)**2])**1 +
                                          2 * max([0, ((sp.mueff - 1) / (N + 1))**mueff_exponent - 1]) + sp.cs
                                          )
        sp.cmean = float(opts['CMA_cmean'])
        # sp.kappa = 1  # 4-D, lam=16, rank1, kappa < 4 does not influence convergence rate
                        # in larger dim it does, 15-D with defaults, kappa=8 factor 2
        if verbose:
            if not sp.CMA_on:
                print('covariance matrix adaptation turned off')
            if opts['CMA_mu'] != None:
                print('mu = %f' % (sp.mu_f))
39917
def _basis_spline_factory(coef, degree, knots, der, ext):
        """"""
        return functools.partial(interpolate.splev, tck=(knots, coef, degree), der=der, ext=ext)
1345552
def kernel(x1, x2, method='gaussian', sigma=1, **kwargs):
    """"""
    if method.lower() in ['gaussian', 'gauss', 'rbf']:
        K = np.exp(-dist(x1, x2) / (2 * sigma**2))
    return K
715485
def stSpectralCentroidAndSpread(X, fs):
    """"""
    ind = (numpy.arange(1, len(X) + 1)) * (fs/(2.0 * len(X)))

    Xt = X.copy()
    Xt = Xt / Xt.max()
    NUM = numpy.sum(ind * Xt)
    DEN = numpy.sum(Xt) + eps

    # Centroid:
    C = (NUM / DEN)

    # Spread:
    S = numpy.sqrt(numpy.sum(((ind - C) ** 2) * Xt) / DEN)

    # Normalize:
    C = C / (fs / 2.0)
    S = S / (fs / 2.0)

    return (C, S)
111797
def _init_random_gaussians(self, X):
        """"""

        n_samples = np.shape(X)[0]
        self.priors = (1 / self.k) * np.ones(self.k)
        for _ in range(self.k):
            params = {}
            params["mean"] = X[np.random.choice(range(n_samples))]
            params["cov"] = calculate_covariance_matrix(X)
            self.parameters.append(params)
1069790
def _fix_pooling(pool_type, inputs, new_attr):
    """"""
    stride = new_attr.get('stride')
    kernel = new_attr.get('kernel')
    padding = new_attr.get('pad')
    p_value = new_attr.get('p_value')

    # Adding default stride.
    if stride is None:
        stride = (1,) * len(kernel)

    # Add padding attr if not provided.
    if padding is None:
        padding = (0,) * len(kernel) * 2

    # Mxnet Pad operator supports only 4D/5D tensors.
    # For 1D case, these are the steps:
    #    Step 1. Add extra dummy dimension to make it 4D. Adding to  axis = 2
    #    Step 2. Apply padding to this changed tensor
    #    Step 3. Remove the extra dimension added in step 1.
    if len(kernel) == 1:
        dummy_axis = 2
        # setting 0 padding to the new dim to be added.
        padding = (0, padding[0], 0, padding[1])
        pad_width = (0, 0, 0, 0) + _pad_sequence_fix(padding, kernel_dim=2)

        # Step 1.
        curr_sym = symbol.expand_dims(inputs[0], axis=dummy_axis)

        # Step 2. Common for all tensor sizes
        new_pad_op = symbol.pad(curr_sym, mode='edge', pad_width=pad_width)

        # Step 3: Removing extra dim added.
        new_pad_op = symbol.split(new_pad_op, axis=dummy_axis, num_outputs=1, squeeze_axis=1)
    else:
        # For 2D/3D cases:
        # Apply padding
        pad_width = (0, 0, 0, 0) + _pad_sequence_fix(padding, kernel_dim=len(kernel))
        curr_sym = inputs[0]

        if pool_type == 'max':
            # For max pool : mode = 'edge', we should replicate the
            # edge values to pad, so that we only include  input data values
            # for calculating 'max'
            new_pad_op = symbol.pad(curr_sym, mode='edge', pad_width=pad_width)
        else:
            # For avg pool, we should add 'zeros' for padding  so mode='constant'
            new_pad_op = symbol.pad(curr_sym, mode='constant', pad_width=pad_width)

    # Apply pooling without pads.
    if pool_type == 'lp':
        new_pooling_op = symbol.Pooling(new_pad_op, pool_type=pool_type, stride=stride, kernel=kernel, p_value=p_value)
    else:
        new_pooling_op = symbol.Pooling(new_pad_op, pool_type=pool_type, stride=stride, kernel=kernel)
    return new_pooling_op
1080716
def gaussianize(Y):
    
    N,P = Y.shape

    YY=toRanks(Y)
    quantiles=(sp.arange(N)+0.5)/N
    gauss = st.norm.isf(quantiles)
    Y_gauss=sp.zeros((N,P))
    for i in range(P):
        Y_gauss[:,i] = gauss[YY[:,i]]
    Y_gauss *= -1
    return Y_gauss
455576
function gaussianKernel1D(s, n) {
    if (typeof n !== 'number') n = 7;
    var kernel = new Float64Array(n),
        halfN = Math.floor(n * 0.5),
        odd = n % 2,
        i;
    if (!s || !n) return kernel;
    for (i = 0; i <= halfN; i++) {
        kernel[i] = gauss(s * (i - halfN - odd * 0.5), s);
    }
    for (; i < n; i++) {
        kernel[i] = kernel[n - 1 - i];
    }
    return kernel;
}