Base calculation for MORAN_BV_Matrix
333910
def _Moran_BV_Matrix_array(variables, w, permutations=0, varnames=None):
    
    if varnames is None:
        varnames = ['x{}'.format(i) for i in range(k)]

    k = len(variables)
    rk = list(range(0, k - 1))
    results = {}
    for i in rk:
        for j in range(i + 1, k):
            y1 = variables[i]
            y2 = variables[j]
            results[i, j] = Moran_BV(y1, y2, w, permutations=permutations)
            results[j, i] = Moran_BV(y2, y1, w, permutations=permutations)
            results[i, j].varnames = {'x': varnames[i], 'y': varnames[j]}
            results[j, i].varnames = {'x': varnames[j], 'y': varnames[i]}
    return results
962851
def _get_a2_value(bbar, dbar, slip, beta, mmax):
        
        return ((dbar - bbar) / bbar) * (slip / beta) *\
            np.exp(-(dbar / 2.) * mmax)
1083012
def cublasSsbmv(handle, uplo, n, k, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasSsbmv_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], n, k,
                                       ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(ctypes.c_float(beta)),
                                       int(y), incy)
    cublasCheckStatus(status)
413672
def compute_n_splits(cv, X, y=None, groups=None):
    
    if not any(is_dask_collection(i) for i in (X, y, groups)):
        return cv.get_n_splits(X, y, groups)

    if isinstance(cv, (_BaseKFold, BaseShuffleSplit)):
        return cv.n_splits

    elif isinstance(cv, PredefinedSplit):
        return len(cv.unique_folds)

    elif isinstance(cv, _CVIterableWrapper):
        return len(cv.cv)

    elif isinstance(cv, (LeaveOneOut, LeavePOut)) and not is_dask_collection(X):
        # Only `X` is referenced for these classes
        return cv.get_n_splits(X, None, None)

    elif isinstance(cv, (LeaveOneGroupOut, LeavePGroupsOut)) and not is_dask_collection(
        groups
    ):
        # Only `groups` is referenced for these classes
        return cv.get_n_splits(None, None, groups)

    else:
        return delayed(cv).get_n_splits(X, y, groups).compute()
1035019
private static int[] elmhes(DenseMatrix A) {
        int n = A.nrows();
        int[] perm = new int[n];

        for (int m = 1; m < n - 1; m++) {
            double x = 0.0;
            int i = m;
            for (int j = m; j < n; j++) {
                if (Math.abs(A.get(j, m - 1)) > Math.abs(x)) {
                    x = A.get(j, m - 1);
                    i = j;
                }
            }
            perm[m] = i;
            if (i != m) {
                for (int j = m - 1; j < n; j++) {
                    double swap = A.get(i, j);
                    A.set(i, j, A.get(m, j));
                    A.set(m, j, swap);
                }
                for (int j = 0; j < n; j++) {
                    double swap = A.get(j, i);
                    A.set(j, i, A.get(j, m));
                    A.set(j, m, swap);
                }
            }
            if (x != 0.0) {
                for (i = m + 1; i < n; i++) {
                    double y = A.get(i, m - 1);
                    if (y != 0.0) {
                        y /= x;
                        A.set(i, m - 1, y);
                        for (int j = m; j < n; j++) {
                            A.sub(i, j, y * A.get(m, j));
                        }
                        for (int j = 0; j < n; j++) {
                            A.add(j, m, y * A.get(j, i));
                        }
                    }
                }
            }
        }

        return perm;
    }
552652
def MarginalBeta(self, i):
        
        alpha0 = self.params.sum()
        alpha = self.params[i]
        return Beta(alpha, alpha0 - alpha)
1285240
def spline_base3d( width, height, depth, nr_knots_x = 10.0, nr_knots_y = 10.0,
        nr_knots_z=10, spline_order = 3, marginal_x = None, marginal_y = None, 
        marginal_z = None):
      
    if not nr_knots_z < depth:
        raise RuntimeError("Too many knots for size of the base")
    basis2d, (knots_x, knots_y) = spline_base2d(height, width, nr_knots_x, 
            nr_knots_y, spline_order, marginal_x, marginal_y)
    if marginal_z is not None:
        knots_z = knots_from_marginal(marginal_z, nr_knots_z, spline_order)
    else:
        knots_z = augknt(np.linspace(0,depth+1, nr_knots_z), spline_order)
    z_eval = np.arange(1,depth+1).astype(float)
    spline_setz = spcol(z_eval, knots_z, spline_order)
    bspline = np.zeros((basis2d.shape[0]*len(z_eval), height*width*depth))
    basis_nr = 0
    for spline_a in spline_setz.T:
        for spline_b in basis2d:
            spline_b = spline_b.reshape((height, width))
            bspline[basis_nr, :] = (spline_b[:,:,np.newaxis] * spline_a[:]).flat
            basis_nr +=1
    return bspline, (knots_x, knots_y, knots_z)
539281
function bfrac(a, b, x, y) {
      var lambda;

      lambda = a - (a + b) * x;
      return brcomp(a, b, x, y) + Math.log(contFrac(
         function(i) { // betas from DiDonato
            var n;

            n = i - 1;
            return i === 0 ? 0
                           : n + n * (b - n) * x / (a + 2 * n - 1) +
                             (a + n) / (a + 2 * n + 1) *
                                (lambda + 1 + n * (1 + y));
         },
         function(i) { // alphas from DiDonato
            var n;

            n = i - 1;
            return i === 1 ? 1
                           : (a + n - 1) * (a + b + n - 1) /
                             (a + 2 * n - 1) / (a + 2 * n - 1) *
                             n * (b - n) * x * x;
         }
      ));
   }
241802
def tensor_product_basis(*bases):
    
    dim = np.prod([basis.data.shape[1] for basis in bases])
    tp_basis = np.zeros((dim**2, dim, dim), dtype=complex)

    for idx_factors, factors in enumerate(it.product(*[basis.data for basis in bases])):
        tp_basis[idx_factors, :, :] = reduce(np.kron, factors)

    return TomographyBasis(tp_basis,
        sum((
            factor.dims for factor in bases
        ), []),
        list(map(
        r"\otimes".join,
        it.product(*[
            basis.labels for basis in bases
        ])
    )))
64137
function computeBaseline(data) {
    var layerNum = data.length;
    var pointNum = data[0].length;
    var sums = [];
    var y0 = [];
    var max = 0;
    var temp;
    var base = {};

    for (var i = 0; i < pointNum; ++i) {
        for (var j = 0, temp = 0; j < layerNum; ++j) {
            temp += data[j][i][1];
        }
        if (temp > max) {
            max = temp;
        }
        sums.push(temp);
    }

    for (var k = 0; k < pointNum; ++k) {
        y0[k] = (max - sums[k]) / 2;
    }
    max = 0;

    for (var l = 0; l < pointNum; ++l) {
        var sum = sums[l] + y0[l];
        if (sum > max) {
            max = sum;
        }
    }
    base.y0 = y0;
    base.max = max;

    return base;
}
948400
def sim_combi(self):
        
        n_samples = 500
        sigma_glob = 1.8

        X = np.zeros((n_samples,3))

        X[:,0] = np.random.uniform(-sigma_glob,sigma_glob,n_samples)
        X[:,1] = np.random.uniform(-sigma_glob,sigma_glob,n_samples)

        func = self.funcs['tanh']

        # XOR type
#         X[:,2] = (func(X[:,0])*sp.stats.norm.pdf(X[:,1],0,0.2)
#                   + func(X[:,1])*sp.stats.norm.pdf(X[:,0],0,0.2))
        # AND type / diagonal
#         X[:,2] = (func(X[:,0]+X[:,1])*sp.stats.norm.pdf(X[:,1]-X[:,0],0,0.2))
        # AND type / horizontal
        X[:,2] = (func(X[:,0])*sp.stats.norm.cdf(X[:,1],1,0.2))

        pl.scatter(X[:,0],X[:,1],c=X[:,2],edgecolor='face')
        pl.show()

        pl.plot(X[:,1],X[:,2],'.')
        pl.show()

        return X
824065
def dImage_wrt_2dVerts_bnd(observed, visible, visibility, barycentric, image_width, image_height, num_verts, f, bnd_bool):
    """"""

    n_channels = np.atleast_3d(observed).shape[2]
    shape = visibility.shape

    # Step 1: get the structure ready, ie the IS and the JS
    IS = np.tile(col(visible), (1, 2*f.shape[1])).ravel()
    JS = col(f[visibility.ravel()[visible]].ravel())
    JS = np.hstack((JS*2, JS*2+1)).ravel()

    pxs = np.asarray(visible % shape[1], np.int32)
    pys = np.asarray(np.floor(np.floor(visible) / shape[1]), np.int32)

    if n_channels > 1:
        IS = np.concatenate([IS*n_channels+i for i in range(n_channels)])
        JS = np.concatenate([JS for i in range(n_channels)])

    # Step 2: get the data ready, ie the actual values of the derivatives
    ksize = 1
    bndf = bnd_bool.astype(np.float64)
    nbndf = np.logical_not(bnd_bool).astype(np.float64)
    sobel_normalizer = cv2.Sobel(np.asarray(np.tile(row(np.arange(10)), (10, 1)), np.float64), cv2.CV_64F, dx=1, dy=0, ksize=ksize)[5,5]

    bnd_nan = bndf.reshape((observed.shape[0], observed.shape[1], -1)).copy()
    bnd_nan.ravel()[bnd_nan.ravel()>0] = np.nan
    bnd_nan += 1
    obs_nonbnd = np.atleast_3d(observed) * bnd_nan

    ydiffnb, xdiffnb = nangradients(obs_nonbnd)

    observed = np.atleast_3d(observed)

    if observed.shape[2] > 1:
        ydiffbnd, xdiffbnd, _ = np.gradient(observed)
    else:
        ydiffbnd, xdiffbnd = np.gradient(observed.squeeze())
        ydiffbnd = np.atleast_3d(ydiffbnd)
        xdiffbnd = np.atleast_3d(xdiffbnd)

    # This corrects for a bias imposed boundary differences begin spread over two pixels
    # (by np.gradients or similar) but only counted once (since OpenGL's line
    # drawing spans 1 pixel)
    xdiffbnd *= 2.0
    ydiffbnd *= 2.0

    xdiffnb = -xdiffnb
    ydiffnb = -ydiffnb
    xdiffbnd = -xdiffbnd
    ydiffbnd = -ydiffbnd
    # ydiffnb *= 0
    # xdiffnb *= 0

    if False:
        import matplotlib.pyplot as plt
        plt.figure()
        plt.subplot(121)
        plt.imshow(xdiffnb)
        plt.title('xdiffnb')
        plt.subplot(122)
        plt.imshow(xdiffbnd)
        plt.title('xdiffbnd')
        import pdb; pdb.set_trace()

    idxs = np.isnan(xdiffnb.ravel())
    xdiffnb.ravel()[idxs] = xdiffbnd.ravel()[idxs]

    idxs = np.isnan(ydiffnb.ravel())
    ydiffnb.ravel()[idxs] = ydiffbnd.ravel()[idxs]

    if True: # should be right thing
        xdiff = xdiffnb
        ydiff = ydiffnb
    else:  #should be old way
        xdiff = xdiffbnd
        ydiff = ydiffbnd


    # TODO: NORMALIZER IS WRONG HERE
    # xdiffnb = -cv2.Sobel(obs_nonbnd, cv2.CV_64F, dx=1, dy=0, ksize=ksize) / np.atleast_3d(cv2.Sobel(row(np.arange(obs_nonbnd.shape[1])).astype(np.float64), cv2.CV_64F, dx=1, dy=0, ksize=ksize))
    # ydiffnb = -cv2.Sobel(obs_nonbnd, cv2.CV_64F, dx=0, dy=1, ksize=ksize) / np.atleast_3d(cv2.Sobel(col(np.arange(obs_nonbnd.shape[0])).astype(np.float64), cv2.CV_64F, dx=0, dy=1, ksize=ksize))
    #
    # xdiffnb.ravel()[np.isnan(xdiffnb.ravel())] = 0.
    # ydiffnb.ravel()[np.isnan(ydiffnb.ravel())] = 0.
    # xdiffnb.ravel()[np.isinf(xdiffnb.ravel())] = 0.
    # ydiffnb.ravel()[np.isinf(ydiffnb.ravel())] = 0.

    # xdiffnb = np.atleast_3d(xdiffnb)
    # ydiffnb = np.atleast_3d(ydiffnb)
    #
    # xdiffbnd = -cv2.Sobel(observed, cv2.CV_64F, dx=1, dy=0, ksize=ksize) / sobel_normalizer
    # ydiffbnd = -cv2.Sobel(observed, cv2.CV_64F, dx=0, dy=1, ksize=ksize) / sobel_normalizer
    #
    # xdiff = xdiffnb * np.atleast_3d(nbndf)
    # xdiff.ravel()[np.isnan(xdiff.ravel())] = 0
    # xdiff += xdiffbnd*np.atleast_3d(bndf)
    #
    # ydiff = ydiffnb * np.atleast_3d(nbndf)
    # ydiff.ravel()[np.isnan(ydiff.ravel())] = 0
    # ydiff += ydiffbnd*np.atleast_3d(bndf)

    #import pdb; pdb.set_trace()

    #xdiff = xdiffnb
    #ydiff = ydiffnb

    #import pdb; pdb.set_trace()

    datas = []

    # The data is weighted according to barycentric coordinates
    bc0 = col(barycentric[pys, pxs, 0])
    bc1 = col(barycentric[pys, pxs, 1])
    bc2 = col(barycentric[pys, pxs, 2])
    for k in range(n_channels):
        dxs = xdiff[pys, pxs, k]
        dys = ydiff[pys, pxs, k]
        if f.shape[1] == 3:
            datas.append(np.hstack((col(dxs)*bc0,col(dys)*bc0,col(dxs)*bc1,col(dys)*bc1,col(dxs)*bc2,col(dys)*bc2)).ravel())
        else:
            datas.append(np.hstack((col(dxs)*bc0,col(dys)*bc0,col(dxs)*bc1,col(dys)*bc1)).ravel())

    data = np.concatenate(datas)

    ij = np.vstack((IS.ravel(), JS.ravel()))
    result = sp.csc_matrix((data, ij), shape=(image_width*image_height*n_channels, num_verts*2))

    return result
241475
def d2Ibr_dV2(Ybr, V, lam):
    
    nb = len(V)
    diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = spdiag(mul(-(Ybr.T * lam), V))
    Hva = -1j * Haa * diaginvVm
    Hav = Hva
    Hvv = spmatrix([], [], [], (nb, nb))

    return Haa, Hav, Hva, Hvv
942549
def d2Sbr_dV2(Cbr, Ybr, V, lam):
    
    nb = len(V)

    diaglam = spdiag(lam)
    diagV = spdiag(V)

    A = Ybr.H * diaglam * Cbr
    B = conj(diagV) * A * diagV
    D = spdiag(mul((A*V), conj(V)))
    E = spdiag(mul((A.T * conj(V)), V))
    F = B + B.T
    G = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = F - D - E
    Hva = 1j * G * (B - B.T - D + E)
    Hav = Hva.T
    Hvv = G * F * G

    return Haa, Hav, Hva, Hvv
370157
private static void hqr(DenseMatrix A, double[] d, double[] e) {
        int n = A.nrows();
        int nn, m, l, k, j, its, i, mmin;
        double z, y, x, w, v, u, t, s, r = 0.0, q = 0.0, p = 0.0, anorm = 0.0;

        for (i = 0; i < n; i++) {
            for (j = Math.max(i - 1, 0); j < n; j++) {
                anorm += Math.abs(A.get(i, j));
            }
        }
        nn = n - 1;
        t = 0.0;
        while (nn >= 0) {
            its = 0;
            do {
                for (l = nn; l > 0; l--) {
                    s = Math.abs(A.get(l - 1, l - 1) + Math.abs(A.get(l, l)));
                    if (s == 0.0) {
                        s = anorm;
                    }
                    if (Math.abs(A.get(l, l - 1)) <= Math.EPSILON * s) {
                        A.set(l, l - 1, 0.0);
                        break;
                    }
                }
                x = A.get(nn, nn);
                if (l == nn) {
                    d[nn--] = x + t;
                } else {
                    y = A.get(nn - 1, nn - 1);
                    w = A.get(nn, nn - 1) * A.get(nn - 1, nn);
                    if (l == nn - 1) {
                        p = 0.5 * (y - x);
                        q = p * p + w;
                        z = Math.sqrt(Math.abs(q));
                        x += t;
                        if (q >= 0.0) {
                            z = p + Math.copySign(z, p);
                            d[nn - 1] = d[nn] = x + z;
                            if (z != 0.0) {
                                d[nn] = x - w / z;
                            }
                        } else {
                            d[nn] = x + p;
                            e[nn] = -z;
                            d[nn - 1] = d[nn];
                            e[nn - 1] = -e[nn];
                        }
                        nn -= 2;
                    } else {
                        if (its == 30) {
                            throw new IllegalStateException("Too many iterations in hqr");
                        }
                        if (its == 10 || its == 20) {
                            t += x;
                            for (i = 0; i < nn + 1; i++) {
                                A.sub(i, i, x);
                            }
                            s = Math.abs(A.get(nn, nn - 1)) + Math.abs(A.get(nn - 1, nn - 2));
                            y = x = 0.75 * s;
                            w = -0.4375 * s * s;
                        }
                        ++its;
                        for (m = nn - 2; m >= l; m--) {
                            z = A.get(m, m);
                            r = x - z;
                            s = y - z;
                            p = (r * s - w) / A.get(m + 1, m) + A.get(m, m + 1);
                            q = A.get(m + 1, m + 1) - z - r - s;
                            r = A.get(m + 2, m + 1);
                            s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                            p /= s;
                            q /= s;
                            r /= s;
                            if (m == l) {
                                break;
                            }
                            u = Math.abs(A.get(m, m - 1)) * (Math.abs(q) + Math.abs(r));
                            v = Math.abs(p) * (Math.abs(A.get(m - 1, m - 1)) + Math.abs(z) + Math.abs(A.get(m + 1, m + 1)));
                            if (u <= Math.EPSILON * v) {
                                break;
                            }
                        }
                        for (i = m; i < nn - 1; i++) {
                            A.set(i + 2, i, 0.0);
                            if (i != m) {
                                A.set(i + 2, i - 1, 0.0);
                            }
                        }
                        for (k = m; k < nn; k++) {
                            if (k != m) {
                                p = A.get(k, k - 1);
                                q = A.get(k + 1, k - 1);
                                r = 0.0;
                                if (k + 1 != nn) {
                                    r = A.get(k + 2, k - 1);
                                }
                                if ((x = Math.abs(p) + Math.abs(q) + Math.abs(r)) != 0.0) {
                                    p /= x;
                                    q /= x;
                                    r /= x;
                                }
                            }
                            if ((s = Math.copySign(Math.sqrt(p * p + q * q + r * r), p)) != 0.0) {
                                if (k == m) {
                                    if (l != m) {
                                        A.set(k, k - 1, -A.get(k, k - 1));
                                    }
                                } else {
                                    A.set(k, k - 1, -s * x);
                                }
                                p += s;
                                x = p / s;
                                y = q / s;
                                z = r / s;
                                q /= p;
                                r /= p;
                                for (j = k; j < nn + 1; j++) {
                                    p = A.get(k, j) + q * A.get(k + 1, j);
                                    if (k + 1 != nn) {
                                        p += r * A.get(k + 2, j);
                                        A.sub(k + 2, j, p * z);
                                    }
                                    A.sub(k + 1, j, p * y);
                                    A.sub(k, j, p * x);
                                }
                                mmin = nn < k + 3 ? nn : k + 3;
                                for (i = l; i < mmin + 1; i++) {
                                    p = x * A.get(i, k) + y * A.get(i, k + 1);
                                    if (k + 1 != nn) {
                                        p += z * A.get(i, k + 2);
                                        A.sub(i, k + 2, p * r);
                                    }
                                    A.sub(i, k + 1, p * q);
                                    A.sub(i, k, p);
                                }
                            }
                        }
                    }
                }
            } while (l + 1 < nn);
        }
    }
744886
def dSbus_dV(Y, V):
    
    I = Y * V

    diagV = spdiag(V)
    diagIbus = spdiag(I)
    diagVnorm = spdiag(div(V, abs(V))) # Element-wise division.

    dS_dVm = diagV * conj(Y * diagVnorm) + conj(diagIbus) * diagVnorm
    dS_dVa = 1j * diagV * conj(diagIbus - Y * diagV)

    return dS_dVm, dS_dVa
1109909
function computeDimension(topLeft, topRight, bottomLeft, matrix) {
        const moduleSize = (sum(countBlackWhiteRun(topLeft, bottomLeft, matrix, 5)) / 7 + // Divide by 7 since the ratio is 1:1:3:1:1
            sum(countBlackWhiteRun(topLeft, topRight, matrix, 5)) / 7 +
            sum(countBlackWhiteRun(bottomLeft, topLeft, matrix, 5)) / 7 +
            sum(countBlackWhiteRun(topRight, topLeft, matrix, 5)) / 7) / 4;
        if (moduleSize < 1) {
            throw new Error("Invalid module size");
        }
        const topDimension = Math.round(distance(topLeft, topRight) / moduleSize);
        const sideDimension = Math.round(distance(topLeft, bottomLeft) / moduleSize);
        let dimension = Math.floor((topDimension + sideDimension) / 2) + 7;
        switch (dimension % 4) {
            case 0:
                dimension++;
                break;
            case 2:
                dimension--;
                break;
        }
        return { dimension, moduleSize };
    }