Kowalik & Osborne function (lmder1 test #9)
480064
def _lmder1_kowalik_osborne():
    """"""
    v = np.asfarray([4, 2, 1, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625])
    y2 = np.asfarray([0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456,
                      0.0342, 0.0323, 0.0235, 0.0246])

    def func(params, vec):
        tmp1 = v * (v + params[1])
        tmp2 = v * (v + params[2]) + params[3]
        vec[:] = y2 - params[0] * tmp1 / tmp2

    def jac(params, jac):
        tmp1 = v * (v + params[1])
        tmp2 = v * (v + params[2]) + params[3]
        jac[0] = -tmp1 / tmp2
        jac[1] = -v * params[0] / tmp2
        jac[2] = jac[0] * jac[1]
        jac[3] = jac[2] / v

    guess = np.asfarray([0.25, 0.39, 0.415, 0.39])

    _lmder1_driver(11, func, jac, guess,
                   0.7289151028829448e-01, 0.1753583772112895e-01,
                   [0.1928078104762493e+00, 0.1912626533540709e+00,
                    0.1230528010469309e+00, 0.1360532211505167e+00])
    _lmder1_driver(11, func, jac, guess * 10,
                   0.2979370075552020e+01, 0.3205219291793696e-01,
                   [0.7286754737686598e+06, -0.1407588031293926e+02,
                    -0.3297779778419661e+08, -0.2057159419780170e+08])
412702
def demo_err():
    
    max_order = 7
    n = 20
    l = 0.25
    fmt1 = '{0: <5s}\t{1: <21s}\t{2: >21s}\t{3: >21s}\t{4: >21s}'
    fmt2 = '{0: <5d}\t{1:20.18f}\t{2: >21.18f}\t{3: >21.18f}\t{4: >21.18f}'
    x = np.cumsum(np.random.rand(n)*l)
    x = np.concatenate((x[::-1]*-1, x))
    lst = []
    derivs = np.zeros(n)
    for order in range(max_order+1):
        print('Order', order)
        for m in range(1+order//2, n+1):
            sub_x = x[n-m:n+m]
            derivs[m-1] = derivatives_at_point_by_finite_diff(
                sub_x, np.exp(sub_x), 0, order)[order]
        print(fmt1.format('m', 'val', 'diff', 'analytical error',
                          'diff/analytical'))
        for m in range(1, n):
            print(fmt2.format(
                (m+1)*2, derivs[m], derivs[m]-derivs[m-1],
                derivs[m]-1, (derivs[m]-derivs[m-1])/(derivs[m]-1)))
        lst.append((derivs[-1], abs(derivs[-1]-derivs[-2])))
    print(np.array(lst))
344534
def mul_at(self, other, k):
        ''''''
        if k > (self.degree + other.degree) or k > self.degree: return 0 # optimization: if the required coefficient is above the maximum coefficient of the resulting polynomial, we can already predict that and just return 0

        term = 0

        for i in _range(min(len(self), len(other))):
            coef1 = self.coefficients[-(k-i+1)]
            coef2 = other.coefficients[-(i+1)]
            if coef1 == 0 or coef2 == 0: continue # log(0) is undefined, skip (and in addition it's a nice optimization)
            term += coef1 * coef2
        return term
284408
function jenksMatrices(data, n_classes) {

        // in the original implementation, these matrices are referred to
        // as `LC` and `OP`
        //
        // * lower_class_limits (LC): optimal lower class limits
        // * variance_combinations (OP): optimal variance combinations for all classes
        var lower_class_limits = [],
            variance_combinations = [],
            // loop counters
            i, j,
            // the variance, as computed at each step in the calculation
            variance = 0;

        // Initialize and fill each matrix with zeroes
        for (i = 0; i < data.length + 1; i++) {
            var tmp1 = [], tmp2 = [];
            // despite these arrays having the same values, we need
            // to keep them separate so that changing one does not change
            // the other
            for (j = 0; j < n_classes + 1; j++) {
                tmp1.push(0);
                tmp2.push(0);
            }
            lower_class_limits.push(tmp1);
            variance_combinations.push(tmp2);
        }

        for (i = 1; i < n_classes + 1; i++) {
            lower_class_limits[1][i] = 1;
            variance_combinations[1][i] = 0;
            // in the original implementation, 9999999 is used but
            // since Javascript has `Infinity`, we use that.
            for (j = 2; j < data.length + 1; j++) {
                variance_combinations[j][i] = Infinity;
            }
        }

        for (var l = 2; l < data.length + 1; l++) {

            // `SZ` originally. this is the sum of the values seen thus
            // far when calculating variance.
            var sum = 0,
                // `ZSQ` originally. the sum of squares of values seen
                // thus far
                sum_squares = 0,
                // `WT` originally. This is the number of
                w = 0,
                // `IV` originally
                i4 = 0;

            // in several instances, you could say `Math.pow(x, 2)`
            // instead of `x * x`, but this is slower in some browsers
            // introduces an unnecessary concept.
            for (var m = 1; m < l + 1; m++) {

                // `III` originally
                var lower_class_limit = l - m + 1,
                    val = data[lower_class_limit - 1];

                // here we're estimating variance for each potential classing
                // of the data, for each potential number of classes. `w`
                // is the number of data points considered so far.
                w++;

                // increase the current sum and sum-of-squares
                sum += val;
                sum_squares += val * val;

                // the variance at this point in the sequence is the difference
                // between the sum of squares and the total x 2, over the number
                // of samples.
                variance = sum_squares - (sum * sum) / w;

                i4 = lower_class_limit - 1;

                if (i4 !== 0) {
                    for (j = 2; j < n_classes + 1; j++) {
                        // if adding this element to an existing class
                        // will increase its variance beyond the limit, break
                        // the class at this point, setting the `lower_class_limit`
                        // at this point.
                        if (variance_combinations[l][j] >=
                            (variance + variance_combinations[i4][j - 1])) {
                            lower_class_limits[l][j] = lower_class_limit;
                            variance_combinations[l][j] = variance +
                                variance_combinations[i4][j - 1];
                        }
                    }
                }
            }

            lower_class_limits[l][1] = 1;
            variance_combinations[l][1] = variance;
        }

        // return the two matrices. for just providing breaks, only
        // `lower_class_limits` is needed, but variances can be useful to
        // evaluate goodness of fit.
        return {
            lower_class_limits: lower_class_limits,
            variance_combinations: variance_combinations
        };
    }
67032
def _param_grad_helper(self, dL_dK, X, X2, target):
        """"""
        if X2 is None: X2 = X
        dist = np.abs(X - X2.T)

        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU

        dk1theta1 = np.exp(-ly*dist)*2*(-lu)/(lu+ly)**3
        #c=np.sqrt(3)
        #t1=c/lu
        #t2=1/ly
        #dk1theta1=np.exp(-dist*ly)*t2*( (2*c*t2+2*t1)/(c*t2+t1)**2 -2*(2*c*t2*t1+t1**2)/(c*t2+t1)**3   )
        
        dk2theta1 = 1*( 
            np.exp(-lu*dist)*dist*(-ly+2*lu-lu*ly*dist+dist*lu**2)*(ly-lu)**(-2) + np.exp(-lu*dist)*(-2+ly*dist-2*dist*lu)*(ly-lu)**(-2) 
            +np.exp(-dist*lu)*(ly-2*lu+ly*lu*dist-dist*lu**2)*2*(ly-lu)**(-3) 
            +np.exp(-dist*ly)*2*(ly-lu)**(-2)
            +np.exp(-dist*ly)*2*(2*lu-ly)*(ly-lu)**(-3)
            )
      
        dk3theta1 = np.exp(-dist*lu)*(lu+ly)**(-2)*((2*lu+ly+dist*lu**2+lu*ly*dist)*(-dist-2/(lu+ly))+2+2*lu*dist+ly*dist)

        dktheta1 = self.varianceU*self.varianceY*(dk1theta1+dk2theta1+dk3theta1)




        dk1theta2 = np.exp(-ly*dist) * ((lu+ly)**(-2)) * (  (-dist)*(2*lu+ly)  +  1  +  (-2)*(2*lu+ly)/(lu+ly)  )

        dk2theta2 = 1*(
            np.exp(-dist*lu)*(ly-lu)**(-2) * ( 1+lu*dist+(-2)*(ly-2*lu+lu*ly*dist-dist*lu**2)*(ly-lu)**(-1) )
            +np.exp(-dist*ly)*(ly-lu)**(-2) * ( (-dist)*(2*lu-ly) -1+(2*lu-ly)*(-2)*(ly-lu)**(-1) )
            )

        dk3theta2 = np.exp(-dist*lu) * (-3*lu-ly-dist*lu**2-lu*ly*dist)/(lu+ly)**3

        dktheta2 = self.varianceU*self.varianceY*(dk1theta2 + dk2theta2 +dk3theta2)



        k1 = np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
        k2 = (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
        k3 = np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
        dkdvar = k1+k2+k3

        target[0] += np.sum(self.varianceY*dkdvar * dL_dK)
        target[1] += np.sum(self.varianceU*dkdvar * dL_dK)
        target[2] += np.sum(dktheta1*(-np.sqrt(3)*self.lengthscaleU**(-2)) * dL_dK)
        target[3] += np.sum(dktheta2*(-self.lengthscaleY**(-2)) * dL_dK)
1075066
def tau_davies(Z, fmax, t, y, dy, normalization='standard', dH=1, dK=3):
    """"""
    N = len(t)
    NH = N - dH  # DOF for null hypothesis
    NK = N - dK  # DOF for periodic hypothesis
    Dt = _weighted_var(t, dy)
    Teff = np.sqrt(4 * np.pi * Dt)
    W = fmax * Teff
    if normalization == 'psd':
        # 'psd' normalization is same as Baluev's z
        return W * np.exp(-Z) * np.sqrt(Z)
    elif normalization == 'standard':
        # 'standard' normalization is Z = 2/NH * z_1
        return (_gamma(NH) * W * (1 - Z) ** (0.5 * (NK - 1))
                * np.sqrt(0.5 * NH * Z))
    elif normalization == 'model':
        # 'model' normalization is Z = 2/NK * z_2
        return (_gamma(NK) * W * (1 + Z) ** (-0.5 * NK)
                * np.sqrt(0.5 * NK * Z))
    elif normalization == 'log':
        # 'log' normalization is Z = 2/NK * z_3
        return (_gamma(NK) * W * np.exp(-0.5 * Z * (NK - 0.5))
                * np.sqrt(NK * np.sinh(0.5 * Z)))
    else:
        raise NotImplementedError("normalization={0}".format(normalization))
270408
def sde(self): 
         
        
        variance = float(self.variances.values) # this is initial variancve in Bayesian linear regression
        t0 = float(self.t0)
        
        F = np.array( ((0,1.0),(0,0) ))
        L = np.array( ((0,),(1.0,)) )
        Qc = np.zeros((1,1))
        H = np.array( ((1.0,0),) )
        
        Pinf   = np.zeros((2,2))
        P0 = np.array( ( (t0**2, t0), (t0, 1) ) ) * variance        
        dF = np.zeros((2,2,1))
        dQc    = np.zeros( (1,1,1) )
        
        dPinf = np.zeros((2,2,1))
        dP0 = np.zeros((2,2,1))
        dP0[:,:,0]  = P0 / variance
  
        return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
141532
def _fD(self, neg_pairs, A):
    
    diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
    return np.log(np.sum(np.sqrt(np.sum(np.dot(diff, A) * diff, axis=1))) + 1e-6)
791933
def iterativeEMDS(X0, N, d, C, b, max_it=10, print_out=False, **kwargs):
    
    from pylocus.basics import mse, projection
    KE = kwargs.get('KE', None)
    KE_projected = KE.copy()
    d = len(X0)
    for i in range(max_it):
        # projection on constraints
        KE_projected, cost, __ = projection(KE_projected, C, b)
        rank = np.linalg.matrix_rank(KE_projected)

        # rank 2 decomposition
        Xhat_KE, Vhat_KE = superMDS(X0, N, d, KE=KE_projected)
        KE_projected = Vhat_KE.dot(Vhat_KE.T)

        error = mse(C.dot(KE_projected), b)

        if (print_out):
            print('cost={:2.2e},error={:2.2e}, rank={}'.format(
                cost, error, rank))
        if cost < 1e-20 and error < 1e-20 and rank == d:
            if (print_out):
                print('converged after {} iterations'.format(i))
            return Xhat_KE, Vhat_KE
    print('iterativeMDS did not converge!')
    return None, None
980362
def _lmder1_powell_singular():
    """"""

    def func(params, vec):
        vec[0] = params[0] + 10 * params[1]
        vec[1] = np.sqrt(5) * (params[2] - params[3])
        vec[2] = (params[1] - 2 * params[2])**2
        vec[3] = np.sqrt(10) * (params[0] - params[3])**2

    def jac(params, jac):
        jac.fill(0)
        jac[0,0] = 1
        jac[0,3] = 2 * np.sqrt(10) * (params[0] - params[3])
        jac[1,0] = 10
        jac[1,2] = 2 * (params[1] - 2 * params[2])
        jac[2,1] = np.sqrt(5)
        jac[2,2] = -2 * jac[2,1]
        jac[3,1] = -np.sqrt(5)
        jac[3,3] = -jac[3,0]

    guess = np.asfarray([3, -1, 0, 1])

    _lmder1_test(4, func, jac, guess)
    _lmder1_test(4, func, jac, guess * 10)
    _lmder1_test(4, func, jac, guess * 100)
1106369
def _der(self,x):
        '''
        
        '''
        if _isscalar(x):
            pos = np.searchsorted(self.x_list,x)
            if pos == 0:
                dydx = self.coeffs[0,1]
            elif (pos < self.n):
                alpha = (x - self.x_list[pos-1])/(self.x_list[pos] - self.x_list[pos-1])
                dydx = (self.coeffs[pos,1] + alpha*(2*self.coeffs[pos,2] + alpha*3*self.coeffs[pos,3]))/(self.x_list[pos] - self.x_list[pos-1])
            else:
                alpha = x - self.x_list[self.n-1]
                dydx = self.coeffs[pos,1] - self.coeffs[pos,2]*self.coeffs[pos,3]*np.exp(alpha*self.coeffs[pos,3])
        else:
            m = len(x)
            pos = np.searchsorted(self.x_list,x)
            dydx = np.zeros(m)
            if dydx.size > 0:
                out_bot   = pos == 0
                out_top   = pos == self.n
                in_bnds   = np.logical_not(np.logical_or(out_bot, out_top))

                # Do the "in bounds" evaluation points
                i = pos[in_bnds]
                coeffs_in = self.coeffs[i,:]
                alpha = (x[in_bnds] - self.x_list[i-1])/(self.x_list[i] - self.x_list[i-1])
                dydx[in_bnds] = (coeffs_in[:,1] + alpha*(2*coeffs_in[:,2] + alpha*3*coeffs_in[:,3]))/(self.x_list[i] - self.x_list[i-1])

                # Do the "out of bounds" evaluation points
                dydx[out_bot] = self.coeffs[0,1]
                alpha = x[out_top] - self.x_list[self.n-1]
                dydx[out_top] = self.coeffs[self.n,1] - self.coeffs[self.n,2]*self.coeffs[self.n,3]*np.exp(alpha*self.coeffs[self.n,3])
        return dydx
659844
def hinit(func, x, t, pos_neg, f0, iord, hmax, rtol, atol, args):
    
    sk = atol + rtol * np.fabs(x)
    dnf = np.sum(np.square(f0 / sk), axis=0)
    dny = np.sum(np.square(x / sk), axis=0)

    h = np.sqrt(dny / dnf) * 0.01

    h = np.min([h, np.fabs(hmax)])
    h = custom_sign(h, pos_neg)

    # perform an explicit Euler step
    xx1 = x + h * f0
    f1 = np.array(func(xx1, t[0] + h, *args))

    # estimate the second derivative of the solution
    der2 = np.sum(np.square((f1 - f0) / sk), axis=0)
    der2 = np.sqrt(der2) / h

    # step size is computed such that h ** iord * max_d(norm(f0), norm(der2)) = 0.01
    der12 = np.max([np.fabs(der2), np.sqrt(dnf)])
    h1 = np.power(0.01 / der12, 1.0 / iord)

    h = np.min([100.0 * np.fabs(h), np.min([np.fabs(h1), np.fabs(hmax)])])

    return custom_sign(h, pos_neg), f0, f1, xx1
10617
function f1(d, m, r) {
    const t = m + d;
    const I = (t << r) | (t >>> (32 - r));
    return ((sBox[0][I >>> 24] ^ sBox[1][(I >>> 16) & 255]) - sBox[2][(I >>> 8) & 255]) + sBox[3][I & 255];
  }
179612
def _mse_converged(self):
        

        mse = mean_squared_error(self.local_prior, self.local_posterior_,
                                 multioutput='uniform_average')
        if mse > self.threshold:
            return False, mse
        else:
            return True, mse
66008
def hessian_local_log_likelihood(self, x):
        
        C, D, u, y = self.C, self.D, self.inputs, self.data
        psi = x.dot(C.T) + u.dot(D.T)
        p = 1. / (1 + np.exp(-psi))
        dp_dpsi = p * (1 - p)
        return np.einsum('tn, ni, nj ->tij', -dp_dpsi, self.C, self.C)
198858
function incrmvariance( W ) {
	if ( !isPositiveInteger( W ) ) {
		throw new TypeError( 'incrmvariance()::invalid input argument. Window size must be a positive integer. Value: `' + W + '`.' );
	}
	var arr = new Array( W ),
		n = W - 1,
		M2 = 0,
		mu = 0,
		N = 0,
		i = -1,
		delta,
		tmp,
		d1,
		d2;
	/**
	* FUNCTION: incrmvariance( [value] )
	*	If a `value` is provided, updates and returns the updated sample variance. If no `value` is provided, returns the current sample variance.
	*
	* @param {Number} [value] - value used to update the moving sample variance
	* @returns {Number|Null} sample variance or null
	*/
	return function incrmvariance( x ) {
		if ( !arguments.length ) {
			if ( N === 0 ) {
				return null;
			}
			if ( N === 1 ) {
				return 0;
			}
			if ( N < W ) {
				return M2 / (N-1);
			}
			return M2 / n;
		}
		// Update the index for managing the circular buffer...
		i = (i+1) % W;

		// Fill up the initial window; else, update the existing window...
		if ( N < W ) {
			arr[ i ] = x;
			N += 1;
			delta = x - mu;
			mu += delta / N;
			M2 += delta * (x - mu);
			if ( N === 1 ) {
				return 0;
			}
			return M2 / (N-1);
		}
		if ( N === 1 ) {
			return 0;
		}
		tmp = arr[ i ];
		arr[ i ] = x;
		delta = x - tmp;
		d1 = tmp - mu;
		mu += delta / W;
		d2 = x - mu;
		M2 += delta * (d1 + d2);

		return M2 / n;
	};
}
266969
def estimate_k(X, max_k):
    
    ks = range(1, max_k)
    fs = np.zeros(len(ks))

    # Special case K=1
    fs[0], Sk = _fK(1)

    # Rest of Ks
    for k in ks[1:]:
        fs[k-1], Sk = _fK(k, Skm1=Sk)
    return np.argmin(fs) + 1