python类meshgrid()的实例源码

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

def saturation_index_countour(lab, elem1, elem2, Ks, labels=False):
    plt.figure()
    plt.title('Saturation index %s%s' % (elem1, elem2))
    resoluion = 100
    n = math.ceil(lab.time.size / resoluion)
    plt.xlabel('Time')
    z = np.log10((lab.species[elem1]['concentration'][:, ::n] + 1e-8) * (
        lab.species[elem2]['concentration'][:, ::n] + 1e-8) / lab.constants[Ks])
    lim = np.max(abs(z))
    lim = np.linspace(-lim - 0.1, +lim + 0.1, 51)
    X, Y = np.meshgrid(lab.time[::n], -lab.x)
    plt.xlabel('Time')
    CS = plt.contourf(X, Y, z, 20, cmap=ListedColormap(sns.color_palette(
        "RdBu_r", 101)), origin='lower', levels=lim, extend='both')
    if labels:
        plt.clabel(CS, inline=1, fontsize=10, colors='w')
    # cbar = plt.colorbar(CS)
    if labels:
        plt.clabel(CS, inline=1, fontsize=10, colors='w')
    cbar = plt.colorbar(CS)
    plt.ylabel('Depth')
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False)
    cbar.ax.set_ylabel('Saturation index %s%s' % (elem1, elem2))
    return ax

def contour_plot_of_rates(lab, r, labels=False, last_year=False):
    plt.figure()
    plt.title('{}'.format(r))
    resoluion = 100
    n = math.ceil(lab.time.size / resoluion)
    if last_year:
        k = n - int(1 / lab.dt)
    else:
        k = 1
    z = lab.estimated_rates[r][:, k - 1:-1:n]
    # lim = np.max(np.abs(z))
    # lim = np.linspace(-lim - 0.1, +lim + 0.1, 51)
    X, Y = np.meshgrid(lab.time[k::n], -lab.x)
    plt.xlabel('Time')
    CS = plt.contourf(X, Y, z, 20, cmap=ListedColormap(
        sns.color_palette("Blues", 51)))
    if labels:
        plt.clabel(CS, inline=1, fontsize=10, colors='w')
    cbar = plt.colorbar(CS)
    plt.ylabel('Depth')
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False)
    cbar.ax.set_ylabel('Rate %s [M/V/T]' % r)
    return ax

def contour_plot_of_delta(lab, element, labels=False, last_year=False):
    plt.figure()
    plt.title('Rate of %s consumption/production' % element)
    resoluion = 100
    n = math.ceil(lab.time.size / resoluion)
    if last_year:
        k = n - int(1 / lab.dt)
    else:
        k = 1
    z = lab.species[element]['rates'][:, k - 1:-1:n]
    lim = np.max(np.abs(z))
    lim = np.linspace(-lim - 0.1, +lim + 0.1, 51)
    X, Y = np.meshgrid(lab.time[k:-1:n], -lab.x)
    plt.xlabel('Time')
    CS = plt.contourf(X, Y, z, 20, cmap=ListedColormap(sns.color_palette(
        "RdBu_r", 101)), origin='lower', levels=lim, extend='both')
    if labels:
        plt.clabel(CS, inline=1, fontsize=10, colors='w')
    cbar = plt.colorbar(CS)
    plt.ylabel('Depth')
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False)
    cbar.ax.set_ylabel('Rate of %s change $[\Delta/T]$' % element)
    return ax

def _meshgrid(self):
        with tf.variable_scope('_meshgrid'):
            x_t = tf.matmul(tf.ones(shape=tf.stack([self.out_height, 1])),
                            tf.transpose(tf.expand_dims(tf.linspace(-1.0, 1.0, self.out_width), 1), [1, 0]))
            y_t = tf.matmul(tf.expand_dims(tf.linspace(-1.0, 1.0, self.out_height), 1),
                            tf.ones(shape=tf.stack([1, self.out_width])))
            x_t_flat = tf.reshape(x_t, (1, -1))
            y_t_flat = tf.reshape(y_t, (1, -1))
            px,py = tf.stack([x_t_flat],axis=2),tf.stack([y_t_flat],axis=2)
            #source control points
            x,y = tf.linspace(-1.,1.,self.Column_controlP_number),tf.linspace(-1.,1.,self.Row_controlP_number)
            x,y = tf.meshgrid(x,y)
            xs,ys = tf.transpose(tf.reshape(x,(-1,1))),tf.transpose(tf.reshape(y,(-1,1)))
            cpx,cpy = tf.transpose(tf.stack([xs],axis=2),perm=[1,0,2]),tf.transpose(tf.stack([ys],axis=2),perm=[1,0,2])
            px, cpx = tf.meshgrid(px,cpx);py, cpy = tf.meshgrid(py,cpy)           
            #Compute distance R
            Rx,Ry = tf.square(tf.subtract(px,cpx)),tf.square(tf.subtract(py,cpy))
            R = tf.add(Rx,Ry)          
            R = tf.multiply(R,tf.log(tf.clip_by_value(R,1e-10,1e+10)))
            #Source coordinates
            ones = tf.ones_like(x_t_flat) 
            grid = tf.concat([ones, x_t_flat, y_t_flat,R],0)
            grid = tf.reshape(grid,[-1])
            grid = tf.reshape(grid,[self.Column_controlP_number*self.Row_controlP_number+3,self.out_height*self.out_width])
            return grid

def Affine_test(self,N,sizex,sizey,sizez,times,stop_time,typeofT,colors):
        for i in range(times):
            # Theta
            idx = np.random.uniform(-1, 1);idy = np.random.uniform(-1, 1);idz = np.random.uniform(-1, 1)
            swithx = np.random.uniform(0,1);swithy = np.random.uniform(0,1);swithz = np.random.uniform(0,1)
            rotatex = np.random.uniform(-1, 1);rotatey = np.random.uniform(-1, 1);rotatez = np.random.uniform(-1, 1)
            cx = np.array([idx,rotatey,rotatez,swithx]);cy = np.array([rotatex,idy,rotatez,swithy]);cz = np.array([rotatex,rotatey,idz,swithz])
            # Source Grid
            x = np.linspace(-sizex, sizex, N);y = np.linspace(-sizey, sizey, N);z = np.linspace(-sizez, sizez, N)
            x, y, z = np.meshgrid(x, y, z)
            xgs, ygs, zgs = x.flatten(), y.flatten(),z.flatten()
            gps = np.vstack([xgs, ygs, zgs, np.ones_like(xgs)]).T
            # transform
            xgt = np.dot(gps, cx);ygt = np.dot(gps, cy);zgt = np.dot(gps, cz)
            # display
            showIm = ShowImage()
            showIm.Show_transform(xgs,ygs,zgs,xgt,ygt,zgt,sizex,sizey,sizez,stop_time,typeofT,N,colors)

def TPS_test(self,N,sizex,sizey,sizez,control_num,show_times,stop_time,typeofT,colors):
        for i in range(show_times):
            # source control points
            x, y, z = np.meshgrid(np.linspace(-1, 1, control_num),np.linspace(-1, 1, control_num), np.linspace(-1, 1, control_num))
            xs = x.flatten();ys = y.flatten();zs = z.flatten()
            cps = np.vstack([xs, ys, zs]).T
            # target control points
            xt = xs + np.random.uniform(-0.3, 0.3, size=xs.size);yt = ys + np.random.uniform(-0.3, 0.3, size=ys.size); zt = zs + np.random.uniform(-0.3, 0.3, size=zs.size)
            # construct T
            T = self.makeT(cps)
            # solve cx, cy (coefficients for x and y)
            xtAug = np.concatenate([xt, np.zeros(4)]);ytAug = np.concatenate([yt, np.zeros(4)]);ztAug = np.concatenate([zt, np.zeros(4)])
            cx = nl.solve(T, xtAug); cy = nl.solve(T, ytAug); cz = nl.solve(T, ztAug)
            # dense grid
            x = np.linspace(-sizex, sizex, 2*sizex); y = np.linspace(-sizey, sizey, 2*sizey); z = np.linspace(-sizez, sizez, 2*sizez)
            x, y, z = np.meshgrid(x, y, z)
            xgs, ygs, zgs = x.flatten(), y.flatten(), z.flatten()
            gps = np.vstack([xgs, ygs, zgs]).T
            # transform
            pgLift = self.liftPts(gps, cps) # [N x (K+3)]
            xgt = np.dot(pgLift, cx.T);ygt = np.dot(pgLift, cy.T); zgt = np.dot(pgLift,cz.T)
            # display
            showIm = ShowImage()
            showIm.Show_transform(xgs,ygs,zgs,xgt,ygt,zgt,sizex,sizey,sizez,stop_time,typeofT,N,colors)

def _initialize_grid(self, mmf):
        self.nr, self.nlat, self.nlon\
            = struct.unpack("3i", mmf.read(12))
        self.dr, self.dlat, self.dlon\
            = struct.unpack("3f", mmf.read(12))
        self.r0, self.lat0, self.lon0\
            = struct.unpack("3f", mmf.read(12))
        self.mr = self.r0 + (self.nr - 1) * self.dr
        self.mlat = self.lat0 + (self.nlat - 1) * self.dlat
        self.mlon = self.lon0 + (self.nlon - 1) * self.dlon
        self.dtheta = radians(self.dlat)
        self.dphi = radians(self.dlon)
        self.ntheta, self.nphi = self.nlat, self.nlon
        self.theta0 = radians(90 - self.lat0)
        self.phi0 = radians(self.lon0)
        r = [self.r0 + self.dr * ir for ir in range(self.nr)]
        theta = [radians(90. - self.lat0 - self.dlat * ilat)
                 for ilat in range(self.nlat)]
        phi = [radians((self.lon0 + self.dlon * ilon) % 360.)
               for ilon in range(self.nlon)]
        R, T, P = np.meshgrid(r, theta, phi, indexing='ij')
        self.nodes = {'r': R, 'theta': T, 'phi': P}

def mapinterpolation(data, x=None, y=None, interpx=1, interpy=1):
    """Interpolate a 2D map."""
    # Get data size :
    dim2, dim1 = data.shape
    # Define xticklabel and yticklabel :
    if x is None:
        x = np.arange(0, dim1, interpx)
    if y is None:
        y = np.arange(0, dim2, interpy)
    # Define the meshgrid :
    Xi, Yi = np.meshgrid(
        np.arange(0, dim1 - 1, interpx), np.arange(0, dim2 - 1, interpy))
    # 2D interpolation :
    datainterp = interp2(data, Xi, Yi)
    # Linearly interpolate vectors :
    xveci = np.linspace(x[0], x[-1], datainterp.shape[0])
    yveci = np.linspace(y[0], y[-1], datainterp.shape[1])
    return datainterp, xveci, yveci