Create strength of connection matrixfrom a function applied to relaxation vectors.
337614
def distance_measure_common(A, func, alpha, R, k, epsilon):
    """"""
    # create test vectors
    x = relaxation_vectors(A, R, k, alpha)

    # apply distance measure function to vectors
    d = func(x)

    # drop distances to self
    (rows, cols) = A.nonzero()
    weak = np.where(rows == cols)[0]
    d[weak] = 0
    C = sparse.csr_matrix((d, (rows, cols)), shape=A.shape)
    C.eliminate_zeros()

    # remove weak connections
    # removes entry e from a row if e > theta * min of all entries in the row
    amg_core.apply_distance_filter(C.shape[0], epsilon, C.indptr,
                                   C.indices, C.data)
    C.eliminate_zeros()

    # Standardized strength values require small values be weak and large
    # values be strong.  So, we invert the distances.
    C.data = 1.0 / C.data

    # Put an identity on the diagonal
    C = C + sparse.eye(C.shape[0], C.shape[1], format='csr')

    # Scale C by the largest magnitude entry in each row
    C = scale_rows_by_largest_entry(C)

    return C
188936
def fitlin(imgarr,refarr):
    
    # Initialize variables
    _mat = np.zeros((3,3),dtype=np.float64)
    _xorg = imgarr[0][0]
    _yorg = imgarr[0][1]
    _xoorg = refarr[0][0]
    _yoorg = refarr[0][1]
    _sigxox = 0.
    _sigxoy = 0.
    _sigxo = 0.
    _sigyox = 0.
    _sigyoy = 0.
    _sigyo = 0.

    _npos = len(imgarr)
    # Populate matrices
    for i in range(_npos):
        _mat[0][0] += np.power((imgarr[i][0] - _xorg),2)
        _mat[0][1] += (imgarr[i][0] - _xorg) * (imgarr[i][1] - _yorg)
        _mat[0][2] += (imgarr[i][0] - _xorg)
        _mat[1][1] += np.power((imgarr[i][1] - _yorg),2)
        _mat[1][2] += imgarr[i][1] - _yorg

        _sigxox += (refarr[i][0] - _xoorg)*(imgarr[i][0] - _xorg)
        _sigxoy += (refarr[i][0] - _xoorg)*(imgarr[i][1] - _yorg)
        _sigxo += refarr[i][0] - _xoorg
        _sigyox += (refarr[i][1] - _yoorg)*(imgarr[i][0] -_xorg)
        _sigyoy += (refarr[i][1] - _yoorg)*(imgarr[i][1] - _yorg)
        _sigyo += refarr[i][1] - _yoorg

    _mat[2][2] = _npos
    _mat[1][0] = _mat[0][1]
    _mat[2][0] = _mat[0][2]
    _mat[2][1] = _mat[1][2]

    # Now invert this matrix
    _mat = linalg.inv(_mat)

    _a  = _sigxox*_mat[0][0]+_sigxoy*_mat[0][1]+_sigxo*_mat[0][2]
    _b  = -1*(_sigxox*_mat[1][0]+_sigxoy*_mat[1][1]+_sigxo*_mat[1][2])
    #_x0 = _sigxox*_mat[2][0]+_sigxoy*_mat[2][1]+_sigxo*_mat[2][2]

    _c  = _sigyox*_mat[1][0]+_sigyoy*_mat[1][1]+_sigyo*_mat[1][2]
    _d  = _sigyox*_mat[0][0]+_sigyoy*_mat[0][1]+_sigyo*_mat[0][2]
    #_y0 = _sigyox*_mat[2][0]+_sigyoy*_mat[2][1]+_sigyo*_mat[2][2]

    _xt = _xoorg - _a*_xorg+_b*_yorg
    _yt = _yoorg - _d*_xorg-_c*_yorg

    return [_a,_b,_xt],[_c,_d,_yt]
474556
def connectivity(ax):

    ''''''
    masked_array = np.ma.array(params.C_YX, mask=params.C_YX==0)
    # if analysis_params.bw:
    #     cmap = plt.get_cmap(gray, 20)
    #     cmap.set_bad('k', 1.)
    # else:
    #     cmap = plt.get_cmap('hot', 20)
    #     cmap.set_bad('k', 0.5)
    # im = ax.imshow(masked_array, cmap=cmap, vmin=0, interpolation='nearest')
    im = ax.pcolormesh(masked_array, cmap=cmap, vmin=0, ) #interpolation='nearest')
    ax.axis(ax.axis('tight'))
    ax.invert_yaxis()
    ax.xaxis.set_ticks_position('top')
    ax.set_xticks(np.arange(9)+0.5)
    ax.set_yticks(np.arange(8)+0.5)
    ax.set_xticklabels(params.X, rotation=270)
    ax.set_yticklabels(params.Y, )
    ax.xaxis.set_label_position('top')
    ax.set_xlabel(r'$X$', labelpad=-1,fontsize=8)
    ax.set_ylabel(r'$Y$', labelpad=0, rotation=0,fontsize=8)

    rect = np.array(ax.get_position().bounds)
    rect[0] += rect[2] + 0.01
    rect[2] = 0.01
    fig = plt.gcf()
    cax = fig.add_axes(rect)

    cbar = plt.colorbar(im, cax=cax)
    #cbar.set_label(r'$C_{YX}$', ha='center')
    cbar.set_label(r'$C_{YX}$', labelpad=0)
1049357
def compute_laplacian_matrix(affinity_matrix, method='auto', **kwargs):
    """"""
    if method == 'auto':
        method = 'geometric'
    return Laplacian.init(method, **kwargs).laplacian_matrix(affinity_matrix)
396245
def build_grad_matrices(V, points):
    
    # See <https://www.allanswered.com/post/lkbkm/#zxqgk>
    mesh = V.mesh()

    bbt = BoundingBoxTree()
    bbt.build(mesh)
    dofmap = V.dofmap()
    el = V.element()
    rows = []
    cols = []
    datax = []
    datay = []
    for i, xy in enumerate(points):
        cell_id = bbt.compute_first_entity_collision(Point(*xy))
        cell = Cell(mesh, cell_id)
        coordinate_dofs = cell.get_vertex_coordinates()

        rows.append([i, i, i])
        cols.append(dofmap.cell_dofs(cell_id))

        v = el.evaluate_basis_derivatives_all(1, xy, coordinate_dofs, cell_id)
        v = v.reshape(3, 2)
        datax.append(v[:, 0])
        datay.append(v[:, 1])

    rows = numpy.concatenate(rows)
    cols = numpy.concatenate(cols)
    datax = numpy.concatenate(datax)
    datay = numpy.concatenate(datay)

    m = len(points)
    n = V.dim()
    dx_matrix = sparse.csr_matrix((datax, (rows, cols)), shape=(m, n))
    dy_matrix = sparse.csr_matrix((datay, (rows, cols)), shape=(m, n))
    return dx_matrix, dy_matrix
163638
def mkron(a, *args):
    """"""
    if not isinstance(a, list):
        a = [a]
    a = list(a)  # copy list
    for i in args:
        if isinstance(i, list):
            a.extend(i)
        else:
            a.append(i)

    c = _vector.vector()
    c.d = 0
    c.n = _np.array([], dtype=_np.int32)
    c.r = _np.array([], dtype=_np.int32)
    c.core = []

    for t in a:
        thetensor = t.tt if isinstance(t, _matrix.matrix) else t
        c.d += thetensor.d
        c.n = _np.concatenate((c.n, thetensor.n))
        c.r = _np.concatenate((c.r[:-1], thetensor.r))
        c.core = _np.concatenate((c.core, thetensor.core))

    c.get_ps()
    return c
263612
def _fully_connected(self, x, out_dim):
        """"""
        x = tf.reshape(x, [self.hps.batch_size, -1])
        w = tf.get_variable(
            "DW", [x.get_shape()[1], out_dim],
            initializer=tf.uniform_unit_scaling_initializer(factor=1.0))
        b = tf.get_variable(
            "biases", [out_dim], initializer=tf.constant_initializer())
        return tf.nn.xw_plus_b(x, w, b)
169894
def relaxation_operators(p):
    
    k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - p)]])
    k1 = np.array([[0.0, np.sqrt(p)], [0.0, 0.0]])
    return k0, k1
449152
def multifuncrs(X, funs, eps=1E-6,
                nswp=10,
                kickrank=5,
                y0=None,
                rmax=999999,  # TODO:infinity \
                kicktype='amr-two',        \
                pcatype='svd',             \
                trunctype='fro',           \
                d2=1,                      \
                do_qr=False,               \
                verb=1):
    

    dtype = np.float64
    if len([x for x in X if x.is_complex]) > 0:
        dtype = np.complex128

    y = y0
    wasrand = False

    nx = len(X)
    d = X[0].d
    n = X[0].n
    rx = np.transpose(np.array([ttx.r for ttx in X]))
    #crx = [tt.tensor.to_list(ttx) for x in X]
    #crx = zip(*crx)
    crx = np.transpose(np.array([tt.tensor.to_list(ttx)
                                 for ttx in X], dtype=np.object))
    crx = np.empty((nx, d), dtype=np.object)
    i = 0
    for ttx in X:
        v = tt.tensor.to_list(ttx)
        j = 0
        for w in v:
            crx[i, j] = w
            j = j + 1
        i = i + 1
    crx = crx.T
    if y is None:
        ry = d2 * np.ones((d + 1,), dtype=np.int32)
        ry[0] = 1
        y = tt.rand(n, d, ry)
        wasrand = True

    ry = y.r
    cry = tt.tensor.to_list(y)

    Ry = np.zeros((d + 1, ), dtype=np.object)
    Ry[0] = np.array([[1.0]], dtype=dtype)
    Ry[d] = np.array([[1.0]], dtype=dtype)
    Rx = np.zeros((d + 1, nx), dtype=np.object)
    Rx[0, :] = np.ones(nx, dtype=dtype)
    Rx[d, :] = np.ones(nx, dtype=dtype)

    block_order = [+d, -d]

    # orth
    for i in range(0, d - 1):
        cr = cry[i]
        cr = reshape(cr, (ry[i] * n[i], ry[i + 1]))
        cr, rv = np.linalg.qr(cr)
        cr2 = cry[i + 1]
        cr2 = reshape(cr2, (ry[i + 1], n[i + 1] * ry[i + 2]))
        cr2 = np.dot(rv, cr2)  # matrix multiplication
        ry[i + 1] = cr.shape[1]
        cr = reshape(cr, (ry[i], n[i], ry[i + 1]))
        cry[i + 1] = reshape(cr2, (ry[i + 1], n[i + 1], ry[i + 2]))
        cry[i] = cr

        Ry[i + 1] = np.dot(Ry[i], reshape(cr, (ry[i], n[i] * ry[i + 1])))
        Ry[i + 1] = reshape(Ry[i + 1], (ry[i] * n[i], ry[i + 1]))
        curind = []
        if wasrand:
            # EVERY DAY I'M SHUFFLIN'
            curind = np.random.permutation(n[i] * ry[i])[:ry[i + 1]]
        else:
            curind = maxvol(Ry[i + 1])
        Ry[i + 1] = Ry[i + 1][curind, :]
        for j in range(0, nx):
            try:
                Rx[i + 1, j] = reshape(crx[i, j],
                                       (rx[i, j], n[i] * rx[i + 1, j]))
            except:
                pass
            Rx[i + 1, j] = np.dot(Rx[i, j], Rx[i + 1, j])
            Rx[i + 1, j] = reshape(Rx[i + 1, j], (ry[i] * n[i], rx[i + 1, j]))
            Rx[i + 1, j] = Rx[i + 1, j][curind, :]

    d2 = ry[d]
    ry[d] = 1
    cry[d - 1] = np.transpose(cry[d - 1], [2, 0, 1])  # permute

    last_sweep = False
    swp = 1

    dy = np.zeros((d, ))
    max_dy = 0

    cur_order = copy.copy(block_order)
    order_index = 1
    i = d - 1
    # can't use 'dir' identifier in python
    dirn = int(math.copysign(1, cur_order[order_index]))

    # DMRG sweeps
    while swp <= nswp or dirn > 0:

        oldy = reshape(cry[i], (d2 * ry[i] * n[i] * ry[i + 1],))

        if not last_sweep:
            # compute the X superblocks
            curbl = np.zeros((ry[i] * n[i] * ry[i + 1], nx), dtype=dtype)
            for j in range(0, nx):
                cr = reshape(crx[i, j], (rx[i, j], n[i] * rx[i + 1, j]))
                cr = np.dot(Rx[i, j], cr)
                cr = reshape(cr, (ry[i] * n[i], rx[i + 1, j]))
                cr = np.dot(cr, Rx[i + 1, j])
                curbl[:, j] = cr.flatten('F')
            # call the function
            newy = funs(curbl)
            # multiply with inverted Ry
            newy = reshape(newy, (ry[i], n[i] * ry[i + 1] * d2))
            newy = np.linalg.solve(Ry[i], newy)  # y = R \ y
            newy = reshape(newy, (ry[i] * n[i] * ry[i + 1], d2))
            newy = reshape(np.transpose(newy), (d2 * ry[i] * n[i], ry[i + 1]))
            newy = np.transpose(np.linalg.solve(
                np.transpose(Ry[i + 1]), np.transpose(newy)))  # y=y/R
            newy = reshape(newy, (d2 * ry[i] * n[i] * ry[i + 1],))
        else:
            newy = oldy

        dy[i] = np.linalg.norm(newy - oldy) / np.linalg.norm(newy)
        max_dy = max(max_dy, dy[i])

        # truncation
        if dirn > 0:  # left-to-right
            newy = reshape(newy, (d2, ry[i] * n[i] * ry[i + 1]))
            newy = reshape(np.transpose(newy), (ry[i] * n[i], ry[i + 1] * d2))
        else:
            newy = reshape(newy, (d2 * ry[i], n[i] * ry[i + 1]))

        r = 0  # defines a variable in global scope

        if kickrank >= 0:
            u, s, v = np.linalg.svd(newy, full_matrices=False)
            v = np.conj(np.transpose(v))
            if trunctype == "fro" or last_sweep:
                r = my_chop2(s, eps / math.sqrt(d) * np.linalg.norm(s))
            else:
                # truncate taking into account the (r+1) overhead in the cross
                # (T.S.: what?)
                cums = abs(s * np.arange(2, len(s) + 2)) ** 2
                cums = np.cumsum(cums[::-1])[::-1]
                cums = cums / cums[0]
                ff = [i for i in range(len(cums)) if cums[i] < eps ** 2 / d]
                if len(ff) == 0:
                    r = len(s)
                else:
                    r = np.amin(ff)
            r = min(r, rmax, len(s))
        else:
            if dirn > 0:
                u, v = np.linalg.qr(newy)
                v = np.conj(np.transpose(v))
                r = u.shape[1]
                s = np.ones((r, ))
            else:
                v, u = np.linalg.qr(np.transpose(newy))
                v = np.conj(v)
                u = np.transpose(u)
                r = u.shape[1]
                s = np.ones((r, ))

        if verb > 1:
            print('=multifuncrs=   block %d{%d}, dy: %3.3e, r: %d' % (i, dirn, dy[i], r))

        # kicks and interfaces
        if dirn > 0 and i < d - 1:
            u = u[:, :r]
            v = np.dot(v[:, :r], np.diag(s[:r]))

            # kick
            radd = 0
            rv = 1
            if not last_sweep and kickrank > 0:
                uk = None
                if kicktype == 'amr-two':
                    # AMR(two)-like kick.

                    # compute the X superblocks
                    ind2 = np.unique(np.random.randint(
                        0, ry[i + 2] * n[i + 1], ry[i + 1]))
                    #ind2 = np.unique(np.floor(np.random.rand(ry[i + 1]) * (ry[i + 2] * n[i + 1])))
                    rkick = len(ind2)
                    curbl = np.zeros((ry[i] * n[i] * rkick, nx), dtype=dtype)
                    for j in range(nx):
                        cr1 = reshape(
                            crx[i, j], (rx[i, j], n[i] * rx[i + 1, j]))
                        cr1 = np.dot(Rx[i, j], cr1)
                        cr1 = reshape(cr1, (ry[i] * n[i], rx[i + 1, j]))
                        cr2 = reshape(
                            crx[i + 1, j], (rx[i + 1, j] * n[i + 1], rx[i + 2, j]))
                        cr2 = np.dot(cr2, Rx[i + 2, j])
                        cr2 = reshape(
                            cr2, (rx[i + 1, j], n[i + 1] * ry[i + 2]))
                        cr2 = cr2[:, ind2]
                        curbl[:, j] = reshape(
                            np.dot(cr1, cr2), (ry[i] * n[i] * rkick,))
                    # call the function
                    uk = funs(curbl)
                    uk = reshape(uk, (ry[i], n[i] * rkick * d2))
                    uk = np.linalg.solve(Ry[i], uk)
                    uk = reshape(uk, (ry[i] * n[i], rkick * d2))
                    if pcatype == 'svd':
                        uk, sk, vk = np.linalg.svd(uk, full_matrices=False)
                        vk = np.conj(np.transpose(vk))
                        uk = uk[:, :min(kickrank, uk.shape[1])]
                    else:
                        # uk = uchol(np.transpose(uk), kickrank + 1) # TODO
                        uk = uk[:, :max(uk.shape[1] - kickrank + 1, 1):-1]
                else:
                    uk = np.random.rand(ry[i] * n[i], kickrank)
                u, rv = np.linalg.qr(np.concatenate((u, uk), axis=1))
                radd = uk.shape[1]
            v = np.concatenate(
                (v, np.zeros((ry[i + 1] * d2, radd), dtype=dtype)), axis=1)
            v = np.dot(rv, np.conj(np.transpose(v)))
            r = u.shape[1]

            cr2 = cry[i + 1]
            cr2 = reshape(cr2, (ry[i + 1], n[i + 1] * ry[i + 2]))
            v = reshape(v, (r * ry[i + 1], d2))
            v = reshape(np.transpose(v), (d2 * r, ry[i + 1]))
            v = np.dot(v, cr2)

            ry[i + 1] = r

            u = reshape(u, (ry[i], n[i], r))
            v = reshape(v, (d2, r, n[i + 1], ry[i + 2]))

            cry[i] = u
            cry[i + 1] = v

            Ry[i + 1] = np.dot(Ry[i], reshape(u, (ry[i], n[i] * ry[i + 1])))
            Ry[i + 1] = reshape(Ry[i + 1], (ry[i] * n[i], ry[i + 1]))
            curind = maxvol(Ry[i + 1])
            Ry[i + 1] = Ry[i + 1][curind, :]
            for j in range(nx):
                Rx[i + 1, j] = reshape(crx[i, j],
                                       (rx[i, j], n[i] * rx[i + 1, j]))
                Rx[i + 1, j] = np.dot(Rx[i, j], Rx[i + 1, j])
                Rx[i + 1, j] = reshape(Rx[i + 1, j],
                                       (ry[i] * n[i], rx[i + 1, j]))
                Rx[i + 1, j] = Rx[i + 1, j][curind, :]
        elif dirn < 0 and i > 0:
            u = np.dot(u[:, :r], np.diag(s[:r]))
            v = np.conj(v[:, :r])

            radd = 0
            rv = 1
            if not last_sweep and kickrank > 0:
                if kicktype == 'amr-two':
                    # compute the X superblocks
                    ind2 = np.unique(np.random.randint(
                        0, ry[i - 1] * n[i - 1], ry[i]))
                    rkick = len(ind2)
                    curbl = np.zeros(
                        (rkick * n[i] * ry[i + 1], nx), dtype=dtype)
                    for j in range(nx):
                        cr1 = reshape(
                            crx[i, j], (rx[i, j] * n[i], rx[i + 1, j]))
                        cr1 = np.dot(cr1, Rx[i + 1, j])
                        cr1 = reshape(cr1, (rx[i, j], n[i] * ry[i + 1]))
                        cr2 = reshape(
                            crx[i - 1, j], (rx[i - 1, j], n[i - 1] * rx[i, j]))
                        cr2 = np.dot(Rx[i - 1, j], cr2)
                        cr2 = reshape(cr2, (ry[i - 1] * n[i - 1], rx[i, j]))
                        cr2 = cr2[ind2, :]
                        curbl[:, j] = reshape(
                            np.dot(cr2, cr1), (rkick * n[i] * ry[i + 1],))
                    # calling the function
                    uk = funs(curbl)
                    uk = reshape(uk, (rkick * n[i] * ry[i + 1], d2))
                    uk = reshape(np.transpose(
                        uk), (d2 * rkick * n[i], ry[i + 1]))
                    uk = np.transpose(np.linalg.solve(
                        np.transpose(Ry[i + 1]), np.transpose(uk)))
                    uk = reshape(uk, (d2 * rkick, n[i] * ry[i + 1]))
                    if pcatype == 'svd':
                        vk, sk, uk = np.linalg.svd(uk, full_matrices=False)
                        uk = np.conj(np.transpose(uk))
                        # TODO: refactor
                        uk = uk[:, :min(kickrank, uk.shape[1])]
                    else:
                        # uk = uchol(uk, kickrank + 1) # TODO
                        uk = uk[:, :max(uk.shape[1] - kickrank + 1, 1):-1]
                else:
                    uk = np.random.rand(n[i] * ry[i + 1], kickrank)
                v, rv = np.linalg.qr(np.concatenate((v, uk), axis=1))
                radd = uk.shape[1]
            u = np.concatenate(
                (u, np.zeros((d2 * ry[i], radd), dtype=dtype)), axis=1)
            u = np.dot(u, np.transpose(rv))
            r = v.shape[1]
            cr2 = cry[i - 1]
            cr2 = reshape(cr2, (ry[i - 1] * n[i - 1], ry[i]))
            u = reshape(u, (d2, ry[i] * r))
            u = reshape(np.transpose(u), (ry[i], r * d2))
            u = np.dot(cr2, u)

            u = reshape(u, (ry[i - 1] * n[i - 1] * r, d2))
            u = reshape(np.transpose(u), (d2, ry[i - 1], n[i - 1], r))
            v = reshape(np.transpose(v), (r, n[i], ry[i + 1]))

            ry[i] = r
            cry[i - 1] = u
            cry[i] = v

            Ry[i] = np.dot(reshape(v, (ry[i] * n[i], ry[i + 1])), Ry[i + 1])
            Ry[i] = reshape(Ry[i], (ry[i], n[i] * ry[i + 1]))
            curind = maxvol(np.transpose(Ry[i]))
            Ry[i] = Ry[i][:, curind]
            for j in range(nx):
                Rx[i, j] = reshape(crx[i, j], (rx[i, j] * n[i], rx[i + 1, j]))
                Rx[i, j] = np.dot(Rx[i, j], Rx[i + 1, j])
                Rx[i, j] = reshape(Rx[i, j], (rx[i, j], n[i] * ry[i + 1]))
                Rx[i, j] = Rx[i, j][:, curind]
        elif dirn > 0 and i == d - 1:
            newy = np.dot(np.dot(u[:, :r], np.diag(s[:r])),
                          np.conj(np.transpose(v[:, :r])))
            newy = reshape(newy, (ry[i] * n[i] * ry[i + 1], d2))
            cry[i] = reshape(np.transpose(newy), (d2, ry[i], n[i], ry[i + 1]))
        elif dirn < 0 and i == 0:
            newy = np.dot(np.dot(u[:, :r], np.diag(s[:r])),
                          np.conj(np.transpose(v[:, :r])))
            newy = reshape(newy, (d2, ry[i], n[i], ry[i + 1]))
            cry[i] = newy

        i = i + dirn
        cur_order[order_index] = cur_order[order_index] - dirn
        if cur_order[order_index] == 0:
            order_index = order_index + 1
            if verb > 0:
                print('=multifuncrs= sweep %d{%d}, max_dy: %3.3e, erank: %g' % (swp, order_index, max_dy,
                                                                                math.sqrt(np.dot(ry[:d], n * ry[1:]) / np.sum(n))))

            if last_sweep:
                break
            if max_dy < eps and dirn < 0:
                last_sweep = True
                kickrank = 0

            if order_index >= len(cur_order):
                cur_order = copy.copy(block_order)
                order_index = 0
                if last_sweep:
                    cur_order = [d - 1]

                max_dy = 0
                swp = swp + 1

            dirn = int(math.copysign(1, cur_order[order_index]))
            i = i + dirn

    cry[d - 1] = np.transpose(cry[d - 1][:, :, :, 0], [1, 2, 0])
    y = tt.tensor.from_list(cry)
    return y
269776
def _construct_functions(self, coefs, **kwargs):
        """"""
        return [self.basis_functions.functions_factory(coef, **kwargs) for coef in coefs]
470228
def min_tens(tens, rmax=10, nswp=10, verb=True, smooth_fun=None):
    """"""
    if smooth_fun is None:
        smooth_fun = lambda p, lam: (math.pi / 2 - np.arctan(p - lam))
    d = tens.d
    Rx = [[]] * (d + 1)  # Python list for the interfaces
    Rx[0] = np.ones((1, 1))
    Rx[d] = np.ones((1, 1))
    Jy = [np.empty(0, dtype=np.int)] * (d + 1)
    ry = rmax * np.ones(d + 1, dtype=np.int)
    ry[0] = 1
    ry[d] = 1
    n = tens.n
    elements_seen = 0
    phi_left = [np.empty(0)] * (d + 1)
    phi_left[0] = np.array([1])
    phi_right = [np.empty(0)] * (d + 1)
    phi_right[d] = np.array([1])
    cores = tt.tensor.to_list(tens)

    # Fill initial multiindex J randomly.
    grid = [np.reshape(range(n[i]), (n[i], 1)) for i in xrange(d)]
    for i in xrange(d - 1):
        ry[i + 1] = min(ry[i + 1], n[i] * ry[i])
        ind = sorted(np.random.permutation(ry[i] * n[i])[0:ry[i + 1]])
        w1 = mkron(np.ones((n[i], 1), dtype=np.int), Jy[i])
        w2 = mkron(grid[i], np.ones((ry[i], 1), dtype=np.int))
        Jy[i + 1] = np.hstack((w1, w2))
        Jy[i + 1] = reshape(Jy[i + 1], (ry[i] * n[i], -1))
        Jy[i + 1] = Jy[i + 1][ind, :]
        phi_left[i + 1] = np.tensordot(phi_left[i], cores[i], 1)
        phi_left[i + 1] = reshape(phi_left[i + 1], (ry[i] * n[i], -1))
        phi_left[i + 1] = phi_left[i + 1][ind, :]
    swp = 0
    dirn = -1
    i = d - 1
    lm = float('Inf')
    while swp < nswp:
        # Right-to-left sweep
        # The idea: compute the current core; compute the function of it;
        # Shift locally or globally? Local shift would be the first try
        # Compute the current core

        if np.size(Jy[i]) == 0:
            w1 = np.zeros((ry[i] * n[i] * ry[i + 1], 0), dtype=np.int)
        else:
            w1 = mkron(np.ones((n[i] * ry[i + 1], 1), dtype=np.int), Jy[i])
        w2 = mkron(mkron(np.ones((ry[i + 1], 1), dtype=np.int),
                         grid[i]), np.ones((ry[i], 1), dtype=np.int))
        if np.size(Jy[i + 1]) == 0:
            w3 = np.zeros((ry[i] * n[i] * ry[i + 1], 0), dtype=np.int)
        else:
            w3 = mkron(Jy[i + 1], np.ones((ry[i] * n[i], 1), dtype=np.int))
        J = np.hstack((w1, w2, w3))

        phi_right[i] = np.tensordot(cores[i], phi_right[i + 1], 1)
        phi_right[i] = reshape(phi_right[i], (-1, n[i] * ry[i + 1]))

        cry = np.tensordot(
            phi_left[i], np.tensordot(
                cores[i], phi_right[
                    i + 1], 1), 1)
        elements_seen += cry.size
        cry = reshape(cry, (ry[i], n[i], ry[i + 1]))
        min_cur = np.min(cry.flatten("F"))
        ind_cur = np.argmin(cry.flatten("F"))
        if lm > min_cur:
            lm = min_cur
            x_full = J[ind_cur, :]
            val = tens[x_full]
            if verb:
                print('New record:', val, 'Point:', x_full, 'elements seen:', elements_seen)
        cry = smooth_fun(cry, lm)
        if dirn < 0 and i > 0:
            cry = reshape(cry, (ry[i], n[i] * ry[i + 1]))
            cry = cry.T
            #q, r = np.linalg.qr(cry)
            u, s, v = mysvd(cry, full_matrices=False)
            ry[i] = min(ry[i], rmax)
            q = u[:, :ry[i]]
            ind = rect_maxvol(q)[0]  # maxvol(q)
            ry[i] = ind.size
            w1 = mkron(np.ones((ry[i + 1], 1), dtype=np.int), grid[i])
            if np.size(Jy[i + 1]) == 0:
                w2 = np.zeros((n[i] * ry[i + 1], 0), dtype=np.int)
            else:
                w2 = mkron(Jy[i + 1], np.ones((n[i], 1), dtype=np.int))
            Jy[i] = np.hstack((w1, w2))
            Jy[i] = reshape(Jy[i], (n[i] * ry[i + 1], -1))
            Jy[i] = Jy[i][ind, :]
            phi_right[i] = np.tensordot(cores[i], phi_right[i + 1], 1)
            phi_right[i] = reshape(phi_right[i], (-1, n[i] * ry[i + 1]))
            phi_right[i] = phi_right[i][:, ind]

        if dirn > 0 and i < d - 1:
            cry = reshape(cry, (ry[i] * n[i], ry[i + 1]))
            q, r = np.linalg.qr(cry)
            #ind = maxvol(q)
            ind = rect_maxvol(q)[0]
            ry[i + 1] = ind.size
            phi_left[i + 1] = np.tensordot(phi_left[i], cores[i], 1)
            phi_left[i + 1] = reshape(phi_left[i + 1], (ry[i] * n[i], -1))
            phi_left[i + 1] = phi_left[i + 1][ind, :]
            w1 = mkron(np.ones((n[i], 1), dtype=np.int), Jy[i])
            w2 = mkron(grid[i], np.ones((ry[i], 1), dtype=np.int))
            Jy[i + 1] = np.hstack((w1, w2))
            Jy[i + 1] = reshape(Jy[i + 1], (ry[i] * n[i], -1))
            Jy[i + 1] = Jy[i + 1][ind, :]

        i += dirn
        if i == d or i == -1:
            dirn = -dirn
            i += dirn
            swp = swp + 1
    return val, x_full
1123496
def _create_gate_variables(self, input_shape, dtype):
    """"""
    if len(input_shape) != 2:
      raise ValueError(
          "Rank of shape must be {} not: {}".format(2, len(input_shape)))
    input_size = input_shape.dims[1].value

    b_shape = [4 * self._hidden_size]

    equiv_input_size = self._hidden_size + input_size
    initializer = basic.create_linear_initializer(equiv_input_size)

    if self._use_batch_norm_h or self._use_batch_norm_x:
      self._w_h = tf.get_variable(
          self.W_GATES + "_H",
          shape=[self._hidden_size, 4 * self._hidden_size],
          dtype=dtype,
          initializer=self._initializers.get(self.W_GATES, initializer),
          partitioner=self._partitioners.get(self.W_GATES),
          regularizer=self._regularizers.get(self.W_GATES))
      self._w_x = tf.get_variable(
          self.W_GATES + "_X",
          shape=[input_size, 4 * self._hidden_size],
          dtype=dtype,
          initializer=self._initializers.get(self.W_GATES, initializer),
          partitioner=self._partitioners.get(self.W_GATES),
          regularizer=self._regularizers.get(self.W_GATES))
    else:
      self._w_xh = tf.get_variable(
          self.W_GATES,
          shape=[self._hidden_size + input_size, 4 * self._hidden_size],
          dtype=dtype,
          initializer=self._initializers.get(self.W_GATES, initializer),
          partitioner=self._partitioners.get(self.W_GATES),
          regularizer=self._regularizers.get(self.W_GATES))
    self._b = tf.get_variable(
        self.B_GATES,
        shape=b_shape,
        dtype=dtype,
        initializer=self._initializers.get(self.B_GATES, initializer),
        partitioner=self._partitioners.get(self.B_GATES),
        regularizer=self._regularizers.get(self.B_GATES))
410863
def glsfit(psr,renormalize=True):
    """"""
    
    mask = psr.deleted == 0
    res, err = psr.residuals(removemean=False)[mask], psr.toaerrs[mask]
    M = psr.designmatrix(updatebats=False,fixunits=True,fixsigns=True)[mask,:]
        
    C = numpy.diag((err * 1e-6)**2)
    
    if renormalize:
        norm = numpy.sqrt(numpy.sum(M**2,axis=0))
        M /= norm
    else:
        norm = np.ones_like(M[0,:])
    
    mtcm = numpy.dot(M.T,numpy.dot(numpy.linalg.inv(C),M))
    mtcy = numpy.dot(M.T,numpy.dot(numpy.linalg.inv(C),res))
    
    xvar = numpy.linalg.inv(mtcm)
        
    c = scipy.linalg.cho_factor(mtcm)
    xhat = scipy.linalg.cho_solve(c,mtcy)

    sol = psr.vals()
    psr.vals(sol + xhat[1:] / norm[1:])
    psr.errs(numpy.sqrt(numpy.diag(xvar)[1:]) / norm[1:])
    
    return chisq(psr)
1183342
def make_nn_funs(layer_sizes, L2_reg, noise_variance, nonlinearity=np.tanh):
    """"""
    shapes = list(zip(layer_sizes[:-1], layer_sizes[1:]))
    num_weights = sum((m+1)*n for m, n in shapes)

    def unpack_layers(weights):
        num_weight_sets = len(weights)
        for m, n in shapes:
            yield weights[:, :m*n]     .reshape((num_weight_sets, m, n)),\
                  weights[:, m*n:m*n+n].reshape((num_weight_sets, 1, n))
            weights = weights[:, (m+1)*n:]

    def predictions(weights, inputs):
        """weights is shape (num_weight_samples x num_weights)
           inputs  is shape (num_datapoints x D)"""
        inputs = np.expand_dims(inputs, 0)
        for W, b in unpack_layers(weights):
            outputs = np.einsum('mnd,mdo->mno', inputs, W) + b
            inputs = nonlinearity(outputs)
        return outputs

    def logprob(weights, inputs, targets):
        log_prior = -L2_reg * np.sum(weights**2, axis=1)
        preds = predictions(weights, inputs)
        log_lik = -np.sum((preds - targets)**2, axis=1)[:, 0] / noise_variance
        return log_prior + log_lik

    return num_weights, predictions, logprob
775106
def make_fc(dim_in, hidden_dim, use_gn=False):
    '''
        
    '''
    if use_gn:
        fc = nn.Linear(dim_in, hidden_dim, bias=False)
        nn.init.kaiming_uniform_(fc.weight, a=1)
        return nn.Sequential(fc, group_norm(hidden_dim))
    fc = nn.Linear(dim_in, hidden_dim)
    nn.init.kaiming_uniform_(fc.weight, a=1)
    nn.init.constant_(fc.bias, 0)
    return fc
1216883
def dense_dropconnect(inputs,
                      output_size,
                      dropconnect_dropout=0.0,
                      name="dense_dropconnect",
                      **kwargs):
  """"""

  if dropconnect_dropout != 0.0:
    tf.logging.info("Applying dropconnect as the kernel regularization.")
    kwargs["kernel_regularizer"] = functools.partial(
        tf.nn.dropout, keep_prob=1.0 - dropconnect_dropout)

  return dense(inputs, output_size, use_bias=True, name=name, **kwargs)
613024
def _fill_diagonals(m, diag_indices):
    """"""
    assert m.ndim == 3, "M must have 3 dimensions"
    assert m.shape[1] == m.shape[2], "M must contain square matrices"
    for r in range(m.shape[0]):
        scipy.fill_diagonal(m[r], 0)
        m[r][diag_indices] -= scipy.sum(m[r], axis=1)