Python numpy 模块，meshgrid() 实例源码

def get_local_wavenumbermesh(self, scaled=True, broadcast=False,
                                 eliminate_highest_freq=False):
        kx = fftfreq(self.N[0], 1./self.N[0])
        ky = rfftfreq(self.N[1], 1./self.N[1])
        if eliminate_highest_freq:
            for i, k in enumerate((kx, ky)):
                if self.N[i] % 2 == 0:
                    k[self.N[i]//2] = 0

        Ks = np.meshgrid(kx, ky[self.rank*self.Np[1]//2:(self.rank*self.Np[1]//2+self.Npf)], indexing='ij', sparse=True)
        if scaled is True:
            Lp = 2*np.pi/self.L
            Ks[0] *= Lp[0]
            Ks[1] *= Lp[1]
        K = Ks
        if broadcast is True:
            K = [np.broadcast_to(k, self.complex_shape()) for k in Ks]
        return K

def normvectorfield(xs,ys,fs,**kw):
    """
    plot normalized vector field

    kwargs
    ======

    - length is a desired length of the lines (default: 1)
    - the rest of kwards are passed to plot
    """
    length = kw.pop('length') if 'length' in kw else 1
    x, y = np.meshgrid(xs, ys)
    # calculate vector field
    vx,vy = fs(x,y)
    # plot vecor field
    norm = length /np.sqrt(vx**2+vy**2)
    plt.quiver(x, y, vx * norm, vy * norm, angles='xy',**kw)

def plot_interpolation(orderx,ordery):
    s = PseudoSpectralDiscretization2D(orderx,XMIN,XMAX,
                                ordery,YMIN,YMAX)
    Xc,Yc = s.get_x2d()
    x = np.linspace(XMIN,XMAX,100)
    y = np.linspace(YMIN,YMAX,100)
    Xf,Yf = np.meshgrid(x,y,indexing='ij')
    f_coarse = f(Xc,Yc)
    f_interpolator = s.to_continuum(f_coarse)
    f_num = f_interpolator(Xf,Yf)
    plt.pcolor(Xf,Yf,f_num)
    cb = plt.colorbar()
    cb.set_label('interpolated function',fontsize=16)
    plt.xlabel('x')
    plt.ylabel('y')
    for postfix in ['.png','.pdf']:
        name = 'orthopoly_interpolated_function'+postfix
        if USE_FIGS_DIR:
            name = 'figs/' + name
        plt.savefig(name,
                    bbox_inches='tight')
    plt.clf()

def __init__(self,
                 orderx,xmin,xmax,
                 ordery,ymin,ymax):
        "Constructor. Needs order and domain in x and y"
        self.orderx,self.ordery = orderx,ordery
        self.stencils = [PseudoSpectralDiscretization1D(orderx,xmin,xmax),
                         PseudoSpectralDiscretization1D(ordery,ymin,ymax)]
        self.stencil_x,self.stencil_y = self.stencils
        self.quads = [s.quads for s in self.stencils]
        self.colocs = [s.colocation_points for s in self.stencils]
        self.x,self.y = self.colocs
        self.colocs2d = np.meshgrid(*self.colocs,indexing='ij')
        self.X,self.Y = self.colocs2d
        self.weights = [s.weights for s in self.stencils]
        self.weights_x,self.weights_y = self.weights
        self.weights2D = np.tensordot(*self.weights,axes=0)

def vectorfield(xs,ys,fs,**kw):
    """
    plot vector field (no normalization!)

    args
    ====
    fs is a function that returns tuple (vx,vy)

    kwargs
    ======

    - length is a desired length of the lines (default: 1)
    - the rest of kwards are passed to plot
    """
    length= kw.pop('length') if 'length' in kw else 1
    x, y=np.meshgrid(xs, ys)
    # calculate vector field
    vx,vy=fs(x,y)
    # plot vecor field
    norm = length 
    plt.quiver(x, y, vx * norm, vy * norm, angles='xy',**kw)

def create_reference_image(size, x0=10., y0=-3., sigma_x=50., sigma_y=30., dtype='float64',
                           reverse_xaxis=False, correct_axes=True, sizey=None, **kwargs):
    """
    Creates a reference image: a gaussian brightness with elliptical
    """
    inc_cos = np.cos(0./180.*np.pi)

    delta_x = 1.
    x = (np.linspace(0., size - 1, size) - size / 2.) * delta_x

    if sizey:
        y = (np.linspace(0., sizey-1, sizey) - sizey/2.) * delta_x
    else:
        y = x.copy()

    if reverse_xaxis:
        xx, yy = np.meshgrid(-x, y/inc_cos)
    elif correct_axes:
        xx, yy = np.meshgrid(-x, -y/inc_cos)
    else:
        xx, yy = np.meshgrid(x, y/inc_cos)

    image = np.exp(-(xx-x0)**2./sigma_x - (yy-y0)**2./sigma_y)

    return image.astype(dtype)

def generate_anchors_pre(height, width, feat_stride, anchor_scales=(8,16,32), anchor_ratios=(0.5,1,2)):
  """ A wrapper function to generate anchors given different scales
    Also return the number of anchors in variable 'length'
  """
  anchors = generate_anchors(ratios=np.array(anchor_ratios), scales=np.array(anchor_scales))
  A = anchors.shape[0]
  shift_x = np.arange(0, width) * feat_stride
  shift_y = np.arange(0, height) * feat_stride
  shift_x, shift_y = np.meshgrid(shift_x, shift_y)
  shifts = np.vstack((shift_x.ravel(), shift_y.ravel(), shift_x.ravel(), shift_y.ravel())).transpose()
  K = shifts.shape[0]
  # width changes faster, so here it is H, W, C
  anchors = anchors.reshape((1, A, 4)) + shifts.reshape((1, K, 4)).transpose((1, 0, 2))
  anchors = anchors.reshape((K * A, 4)).astype(np.float32, copy=False)
  length = np.int32(anchors.shape[0])

  return anchors, length

def plot2d_simplex(simplex, ind):
    fig_dir = "./"
    plt.cla()
    n = 1000
    x1 = np.linspace(-256, 1024, n)
    x2 = np.linspace(-256, 1024, n)
    X, Y = np.meshgrid(x1, x2)
    Z = np.sqrt(X ** 2 + Y ** 2)
    plt.contour(X, Y, Z, levels=list(np.arange(0, 1200, 10)))
    plt.gca().set_aspect("equal")
    plt.xlim((-256, 768))
    plt.ylim((-256, 768))

    plt.plot([simplex[0].x[0], simplex[1].x[0]],
             [simplex[0].x[1], simplex[1].x[1]], color="#000000")
    plt.plot([simplex[1].x[0], simplex[2].x[0]],
             [simplex[1].x[1], simplex[2].x[1]], color="#000000")
    plt.plot([simplex[2].x[0], simplex[0].x[0]],
             [simplex[2].x[1], simplex[0].x[1]], color="#000000")
    plt.savefig(os.path.join(fig_dir, "{:03d}.png".format(ind)))

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]
#  ????

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]
#  ????

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]
#  ????

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]
#  ????

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]
#  ????

def draw_hill(x,y):
    a = np.linspace(-20,20,100)
    print(a)
    b = np.linspace(-20,20,100)
    x = np.array(x)
    y = np.array(y)

    allSSE = np.zeros(shape=(len(a),len(b)))
    for ai in range(0,len(a)):
        for bi in range(0,len(b)):
            a0 = a[ai]
            b0 = b[bi]
            SSE = calc_loss(a=a0,b=b0,x=x,y=y)
            allSSE[ai][bi] = SSE

    a,b = np.meshgrid(a, b)

    return [a,b,allSSE]
#  ????

def plot_grid2D(lons, lats, tec_grid2D, datetime, title_label = ''):

    LATS, LONS = np.meshgrid(lats, lons)

    m = Basemap(llcrnrlon=-180,
                llcrnrlat=-55,
                urcrnrlon=180,
                urcrnrlat=75,
                projection='merc',
                area_thresh=1000,
                resolution='i')

    m.drawstates()
    m.drawcountries()
    m.drawcoastlines()

    parallels = np.arange(-90,90,20)
    m.drawparallels(parallels,labels=[True,False,False,True])
    meridians = np.arange(0,360,40)
    m.drawmeridians(meridians,labels=[True,False,False,True])

    m.scatter(LONS, LATS, c=tec_grid2D, latlon = True, linewidths=0, s=5)
    m.colorbar()

    plt.title('%s\n%s' % (title_label, datetime.isoformat(' ')))

def _meshgrid(self, height, width):
    with tf.variable_scope('_meshgrid'):
      # This should be equivalent to:
      #  x_t, y_t = np.meshgrid(np.linspace(-1, 1, width),
      #                         np.linspace(-1, 1, height))
      #  ones = np.ones(np.prod(x_t.shape))
      #  grid = np.vstack([x_t.flatten(), y_t.flatten(), ones])
      x_t = tf.matmul(tf.ones(shape=tf.pack([height, 1])),
                        tf.transpose(tf.expand_dims(tf.linspace(-1.0, 1.0, width), 1), [1, 0]))
      y_t = tf.matmul(tf.expand_dims(tf.linspace(-1.0, 1.0, height), 1),
                        tf.ones(shape=tf.pack([1, width])))

      x_t_flat = tf.reshape(x_t, (1, -1))
      y_t_flat = tf.reshape(y_t, (1, -1))

      ones = tf.ones_like(x_t_flat)
      grid = tf.concat(0, [x_t_flat, y_t_flat, ones])
      return grid

def make_3d_mask(img_shape, center, radius, shape='sphere'):
    mask = np.zeros(img_shape)
    radius = np.rint(radius)
    center = np.rint(center)
    sz = np.arange(int(max(center[0] - radius, 0)), int(max(min(center[0] + radius + 1, img_shape[0]), 0)))
    sy = np.arange(int(max(center[1] - radius, 0)), int(max(min(center[1] + radius + 1, img_shape[1]), 0)))
    sx = np.arange(int(max(center[2] - radius, 0)), int(max(min(center[2] + radius + 1, img_shape[2]), 0)))
    sz, sy, sx = np.meshgrid(sz, sy, sx)
    if shape == 'cube':
        mask[sz, sy, sx] = 1.
    elif shape == 'sphere':
        distance2 = ((center[0] - sz) ** 2
                     + (center[1] - sy) ** 2
                     + (center[2] - sx) ** 2)
        distance_matrix = np.ones_like(mask) * np.inf
        distance_matrix[sz, sy, sx] = distance2
        mask[(distance_matrix <= radius ** 2)] = 1
    elif shape == 'gauss':
        z, y, x = np.ogrid[:mask.shape[0], :mask.shape[1], :mask.shape[2]]
        distance = ((z - center[0]) ** 2 + (y - center[1]) ** 2 + (x - center[2]) ** 2)
        mask = np.exp(- 1. * distance / (2 * radius ** 2))
        mask[(distance > 3 * radius ** 2)] = 0
    return mask

def run_numerical_test(function, n, k, results):
    n_grid, k_grid = np.meshgrid(n, k)

    if results.shape != n_grid.shape:
        raise ValueError('results passed do not match the shape of n/k')

    # Test both as scalars
    for i, cur_n in enumerate(n):
        for j, cur_k in enumerate(k):
            np.testing.assert_allclose(function(cur_n, cur_k), results[i, j], err_msg=('n: ' + str(cur_n) + ', k: ' + str(cur_k)))

    # Test first as scalar
    for i, cur_n in enumerate(n):
        np.testing.assert_allclose(function(cur_n, k_grid[:, i]), results[i, :], err_msg=('n: ' + str(cur_n) + ', k: ' + str(k_grid[:, i])))

    # Test two as scalar
    for j, cur_k in enumerate(k):
        np.testing.assert_allclose(function(n_grid[j, :], cur_k), results[:, j], err_msg=('n: ' + str(n_grid[j, :]) + ', k: ' + str(cur_k)))

    # Test matrix
    np.testing.assert_allclose(function(n_grid, k_grid), results.T, err_msg=('Matrix computation failed!'))

def genSphCoords():
    """ Generates cartesian (x,y,z) and spherical (theta, phi) coordinates of a sphere
    Returns
    -------
    coords : named tuple
        holds cartesian (x,y,z) and spherical (theta, phi) coordinates
    """
    coords = namedtuple('coords', ['x', 'y', 'z', 'az', 'el'])
    az = _np.linspace(0, 2 * _np.pi, 360)
    el = _np.linspace(0, _np.pi, 181)
    coords.x = _np.outer(_np.cos(az), _np.sin(el))
    coords.y = _np.outer(_np.sin(az), _np.sin(el))
    coords.z = _np.outer(_np.ones(360), _np.cos(el))

    coords.el, coords.az = _np.meshgrid(_np.linspace(0, _np.pi, 181),
                                        _np.linspace(0, 2 * _np.pi, 360))
    return coords

def sph_harm_all(nMax, az, el, type='complex'):
    '''Compute all sphercial harmonic coefficients up to degree nMax.

    Parameters
    ----------
    nMax : (int)
        Maximum degree of coefficients to be returned. n >= 0

    az: (float), array_like
        Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.

    el : (float), array_like
        Elevation (colatitudinal) coordinate [0, pi], also called Phi.

    Returns
    -------
    y_mn : (complex float), array_like
        Complex spherical harmonics of degrees n [0 ... nMax] and all corresponding
        orders m [-n ... n], sampled at [az, el]. dim1 corresponds to az/el pairs,
        dim2 to oder/degree (m, n) pairs like 0/0, -1/1, 0/1, 1/1, -2/2, -1/2 ...
    '''
    m, n = mnArrays(nMax)
    mA, azA = _np.meshgrid(m, az)
    nA, elA = _np.meshgrid(n, el)
    return sph_harm(mA, nA, azA, elA, type=type)

def draw2dsurface(X, Y, zf):
    fig = plt.figure()
    ax = fig.gca(projection='3d')

    X, Y = np.meshgrid(X, Y)
    Z = X*0
    for i in range(len(X)):
        for j in range(len(X[0])):
            Z[i][j] = zf([X[i][j], Y[i][j]])

    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False)

    ax.set_zlim(np.min(Z.flatten()), np.max(Z.flatten()))

    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink=0.5, aspect=5)

    # plt.show()

def make_mask(self, shape):
        try:
            nrow, ncol = shape
        except TypeError:
            nrow = ncol = shape
        # HACK: to make the masks consistent with ``rand_position()``
        ix = self.xc - np.arange(ncol)
        iy = self.yc - np.arange(nrow)
        mx, my = np.meshgrid(ix, iy)
        rho, phi = self.cart2pol(mx, my)
        mask_rho = (rho >= self.rin) & (rho <= self.rout)
        mask_phi = (phi >= self.abegin) & (phi <= self.aend)
        if self.aend > 360:
            mask_phi |= (phi <= (self.aend-360))
        mask = mask_rho & mask_phi
        return mask

def mapFunction( x , y , func , ax = None, arrayInput = False, n = 10, title = None, **kwargs ) :
   """
      Plot function on a regular grid
        x : 1d array
        y : 1d array
        func : function to map
        arrayInput : False if func(x,y) , True if func( [x,y] )
   """

   if ax is None :
      fig , ax = plt.subplots()

   X , Y = np.meshgrid( x , y )

   if not arrayInput :
      Z = func( X.flatten() , Y.flatten() ).reshape(X.shape)
   else :
      Z = func( np.stack( [ X.flatten() , Y.flatten() ]) )

   ax.contourf( X , Y , Z , n , **kwargs)

   if title is not None : ax.set_title(title)

   return ax

def __init__(self, width, sample_spacing=0.025, samples=384):
        ''' Produces a Pinhole.

        Args:
            width (`float`): the width of the pinhole.

            sample_spacing (`float`): spacing of samples in the synthetic image.

            samples (`int`): number of samples per dimension in the synthetic image.

        '''
        self.width = width

        # produce coordinate arrays
        ext = samples / 2 * sample_spacing
        x, y = np.linspace(-ext, ext, samples), np.linspace(-ext, ext, samples)
        xv, yv = np.meshgrid(x, y)
        w = width / 2

        # paint a circle on a black background
        arr = np.zeros((samples, samples))
        arr[sqrt(xv**2 + yv**2) < w] = 1
        super().__init__(data=arr, sample_spacing=sample_spacing, synthetic=True)

def defineSensorLoc(self,XYZ):
        #############################################
        # DEFINE TRANSMITTER AND RECEIVER LOCATIONS
        #
        # XYZ: N X 3 array containing transmitter center locations
        # **NOTE** for this sensor, we know where the receivers are relative to transmitters
        self.TxLoc = XYZ

        dx,dy = np.meshgrid([-0.8,-0.4,0.,0.4,0.8],[-0.8,-0.4,0.,0.4,0.8])
        dx = mkvc(dx)
        dy = mkvc(dy)

        N = np.shape(XYZ)[0]

        X = np.kron(XYZ[:,0],np.ones((25))) + np.kron(np.ones((N)),dx)
        Y = np.kron(XYZ[:,1],np.ones((25))) + np.kron(np.ones((N)),dy)
        Z = np.kron(XYZ[:,2],np.ones((25)))

        self.RxLoc = np.c_[X,Y,Z]

        self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((25)))
        self.RxComp = np.kron(3*np.ones(np.shape(XYZ)[0]),np.ones((25)))

def videoSeq(number_cells,inputIm,simtime,resolution,video_step):

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    X = np.arange(number_cells)
    Y = X
    X, Y = np.meshgrid(X, Y)
    Z = inputIm[:,:,0]
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='coolwarm',
                           linewidth=0, antialiased=False)
    ax.set_title('Input stimulus. Time = 0.0 ms',y=1.08)
    fig.show()
    ax.axes.figure.canvas.draw()

    for t in np.arange(0.0,int(simtime/resolution),video_step/resolution):
        surf.remove()
        Z = inputIm[:,:,t]
        surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='coolwarm',
        linewidth=0, antialiased=False)
        ax.set_title('Input stimulus. Time = %s ms'%str(t),y=1.08)
        ax.axes.figure.canvas.draw()

# Estimation of Local Field Potential(LFP)

def build_random_variables(self, **kwargs):
        # All this is done just once per batch (i.e. until `clear_random_variables` is called)
        np.random.seed()
        imshape = kwargs.get('imshape')
        # Build and scale random fields
        random_field_x = np.random.uniform(-1, 1, imshape) * self.alpha
        random_field_y = np.random.uniform(-1, 1, imshape) * self.alpha
        # Smooth random field (this has to be done just once per reset)
        sdx = gaussian_filter(random_field_x, self.sigma, mode='reflect')
        sdy = gaussian_filter(random_field_y, self.sigma, mode='reflect')
        # Make meshgrid
        x, y = np.meshgrid(np.arange(imshape[1]), np.arange(imshape[0]))
        # Make inversion coefficient
        _inverter = 1. if not self.invert else -1.
        # Distort meshgrid indices (invert if required)
        flow_y, flow_x = (y + _inverter * sdy).reshape(-1, 1), (x + _inverter * sdx).reshape(-1, 1)
        # Set random states
        self.set_random_variable('flow_x', flow_x)
        self.set_random_variable('flow_y', flow_y)

def test_bowtie(self):
        # Test the bowtie filter with sine waves
        wDeg = 3
        nPix = 257
        sf = 1
        orientation = 0

        x = y = np.linspace(-wDeg // 2, wDeg // 2, nPix + 1)
        u, v = np.meshgrid(x, y)

        ramp = np.sin(orientation * np.pi / 180) * u
        ramp -= np.cos(orientation * np.pi / 180) * v
        img = np.sin(2 * np.pi * sf * ramp)
        fimg = fft2(img)

        orientation_widths = [1, 10, 20, 40, 80, 100]
        for x in orientation_widths:
            filt = OrientationFilter('bowtie', 90, x, nPix, .2, nPix + 1,
                                     'triangle')
            filt = filt.filter
            filt = 1 - filt
            filt = fftshift(filt)
            out = ifft2(fimg * filt).real.astype(int)

            self.assertEqual(np.sum(out), 0)

def test_plot_bounds_2d(self):
        x = np.arange(1, 5)
        y = np.arange(5, 10)
        x2d, y2d = np.meshgrid(x, y)
        x_bnds = np.arange(0.5, 4.51, 1.0)
        y_bnds = np.arange(4.5, 9.51, 1.0)
        # the borders are not modified
        x_bnds[0] = 1.0
        x_bnds[-1] = 4.0
        y_bnds[0] = 5.0
        y_bnds[-1] = 9.0
        x2d_bnds, y2d_bnds = np.meshgrid(x_bnds, y_bnds)
        d = psyd.CFDecoder()
        # test x bounds
        bounds = d.get_plotbounds(xr.Variable(('y', 'x'), x2d))
        self.assertAlmostArrayEqual(bounds, x2d_bnds)

        # test y bounds
        bounds = d.get_plotbounds(xr.Variable(('y', 'x'), y2d))
        self.assertAlmostArrayEqual(bounds, y2d_bnds)

def draw(X, pred, means, covariances, output):
    xp = cupy.get_array_module(X)
    for i in six.moves.range(2):
        labels = X[pred == i]
        if xp is cupy:
            labels = labels.get()
        plt.scatter(labels[:, 0], labels[:, 1], c=np.random.rand(3))
    if xp is cupy:
        means = means.get()
        covariances = covariances.get()
    plt.scatter(means[:, 0], means[:, 1], s=120, marker='s', facecolors='y',
                edgecolors='k')
    x = np.linspace(-5, 5, 1000)
    y = np.linspace(-5, 5, 1000)
    X, Y = np.meshgrid(x, y)
    for i in six.moves.range(2):
        Z = mlab.bivariate_normal(X, Y, np.sqrt(covariances[i][0]),
                                  np.sqrt(covariances[i][1]),
                                  means[i][0], means[i][1])
        plt.contour(X, Y, Z)
    plt.savefig(output)

def plot(self, nmin=-3.5, nmax=1.5):
        """Plots the field magnitude."""
        x, y = meshgrid(
            linspace(XMIN/ZOOM+XOFFSET, XMAX/ZOOM+XOFFSET, 200),
            linspace(YMIN/ZOOM, YMAX/ZOOM, 200))
        z = zeros_like(x)
        for i in range(x.shape[0]):
            for j in range(x.shape[1]):
                z[i, j] = log10(self.magnitude([x[i, j], y[i, j]]))
        levels = arange(nmin, nmax+0.2, 0.2)
        cmap = pyplot.cm.get_cmap('plasma')
        pyplot.contourf(x, y, numpy.clip(z, nmin, nmax),
                        10, cmap=cmap, levels=levels, extend='both')


# pylint: disable=too-few-public-methods

def checkDiamondLattice(JJ, II, L, d, offset, CSmooth, doPlot = False):
    arr = np.arange(0, L+1, 2*d)
    narr = -arr
    arr = narr.tolist()[::-1] + arr.tolist()[1::]
    X, Y = np.meshgrid(arr, arr)
    Q1 = np.zeros((X.size, 2))
    Q1[:, 0] = Y.flatten()
    Q1[:, 1] = X.flatten()
    arr2 = np.arange(d, L+1, 2*d)
    narr = -arr2
    arr2 = narr.tolist()[::-1] + arr2.tolist()
    X, Y = np.meshgrid(arr2, arr2)
    Q2 = np.zeros((X.size, 2))
    Q2[:, 0] = Y.flatten()
    Q2[:, 1] = X.flatten()
    Q = np.concatenate((Q1, Q2), 0)
    return checkLattice(Q, JJ, II, L, d, offset, CSmooth, doPlot)

def make2ShakingCircles(NFrames = 200, T1 = 20, T2 = 20*np.pi, A1 = 10, A2 = 10, ydim = 50):
    print("T1 = ", T1, ", T2 = ", T2)
    IDims = (ydim, ydim*2, 3)
    I = np.zeros((NFrames, ydim*ydim*2*3))

    [X, Y] = np.meshgrid(np.arange(ydim*2), np.arange(ydim))
    yc = float(ydim/2)
    R = float(ydim/8)
    for i in range(NFrames):
        x1c = float(ydim/2) - A1*np.cos(2*np.pi*i/T1)
        x2c = 3*float(ydim/2) - A2*np.cos(2*np.pi*i/T2)
        f = np.zeros((X.shape[0], X.shape[1], 3))
        for k in range(3):
            f[:, :, k] = ((X-x1c)**2 + (Y-yc)**2 < R**2) + ((X-x2c)**2 + (Y-yc)**2 < R**2)
        I[i, :] = f.flatten()
    I = 0.5*(I + 1)
    return (I, IDims)

#############################################################
####                 OTHER VIDEO TOOLS                  #####
#############################################################

def doTotalOrderExperiment(N, binaryWeights = False):
    I, J = np.meshgrid(np.arange(N), np.arange(N))
    I = I[np.triu_indices(N, 1)]
    J = J[np.triu_indices(N, 1)]
    #[I, J] = [I[0:N-1], J[0:N-1]]
    NEdges = len(I)
    R = np.zeros((NEdges, 2))
    R[:, 0] = J
    R[:, 1] = I

    #W = np.random.rand(NEdges)
    W = np.ones(NEdges)

    #Note: When using binary weights, Y is not necessarily a cocycle
    Y = I - J
    if binaryWeights:
        Y = np.ones(NEdges)

    (s, I, H) = doHodge(R, W, Y, verbose = True)
    printConsistencyRatios(Y, I, H, W)

#Random flip experiment with linear order

def plot_cost_to_go_mountain_car(env, estimator, num_tiles=20):
    x = np.linspace(env.observation_space.low[0], env.observation_space.high[0], num=num_tiles)
    y = np.linspace(env.observation_space.low[1], env.observation_space.high[1], num=num_tiles)
    X, Y = np.meshgrid(x, y)
    Z = np.apply_along_axis(lambda _: -np.max(estimator.predict(_)), 2, np.dstack([X, Y]))

    fig = plt.figure(figsize=(10, 5))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
                           cmap=matplotlib.cm.coolwarm, vmin=-1.0, vmax=1.0)
    ax.set_xlabel('Position')
    ax.set_ylabel('Velocity')
    ax.set_zlabel('Value')
    ax.set_title("Mountain \"Cost To Go\" Function")
    fig.colorbar(surf)
    plt.show()

def test_megkde_2d_basic():
    # Draw from normal, fit KDE, see if sampling from kde's pdf recovers norm
    np.random.seed(1)
    data = np.random.multivariate_normal([0, 1], [[1.0, 0.], [0., 0.75 ** 2]], size=10000)
    xs, ys = np.linspace(-4, 4, 50), np.linspace(-4, 4, 50)
    xx, yy = np.meshgrid(xs, ys, indexing='ij')
    samps = np.vstack((xx.flatten(), yy.flatten())).T
    zs = MegKDE(data).evaluate(samps).reshape(xx.shape)
    zs_x = zs.sum(axis=1)
    zs_y = zs.sum(axis=0)
    cs_x = zs_x.cumsum()
    cs_x /= cs_x[-1]
    cs_x[0] = 0
    cs_y = zs_y.cumsum()
    cs_y /= cs_y[-1]
    cs_y[0] = 0
    samps_x = interp1d(cs_x, xs)(np.random.uniform(size=10000))
    samps_y = interp1d(cs_y, ys)(np.random.uniform(size=10000))
    mu_x, std_x = norm.fit(samps_x)
    mu_y, std_y = norm.fit(samps_y)
    assert np.isclose(mu_x, 0, atol=0.1)
    assert np.isclose(std_x, 1.0, atol=0.1)
    assert np.isclose(mu_y, 1, atol=0.1)
    assert np.isclose(std_y, 0.75, atol=0.1)

def test_grid_data(self):
        x, y = np.linspace(-3, 3, 200), np.linspace(-5, 5, 200)
        xx, yy = np.meshgrid(x, y, indexing='ij')
        xs, ys = xx.flatten(), yy.flatten()
        chain = np.vstack((xs, ys)).T
        pdf = (1 / (2 * np.pi)) * np.exp(-0.5 * (xs * xs + ys * ys / 4))
        c = ChainConsumer()
        c.add_chain(chain, parameters=['x', 'y'], weights=pdf, grid=True)
        summary = c.analysis.get_summary()
        x_sum = summary['x']
        y_sum = summary['y']
        expected_x = np.array([-1.0, 0.0, 1.0])
        expected_y = np.array([-2.0, 0.0, 2.0])
        threshold = 0.1
        assert np.all(np.abs(expected_x - x_sum) < threshold)
        assert np.all(np.abs(expected_y - y_sum) < threshold)

def surface_bounds(f, bounds, cmap=cm.magma_r, resolution=[100, 100], type='2D'):
    """
    Plot a function over a grid
    :param f: function to plot in Z
    :param bounds:
    :param cmap: cmap
    :param resolution:
    :param type: '2D' or '3D'
    :return:
    """
    resolution = np.array(resolution)
    assert bounds.get_n_dim() == 2, 'Bounds are not 2D'

    x = np.linspace(start=bounds.get_min(0), stop=bounds.get_max(0), num=resolution[0])
    y = np.linspace(start=bounds.get_min(1), stop=bounds.get_max(1), num=resolution[1])
    X, Y = np.meshgrid(x, y)
    Z = f(np.vstack((X.flatten(), Y.flatten())).T).reshape(X.shape)  # Evaluate grid on the function
    surface(X=X, Y=Y, Z=Z, type=type, cmap=cmap)


# def surface_grid(f, x, y)

def visualize(config, vae):
    if(config['n_z'] != 2):
        print("Skipping visuals since n_z is not 2")
        return
    nx = ny = 20
    x_values = np.linspace(-3, 3, nx)
    y_values = np.linspace(-3, 3, ny)

    canvas = np.empty((28*ny, 28*nx))
    for i, yi in enumerate(x_values):
        for j, xi in enumerate(y_values):
            z_mu = np.array([[xi, yi]])
            x_mean = vae.generate(np.tile(z_mu, [config['batch_size'], 1]))
            canvas[(nx-i-1)*28:(nx-i)*28, j*28:(j+1)*28] = x_mean[0].reshape(28, 28)

    plt.figure(figsize=(8, 10))        
    Xi, Yi = np.meshgrid(x_values, y_values)
    plt.imshow(canvas, origin="upper")
    plt.tight_layout()

    img = "samples/2d-visualization.png"
    plt.savefig(img)
    hc.io.sample(config, [{"label": "2d visualization", "image": img}])

def dihedral_transform_batch(x):

    g = np.random.randint(low=0, high=8, size=x.shape[0])

    h, w = x.shape[-2:]
    hh = (h - 1) / 2.
    hw = (w - 1) / 2.

    I, J = np.meshgrid(np.linspace(-hh, hh, x.shape[-2]), np.linspace(-hw, hw, x.shape[-1]))
    C = np.r_[[I, J]]
    D4C = np.einsum('...ij,jkl->...ikl', D4, C)
    D4C[:, 0] += hh
    D4C[:, 1] += hw
    D4C = D4C.astype(int)

    x_out = np.empty_like(x)
    for i in range(x.shape[0]):
        I, J = D4C[g[i]]
        x_out[i, :] = x[i][:, J, I]

    return x_out