Python pylab 模块，draw() 实例源码

def show(self):
        r"""Display an image of the transfer function

        This function loads up matplotlib and displays the current transfer function.

        Parameters
        ----------

        Examples
        --------

        >>> tf = TransferFunction( (-10.0, -5.0) )
        >>> tf.add_gaussian(-9.0, 0.01, 1.0)
        >>> tf.show()
        """
        import pylab
        pylab.clf()
        pylab.plot(self.x, self.y, 'xk-')
        pylab.xlim(*self.x_bounds)
        pylab.ylim(0.0, 1.0)
        pylab.draw()

def plot_2D_heat_map(states,p,labels, inter=False):
    import pylab as pl
    X = np.unique(states[0,:])
    Y = np.unique(states[1,:])
    X_len = len(X)
    Y_len = len(Y)
    Z = np.zeros((X.max()+1,Y.max()+1))
    for i in range(len(p)):
        Z[states[0,i],states[1,i]] = p[i]
    pl.clf()    
    pl.imshow(Z.T, origin='lower')
    pl.xlabel(labels[0])
    pl.ylabel(labels[1])
    if inter== True:
        pl.draw()
    else:
        pl.show()

def plot_2D_contour(states,p,labels,inter=False):
    import pylab as pl

    from pyme.statistics import expectation as EXP
    exp = EXP((states,p)) 
    X = np.unique(states[0,:])
    Y = np.unique(states[1,:])
    X_len = len(X)
    Y_len = len(Y)
    Z = np.zeros((X.max()+1,Y.max()+1))
    for i in range(len(p)):
        Z[states[0,i],states[1,i]] = p[i]

    Z = np.where(Z < 1e-8,0.0,Z)
    pl.clf()
    XX, YY = np.meshgrid(X,Y)   
    pl.contour(range(X.max()+1),range(Y.max()+1),Z.T)
    pl.axhline(y=exp[1])
    pl.axvline(x=exp[0])
    pl.xlabel(labels[0])
    pl.ylabel(labels[1])
    if inter == True:
        pl.draw()
    else:
        pl.show()

def _plot(self):
        # Called from the main thread
        pylab.ion()

        if not getattr(self, 'data_available', False):
            return

        if self.peaks is not None:

            for key in self.sign_peaks:
                for channel in self.peaks[key].keys():
                    self.rates[key][int(channel)] += len(self.peaks[key][channel])

            pylab.scatter(self.positions[0, :], self.positions[1, :], c=self.rates[key])

        pylab.gca().set_title('Buffer %d' %self.counter)
        pylab.draw()
        return

def show(X, C, centroids, keep = False):
    import time
    time.sleep(0.5)
    plt.cla()
    plt.plot(X[C == 0, 0], X[C == 0, 1], '*b',
         X[C == 1, 0], X[C == 1, 1], '*r',
         X[C == 2, 0], X[C == 2, 1], '*g')
    plt.plot(centroids[:,0],centroids[:,1],'*m',markersize=20)
    plt.draw()
    if keep :
        plt.ioff()
        plt.show()

# generate 3 cluster data
# data = np.genfromtxt('data1.csv', delimiter=',')

def wrcontour(dir, var, **kwargs):
    fig = plt.figure()
    rect = [0.1, 0.1, 0.8, 0.8]
    ax = WindroseAxes(fig, rect)
    fig.add_axes(ax)
    ax.contour(dir, var, **kwargs)
    l = ax.legend(axespad=-0.10)
    plt.setp(l.get_texts(), fontsize=8)
    plt.draw()
    plt.show()
    return ax

def wrbox(dir, var, **kwargs):
    fig = plt.figure()
    rect = [0.1, 0.1, 0.8, 0.8]
    ax = WindroseAxes(fig, rect)
    fig.add_axes(ax)
    ax.box(dir, var, **kwargs)
    l = ax.legend(axespad=-0.10)
    plt.setp(l.get_texts(), fontsize=8)
    plt.draw()
    plt.show()
    return ax

def wrbar(dir, var, **kwargs):
    fig = plt.figure()
    rect = [0.1, 0.1, 0.8, 0.8]
    ax = WindroseAxes(fig, rect)
    fig.add_axes(ax)
    ax.bar(dir, var, **kwargs)
    l = ax.legend(axespad=-0.10)
    plt.setp(l.get_texts(), fontsize=8)
    plt.draw()
    plt.show()
    return ax

def add_gaussian(self, location, width, height):
        r"""Add a Gaussian distribution to the transfer function.

        Typically, when rendering isocontours, a Gaussian distribution is the
        easiest way to draw out features.  The spread provides a softness.
        The values are calculated as :math:`f(x) = h \exp{-(x-x_0)^2 / w}`.

        Parameters
        ----------
        location : float
            The centroid of the Gaussian (:math:`x_0` in the above equation.)
        width : float
            The relative width (:math:`w` in the above equation.)
        height : float
            The peak height (:math:`h` in the above equation.)  Note that while
            values greater 1.0 will be accepted, the values of the transmission
            function are clipped at 1.0.

        Examples
        --------

        >>> tf = TransferFunction( (-10.0, -5.0) )
        >>> tf.add_gaussian(-9.0, 0.01, 1.0)
        """
        vals = height * np.exp(-(self.x - location)**2.0/width)
        self.y = np.clip(np.maximum(vals, self.y), 0.0, np.inf)
        self.features.append(('gaussian', "location(x):%3.2g" % location, 
                              "width(x):%3.2g" % width, "height(y):%3.2g" % height))

def add_gaussian(self, location, width, height):
        r"""Add a Gaussian distribution to the transfer function.

        Typically, when rendering isocontours, a Guassian distribution is the
        easiest way to draw out features.  The spread provides a softness.
        The values are calculated as :math:`f(x) = h \exp{-(x-x_0)^2 / w}`.

        Parameters
        ----------
        location : float
            The centroid of the Gaussian (:math:`x_0` in the above equation.)
        width : float
            The relative width (:math:`w` in the above equation.)
        height : list of 4 float
            The peak height (:math:`h` in the above equation.)  Note that while
            values greater 1.0 will be accepted, the values of the transmission
            function are clipped at 1.0.  This must be a list, and it is in the
            order of (red, green, blue, alpha).

        Examples
        --------
        This adds a red spike.

        >>> tf = ColorTransferFunction( (-10.0, -5.0) )
        >>> tf.add_gaussian(-9.0, 0.01, [1.0, 0.0, 0.0, 1.0])
        """
        for tf, v in zip(self.funcs, height):
            tf.add_gaussian(location, width, v)
        self.features.append(('gaussian', "location(x):%3.2g" % location, \
                              "width(x):%3.2g" % width, \
                              "height(y):(%3.2g, %3.2g, %3.2g, %3.2g)" % 
                              (height[0], height[1], height[2], height[3])))

def show_tf(self):
        if self._pylab is None: 
            import pylab
            self._pylab = pylab
        if self._tf_figure is None:
            self._tf_figure = self._pylab.figure(2)
            self.transfer_function.show(ax=self._tf_figure.axes)
        self._pylab.draw()

def draw(self):
        self._pylab.draw()

def show_results_plot(self):
        pl.figure(1)
        pl.plot(self.x, self.y,label = "firing angle = %.1f °"%self.theta)
        pl.draw()
        pl.legend(loc='upper right', shadow=True, fontsize='small')
        print("\ninitial velocity:", _input.initial_velocity, "m/s")
        print("time step:", _input.time_step, "s")
        print("firing angle:", self.theta, "°")
        print("falling fange:%.4f km"%self.x[-1], "\n")

#Users input the initial values

def plot_2D_heat_map(states,p,labels):
    import pylab as pl
    X = np.unique(states[0,:])
    Y = np.unique(states[1,:])
    X_len = len(X)
    Y_len = len(Y)
    Z = np.zeros((X.max()+1,Y.max()+1))
    for i in range(len(p)):
        Z[states[0,i],states[1,i]] = p[i]
    pl.clf()    
    pl.imshow(Z.T, origin='lower')
    pl.xlabel(labels[0])
    pl.ylabel(labels[1])
    pl.draw()
    #pl.show()

def plot_marginals(state_space,p,name,t,labels = False,interactive = False):
    import matplotlib

    import matplotlib.pyplot as pl
    if interactive == True: 
        pl.ion()
    pl.clf()
    pl.suptitle("time: "+ str(t)+" units")
    #print("time : "+ str(t))
    D = state_space.shape[1]

    for i in range(D):
        marg_X = np.unique(state_space[:,i])
        A = np.where(marg_X[:,np.newaxis] == state_space[:,i].T[np.newaxis,:],1,0)
        marg_p = np.dot(A,p)
        pl.subplot(int(D/2)+1,2,i+1)
        pl.plot(marg_X,marg_p)
        pl.yticks(np.linspace(np.amin(marg_p), np.amax(marg_p), num=3))
        pl.axvline(np.sum(marg_X*marg_p),color= 'r')
        pl.axvline(marg_X[np.argmax(marg_p)],color='g')
        if labels == False:
            pl.xlabel("Specie: " + str(i+1))
        else:
            pl.xlabel(labels[i])
    if interactive == True:
        pl.draw()
    else:
        pl.tight_layout()
        pl.show()

def sample_cond(self, X):
        """ActInfHebbianSOM.sample_cond: draw single sample from model conditioned on X"""
        # print("%s.sample_cond X.shape = %s, %d" % (self.__class__.__name__, X.shape, 0))

        # fix the SOMs with learning rate constant 0
        self.filter_e_lr = self.filter_e.map._learning_rate
        self.filter_p_lr = self.filter_p.map._learning_rate
        # print("fit_hebb", self.filter_e.map._learning_rate)
        self.filter_e.map._learning_rate = self.CT(0.0)
        self.filter_p.map._learning_rate = self.CT(0.0)

        e_shape = (np.prod(self.filter_e.map._shape), 1)
        p_shape = (np.prod(self.filter_p.map._shape), 1)

        # activate input network
        self.filter_e.learn(X)

        # pl.plot(self.filter_e.

        # propagate activation via hebbian associative links
        if self.hebblink_use_activity:
            e_ = self.filter_e.activity.reshape((np.prod(self.filter_e.map._shape), 1))
            e_ = (e_ == np.max(e_)) * 1.0
            e2p_activation = np.dot(self.hebblink_filter.T, e_)
            # print("e2p_activation", e2p_activation)
            self.filter_p.activity = np.clip((e2p_activation / (np.sum(e2p_activation) + 1e-9)).reshape(self.filter_p.map._shape), 0, np.inf)
        else:
            e2p_activation = np.dot(self.hebblink_filter.T, self.filter_e.distances(e).flatten().reshape(e_shape))

        # sample the output network with
        sidx = self.filter_p.sample(1)[0]
        e2p_w_p_weights = self.filter_p.neuron(self.filter_p.flat_to_coords(sidx))
        # e2p_w_p_weights = self.filter_p.neuron(self.filter_p.flat_to_coords(np.argmax(self.filter_p.activity)))

        return e2p_w_p_weights.reshape((1, self.odim))
        # ret = np.random.normal(e2p_w_p_weights, self.filter_p.sigmas[sidx] * 0.001, (1, self.odim))
        # ret = np.random.normal(e2p_w_p_weights, 0.01, (1, self.odim))
        # return ret

def show(self, genes, color):
        if self.vis:
            self.get_fig(genes, color)
            pylab.draw()
            pause(0.05)
            pylab.clf()
            self.reset()

def redraw(self):
        global hist_box, HaveTix, marker_size
        bshot.label.set_text(str(self.shot))
        status=call_spec()
        if HaveTix:  # update shot field in either case, only update history if good
            # this updates hist_box if the shot was changed by the other (matplotlib) widgets
            hist_box.set_silent(str(self.shot))
            if status==True:
                hist_box.add_history(str(self.shot))

        print("marker_size", marker_size)
#        if marker_size>0: plot_flucstrucs_for_shot(self.shot, size_factor=marker_size, savefile='')
#        pl.draw()  # what does this do?
        return(status) # False if no data

def remove_last_plot_xy(self, name="xy"):
        """Remove last individual x-y trajectory generated with plot_xy()."""
        import pylab
        l = self._lines[name].pop()
        l.remove()    
        trajectory = self._lines_annotation[name].pop()
        pylab.draw()
        print "remove_last_plot_xy(): Removed trajectory path %(trajectory)r." % vars()

    #
    #------------------------------------------------------------

def remove_last_transition(db):
    import pylab
    l = db._lines['transitions'].pop()
    l.remove()    
    trajectory = db._lines_annotation['transitions'].pop()
    pylab.draw()
    print "Removed trajectory path %(trajectory)r." % vars()

def onpick(self, event):
        for dataind in event.ind:
            site = self.dataObject.sites[dataind]
            if len(str(site)) == 1:
                site = '0'+str(site)
            x1, y1 = self.dataObject.extracts[dataind]
            t = self.dataObject.times[dataind]
            f = glob( '{0}ANGLEsite{1}_extract{2}&{3}_t{4}.png'.format(self.imagespath, int(site), x1, y1, t))
            if len(f) >1 :
                print f
            elif len(f) == 0:
                print '{0}ANGLEsite{1}_extract{2}&{3}_t{4}.png'.format(self.imagespath, int(site), x1, y1, t)
            subprocess.call(['eog', f[0]], stdout=open(os.devnull, 'wb'), stderr=open(os.devnull, 'wb'))
            val = input('Good (1) or not (0) ? If you want the stats, enter a negative number. ')
            self.stats[dataind] = val
            while val < 0 :
                print "True : {0}, False {1}".format(len(self.stats[self.stats==1]), len(self.stats[self.stats==0]))
                val = input('Good or not ? ')
            if val == 0:
                self.ax.scatter(self.X[dataind], self.Y[dataind], c= 'k', marker = 'v')
                print "!!! BUG : ", f[0]
            else:
                self.ax.scatter(self.X[dataind], self.Y[dataind], c = 'r', marker = 'v')

            pylab.draw()
        return True

def plot_rels(data, labels=None, colors=None, outfile="rels", latent=None, alpha=0.8, title=''):
    ns, n = data.shape
    if labels is None:
        labels = map(str, range(n))
    ncol = 5
    nrow = int(np.ceil(float(n * (n - 1) / 2) / ncol))

    fig, axs = pylab.subplots(nrow, ncol)
    fig.set_size_inches(5 * ncol, 5 * nrow)
    pairs = list(combinations(range(n), 2))
    if colors is not None:
        colors = (colors - np.min(colors)) / (np.max(colors) - np.min(colors))

    for ax, pair in zip(axs.flat, pairs):
        diff_x = max(data[:, pair[0]]) - min(data[:, pair[0]])
        diff_y = max(data[:, pair[1]]) - min(data[:, pair[1]])
        ax.set_xlim([min(data[:, pair[0]]) - 0.05 * diff_x, max(data[:, pair[0]]) + 0.05 * diff_x])
        ax.set_ylim([min(data[:, pair[1]]) - 0.05 * diff_y, max(data[:, pair[1]]) + 0.05 * diff_y])
        ax.scatter(data[:, pair[0]], data[:, pair[1]], c=colors, cmap=pylab.get_cmap("jet"),
                       marker='.', alpha=alpha, edgecolors='none', vmin=0, vmax=1)

        ax.set_xlabel(shorten(labels[pair[0]]))
        ax.set_ylabel(shorten(labels[pair[1]]))

    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.scatter(data[:, 0], data[:, 1], marker='.')

    fig.suptitle(title, fontsize=16)
    pylab.rcParams['font.size'] = 12  #6
    # pylab.draw()
    # fig.set_tight_layout(True)
    pylab.tight_layout()
    pylab.subplots_adjust(top=0.95)
    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.set_visible(False)
    filename = outfile + '.png'
    if not os.path.exists(os.path.dirname(filename)):
        os.makedirs(os.path.dirname(filename))
    fig.savefig(outfile + '.png')
    pylab.close('all')
    return True


# Hierarchical graph visualization utilities

def contour(self, dir, var, **kwargs):
        """
        Plot a windrose in linear mode. For each var bins, a line will be
        draw on the axes, a segment between each sector (center to center).
        Each line can be formated (color, width, ...) like with standard plot
        pylab command.

        Mandatory:
        * dir : 1D array - directions the wind blows from, North centred
        * var : 1D array - values of the variable to compute. Typically the wind
        speeds
        Optional:
        * nsector: integer - number of sectors used to compute the windrose
        table. If not set, nsectors=16, then each sector will be 360/16=22.5°,
        and the resulting computed table will be aligned with the cardinals
        points.
        * bins : 1D array or integer- number of bins, or a sequence of
        bins variable. If not set, bins=6, then
            bins=linspace(min(var), max(var), 6)
        * blowto : bool. If True, the windrose will be pi rotated,
        to show where the wind blow to (usefull for pollutant rose).
        * colors : string or tuple - one string color ('k' or 'black'), in this
        case all bins will be plotted in this color; a tuple of matplotlib
        color args (string, float, rgb, etc), different levels will be plotted
        in different colors in the order specified.
        * cmap : a cm Colormap instance from matplotlib.cm.
          - if cmap == None and colors == None, a default Colormap is used.

        others kwargs : see help(pylab.plot)

        """

        bins, nbins, nsector, colors, angles, kwargs = self._init_plot(dir, var,
                                                                       **kwargs)

        #closing lines
        angles = np.hstack((angles, angles[-1]-2*np.pi/nsector))
        vals = np.hstack((self._info['table'],
                         np.reshape(self._info['table'][:,0],
                                   (self._info['table'].shape[0], 1))))

        offset = 0
        for i in range(nbins):
            val = vals[i,:] + offset
            offset += vals[i, :]
            zorder = ZBASE + nbins - i
            patch = self.plot(angles, val, color=colors[i], zorder=zorder,
                              **kwargs)
            self.patches_list.extend(patch)
        self._update()

def contourf(self, dir, var, **kwargs):
        """
        Plot a windrose in filled mode. For each var bins, a line will be
        draw on the axes, a segment between each sector (center to center).
        Each line can be formated (color, width, ...) like with standard plot
        pylab command.

        Mandatory:
        * dir : 1D array - directions the wind blows from, North centred
        * var : 1D array - values of the variable to compute. Typically the wind
        speeds
        Optional:
        * nsector: integer - number of sectors used to compute the windrose
        table. If not set, nsectors=16, then each sector will be 360/16=22.5°,
        and the resulting computed table will be aligned with the cardinals
        points.
        * bins : 1D array or integer- number of bins, or a sequence of
        bins variable. If not set, bins=6, then
            bins=linspace(min(var), max(var), 6)
        * blowto : bool. If True, the windrose will be pi rotated,
        to show where the wind blow to (usefull for pollutant rose).
        * colors : string or tuple - one string color ('k' or 'black'), in this
        case all bins will be plotted in this color; a tuple of matplotlib
        color args (string, float, rgb, etc), different levels will be plotted
        in different colors in the order specified.
        * cmap : a cm Colormap instance from matplotlib.cm.
          - if cmap == None and colors == None, a default Colormap is used.

        others kwargs : see help(pylab.plot)

        """

        bins, nbins, nsector, colors, angles, kwargs = self._init_plot(dir, var,
                                                                       **kwargs)
        null = kwargs.pop('facecolor', None)
        null = kwargs.pop('edgecolor', None)

        #closing lines
        angles = np.hstack((angles, angles[-1]-2*np.pi/nsector))
        vals = np.hstack((self._info['table'],
                         np.reshape(self._info['table'][:,0],
                                   (self._info['table'].shape[0], 1))))
        offset = 0
        for i in range(nbins):
            val = vals[i,:] + offset
            offset += vals[i, :]
            zorder = ZBASE + nbins - i
            xs, ys = poly_between(angles, 0, val)
            patch = self.fill(xs, ys, facecolor=colors[i],
                              edgecolor=colors[i], zorder=zorder, **kwargs)
            self.patches_list.extend(patch)

def bar(self, dir, var, **kwargs):
        """
        Plot a windrose in bar mode. For each var bins and for each sector,
        a colored bar will be draw on the axes.

        Mandatory:
        * dir : 1D array - directions the wind blows from, North centred
        * var : 1D array - values of the variable to compute. Typically the wind
        speeds
        Optional:
        * nsector: integer - number of sectors used to compute the windrose
        table. If not set, nsectors=16, then each sector will be 360/16=22.5°,
        and the resulting computed table will be aligned with the cardinals
        points.
        * bins : 1D array or integer- number of bins, or a sequence of
        bins variable. If not set, bins=6 between min(var) and max(var).
        * blowto : bool. If True, the windrose will be pi rotated,
        to show where the wind blow to (usefull for pollutant rose).
        * colors : string or tuple - one string color ('k' or 'black'), in this
        case all bins will be plotted in this color; a tuple of matplotlib
        color args (string, float, rgb, etc), different levels will be plotted
        in different colors in the order specified.
        * cmap : a cm Colormap instance from matplotlib.cm.
          - if cmap == None and colors == None, a default Colormap is used.
        edgecolor : string - The string color each edge bar will be plotted.
        Default : no edgecolor
        * opening : float - between 0.0 and 1.0, to control the space between
        each sector (1.0 for no space)

        """

        bins, nbins, nsector, colors, angles, kwargs = self._init_plot(dir, var,
                                                                       **kwargs)
        null = kwargs.pop('facecolor', None)
        edgecolor = kwargs.pop('edgecolor', None)
        if edgecolor is not None:
            if not isinstance(edgecolor, str):
                raise ValueError('edgecolor must be a string color')
        opening = kwargs.pop('opening', None)
        if opening is None:
            opening = 0.8
        dtheta = 2*np.pi/nsector
        opening = dtheta*opening

        for j in range(nsector):
            offset = 0
            for i in range(nbins):
                if i > 0:
                    offset += self._info['table'][i-1, j]
                val = self._info['table'][i, j]
                zorder = ZBASE + nbins - i
                patch = Rectangle((angles[j]-opening/2, offset), opening, val,
                    facecolor=colors[i], edgecolor=edgecolor, zorder=zorder,
                    **kwargs)
                self.add_patch(patch)
                if j == 0:
                    self.patches_list.append(patch)
        self._update()

def box(self, dir, var, **kwargs):
        """
        Plot a windrose in proportional bar mode. For each var bins and for each
        sector, a colored bar will be draw on the axes.

        Mandatory:
        * dir : 1D array - directions the wind blows from, North centred
        * var : 1D array - values of the variable to compute. Typically the wind
        speeds
        Optional:
        * nsector: integer - number of sectors used to compute the windrose
        table. If not set, nsectors=16, then each sector will be 360/16=22.5°,
        and the resulting computed table will be aligned with the cardinals
        points.
        * bins : 1D array or integer- number of bins, or a sequence of
        bins variable. If not set, bins=6 between min(var) and max(var).
        * blowto : bool. If True, the windrose will be pi rotated,
        to show where the wind blow to (usefull for pollutant rose).
        * colors : string or tuple - one string color ('k' or 'black'), in this
        case all bins will be plotted in this color; a tuple of matplotlib
        color args (string, float, rgb, etc), different levels will be plotted
        in different colors in the order specified.
        * cmap : a cm Colormap instance from matplotlib.cm.
          - if cmap == None and colors == None, a default Colormap is used.
        edgecolor : string - The string color each edge bar will be plotted.
        Default : no edgecolor

        """

        bins, nbins, nsector, colors, angles, kwargs = self._init_plot(dir, var,
                                                                       **kwargs)
        null = kwargs.pop('facecolor', None)
        edgecolor = kwargs.pop('edgecolor', None)
        if edgecolor is not None:
            if not isinstance(edgecolor, str):
                raise ValueError('edgecolor must be a string color')
        opening = np.linspace(0.0, np.pi/16, nbins)

        for j in range(nsector):
            offset = 0
            for i in range(nbins):
                if i > 0:
                    offset += self._info['table'][i-1, j]
                val = self._info['table'][i, j]
                zorder = ZBASE + nbins - i
                patch = Rectangle((angles[j]-opening[i]/2, offset), opening[i],
                    val, facecolor=colors[i], edgecolor=edgecolor,
                    zorder=zorder, **kwargs)
                self.add_patch(patch)
                if j == 0:
                    self.patches_list.append(patch)
        self._update()

def allplot(xb,yb,bins=30,fig=1,xlabel='x',ylabel='y'):
    """
Input:
X,Y : objects referring to the variables produced by PyMC that you want
to analyze. Example: X=M.theta, Y=M.slope.

Inherited from Tommy LE BLANC's code at astroplotlib|STSCI.
    """
    #X,Y=xb.trace(),yb.trace()
    X,Y=xb,yb

    #pylab.rcParams.update({'font.size': fontsize})
    fig=pylab.figure(fig)
    pylab.clf()

    gs = pylab.GridSpec(2, 2, width_ratios=[3,1], height_ratios=[1,3], wspace=0.07, hspace=0.07)
    scat=pylab.subplot(gs[2])
    histx=pylab.subplot(gs[0], sharex=scat)
    histy=pylab.subplot(gs[3], sharey=scat)
    #scat=fig.add_subplot(2,2,3)
    #histx=fig.add_subplot(2,2,1, sharex=scat)
    #histy=fig.add_subplot(2,2,4, sharey=scat)

    # Scatter plot
    scat.plot(X, Y,linestyle='none', marker='o', color='green', mec='green',alpha=.2, zorder=-99)

    gkde = scipy.stats.gaussian_kde([X, Y])
    x,y = numpy.mgrid[X.min():X.max():(X.max()-X.min())/25.,Y.min():Y.max():(Y.max()-Y.min())/25.]
    z = numpy.array(gkde.evaluate([x.flatten(), y.flatten()])).reshape(x.shape)
    scat.contour(x, y, z, linewidths=2)
    scat.set_xlabel(xlabel)
    scat.set_ylabel(ylabel)

    # X-axis histogram
    histx.hist(X, bins, histtype='stepfilled')
    pylab.setp(histx.get_xticklabels(), visible=False)  # no X label
    #histx.xaxis.set_major_formatter(pylab.NullFormatter()) # no X label

    # Y-axis histogram
    histy.hist(Y, bins, histtype='stepfilled', orientation='horizontal')
    pylab.setp(histy.get_yticklabels(), visible=False)  # no Y label
    #histy.yaxis.set_major_formatter(pylab.NullFormatter()) # no Y label

    #pylab.minorticks_on()
    #pylab.subplots_adjust(hspace=0.1)
    #pylab.subplots_adjust(wspace=0.1)
    pylab.draw()
    pylab.show()

def plot_rels(data, labels=None, colors=None, outfile="rels", latent=None, alpha=0.8):
    ns, n = data.shape
    if labels is None:
        labels = list(map(str, range(n)))
    ncol = 5
    # ncol = 4
    nrow = int(np.ceil(float(n * (n - 1) / 2) / ncol))
    #nrow=1
    #pylab.rcParams.update({'figure.autolayout': True})
    fig, axs = pylab.subplots(nrow, ncol)
    fig.set_size_inches(5 * ncol, 5 * nrow)
    #fig.set_canvas(pylab.gcf().canvas)
    pairs = list(combinations(range(n), 2))  #[:4]
    pairs = sorted(pairs, key=lambda q: q[0]**2+q[1]**2)  # Puts stronger relationships first
    if colors is not None:
        colors = (colors - np.min(colors)) / (np.max(colors) - np.min(colors)).clip(1e-7)

    for ax, pair in zip(axs.flat, pairs):
        if latent is None:
            ax.scatter(data[:, pair[0]], data[:, pair[1]], marker='.', edgecolors='none', alpha=alpha)
        else:
            # cs = 'rgbcmykrgbcmyk'
            markers = 'x+.o,<>^^<>,+x.'
            for j, ind in enumerate(np.unique(latent)):
                inds = (latent == ind)
                ax.scatter(data[inds, pair[0]], data[inds, pair[1]], c=colors[inds], cmap=pylab.get_cmap("jet"),
                           marker=markers[j], alpha=0.5, edgecolors='none', vmin=0, vmax=1)

        ax.set_xlabel(shorten(labels[pair[0]]))
        ax.set_ylabel(shorten(labels[pair[1]]))

    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.scatter(data[:, 0], data[:, 1], marker='.')

    pylab.rcParams['font.size'] = 12  #6
    pylab.draw()
    #fig.set_tight_layout(True)
    fig.tight_layout()
    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.set_visible(False)
    filename = outfile + '.png'
    if not os.path.exists(os.path.dirname(filename)):
        os.makedirs(os.path.dirname(filename))
    fig.savefig(outfile + '.png')  #df')
    pylab.close('all')
    return True

def sample_cond(self, X):
        """smpHebbianSOM.sample_cond: draw single sample from model conditioned on X"""
        # print("%s.sample_cond X.shape = %s, %d" % (self.__class__.__name__, X.shape, 0))

        # fix the SOMs with learning rate constant 0
        self.filter_e_lr = self.filter_e.map._learning_rate
        self.filter_p_lr = self.filter_p.map._learning_rate
        # print("fit_hebb", self.filter_e.map._learning_rate)
        self.filter_e.map._learning_rate = self.CT(0.0)
        self.filter_p.map._learning_rate = self.CT(0.0)

        e_shape = (np.prod(self.filter_e.map._shape), 1)
        p_shape = (np.prod(self.filter_p.map._shape), 1)

        # activate input network
        self.filter_e.learn(X)

        # pl.plot(self.filter_e.

        # propagate activation via hebbian associative links
        if self.hebblink_use_activity:
            e_ = self.filter_e.activity.reshape((np.prod(self.filter_e.map._shape), 1))
            e_ = (e_ == np.max(e_)) * 1.0
            e2p_activation = np.dot(self.hebblink_filter.T, e_)
            # print("e2p_activation", e2p_activation)
            self.filter_p.activity = np.clip((e2p_activation / (np.sum(e2p_activation) + 1e-9)).reshape(self.filter_p.map._shape), 0, np.inf)
        else:
            e2p_activation = np.dot(self.hebblink_filter.T, self.filter_e.distances(e).flatten().reshape(e_shape))

        # sample the output network with
        sidxs = self.filter_p.sample(100)
        # print("sidxs", stats.mode(sidxs)[0], sidxs)
        # sidx = self.filter_p.sample(1)[0]
        # find the mode (most frequent realization) of distribution
        sidx = stats.mode(sidxs)[0][0]
        e2p_w_p_weights = self.filter_p.neuron(self.filter_p.flat_to_coords(sidx))
        # e2p_w_p_weights = self.filter_p.neuron(self.filter_p.flat_to_coords(np.argmax(self.filter_p.activity)))

        # ret = np.random.normal(e2p_w_p_weights, self.filter_p.sigmas[sidx], (1, self.odim))
        ret = np.random.normal(e2p_w_p_weights, np.sqrt(self.filter_p.sigmas[sidx]), (1, self.odim))
        # ret = np.random.normal(e2p_w_p_weights, 0.01, (1, self.odim))
        # print("hebbsom sample", sidx, e2p_w_p_weights) # , sidxs) # , self.filter_p.sigmas[sidx])
        # ret = e2p_w_p_weights.reshape((1, self.odim))
        return ret

def add_markers(self, ax=None, where=0.0, orientation='horizontal',
                    jitter=0, **kwargs):

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

        # draw the positions
        if 'marker' not in kwargs:
            if orientation == 'horizontal':
                kwargs['marker'] = '|'
            else:
                kwargs['marker'] = '_'

        if ('facecolor' not in kwargs.keys()) | ('fc' not in kwargs.keys()) | \
           ('markerfacecolor' not in kwargs.keys()) | ('mfc' not in kwargs.keys()):
            kwargs['markerfacecolor'] = 'None'
        if ('edgecolor' not in kwargs.keys()) | ('ec' not in kwargs.keys()) | \
           ('markeredgecolor' not in kwargs.keys()) | ('mec' not in kwargs.keys()):
            kwargs['markeredgecolor'] = 'k'
        if ('linestyle' not in kwargs.keys()) | ('ls' not in kwargs.keys()):
            kwargs['linestyle'] = 'None'
        if ('size' not in kwargs.keys()) | ('markersize' not in kwargs.keys()):
            kwargs['markersize'] = 3

        if orientation == 'horizontal':
            # Draw the lines
            if jitter > 0:
                pos = np.random.uniform(low=float(where - jitter),
                                        high=float(where + jitter),
                                        size=len(self.x))
                ax.plot(self.x, pos, **kwargs)
            else:
                ax.plot(self.x, float(where) * np.ones(len(self.x)), **kwargs)

            plt.draw_if_interactive()

        elif orientation == 'vertical':
            # Draw the lines
            if jitter > 0.:
                pos = np.random.uniform(low=float(where - jitter),
                                        high=float(where + jitter),
                                        size=len(self.x))
                ax.plot(pos, self.x, **kwargs)
            else:
                ax.plot(float(where) * np.ones(len(self.x)), self.x, marker='_',
                        **kwargs)

        plt.draw_if_interactive()

def rh_model_sweep_generate_input_grid_a(self):
        sweepsteps = 11 # 21
        # extero config
        dim_axes = [np.linspace(self.environment.conf.s_mins[i], self.environment.conf.s_maxs[i], sweepsteps) for i in range(self.environment.conf.s_ndims)]
        # dim_axes = [np.linspace(self.environment.conf.s_mins[i], self.environment.conf.s_maxs[i], sweepsteps) for i in range(self.mdl.idim)]
        print "rh_model_sweep_generate_input_grid: s_ndims = %d, dim_axes = %s" % (self.environment.conf.s_ndims, dim_axes,)
        full_axes = np.meshgrid(*tuple(dim_axes), indexing='ij')
        print "rh_model_sweep_generate_input_grid: full_axes = %s, %s" % (len(full_axes), full_axes,)

        for i in range(len(full_axes)):
            print i, full_axes[i].shape
            print i, full_axes[i].flatten()

        # return proxy
        error_grid = np.vstack([full_axes[i].flatten() for i in range(len(full_axes))])
        print "error_grid", error_grid.shape

        # draw state / goal configurations
        X_accum = []
        states = np.linspace(-1, 1, sweepsteps)
        # for state in range(1): # sweepsteps):
        for state in states:
            # randomize initial position
            # self.M_prop_pred = self.environment.compute_motor_command(np.random.uniform(-1.0, 1.0, (1, self.odim)))
            # draw random goal and keep it fixed
            # self.goal_prop = self.environment.compute_motor_command(np.random.uniform(-1.0, 1.0, (1, self.odim)))
            self.goal_prop = self.environment.compute_motor_command(np.ones((1, self.odim)) * state)
            # self.goal_prop = np.random.uniform(self.environment.conf.m_mins, self.environment.conf.m_maxs, (1, self.odim))
            GOALS = np.repeat(self.goal_prop, error_grid.shape[1], axis = 0) # as many goals as error components
            # FIXME: hacks for M1/M2
            if self.mdl.idim == 3:
                X = GOALS
            elif self.mdl.idim == 6:
                X = np.hstack((GOALS, error_grid.T))
            else:
                X = np.hstack((GOALS, error_grid.T))
            X_accum.append(X)

        X_accum = np.array(X_accum)

        # don't need this?
        # X_accum = X_accum.reshape((X_accum.shape[0] * X_accum.shape[1], X_accum.shape[2]))

        print "X_accum.shape = %s, mdl.idim = %d, mdl.odim = %d" % (X_accum.shape, self.mdl.idim, self.mdl.odim)
        # print X_accum
        X = X_accum
        # X's and pred's indices now mean: slowest: goal, e1, e2, fastest: e3
        self.X_model_sweep = X.copy()
        return X
        # print "self.X_model_sweep.shape", self.X_model_sweep.shape

def check_matplotlib_backends():
    # from http://stackoverflow.com/questions/5091993/list-of-all-available-matplotlib-backends
    # get the directory where the backends live
    backends_dir = os.path.dirname(matplotlib.backends.__file__)

    # filter all files in that directory to identify all files which provide a backend
    backend_fnames = filter(is_backend_module, os.listdir(backends_dir))

    backends = [backend_fname_formatter(fname) for fname in backend_fnames]

    print("supported backends: \t" + str(backends))

    # validate backends
    backends_valid = []
    for b in backends:
        try:
            plt.switch_backend(b)
            backends_valid += [b]
        except:
            continue

    print("valid backends: \t" + str(backends_valid))


    # try backends performance
    for b in backends_valid:

        pylab.ion()
        try:
            plt.switch_backend(b)


            pylab.clf()
            tstart = time.time()               # for profiling
            x = range(0,2*pylab.pi,0.01)            # x-array
            line, = pylab.plot(x,pylab.sin(x))
            for i in range(1,200):
                line.set_ydata(pylab.sin(x+i/10.0))  # update the data
                pylab.draw()                         # redraw the canvas

            print(b + ' FPS: \t' , 200/(time.time()-tstart))
            pylab.ioff()

        except:
            print(b + " error :(")