def test_plot_surf_stat_map():
    mesh = _generate_surf()
    rng = np.random.RandomState(0)
    bg = rng.randn(mesh[0].shape[0], )
    data = 10 * rng.randn(mesh[0].shape[0], )

    # Plot mesh with stat map
    plot_surf_stat_map(mesh, stat_map=data)
    plot_surf_stat_map(mesh, stat_map=data, colorbar=True)
    plot_surf_stat_map(mesh, stat_map=data, alpha=1)

    # Plot mesh with background and stat map
    plot_surf_stat_map(mesh, stat_map=data, bg_map=bg)
    plot_surf_stat_map(mesh, stat_map=data, bg_map=bg,
                       bg_on_data=True, darkness=0.5)
    plot_surf_stat_map(mesh, stat_map=data, bg_map=bg, colorbar=True,
                       bg_on_data=True, darkness=0.5)

    # Apply threshold
    plot_surf_stat_map(mesh, stat_map=data, bg_map=bg,
                       bg_on_data=True, darkness=0.5,
                       threshold=0.3)
    plot_surf_stat_map(mesh, stat_map=data, bg_map=bg, colorbar=True,
                       bg_on_data=True, darkness=0.5,
                       threshold=0.3)

    # Change vmax
    plot_surf_stat_map(mesh, stat_map=data, vmax=5)
    plot_surf_stat_map(mesh, stat_map=data, vmax=5, colorbar=True)

    # Change colormap
    plot_surf_stat_map(mesh, stat_map=data, cmap='cubehelix')
    plot_surf_stat_map(mesh, stat_map=data, cmap='cubehelix', colorbar=True)

    # Plot to axes
    axes = plt.subplots(ncols=2, subplot_kw={'projection': '3d'})[1]
    for ax in axes.flatten():
        plot_surf_stat_map(mesh, stat_map=data, ax=ax)
    axes = plt.subplots(ncols=2, subplot_kw={'projection': '3d'})[1]
    for ax in axes.flatten():
        plot_surf_stat_map(mesh, stat_map=data, ax=ax, colorbar=True)

    fig = plot_surf_stat_map(mesh, stat_map=data, colorbar=False)
    assert len(fig.axes) == 1
    # symmetric_cbar
    fig = plot_surf_stat_map(
        mesh, stat_map=data, colorbar=True, symmetric_cbar=True)
    assert len(fig.axes) == 2
    yticklabels = fig.axes[1].get_yticklabels()
    first, last = yticklabels[0].get_text(), yticklabels[-1].get_text()
    assert float(first) == - float(last)
    # no symmetric_cbar
    fig = plot_surf_stat_map(
        mesh, stat_map=data, colorbar=True, symmetric_cbar=False)
    assert len(fig.axes) == 2
    yticklabels = fig.axes[1].get_yticklabels()
    first, last = yticklabels[0].get_text(), yticklabels[-1].get_text()
    assert float(first) != - float(last)
    # Save execution time and memory
    plt.close()

def test_cursor_data():
    from matplotlib.backend_bases import MouseEvent

    fig, ax = plt.subplots()
    im = ax.imshow(np.arange(100).reshape(10, 10), origin='upper')

    x, y = 4, 4
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) == 44

    # Now try for a point outside the image
    # Tests issue #4957
    x, y = 10.1, 4
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) is None

    # Hmm, something is wrong here... I get 0, not None...
    # But, this works further down in the tests with extents flipped
    #x, y = 0.1, -0.1
    #xdisp, ydisp = ax.transData.transform_point([x, y])
    #event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    #z = im.get_cursor_data(event)
    #assert z is None, "Did not get None, got %d" % z

    ax.clear()
    # Now try with the extents flipped.
    im = ax.imshow(np.arange(100).reshape(10, 10), origin='lower')

    x, y = 4, 4
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) == 44

    fig, ax = plt.subplots()
    im = ax.imshow(np.arange(100).reshape(10, 10), extent=[0, 0.5, 0, 0.5])

    x, y = 0.25, 0.25
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) == 55

    # Now try for a point outside the image
    # Tests issue #4957
    x, y = 0.75, 0.25
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) is None

    x, y = 0.01, -0.01
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) is None

    def test_heatmap_ticklabel_rotation(self):

        f, ax = plt.subplots(figsize=(2, 2))
        mat.heatmap(self.df_norm, ax=ax)

        for t in ax.get_xticklabels():
            nt.assert_equal(t.get_rotation(), 0)

        for t in ax.get_yticklabels():
            nt.assert_equal(t.get_rotation(), 90)

        plt.close(f)

        df = self.df_norm.copy()
        df.columns = [str(c) * 10 for c in df.columns]
        df.index = [i * 10 for i in df.index]

        f, ax = plt.subplots(figsize=(2, 2))
        mat.heatmap(df, ax=ax)

        for t in ax.get_xticklabels():
            nt.assert_equal(t.get_rotation(), 90)

        for t in ax.get_yticklabels():
            nt.assert_equal(t.get_rotation(), 0)

        plt.close(f)

 def decode_uniform_samples_from_latent_space(_):
     fig, ax = plt.subplots()
     nx = ny = 20
     extent_x = extent_y = [-3, 3]
     extent = numpy.array(extent_x + extent_y)
     x_values = numpy.linspace(*(extent_x + [nx]))
     y_values = numpy.linspace(*(extent_y + [nx]))
     full_extent = extent * (nx + 1) / float(nx)
     canvas = numpy.empty((28 * ny, 28 * nx))
     for ii, yi in enumerate(x_values):
         for j, xi in enumerate(y_values):
             n = ii * nx + j + 1
             sys.stdout.write("\rsampling p(X|z), sample %d/%d" % (n, nx*ny))
             sys.stdout.flush()
             np_z = numpy.array([[xi, yi]])
             x_mean = sess.run(prior_model(latent=numpy.reshape(np_z, newshape=(1, LATENT_DIM)))[0])
             canvas[(nx - ii - 1) * 28:(nx - ii) * 28, j * 28:(j + 1) * 28] = x_mean[0].reshape(28, 28)
     with seaborn.axes_style('ticks'):
         seaborn.set_context(context='notebook', font_scale=1.75)
         fig, ax = plt.subplots(figsize=(12, 9))
     ax.imshow(canvas, extent=full_extent)
     ax.xaxis.set_ticks(numpy.linspace(*(extent_x + [nx])))
     ax.yaxis.set_ticks(numpy.linspace(*(extent_y + [ny])))
     ax.set_xlabel('z_1')
     ax.set_ylabel('z_2')
     ax.set_title('P(X|z); decoding latent space; (CONV, BNAE, IND_ERROR) = (%d,%d,%d)' % (CONV, BNAE, IND_ERROR))
     plt.show()
     plt.savefig(os.path.join(FLAGS.viz_dir, 'P(X|z).png'))
     return fig, ax

def accuracy(target, prediction, label="Classifier", c=np.zeros((0,0))):
    correct = (target == prediction)
    correct = np.array((correct, correct))
    compare = np.array((target, prediction))
    
    showC = c != np.zeros((0,0))
    
    if (showC):
        fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, figsize=(6,10))
    else:
        fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(6,8))
    
    dim = [0,compare.shape[1],0,compare.shape[0]]
    ax1.imshow(compare, extent=dim, aspect='auto', interpolation='nearest')
    ax1.set_title(label + ": Prediction vs. Target")
    
    imgPlt = ax2.imshow(correct, extent=dim, aspect='auto', interpolation='nearest')
    imgPlt.set_cmap('RdYlGn')
    ax2.set_title(label + " Prediction Accuracy")
    
    if (showC):
        ax3.plot(c)
        ax3.set_title("Concentration")
        ax3.set_yscale('log')
        ax3.set_ylim(0.02,0.7)
    
    plt.draw()

def main():
    from matplotlib import pyplot as plt
    
    fig, ax = plt.subplots()

    L = 80. # granule cell layer thickness in um
    r_values = np.arange(0.1, 160, 0.02)
    rho = 6.6e-4 # glom density in um^-3

    for l in np.arange(4.65, 4.8, 0.01):
        d_values = np.array([d(r,l,L) for r in r_values])
        mean = (d_values * r_values).sum() / d_values.sum()
        fraction_above_30 = d_values[r_values>30].sum() / d_values.sum()
        fraction_above_40 = d_values[r_values>39.5].sum() / d_values.sum()

        #print l, mean, fraction_above_30, fraction_above_40
        ax.plot(r_values, d_values)
    
    l = 5.
    k_palkovits = 4.17

    n_dendrites_range = np.arange(1.,8.,1.)
    degree_distribution = np.array([poisson(k_palkovits, n) for n in n_dendrites_range])
    fig, ax = plt.subplots()
    print degree_distribution
    ax.bar(n_dendrites_range, degree_distribution)
        
    plt.show()

def check_scaling_distribution(x_unscl, x_scl, y_unscl, y_scl, lat, lev,
                               figpath):
    # For input variables
    fig, ax = plt.subplots(2, 2)
    _plot_distribution(unpack(x_unscl, 'T'), lat, lev, fig, ax[0, 0],
                       './figs/', 'T (unscaled) [K]', '')
    _plot_distribution(unpack(x_scl, 'T'), lat, lev, fig, ax[0, 1],
                       './figs/', 'T (scaled) []', '')
    _plot_distribution(unpack(x_unscl, 'q'), lat, lev, fig, ax[1, 0],
                       './figs/', 'q (unscaled) [g/kg]', '')
    _plot_distribution(unpack(x_scl, 'q'), lat, lev, fig, ax[1, 1],
                       './figs/', 'q (scaled) []', '')
    fig.savefig(figpath + 'input_scaling_check.png', bbox_inches='tight',
                dpi=450)
    plt.close()
    # For output variables
    fig, ax = plt.subplots(2, 2)
    _plot_distribution(unpack(y_unscl, 'T'), lat, lev, fig, ax[0, 0],
                       './figs/', 'T tend (unscaled) [K/day]', '')
    _plot_distribution(unpack(y_scl, 'T'), lat, lev, fig, ax[0, 1],
                       './figs/', 'T tend (scaled) []', '')
    _plot_distribution(unpack(y_unscl, 'q'), lat, lev, fig, ax[1, 0],
                       './figs/', 'q tend (unscaled) [g/kg/day]', '')
    _plot_distribution(unpack(y_scl, 'q'), lat, lev, fig, ax[1, 1],
                       './figs/', 'q tend(scaled) []', '')
    fig.savefig(figpath + 'output_scaling_check.png', bbox_inches='tight',
                dpi=450)
    plt.close()

def plot_toy_data(algo_name, features, labels, thetas):
    """
    Plots the toy data in 2D.
    Arguments:
    * features - an Nx2 ndarray of features (points)
    * labels - a length-N vector of +1/-1 labels
    * thetas - the tuple (theta, theta_0) that is the output of the learning algorithm
    * algorithm - the string name of the learning algorithm used
    """
    # plot the points with labels represented as colors
    plt.subplots()
    colors = ['b' if label == 1 else 'r' for label in labels]
    plt.scatter(features[:, 0], features[:, 1], s=40, c=colors)
    xmin, xmax = plt.axis()[:2]

    # plot the decision boundary
    theta, theta_0 = thetas
    xs = np.linspace(xmin, xmax)
    ys = -(theta[0]*xs + theta_0) / (theta[1] + 1e-16)
    plt.plot(xs, ys, 'k-')

    # show the plot
    algo_name = ' '.join((word.capitalize() for word in algo_name.split(' ')))
    plt.suptitle('Classified Toy Data ({})'.format(algo_name))
    plt.show()

def _get_axes(dim, axes=None, triangle=False, subplots_kwargs=dict()):
    """
    Parameters
    ----------
    dim : int
        Dimensionality of the orbit.
    axes : array_like (optional)
        Array of matplotlib Axes objects.
    triangle : bool (optional)
        Make a triangle plot instead of plotting all projections in a single row.
    subplots_kwargs : dict (optional)
        Dictionary of kwargs passed to :func:`~matplotlib.pyplot.subplots`.
    """

    import matplotlib.pyplot as plt

    if dim == 3:
        if triangle and axes is None:
            figsize = subplots_kwargs.pop('figsize', (8,8))
            sharex = subplots_kwargs.pop('sharex', True)
            sharey = subplots_kwargs.pop('sharey', True)
            fig,axes = plt.subplots(2,2,figsize=figsize, sharex=sharex, sharey=sharey,
                                    **subplots_kwargs)
            axes[0,1].set_visible(False)
            axes = axes.flat
            axes = [axes[0],axes[2],axes[3]]

        elif triangle and axes is not None:
            try:
                axes = axes.flat
            except:
                pass

            if len(axes) == 4:
                axes = [axes[0],axes[2],axes[3]]

        elif not triangle and axes is None:
            figsize = subplots_kwargs.pop('figsize', (10,3.5))
            fig,axes = plt.subplots(1, 3, figsize=figsize, **subplots_kwargs)

    elif dim <= 2:
        if axes is not None:
            try:
                len(axes)
            except TypeError:  # single axes object
                axes = [axes]

        else:
            if dim ==1:
                figsize = subplots_kwargs.pop('figsize', (14,8))
            elif dim == 2:
                figsize = subplots_kwargs.pop('figsize', (8,8))

            fig,axes = plt.subplots(1, 1, figsize=figsize, **subplots_kwargs)
            axes = [axes]

    else:
        raise ValueError("Orbit must have dimensions <= 3.")

    return axes

def save_plot_data(interval,data_1, data_2, data_3, names, style, p_title, s_title, format):
  close='all'
  if ( style == 0): 
    
    fig, ax = plt.subplots()
    ax.plot(interval,data_1, color='black', label=names[1])
    ax.plot(interval,data_2, color='red' , label=names[2])
    ax.plot(interval,data_3, color='blue',linestyle='-', label=names[3])
    ax.set_yscale('symlog')
    plt.title(p_title)
    plt.xlabel('Liquid') #das r hier markiert, dass jetz latex code kommt
    plt.ylabel('Amount of particles')
    plt.xlim(0.5,4)
    #plt.ylim([-0.5,2])
    plt.legend(loc=1, prop={'size':12})   
    
    save(s_title, ext=format, close=False, verbose=True)
  elif (style == 1):
  
    width = 0.05
    opacity=0.7
    fig, ax = plt.subplots()
    ax.bar(interval,data_2,width, color='red', label=names[2])
    ax.bar(interval,data_3,width, color='black', alpha=opacity, label=names[3])
    plt.title('First plot')
    plt.xlabel('Liquid') #das r hier markiert, dass jetz latex code kommt
    plt.ylabel('Amount of particles')
    plt.xlim(0.5,2)
    #plt.ylim([-0.5,2])
    plt.legend(loc=1, prop={'size':12})   
    
    save(s_title, ext='pdf', close=False, verbose=True)

def test_add_background_image():
    """Test adding background image to a figure."""
    rng = np.random.RandomState(0)
    for ii in range(2):
        f, axs = plt.subplots(1, 2)
        x, y = rng.randn(2, 10)
        im = rng.randn(10, 10)
        axs[0].scatter(x, y)
        axs[1].scatter(y, x)
        for ax in axs:
            ax.set_aspect(1)

        # Background without changing aspect
        if ii == 0:
            ax_im = add_background_image(f, im)
            return
            assert (ax_im.get_aspect() == 'auto')
            for ax in axs:
                assert (ax.get_aspect() == 1)
        else:
            # Background with changing aspect
            ax_im_asp = add_background_image(f, im, set_ratios='auto')
            assert (ax_im_asp.get_aspect() == 'auto')
            for ax in axs:
                assert (ax.get_aspect() == 'auto')
        plt.close('all')

    # Make sure passing None as image returns None
    f, axs = plt.subplots(1, 2)
    assert (add_background_image(f, None) is None)
    plt.close('all')

def execute(model, data, savepath, *args, **kwargs):

    fluence_divisions = [3.3E18, 3.3E19, 3.3E20]
    flux_divisions = [5e11,2e11,1e11]

    fig, ax = plt.subplots(1,3, figsize = (30,10))
    for x in range(len(fluence_divisions)):
        model = model
        data.remove_all_filters()
        data.add_inclusive_filter("fluence n/cm2", '<', fluence_divisions[x])
        l_train = len(data.get_y_data())
        model.fit(data.get_x_data(), np.array(data.get_y_data()).ravel())

        data.remove_all_filters()
        data.add_inclusive_filter("fluence n/cm2", '>=', fluence_divisions[x])
        l_test = len(data.get_y_data())
        Ypredict = model.predict(data.get_x_data())
        RMSE = np.sqrt(mean_squared_error(Ypredict, np.array(data.get_y_data()).ravel()))

        matplotlib.rcParams.update({'font.size': 26})
        ax[x].scatter(data.get_y_data(), Ypredict, color='black', s=10)
        ax[x].plot(ax[x].get_ylim(), ax[x].get_ylim(), ls="--", c=".3")
        ax[x].set_xlabel('Measured ∆sigma (Mpa)')
        ax[x].set_ylabel('Predicted ∆sigma (Mpa)')
        ax[x].set_title('Testing Fluence > {}'.format(fluence_divisions[x]))
        ax[x].text(.1, .88, 'RMSE: {:.3f}'.format(RMSE),fontsize = 30, transform=ax[x].transAxes)
        ax[x].text(.1, .83, 'Train: {}, Test: {}'.format(l_train, l_test), transform=ax[x].transAxes)

    fig.tight_layout()
    plt.subplots_adjust(bottom = .2)
    fig.savefig(savepath.format("fluence_extrapolation"), dpi=150, bbox_inches='tight')
    plt.close()

    fig, ax = plt.subplots(1, 3, figsize=(30, 10))
    for x in range(len(flux_divisions)):
        model = model
        data.remove_all_filters()
        data.add_inclusive_filter("flux n/cm2/s", '>', flux_divisions[x])
        l_train = len(data.get_y_data())
        model.fit(data.get_x_data(), np.array(data.get_y_data()).ravel())

        data.remove_all_filters()
        data.add_inclusive_filter("flux n/cm2/s", '<=', flux_divisions[x])
        l_test = len(data.get_y_data())
        Ypredict = model.predict(data.get_x_data())
        RMSE = np.sqrt(mean_squared_error(Ypredict, np.array(data.get_y_data()).ravel()))

        matplotlib.rcParams.update({'font.size': 26})
        ax[x].scatter(data.get_y_data(), Ypredict, color='black', s=10)
        ax[x].plot(ax[x].get_ylim(), ax[x].get_ylim(), ls="--", c=".3")
        ax[x].set_xlabel('Measured ∆sigma (Mpa)')
        ax[x].set_ylabel('Predicted ∆sigma (Mpa)')
        ax[x].set_title('Testing Flux < {:.0e}'.format(flux_divisions[x]))
        ax[x].text(.1, .88, 'RMSE: {:.3f}'.format(RMSE), fontsize=30, transform=ax[x].transAxes)
        ax[x].text(.1, .83, 'Train: {}, Test: {}'.format(l_train, l_test), transform=ax[x].transAxes)

    fig.tight_layout()
    plt.subplots_adjust(bottom=.2)
    fig.savefig(savepath.format("flux_extrapolation"), dpi=150, bbox_inches='tight')
    plt.close()

def create_plots(params, w_holder, rate_holder):
    print("Creating plots..")
    N_inh_neurons = rate_holder.shape[0]
#    # all spikes
#    plt.figure()
#    plt.plot(SpikeMon.t/ms, SpikeMon.i, '.k', markersize=.1)
#    plt.xlabel("Time (ms)")
#    plt.ylabel("Neuron index")
    
    rate_interval = params["rate_interval"]
    rho_0 = params["rho_0"]
    simtime = params["simtime"]
    dt = params["dt"]

    
    avg_w_stream = np.average(w_holder, axis=0)

    r_idxes = np.random.uniform(rate_holder.shape[0], size=plot_n_rates)
    r_idxes = r_idxes.astype(int)    
    r_stream = rate_holder[r_idxes, :]
    avg_r_stream = np.average(rate_holder, axis=0)

    r_times = np.arange(rate_interval/ms, simtime/ms, rate_interval/ms) * ms  
    w_times = np.arange(0, simtime/ms, rate_interval/ms) * ms

    fig, axes = plt.subplots(2, figsize=(15, 10))
    axes[0].plot(r_times/second, r_stream.T, color="red",
                 alpha=.2, linewidth=.3)
    axes[0].plot(r_times/second, avg_r_stream, color="red", linewidth=2,
                 label="firing_rate")
    axes[0].hlines(rho_0, 0, r_times[-1], linestyles="--")
    axes[0].set_xlim([0, r_times[-1]])
    axes[0].set_xlabel("time [s]")
    axes[0].set_ylabel("firing rate [Hz]")
    axes[0].set_title(str(plot_n_rates) + \
                      " randomly selected firing rates estimated every " + \
                      str(rate_interval))
    axes[1].plot(w_times/second, w_holder.T, color="gray", alpha=.2,
                 linewidth=.3)
    axes[1].plot(w_times/second, avg_w_stream, color="black")
    axes[1].hlines(0, 0, w_times[-1], linestyles="--")
    axes[1].set_ylim([-1, np.amax(w_holder)+10])
    axes[1].set_xlim([0, w_times[-1]])
    axes[1].set_xlabel("time [s]")
    axes[1].set_ylabel("Inh to exc weight")
    axes[1].set_title(str(w_holder.shape[0]) + \
                      " randomly selected inh-to-exc weights")

    
    # firing rate plot as matrix
    rate_vector = rate_holder[:, -1]
    matrix_axis = np.floor(np.sqrt(len(rate_vector)))
    rate_vector = rate_vector[:matrix_axis**2]
    rate_mat = np.reshape(rate_vector, (int(np.sqrt(N_inh_neurons)), -1))
    fig, ax = plt.subplots()
    ax.pcolor(rate_mat, cmap="Reds")
    plt.title("Inh firing rate estimated with counting spikes")
    plt.xticks([]); plt.yticks([]);
    
    plt.show()

def main():
    
    Lx = 1.
    Ly = 1.
    V = Lx*Ly
    N = 101
    dk = (2.*np.pi)/Lx
    
    #-----------------------------------------#
    #-- Generate Gaussian random field -------#
    #-----------------------------------------#
    
    kspace_field = gen_2Dgauss(N, Lx, Ly, model0)
    
    config_field = np.fft.ifft2(kspace_field)*kspace_field.size**0.5
    
    fig, ax = pl.subplots()
    im = ax.imshow(config_field.real, cmap=pl.cm.jet)
    fig.colorbar(im, ax=ax)
    pl.title("field.real, config_space")
    pl.show()
    
    fig, ax = pl.subplots()
    im = ax.imshow(config_field.imag, cmap=pl.cm.jet)
    fig.colorbar(im, ax=ax)
    pl.title("field.imag, config_space")
    pl.show()

def plot_gens(images, rowlabels, losses):
    '''
    From great jupyter notebook by Tim Sainburg:
    http://github.com/timsainb/Tensorflow-MultiGPU-VAE-GAN
    '''
    examples = 8
    fig, ax = plt.subplots(nrows=len(images), ncols=examples, figsize=(18, 8))
    for i in range(examples):
        for j in range(len(images)):
            ax[(j, i)].imshow(create_image(images[j][i]), cmap=plt.cm.gray,
                              interpolation='nearest')
            ax[(j, i)].axis('off')
    title = ''
    for i in rowlabels:
        title += ' {}, '.format(i)
    fig.suptitle('Top to Bottom: {}'.format(title))
    plt.show()
    #fig.savefig(''.join(['imgs/test_',str(epoch).zfill(4),'.png']),dpi=100)
    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 10), linewidth = 4)

    D_plt, = plt.semilogy((losses['discriminator']), linewidth=4, ls='-',
                          color='b', alpha=.5, label='D')
    G_plt, = plt.semilogy((losses['generator']), linewidth=4, ls='-',
                          color='k', alpha=.5, label='G')

    plt.gca()
    leg = plt.legend(handles=[D_plt, G_plt],
                     fontsize=20)
    leg.get_frame().set_alpha(0.5)
    plt.show()

def test_continuous_distortion():

	#load unlensed image from png file
	image_unlensed = plt.imread(os.path.join(dataExtern(),"lensing/lens.png"))[:,:,0]
	pos_original = (np.array(np.where(image_unlensed>0)) * 2.0/image_unlensed.shape[0]) * deg
	pos_original = np.roll(pos_original,1,axis=0)
	pos_original[1] *= -1
	pos_original[1] += 2.0*deg

	#Perform forward ray tracing with grid interpolation to compute the image distortion
	pos_apparent = tracer.shootForward(pos_original,z=2.0,save_intermediate=True)

	#Plot some of the distorted images
	fig,ax = plt.subplots(2,5,figsize=(40,16))
	for i in range(2):
		for j in range(5):
			ax[i,j].scatter(pos_apparent[(5*i + j)*4,0],pos_apparent[(5*i + j)*4,1]) 
			ax[i,j].set_xlabel(r"$x$({0})".format(pos_original.unit.to_string()))
			ax[i,j].set_ylabel(r"$y$({0})".format(pos_original.unit.to_string()))
			ax[i,j].set_title("z={0:.2f}".format(tracer.redshift[(5*i + j)*4]))

	#Save the figure
	fig.tight_layout()
	fig.savefig("lens_distortion_continuous.png")

	#Plot each single frame
	fig,ax = plt.subplots()
	scatterplot = ax.scatter(pos_original[0],pos_original[1])
	ax.set_xlabel(r"$x$({0})".format(pos_original.unit.to_string()))
	ax.set_ylabel(r"$y$({0})".format(pos_original.unit.to_string()))
	
	for n in range(pos_apparent.shape[0]):
		scatterplot.set_offsets(pos_apparent[n].transpose())
		ax.set_title("z={:.2f}".format(tracer.redshift[n]))	
		fig.savefig("distortion{0}.png".format(n))	 

def test_cursor_data():
    from matplotlib.backend_bases import MouseEvent

    fig, ax = plt.subplots()
    im = ax.imshow(np.arange(100).reshape(10, 10), origin='upper')

    x, y = 4, 4
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) == 44

    ax.clear()
    im = ax.imshow(np.arange(100).reshape(10, 10), origin='lower')

    x, y = 4, 4
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) == 44

    fig, ax = plt.subplots()
    im = ax.imshow(np.arange(100).reshape(10, 10), extent=[0, 0.5, 0, 0.5])

    x, y = 0.25, 0.25
    xdisp, ydisp = ax.transData.transform_point([x, y])

    event = MouseEvent('motion_notify_event', fig.canvas, xdisp, ydisp)
    assert im.get_cursor_data(event) == 55

def main2():
    print(__file__ + " start!!")

    dl = 1.0  # course tick
    cx, cy, cyaw, ck = get_straight_course3(dl)

    sp = calc_speed_profile(cx, cy, cyaw, TARGET_SPEED)

    initial_state = State(x=cx[0], y=cy[0], yaw=0.0, v=0.0)

    t, x, y, yaw, v, d, a = do_simulation(
        cx, cy, cyaw, ck, sp, dl, initial_state)

    if show_animation:  # pragma: no cover
        plt.close("all")
        plt.subplots()
        plt.plot(cx, cy, "-r", label="spline")
        plt.plot(x, y, "-g", label="tracking")
        plt.grid(True)
        plt.axis("equal")
        plt.xlabel("x[m]")
        plt.ylabel("y[m]")
        plt.legend()

        plt.subplots()
        plt.plot(t, v, "-r", label="speed")
        plt.grid(True)
        plt.xlabel("Time [s]")
        plt.ylabel("Speed [kmh]")

        plt.show()

def get_axis(figsize, size, aspect, ax):
    import matplotlib as mpl
    import matplotlib.pyplot as plt

    if figsize is not None:
        if ax is not None:
            raise ValueError('cannot provide both `figsize` and '
                             '`ax` arguments')
        if size is not None:
            raise ValueError('cannot provide both `figsize` and '
                             '`size` arguments')
        _, ax = plt.subplots(figsize=figsize)
    elif size is not None:
        if ax is not None:
            raise ValueError('cannot provide both `size` and `ax` arguments')
        if aspect is None:
            width, height = mpl.rcParams['figure.figsize']
            aspect = width / height
        figsize = (size * aspect, size)
        _, ax = plt.subplots(figsize=figsize)
    elif aspect is not None:
        raise ValueError('cannot provide `aspect` argument without `size`')

    if ax is None:
        ax = plt.gca()

    return ax

def plot(args):
    from pickle import load
    from sys import stdin
    import matplotlib.pyplot as plt
    from itertools import starmap
    from numpy.random import seed, random
    from coherent_point_drift.geometry import rigidXform, RMSD

    seed(4) # For color choice
    reference = load(stdin.buffer)

    rmsds = []
    fig, ax = plt.subplots(1, 1)
    for degradation, fit in loadAll(stdin.buffer):
        color = random(3)
        degraded = degrade(reference, *degradation)
        ax.scatter(degraded[:, 0], degraded[:, 1], marker='o', color=color, alpha=0.2)
        fitted = rigidXform(degraded, *fit)
        ax.scatter(fitted[:, 0], fitted[:, 1], marker='+', color=color)
        rmsds.append(RMSD(reference, fitted))
    ax.scatter(reference[:, 0], reference[:, 1], marker='D', color='black')
    ax.set_xticks([])
    ax.set_yticks([])

    if len(rmsds) > 1:
        fig, ax = plt.subplots(1, 1)
        ax.violinplot(rmsds)
        ax.set_ylabel("RMSD")

    plt.show()