def plot(shiftRelabel):
   xs, srBias, srError = np.array(shiftRelabel).T
   #xs2, trError, trBias = np.array(thresholdRelabel).T

   #f, (ax1, ax2) = plt.subplots(2,1)
   width = 3

   #plt.plot(figsize=(plt.figaspect(2.0)))

   plt.plot(xs, srError, label='Label error', linewidth=width)
   plt.plot(xs, srBias, label='Bias', linewidth=width)
   plt.title("Shifted Decision Boundary Bias vs Error")
   plt.gca().invert_xaxis()
   plt.axhline(0, color='black')
   plt.figaspect(10.0)
   #handles, labels = plt.get_legend_handles_labels()

   # ax2.plot(xs2, trError, label='Label error', linewidth=width)
   # ax2.plot(xs2, trBias, label='Bias', linewidth=width)
   # ax2.set_title("Margin Threshold Relabeling Bias vs Error")
   # ax2.axhline(0, color='black')

   #plt.subplots_adjust(hspace=.75)
   #plt.figlegend(handles, labels, 'center right')
   plt.legend(loc='center right')
   #plt.savefig("relabeling-msr-tradeoffs.pdf",bbox_inches='tight')
   plt.show()

def test_figaspect():
    w, h = plt.figaspect(np.float64(2) / np.float64(1))
    assert h / w == 2
    w, h = plt.figaspect(2)
    assert h / w == 2
    w, h = plt.figaspect(np.zeros((1, 2)))
    assert h / w == 0.5
    w, h = plt.figaspect(np.zeros((2, 2)))
    assert h / w == 1

def main(argv):
    
    plt.figure(figsize=plt.figaspect(0.55))
    xd=np.linspace(0,4000,100)
    td = getTC(xd)
    plt1, = plt.plot(td,xd,'b--',lw=3)
    plt.grid()
    plt.xlabel(r'$\theta_D \rm[deg]$',fontsize=20)
    plt.ylabel(r'$x_d \rm[m]$',fontsize=20)
    
    plt.text(91.55, 1200, 'Zona permitida',
        bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
    
    plt.legend([plt1],["$x_d^{CUT}$"],loc=2)
    
    plt.savefig("thetaDCut.pdf",format="pdf")
    
    
    thD = np.arange(90.1,95,0.00001)
    cThD = np.cos(thD*kdeg2rad)
    #Xd = np.arange(100,1700,200)
    Xd = np.arange(100,3700,400)
    
    plt.figure(figsize=plt.figaspect(0.55))
    plots = []
    xds = []
    for x in Xd:
        cThE = np.sqrt((R**2 - ((R+x)**2) * (1-cThD**2))/R**2)
        auxP, = plt.plot(thD,(x*(R*cThE-(R+x)*cThD))/((R**2 - (R+x)**2)*cThD),c=cm.gist_rainbow(float(x-5)/Xd[-1],1),lw=2.5)
        plots.append(auxP)
        xds.append(str(x))
    plt.grid()
    plt.legend(plots,xds,loc=4,title=r'$\bf{x_d \rm [m]}$')
    plt.xlabel(r'$\theta_D \rm[deg]$',fontsize=20)
    plt.ylabel(r'$\frac{l_{plano}}{l_{curvo}}$',fontsize=20)
    plt.savefig("lPlane_lCurve.pdf",format="pdf")
    
    
    plt.figure(figsize=plt.figaspect(0.55))
    plots = []
    xds = []
    for x in Xd:
        cThE = -np.sqrt((R**2 - ((R+x)**2) * (1-cThD**2))/R**2)
        thE = np.arccos(cThE)*krad2deg
        auxP, = plt.plot(thD,thE,lw=2.5,c=cm.gist_rainbow(float(x-5)/Xd[-1],1))
        plots.append(auxP)
        xds.append(str(x))
    plt.grid()
    plt.legend(plots,xds,loc=4,title=r'$\bf{x_d \rm [m]}$')
    plt.xlabel(r'$\theta_D \rm[deg]$',fontsize=20)
    plt.ylabel(r'$\theta_E \rm[deg]$',fontsize=20)
    plt.savefig("thE_thD.pdf",format="pdf")

def _plot_orbits(x, y, z):
    from matplotlib import pyplot, rc
    fig = pyplot.figure(figsize=(14,14))
    font = {'size' : 30}
    rc('font', **font)
    pyplot.figaspect(1.0)
    pyplot.scatter(x[0].value_in(units.AU), y[0].value_in(units.AU), 200, lw=0)
    pyplot.plot(x.value_in(units.AU), y.value_in(units.AU), lw=2)
    pyplot.xlabel("$X [AU]$")
    pyplot.ylabel("$Y [AU]$")
    pyplot.xlim(-15000, 15000)
    pyplot.ylim(-30000, 10000)
#    pyplot.show()
    pyplot.savefig("SStars_1Jan2001_orbits")

def showslice2(thedata, thelabel, minval, maxval, colormap):
    # initialize and show a 2D slice from a dataset in greyscale
    plt.figure(figsize=plt.figaspect(1.0))
    theshape = thedata.shape
    numslices = theshape[0]
    ysize = theshape[1]
    xsize = theshape[2]
    slicesqrt = int(np.ceil(np.sqrt(numslices)))
    theslice = np.zeros((ysize * slicesqrt, xsize * slicesqrt))
    for i in range(numslices):
        ypos = int(i / slicesqrt) * ysize
        xpos = int(i % slicesqrt) * xsize
        theslice[ypos:ypos + ysize, xpos:xpos + xsize] = thedata[i, :, :]
    if plt.isinteractive():
        plt.ioff()
    plt.axis('off')
    plt.axis('equal')
    plt.subplots_adjust(hspace=0.0)
    plt.axes([0, 0, 1, 1], frameon=False)
    if colormap == 0:
        thecmap = cm.gray
    else:
        mycmdata1 = {
            'red': ((0., 0., 0.), (0.5, 1.0, 0.0), (1., 1., 1.)),
            'green': ((0., 0., 0.), (0.5, 1.0, 1.0), (1., 0., 0.)),
            'blue': ((0., 0., 0.), (0.5, 1.0, 0.0), (1., 0., 0.))
        }
        thecmap = colors.LinearSegmentedColormap('mycm', mycmdata1)
    plt.imshow(theslice, vmin=minval, vmax=maxval, interpolation='nearest', label=thelabel, aspect='equal',
               cmap=thecmap)

def main():
    import matplotlib.pyplot as plt
    from sklearn.datasets.samples_generator import make_blobs
    n_centers = 3
    X, y = make_blobs(n_samples=1000, centers=n_centers, n_features=2,
                    cluster_std=0.7, random_state=0)

    # Run this K-Means
    import kmeans
    t0 = time.time()
    y_pred, centers, obj_val_seq = kmeans.kmeans(X, n_centers)
    t1 = time.time()
    print("Final obj val: {}".format(obj_val_seq[-1]))
    print("Time taken (this implementation): {}".format(t1 - t0))

    # Run scikit-learn's K-Means
    from sklearn.cluster import k_means
    t0 = time.time()
    centers, y_pred, obj_val = k_means(X, n_centers, random_state=0)
    t1 = time.time()
    print("Final obj val: {}".format(obj_val))
    print("Time taken (Scikit, 1 job): {}".format(t1 - t0))

    # Plot change in objective value over iteration
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(obj_val_seq, 'b-', marker='*')
    fig.suptitle("Change in K-means objective value across iterations")
    ax.set_xlabel("Iteration")
    ax.set_ylabel("Objective value")
    fig.show()

    # Plot data
    from itertools import cycle
    colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
    fig = plt.figure(figsize=plt.figaspect(0.5))  # Make twice as wide to accomodate both plots
    ax = fig.add_subplot(121)
    ax.set_title("Data with true labels and final centers")
    for k, color in zip(range(n_centers), colors):
        ax.plot(X[y==k, 0], X[y==k, 1], color + '.')

    initial_centers = kmeans.init_centers(X, n_centers, 2) # This is valid because we always use the same random seed.
    # Plot initial centers
    for x in initial_centers:
        ax.plot(x[0], x[1], "mo", markeredgecolor="k", markersize=8)

    # Plot final centers
    for x in centers:
        ax.plot(x[0], x[1], "co", markeredgecolor="k", markersize=8)

    # Plot assignments
    colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
    ax = fig.add_subplot(122)
    ax.set_title("Data with final assignments")
    for k, color in zip(range(n_centers), colors):
        ax.plot(X[y_pred==k, 0], X[y_pred==k, 1], color + '.')

    fig.tight_layout()
    fig.gca()
    fig.show()

def draw_trajectories_eps(paths_list):
	"""
	Same as draw_trajectories with additional labels
	Plots trajectories listed in the paths_list tuple.
	Paths tuple contains other tuples of the form
	(x0,x1,...) where xi is the i^th coordinate of
	the diffusion.
	"""

	# set up figures
	fig = plt.figure(figsize=plt.figaspect(2.))
	fig.suptitle('Paths simulated using Euler-Maruyama scheme.')

	for (i,paths) in enumerate(paths_list):
		ax = fig.add_subplot(1, len(paths_list), i+1, projection='3d')

		# plot the path of a realisation of a diffusion
		for path in paths:
			ax.plot(path[0], path[1], path[2])

		# add labels
		ax.grid(True)
		text = "time step: " + str(np.power(10.0, (-i-1)))
		ax.set_title(text)
	plt.show()

def plot_3d(X,Y,Z,plot_type = "surface",ax = None,color = "blue"):
    from scipy.interpolate import griddata
    if ax == None:
        fig = plt.figure(figsize=plt.figaspect(0.5))
        ax = fig.add_subplot(1, 1, 1, projection='3d')

    if plot_type == "scatter":
        ax.scatter(X,Y,Z)
        return
    x = np.array(X)
    y = np.array(Y)
    z = np.array(Z)
    xi = np.linspace(x.min(),x.max(),100)
    yi = np.linspace(y.min(),y.max(),100)
    # VERY IMPORTANT, to tell matplotlib how is your data organized
    zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='linear')
    xig, yig = np.meshgrid(xi, yi)


    if plot_type == "contour":
        CS = plt.contour(xi,yi,zi,15,linewidths=0.5,color='k')
    else:
        surf = ax.plot_surface(xig, yig, zi,linewidth=0,color=color)#,cmap=cm.coolwarm)

    return ax

    def componentes(self):

        fig = plt.figure(figsize=plt.figaspect(0.5))

        ax = fig.add_subplot(2, 3, 1, projection='3d')
        surf = ax.plot_surface(self.X, self.Y,self.SOL[:,:,0],cmap=cm.coolwarm)
        plt.xlabel("longitud (m)")
        plt.ylabel("tiempo (horas)")
        plt.title("ciclohexanol (kmol/h)")

        ax = fig.add_subplot(2, 3, 2, projection='3d')
        surf = ax.plot_surface(self.X, self.Y, self.SOL[:,:,1],cmap=cm.coolwarm)
        plt.xlabel("longitud (m)")
        plt.ylabel("tiempo (horas)")
        plt.title("ciclohexanona (kmol/h)")

        ax = fig.add_subplot(2, 3, 3, projection='3d')
        surf = ax.plot_surface(self.X, self.Y, self.SOL[:,:,2],cmap=cm.coolwarm)
        plt.xlabel("longitud (m)")
        plt.ylabel("tiempo (horas)")
        plt.title("fenol (kmol/h)")

        ax = fig.add_subplot(2, 3, 4, projection='3d')
        surf = ax.plot_surface(self.X, self.Y, self.SOL[:,:,3]*1000,cmap=cm.coolwarm)
        plt.xlabel("longitud (m)")
        plt.ylabel("tiempo (horas)")
        plt.title("ciclohexenona (mol/h)")

        ax = fig.add_subplot(2, 3, 6, projection='3d')
        surf = ax.plot_surface(self.X, self.Y, self.SOL[:,:,4],cmap=cm.coolwarm)
        plt.xlabel("longitud (m)")
        plt.ylabel("tiempo (horas)")
        plt.title("H2 (kmol/h)")

        plt.show()

    def moduloThiele(self):
        mL = np.zeros((self.nt,self.nl))
        rr = kn.CuZn()
        for i in range(self.nt):
            for j in range(self.nl):
                n = self.SOL[i,j,0:5] #kmol/h
                T = self.SOL[i,j,5]+273.15# K
                P = self.SOL[i,j,7] #atm
                a = self.SOL[i,j,8]
                L = datos.Dint/2/2 #m
                r =  -rr.dnjdW(n,T,P,a)[0]/3600/1000 #kmol/(kgs), (tiene que ser positiva)
                if (r<0): r = 0
                C_ol = n[0]/((sum(n))*datos.R*T/P) #kmol/m3
                k = r*datos.ro_l/C_ol # s(-1)
                D_OL = datos.D0_OL*np.exp(-datos.Ed_OL/((datos.R*101325)*T)) #cambiar la R a kJ/molK
                De = D_OL*datos.e_cat**2 #m2/s
                X = 0.712505333333 #grado de conversion de ciclohexanol en equilibrio
                mL[i,j] = L*(k/De/X)**0.5 #adimensional

        fig = plt.figure(figsize=plt.figaspect(0.5))
        ax = fig.add_subplot(1, 1, 1, projection='3d')
        surf = ax.plot_surface(self.X, self.Y, mL, cmap=cm.coolwarm)
        plt.xlabel("longitud (m)")
        plt.ylabel("tiempo (horas)")
        plt.title("modulo de Thiele")
        plt.show()

    def run(self):
        import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d import Axes3D
        import matplotlib.animation as ani

        self.fig = plt.figure(figsize=plt.figaspect(0.5))
        # setup 3d axis
        self.ax3d = self.fig.add_subplot(1, 2, 1, projection="3d")
        self.ax3d.set_xlabel("X")
        self.ax3d.set_ylabel("Y")
        self.ax3d.set_zlabel("Z")
        self.ax3d_lines = [self.ax3d.plot([1, 1], [2, 2], [3, 3], c="b")[0] for i in range(17 * 12)]
        self.ax3d_scatters = [self.ax3d.scatter(0, 0) for i in range(17)]
        # setup 2d axis
        self.ax_flat = self.fig.add_subplot(1, 2, 2)
        self.ax_flat.set_xlim([-17, 17])
        self.ax_flat.set_ylim([-17, 17])
        self.ax_flat.set_xlabel("X")
        self.ax_flat.set_ylabel("Y")
        self.ax_flat.grid()
        self.ax_flat_scatters = [self.ax_flat.scatter(0, 0) for i in range(2)]
        self.ax_flat_lines = [self.ax_flat.plot([], [], c="b")[0] for i in range(5)]
        self.BASE_CUBE_POINTS = np.array(
            [[-1, -1, -1], [-1, -1, 1], [-1, 1, -1], [-1, 1, 1], [1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1]]
        )
        self.BASE_CUBE_POINTS = np.transpose(np.insert(self.BASE_CUBE_POINTS, 3, 1, axis=1))
        a = ani.FuncAnimation(self.fig, self.runani, self.data_gen, blit=False, interval=300)
        plt.show()

def esfera(meia):
    fig = plt.figure(figsize=plt.figaspect(1)) 
    ax = fig.add_subplot(111, projection='3d')  

    vx= eval(input("Raio da esfera(r): "))
    coefs = (vx, vx, vx)  
    rx, ry, rz = [1/np.sqrt(coef) for coef in coefs]

    if meia == 1:
        u = np.linspace(0,1 * np.pi, 100)
    else:
        u = np.linspace(0,2 * np.pi, 100)
    v = np.linspace(0,  np.pi, 100)

    x = rx * np.outer(np.cos(u), np.sin(v))
    y = ry * np.outer(np.sin(u), np.sin(v))
    z = rz * np.outer(np.ones_like(u), np.cos(v))


    ax.plot_surface(x, y, z,  rstride=4, cstride=4, color='b')

    max_radius = max(rx, ry, rz)
    for axis in 'xyz':
        getattr(ax, 'set_{}lim'.format(axis))((-max_radius, max_radius))

    plt.show()

def plotAllTheThings(pl1,pl2,thete,lines, contourCount = 20,simbol="x",figname=""):
	# Twice as wide as it is tall.
	fig = plt.figure(figsize=plt.figaspect(0.4))

	#---- First subplot
	ax = fig.add_subplot(1, 2, 1)
	X,Y,Z = pl1
	for t in thete:
		ax.plot(t[0],t[1],simbol,color="blue")

	levels = np.arange(Z.min(), Z.max(), (Z.max()-Z.min())/contourCount)
	ax.contour(X, Y, Z,levels)
	ax.set_title("Konvergenca theta_0 in theta_1")


	#---- Second subplot
	a,b = pl2
	a = a[1]
	ax = fig.add_subplot(1,2,2)
	ax.plot(a, b, 'o')
	
	lsp = np.linspace(a.min(),a.max(),1+a.max()-a.min())
	print lsp
	for ll in lines:
		l = ll[0]+ll[1]*lsp
		ax.plot(lsp,l)
	ax.set_title("Prikaz tock z prilegajoco premico")
	if figname == "":
		plt.show()
	else:
		plt.savefig(figname)

 def plotMe(self):
     fig = plt.figure(figsize=plt.figaspect(1))
     # Square figure
     ax = fig.add_subplot(111, projection='3d')
     getattr(ax, 'set_{}lim'.format('x'))((-np.pi / 4, np.pi / 4))
     getattr(ax, 'set_{}lim'.format('y'))((-np.pi / 4, np.pi / 4))
     getattr(ax, 'set_{}lim'.format('z'))((-np.pi / 4, 2.25 * np.pi))

def plot_ellipse_mpl(Dxx, Dyy, Dzz, rstride=4, cstride=4, color='b'):
    """
    Plot an ellipse shape (in contrast to a tensor ellipsoid, see
    plot_ellipsoid_mpl for that) in 3d using matplotlib
    
    """
    
    fig = plt.figure(figsize=plt.figaspect(1))  # Square figure
    ax = fig.add_subplot(111, projection='3d')

    coefs = Dxx, Dyy, Dzz # Coefficients in a0/c x**2 + a1/c y**2 + a2/c z**2 =

    # 1 Radii corresponding to the coefficients:
    rx, ry, rz = Dxx, Dyy, Dzz#[1/np.sqrt(coef) for coef in coefs]

    # Set of all spherical angles:
    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)

    # Cartesian coordinates that correspond to the spherical angles:
    # (this is the equation of an ellipsoid):
    x = rx * np.outer(np.cos(u), np.sin(v))
    y = ry * np.outer(np.sin(u), np.sin(v))
    z = rz * np.outer(np.ones_like(u), np.cos(v))

    # Plot:
    ax.plot_surface(x, y, z,  rstride=rstride, cstride=cstride, color=color)

    # Adjustment of the axes, so that they all have the same span:
    max_radius = max(rx, ry, rz)
    for axis in 'xyz':
        getattr(ax, 'set_{}lim'.format(axis))((-max_radius, max_radius))

    return fig

def test():
    # Twice as wide as it is tall.
    fig = plt.figure(figsize=plt.figaspect(0.5))

    #---- First subplot
    ax = fig.add_subplot(1, 2, 1, projection='3d')
    X = np.arange(-5, 5, 0.25)
    Y = np.arange(-5, 5, 0.25)
    X, Y = np.meshgrid(X, Y)
    R = np.sqrt(X**2 + Y**2)
    Z = np.sin(R)
    print type(X)
    print X.shape
    print Y.shape
    print Z.shape
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
            linewidth=0, antialiased=False)
    ax.set_zlim3d(-1.01, 1.01)

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

    #---- Second subplot
    ax = fig.add_subplot(1, 2, 2, projection='3d')
    X, Y, Z = get_test_data(0.05)
    ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)

    plt.show()

def plotFile(file):
    
    data = []
    
    for row in file:
        data.append(map(float,row.strip().split()))
    
    adata = np.array(data)
    adata = adata[np.argsort(adata[:,1])]
    norm = max(adata[:,2])
    adata[:,2] = np.multiply(adata[:,2],1./norm)
    #~ p0 = [1,-1,-80,-0.5,0.1]
    p0 = [1,-40,3660,1./375,13./30,0]
    #~ p0 = [1,-0.5,1,100,1,0]
    pOut = leastsq(regF,p0, args=(adata[:,2],(adata[:,0],adata[:,1])),full_output=1)
    print pOut[0]
    
    xs = adata[:,0]
    ys = adata[:,1]
    zs = adata[:,2]
    
    nPoints = 31
    xd = np.linspace(0, 300, nPoints)
    th = np.linspace(87, 90, nPoints)
    X,Y = np.meshgrid(xd, th)
    Z = fitF((X, Y),pOut[0])
    #~ aspect = 1./(1+np.sqrt(5)/2.)
    fig = plt.figure(figsize=plt.figaspect(0.55))
    ax = Axes3D(fig)
    ax.scatter(xs,ys,zs,c='y',marker='o',s=2000)
    ax.plot_wireframe(X, Y, Z,cmap='hot')

def plot(filename):

	states=[ 2, 10, 18, 26, 42, 50, 58, 74, 90, 98, 114, 122, 138, 162, 178, 194, 202, 218, 226, 242, 258, 274, 290, 298, 322, 338, 354, 370, 386, 394, 426, 442, 450, 466, 482, 498, 506, 522, 554, 570, 586, 602, 610, 634, 650, 666, 682, 698, 714, 730, 746, 754, 770, 802, 810, 842, 858, 874, 882, 914, 930, 946, 962, 978, 994, 1010, 1018, 1034, 1058, 1090, 1106, 1122, 1138, 1154, 1186, 1202, 1218, 1226, 1242, 1266, 1282, 1314, 1330, 1346, 1362, 1394, 1418, 1434, 1450, 1466, 1482, 1498, 1514, 1522, 1538, 1554, 1586, 1594, 1610, 1642, 1658, 1690, 1706, 1722, 1738, 1754, 1770, 1778, 1802, 1834, 1850, 1866, 1882, 1898, 1930, 1946, 1962, 1978, 1994, 2010, 2018, 2066, 2082, 2098, 2114, 2138, 2170, 2186, 2202, 2218, 2234, 2250, 2258, 2274, 2306, 2322, 2354, 2370, 2402, 2418, 2434, 2450, 2458, 2474, 2490, 2514, 2530, 2546, 2562, 2578, 2610, 2626, 2642, 2658, 2706, 2722, 2738, 2746, 2778, 2810, 2826, 2850, 2866, 2882, 2898, 2914, 2930, 2946, 2962, 2978, 3010, 3026, 3034, 3066, 3082, 3098, 3130, 3162, 3194, 3210, 3218, 3234, 3266, 3282, 3298, 3306, 3338, 3370, 3386, 3402, 3418, 3450, 3466, 3482, 3498, 3514, 3530, 3562, 3578, 3586, 3602, 3626, 3658, 3674, 3706, 3722, 3738, 3754, 3770, 3786, 3802, 3834, 3850, 3866, 3882, 3922, 3938, 3954, 3986, 4002, 4018, 4034]

	
	data = loadtxt(filename)

	l = len(data)
	dataStates = ones(l)

	for i in range(l):
		dataStates[i] = states[(int)(data[i,0]-1+0.1)]

	fig = plt.figure(figsize=plt.figaspect(1.))
	ax = fig.gca()

	ax.set_xlabel(r'Basis size',fontsize=16)
	ax.set_ylabel(r'Corr. energy',fontsize=16)

	for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14) 

	for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)

	plt.plot(dataStates[:], data[:,2],'b')
	#ax.set_xticks(dataStates, minor=False)

	#savefig('test.pdf', format='pdf')

	plt.show()

def show3D():
    # imports specific to the plots in this example
    import numpy as np
    from matplotlib import cm
    from mpl_toolkits.mplot3d.axes3d import get_test_data
    
    # Twice as wide as it is tall.
    fig = plt.figure(figsize=plt.figaspect(0.5))
    
    #---- First subplot
    ax = fig.add_subplot(1, 2, 1, projection='3d')
    X = np.arange(-5, 5, 0.1)
    Y = np.arange(-5, 5, 0.1)
    X, Y = np.meshgrid(X, Y)
    R = np.sqrt(X**2 + Y**2)
    Z = np.sin(R)
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,
            linewidth=0, antialiased=False)
    ax.set_zlim3d(-1.01, 1.01)
    
    fig.colorbar(surf, shrink=0.5, aspect=10)
    
    #---- Second subplot
    ax = fig.add_subplot(1, 2, 2, projection='3d')
    X, Y, Z = get_test_data(0.05)
    ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)

    mystyle.printout_plain('3dGraph.png')

def Main():
	myPoints = importPoints('points.csv', 3) # Import the Given Point Cloud

	C = np.array([[0,0,100]]) # View Point, which is well above the point cloud in z direction
	flippedPoints = sphericalFlip(myPoints, C, math.pi) # Reflect the point cloud about a sphere centered at C
	myHull = convexHull(flippedPoints) # Take the convex hull of the center of the sphere and the deformed point cloud


	# Plot
	fig = plt.figure(figsize = plt.figaspect(0.5))
	plt.title('Cloud Points With All Points (Left) vs. Visible Points Viewed from Well Above (Right)')
	
	# First subplot
	ax = fig.add_subplot(1,2,1, projection = '3d')
	ax.scatter(myPoints[:, 0], myPoints[:, 1], myPoints[:, 2], c='r', marker='^') # Plot all points
	ax.set_xlabel('X Axis')
	ax.set_ylabel('Y Axis')
	ax.set_zlabel('Z Axis')

	# Second subplot
	ax = fig.add_subplot(1,2,2, projection = '3d')
	for vertex in myHull.vertices[:-1]: # Exclude Origin, which is the last element
		ax.scatter(myPoints[vertex, 0], myPoints[vertex, 1], myPoints[vertex, 2], c='b', marker='o') # Plot visible points
	ax.set_xlabel('X Axis')
	ax.set_ylabel('Y Axis')
	ax.set_zlabel('Z Axis')

	plt.show()

	return