def show_minimum_bounding_ellipsoids(Xr, with_sampled = False):
    clustellip = Clustered_Sampler(active_samples = Xr,likelihood_constraint=0.0, enlargement= 1.0,no=1)
    ellipsoids = clustellip.ellipsoid_set
    X0 = [i[0] for i in Xr]
    Y0 = [i[1] for i in Xr] 
    plt.figure()
    ax = plt.gca()
    points = []
    #print len(ellipsoids)
    for i in range(len(ellipsoids)):
        if ellipsoids[i]!=None:
            a, b = np.linalg.eig(ellipsoids[i].covariance_matrix)
            c = np.dot(b, np.diag(np.sqrt(a)))
            width = np.sqrt(np.sum(c[:,1]**2)) * 2.
            height = np.sqrt(np.sum(c[:,0]**2)) * 2.
            angle = math.atan(c[1,1] / c[0,1]) * 180./math.pi
            points = ellipsoids[i].sample(n_points = 50)
            ellipse = Ellipse(ellipsoids[i].centroid, width, height, angle)
            ellipse.set_facecolor('None')
            ax.add_patch(ellipse)
    X = [i[0] for i in points]
    Y = [i[1] for i in points]
    if with_sampled == True:
        plt.plot(X, Y, 'bo')
    plt.plot(X0, Y0, 'ro')
    plt.show()

def test_polar_coord_annotations():
    # You can also use polar notation on a catesian axes.  Here the
    # native coordinate system ('data') is cartesian, so you need to
    # specify the xycoords and textcoords as 'polar' if you want to
    # use (theta, radius)
    from matplotlib.patches import Ellipse
    el = Ellipse((0,0), 10, 20, facecolor='r', alpha=0.5)

    fig = plt.figure()
    ax = fig.add_subplot( 111, aspect='equal' )

    ax.add_artist( el )
    el.set_clip_box( ax.bbox )

    ax.annotate('the top',
                xy=(np.pi/2., 10.),      # theta, radius
                xytext=(np.pi/3, 20.),   # theta, radius
                xycoords='polar',
                textcoords='polar',
                arrowprops=dict(facecolor='black', shrink=0.05),
                horizontalalignment='left',
                verticalalignment='baseline',
                clip_on=True, # clip to the axes bounding box
                )

    ax.set_xlim( -20, 20 )
    ax.set_ylim( -20, 20 )

def graficaGrilla(dataGrilla,name,framesNumber,colour,xPixels,yPixels):	
	from matplotlib.patches import Ellipse
	from pylab import figure, show, savefig

	fig = figure()
	ax = fig.add_subplot(111, aspect='equal')
	# Each row of dataGrilla contains 
	# N == "framesNumbers" , signal
	# A radius of the RF ellipse
	# B radius of the RF ellipse
	# Angle of the RF ellipse
	# X coordinate of the RF ellipse
	# Y coordinate of the RF ellipse

	ax = fig.add_subplot(111, aspect='equal')
	for unit in range(dataGrilla.shape[0]):
		eWidth = dataGrilla[unit][framesNumber-1+1]
		eHeight = dataGrilla[unit][framesNumber-1+2]
		eAngle = dataGrilla[unit][framesNumber-1+3]
		eXY = [dataGrilla[unit][framesNumber-1+4],  dataGrilla[unit][framesNumber-1+5]]
		e = Ellipse(xy=eXY, width=eWidth, height=eHeight, angle=eAngle)
		ax.add_artist(e)
		e.set_alpha(0.2)
		e.set_facecolor(colour)
	
	ax.set_xlim(0, xPixels)
	ax.set_ylim(0, yPixels)
	savefig(name, dpi=None, bbox_inches='tight')
	return 0

 def PlotLens(self,CenterXY, thickness, LensHeight, rotation, alpha=0.5):
     # Definition to plot a lens.
     lens = Ellipse(xy=CenterXY, width=thickness, height=LensHeight, angle=-rotation)
     self.ax.add_artist(lens)
     lens.set_clip_box(self.ax.bbox)
     lens.set_alpha(alpha)
     return True

def plot(data, mus=None, sigmas=None, colors='black', figname='fig.png'):
    # Create figure
    fig = plt.figure()

    # Plot data
    x = data[:, 0]
    y = data[:, 1]
    plt.scatter(x, y, 24, c=colors)

    # Plot cluster centers
    if mus is not None:
        x = [float(m[0]) for m in mus]
        y = [float(m[1]) for m in mus]
        plt.scatter(x, y, 99, c='red')

    # Plot ellipses for each cluster
    if sigmas is not None:
        for sig_ix in range(K):
            ax = fig.gca()
            cov = np.array(sigmas[sig_ix])
            lam, v = np.linalg.eig(cov)
            lam = np.sqrt(lam)
            ell = Ellipse(xy=(x[sig_ix], y[sig_ix]),
                          width=lam[0]*4, height=lam[1]*4,
                          angle=np.rad2deg(np.arccos(v[0, 0])),
                          color='blue')
            ell.set_facecolor('none')
            ax.add_artist(ell)

    # Save figure
    fig.savefig(figname)

def drawDef(dfeat,dy,dx,mindef=0.001,distr="father"):
    """
        auxiliary funtion to draw recursive levels of deformation
    """
    from matplotlib.patches import Ellipse
    pylab.ioff()
    if distr=="father":
        py=[0,0,2,2];px=[0,2,0,2]
    if distr=="child":
        py=[0,1,1,2];px=[1,2,0,1]
    ordy=[0,0,1,1];ordx=[0,1,0,1]
    x1=-0.5+dx;x2=2.5+dx
    y1=-0.5+dy;y2=2.5+dy
    if distr=="father":       
        pylab.fill([x1,x1,x2,x2,x1],[y1,y2,y2,y1,y1],"r", alpha=0.15, edgecolor="b",lw=1)    
    for l in range(len(py)):
        aux=dfeat[ordy[l],ordx[l],:].clip(-1,-mindef)
        wh=numpy.exp(-mindef/aux[0])/numpy.exp(1);hh=numpy.exp(-mindef/aux[1])/numpy.exp(1)
        e=Ellipse(xy=[(px[l]+dx),(py[l]+dy)], width=wh, height=hh, alpha=0.35)
        x1=-0.75+dx+px[l];x2=0.75+dx+px[l]
        y1=-0.76+dy+py[l];y2=0.75+dy+py[l]
        col=numpy.array([wh*hh]*3).clip(0,1)
        if distr=="father":
            col[0]=0       
        e.set_facecolor(col)
        pylab.gca().add_artist(e)
        if distr=="father":       
            pylab.fill([x1,x1,x2,x2,x1],[y1,y2,y2,y1,y1],"b", alpha=0.15, edgecolor="b",lw=1)            

    def error_ellipse(self, i, j, nstd=1, space_factor=5.0, clr="b", alpha=0.5, lw=1.5):
        """
        return the plot of the error ellipse from the covariance matrix
        use ideas from error_ellipse.m from
        http://www.mathworks.com/matlabcentral/fileexchange/4705-error-ellipse
        ( the version used here was taken from an earlier version, it looks to have been updated there to a new version.)

        * (i,j) specify the ith and jth parameters to be used and
        * nstd specifies the number of standard deviations to plot, the default is 1 sigma.
        * space_factor specifies the number that divides the width/height, the result of which is then added as extra space to the figure
        """

        def eigsorted(cov):
            vals, vecs = np.linalg.eigh(cov)
            order = vals.argsort()[::1]
            return vals[order], vecs[:, order]

        self.marginalize(param_list=[i, j])
        vals, vecs = eigsorted(self.marginal_covariance_matrix)
        theta = np.degrees(np.arctan2(*vecs[:, 0][::-1]))
        # print theta
        width, height = 2 * nstd * np.sqrt(vals)
        xypos = [self.parameter_values[i], self.parameter_values[j]]
        ellip = Ellipse(xy=xypos, width=width, height=height, angle=theta, color=clr, alpha=alpha)
        # ellip.set_facecolor("white")
        ellip.set_linewidth(lw)
        ellip_vertices = ellip.get_verts()
        xl = [ellip_vertices[k][0] for k in range(len(ellip_vertices))]
        yl = [ellip_vertices[k][1] for k in range(len(ellip_vertices))]
        dx = (max(xl) - min(xl)) / space_factor
        dy = (max(yl) - min(yl)) / space_factor
        xyaxes = [min(xl) - dx, max(xl) + dx, min(yl) - dy, max(yl) + dy]
        return ellip, xyaxes

def visualizer(data, x, y, lengths):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    fig.suptitle("Clustered Results")
    colors = plt.get_cmap("gist_rainbow")
    norm = col.Normalize(vmin = 0, vmax = len(lengths))
    color_map = cm.ScalarMappable(cmap = colors, norm = norm)
    index = 0
    for i in range(len(lengths)):
        plt.scatter(data[index:index + lengths[i], x], data[index:index + lengths[i], y], c = color_map.to_rgba(i))
        if lengths[i] != 1:
            cov = numpy.cov(data[index:index + lengths[i], x], data[index:index + lengths[i], y], ddof = 0)
        else:
            cov = numpy.asarray([[0, 0], [0,0]])
        values, vectors = numpy.linalg.eig(cov)
        angle = numpy.arctan((vectors[0,1] - vectors[1,1])/(vectors[0,0] - vectors[1,0]))
        ellipse = Ellipse(xy = (numpy.mean(data[index:index + lengths[i], x]), numpy.mean(data[index:index + lengths[i], y])), width = 2 * 2 * numpy.sqrt(values[0]),
                          height = 2 * 2 * numpy.sqrt(values[1]), angle = numpy.rad2deg(angle))
        ellipse.set_edgecolor(color_map.to_rgba(i))
        ellipse.set_facecolor("none")
        ax.add_artist(ellipse)
        plt.scatter(numpy.mean(data[index:index + lengths[i], x]), numpy.mean(data[index:index + lengths[i], y]), c = color_map.to_rgba(i), marker = 'x')
        index += lengths[i]
    plt.show()
    plt.close()

def motorplot_group_scatter(fig, ax, data, hz, title, markersz = 2, ellipseht=0.2, ylim=[]):
  for d in data:
    x = d.get('i')
    y = d.get('hz')
    label = d.get('label')
    ax.scatter(x, y, marker='o', s=markersz, label=label)

  xi = ax.get_xlim()[0]
  xf = ax.get_xlim()[1]
  e = Ellipse(((xi + xf)/2,hz),xf-xi, ellipseht)
  e.set_alpha(0.2)
  ax.add_artist(e)
  ax.set_title(title)
  ax.set_xlabel('Sample')
  ax.set_ylabel('Frequency (Hz)')

  box = ax.get_position()
  ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
  plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize='large', markerscale=2)

  ax.title.set_fontsize(16)
  ax.xaxis.label.set_fontsize(16)
  for item in ax.get_xticklabels():
    item.set_fontsize(12)
  ax.yaxis.label.set_fontsize(16)
  for item in ax.get_yticklabels():
    item.set_fontsize(12)
  #ax.yaxis.get_major_formatter().set_powerlimits((0, 1))
  if ylim:
    ax.set_ylim(ylim)

def motorplot_scatter(fig, ax, session, hz, title, markersz = 2, ylim=[], color='r'):
  x = session.get('i')
  y = session.get('hz')
  label = session.get('label')
  #title = 'Actuator Frequency - ' + label 
  ax.scatter(x, y, marker='o', s=markersz, label=label, color=color)
  xi = ax.get_xlim()[0]
  xf = ax.get_xlim()[1]
  e = Ellipse(((xi + xf)/2,hz),xf-xi, 0.1)
  e.set_alpha(0.2)
  ax.add_artist(e)
  ax.set_title(title)
  ax.set_xlabel('Sample')
  ax.set_ylabel('Frequency (Hz)')

  ax.title.set_fontsize(16)
  ax.xaxis.label.set_fontsize(16)
  for item in ax.get_xticklabels():
    item.set_fontsize(12)
  ax.yaxis.label.set_fontsize(16)
  for item in ax.get_yticklabels():
    item.set_fontsize(12)
  #ax.yaxis.get_major_formatter().set_powerlimits((0, 1))
  if ylim:
    ax.set_ylim(ylim)

def plot_ellipse(x, y, ax, color):
	(xy, width, height, angle) = get_ellipse(x, y)
	e = Ellipse(xy, width, height, angle)
	ax.add_artist(e)
	#e.set_clip_box(ax.bbox)
	e.set_alpha(0.2)
	e.set_facecolor(color)

 def plot_galaxies(self, ax, gals=None):
     if gals == None:
         gals = self.galaxies
     for gal in gals:
         e = Ellipse((gal.x, gal.y), gal.a(), gal.b(), angle=gal.theta, linewidth=2, fill=True)
         ax.add_artist(e)
         e.set_clip_box(ax.bbox)

    def calculate_ellipse(center_x, center_y, covariance_matrix, n_std_dev=3):
        values, vectors = np.linalg.eigh(covariance_matrix)
        order = values.argsort()[::-1]
        values = values[order]
        vectors = vectors[:, order]

        theta = np.degrees(np.arctan2(*vectors[:, 0][::-1]))

        # make all angles positive
        if theta < 0:
            theta += 360

        # Width and height are "full" widths, not radius
        width, height = 2 * n_std_dev * np.sqrt(values)

        ellipse = Ellipse(
            xy=[center_x, center_y],
            width=width,
            height=height,
            angle=theta
        )

        ellipse.set_alpha(0.3)
        ellipse.set_facecolor((1, 0, 0))

        return ellipse

def plotCovariance(center, matrix, color):
    eigvalue, eigvector = lng.eig(matrix)
    theta = atan2(eigvector[0][1], eigvector[0][0])
    ellipse = Ellipse([center[0], center[1]], 2*pow(eigvalue[0], 0.5), 2*pow(eigvalue[1], 0.5), center[2])
    ellipse.set_facecolor(color)
    fig = plt.gcf()
    fig.gca().add_artist(ellipse)

def photogeneration(loc):
    arl = 0.10
    sep = 0.05
    arprop = dict(mutation_scale=10,
                  lw=1, arrowstyle="-|>", 
                  zorder=20)
    x, y, z = loc
    
    #electron
    electron, = ax.plot3D([x-sep/2], [y], [z],
                          c='c', marker='$-$', markersize=7, ls='none')
    v = ((x-sep/2,x-sep/2), (y,y), (z+0.02, z+arl))
    ar = Arrow3D(*v, color='c', **arprop)
    ax.add_artist(ar)
    
    #hole
    hole, = ax.plot3D([x+sep/2], [y], [z],
              c='r', marker='$+$', markersize=7, ls='none')
    v = ((x+sep/2,x+sep/2), (y,y), (z-0.02, z-arl))
    ar = Arrow3D(*v, color='r', **arprop)
    ax.add_artist(ar)
    
    #photon
    draw_photon(ax, (0,0,1), location=(x,y,z-0.35), length=0.3, amplitude=0.02)
    
    #encircle
    ell = Ellipse((x,y), 2.5*sep, 2.5*sep, ls=':', ec='y', fc='none', lw=1)
    art3d.pathpatch_2d_to_3d(ell, z, 'z')
    ell._path2d = ell._path
    ax.add_patch(ell)
    
    return electron, hole

class ErrEllipse(object):
    """
    A class for creating and drawing an ellipse in matplotlib.
    """
    def __init__(self, x_cen, y_cen, a, b, pa, color='b'):
        """
        :param x_cen: x coordinate at center of the ellipse
        :param y_cen: y coordinate at center of the ellipse
        :param a: size of semi-major axes of the ellipse
        :param b: size of semi-minor axes of the ellipse
        :param pa: position angle of a to x  (90 ==> a is same orientation as x)
        """

        self.center = (x_cen, y_cen)
        self.a = max(a, 10)
        self.b = max(b, 10)
        self.pa = pa

        angle = self.pa - 90

        self.artist = Ellipse(self.center, self.a, self.b, angle=angle,
                              linewidth=3, edgecolor=color, facecolor='#E47833',
                              alpha=0.1)

    def add_to_axes(self, axes):
        self.artist.set_clip_box(axes.bbox)
        axes.add_patch(self.artist)

    def plot_error_ellipse_xy(cls, figure, x, y, abom, style): 
        """
        plots standard error ellipse
        x, y : position of center of ellipse
        abom : a, b, om parameters of ellipse

        method do not handle axes orientation (self.figure.axesXY)
        """

        from matplotlib.patches import Ellipse
        from math import pi

        a, b, om = abom

        #print "swapXY:", figure.is_swap_xy()
        if figure.is_swap_xy():
            y, x = x, y
            om = pi/2 - om
            #om = om - pi/2

        ell = Ellipse((x, y), 2*a, 2*b, om*180.0/pi) #?
                      #transform=self.axes.transData + self.ellTrans)

        figure.set_style(style, ell)
        
        ell.set_clip_box(figure.gca().bbox) # cut edges outside axes box?
        figure.gca().add_artist(ell)

def plot_data(data0, label0, title, subtitle, *args):
    
    colors = iter(cm.rainbow(np.linspace(0,1,5)))
    plt.clf()
    fig = plt.figure()
    ax = plt.gca()
    if (len(args) == 2):
        means, covmats = args
    elif (len(args) == 4):
        means, covmats, traindata, trainlabel = args
        trainlabel = trainlabel.reshape(trainlabel.size)
    for i in range(1,6):
        cl = next(colors)
        plt.scatter(data0[(label0==i).reshape((label0==i).size),0],
        data0[(label0==i).reshape((label0==i).size),1],color=cl)
        # Draw ellipse with 1 standard deviation
        if (len(args) >= 2):
            # Compute eigen vectors and value for covmat to plot an ellipse
            lambda_, v = np.linalg.eig(covmats[i-1])
            lambda_ = np.sqrt(lambda_)
            ell = Ellipse(xy=(means[0,i-1], means[1,i-1]),
                        width=lambda_[0]*2, height=lambda_[1]*2,
                        angle=np.rad2deg(np.arccos(v[0, 0])))
            ell.set_facecolor('none')
            ax.add_artist(ell)
            #Add mean points
            plt.scatter(means[0,i-1], means[1,i-1],c='black',s=30);
            if (len(args) == 4):
                plt.scatter(traindata[trainlabel==i,0],
                traindata[trainlabel==i,1],color=cl,
                edgecolor='black')
    fig.suptitle(title)
    ax.set_title(subtitle,fontsize='10')
    fig.show()
    return

def plotCovEllipse(cov, pos, nstd=2, ax=None, with_line=False, **kwargs):
    """
    Plots an `nstd` sigma error ellipse based on the specified covariance
    matrix (`cov`). Additional keyword arguments are passed on to the
    ellipse patch artist.
    Parameters
    ----------
        cov : The 2x2 covariance matrix to base the ellipse on
        pos : The location of the center of the ellipse. Expects a 2-element
            sequence of [x0, y0].
        nstd : The radius of the ellipse in numbers of standard deviations.
            Defaults to 2 standard deviations.
        ax : The axis that the ellipse will be plotted on. Defaults to the
            current axis.
        Additional keyword arguments are pass on to the ellipse patch.
    Returns
    -------
        A matplotlib ellipse artist
    """

    def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:, order]

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

    # largest eigenvalue is first
    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:, 0][::-1]))

    # Width and height are "full" widths, not radius
    width, height = 2 * nstd * np.sqrt(vals)
    ellip = Ellipse(xy=pos, width=width, height=height, angle=theta, **kwargs)

    if 'alpha' not in kwargs.keys():
        ellip.set_alpha(0.3)
    if 'color' not in kwargs.keys():# and 'c' not in kwargs.keys():
        ellip.set_facecolor('red')

    ax.add_patch(ellip)

    # THEN just f***ing plot an invisible line across the ellipse.
    if with_line:
        # brute forcing axes limits so they contain ellipse patch
        # maybe a cleaner way of doing this, but I couldn't work it out
        x_extent = 0.5*(abs(width*np.cos(np.radians(theta))) +
                        abs(height*np.sin(np.radians(theta))))
        y_extent = 0.5*(abs(width*np.sin(np.radians(theta))) +
                        abs(height*np.cos(np.radians(theta))))

        lx = pos[0] - x_extent
        ux = pos[0] + x_extent
        ly = pos[1] - y_extent
        uy = pos[1] + y_extent
        ax.plot((lx, ux), (ly, uy), alpha=0.)

    return ellip

def grad(im,Ix,Iy,prozor=5):

    #koordinate oko kojih ce se gledati gradijent
    imshow(im)
    koord = ginput(1)
    x = round(koord[0][0])
    y = round(koord[0][1])
    
    #uzmi gradijente Ix i Iy unutar tog prozora
    IxF = (Ix[x-prozor:x+prozor+1,y-prozor:y+prozor+1]).flatten()
    IyF = (Iy[x-prozor:x+prozor+1,y-prozor:y+prozor+1]).flatten()

    #neka kvazi matematika...
    width = abs(max(IxF)) + abs(min(IxF))
    height = abs(max(IyF)) + abs(min(IyF))
    xS = (max(IxF) + min(IxF))/2
    yS = (max(IyF) + min(IyF))/2

    """x_segment = width/15.0
    y_segment = width/15.0
    broj_pojava = zeros((15,15),dtype='int')

    #procjeni gdje je najgusce
    for i in range(15):
        for j in range(15):
            for k in IxF:
                for m in IyF:
                    if IxF[k] >= i*x_segment and IxF[k] <= (i+1)*x_segment:
                        if IyF[m] >= j*y_segment and IyF[m] <= (j+1)*y_segment:
                            broj_pojava[i][j] += 1

    x_pom,y_pom = unravel_index(broj_pojava.argmax(), broj_pojava.shape)

    x = (x_pom*x_segment + (x_pom+1)*x_segment)/2
    y = (y_pom*y_segment + (y_pom+1)*y_segment)/2

    angle = arctan(y*1.0/x*1.0)*180/pi"""
            
    
    
    mod = sqrt(IxF**2 + IyF**2)
    pom = argmax(mod)
    angle = arctan(IyF[pom]*1.0/IxF[pom]*1.0)*180/pi
    if width < max(mod):
        width = max(mod)
    
    ells = Ellipse((xS,yS),width,height,angle)
    ells.set_alpha(0.3)


    
    fig = figure()
    ax = fig.add_subplot(111,aspect='equal')
    ax.add_artist(ells)
    ax.plot(IxF,IyF,'mo') #crtaj histogram
    gray()
    grid('on')
    show()
    return