def make_chromatin_movies():
    """
    """
    dnachain.MultiSolenoidVolume(voxelheight=1500, separation=400)\
        .to_line_plot()
    fig = mlab.gcf()
    fig.scene.set_size([1024, 768])
    fig.scene.render()
    time.sleep(1)
    mlab.savefig("example/img/chromatin_multi_straight0000.png")
    step = 5
    for ii in range(1, int(720/step)):
        fig.scene.camera.azimuth(step)
        fig.scene.render()
        fname = "example/img/chromatin_multi_straight{:04d}.png".format(ii)
        mlab.savefig(fname)

    dnachain.MultiSolenoidVolume(voxelheight=1500, separation=400, turn=True)\
        .to_line_plot()
    fig = mlab.gcf()
    fig.scene.set_size([1024, 768])
    fig.scene.render()
    time.sleep(1)
    mlab.savefig("example/img/chromatin_multi_turn0000.png")
    step = 5
    for ii in range(1, int(720/step)):
        fig.scene.camera.azimuth(step)
        fig.scene.render()
        fname = "example/img/chromatin_multi_turn{:04d}.png".format(ii)
        mlab.savefig(fname)

    dnachain.TurnedSolenoid().to_line_plot()
    fig = mlab.gcf()
    fig.scene.set_size([1024, 768])
    fig.scene.render()
    time.sleep(1)
    mlab.savefig("example/img/chromatin_single_turn0000.png")
    step = 5
    for ii in range(1, int(720/step)):
        fig.scene.camera.azimuth(step)
        fig.scene.render()
        fname = "example/img/chromatin_single_turn{:04d}.png".format(ii)
        mlab.savefig(fname)

    dnachain.Solenoid().to_line_plot()
    fig = mlab.gcf()
    fig.scene.set_size([1024, 768])
    fig.scene.render()
    time.sleep(1)
    mlab.savefig("example/img/chromatin_single_straight0000.png")
    step = 5
    for ii in range(1, int(720/step)):
        fig.scene.camera.azimuth(step)
        fig.scene.render()
        fname = "example/img/chromatin_single_straight{:04d}.png".format(ii)
        mlab.savefig(fname)
    return None

def lorenz_animation(N=10, res=1000, step=2, t=10, seed_=120, atol=1E-15,
                     rtol=1E-13, delay=10, sigma=10., beta=8./3, rho=28.):
    """ Animate the trajectories given by the Lorenz equations for 'N' starting points.
    Choose random x, y, and z values between -15 and 15.
    Seed the random number generator with 'seed_'.
    Use a resolution of 'res' for the points in the plot.
    Plot the time values between 0 ant 't'.
    When computing the trajectories, pass the tolerances
    'atol' and 'rtol' to the ODE solver.
    At each update, add 'step' points to the plot.
    Use a delay of 'delay' at each update in the animation.
    Use different colors for each trajectory.
    Use the values of 'sigma', 'beta', and 'rho' in the Lorenz ODE. """
    
    # Get initial conditions.
    seed(seed_)
    x0 = -15 + 30 * rand(N, 3)
    
    # Solve for the trajectories.
    t = np.linspace(0, t, res)
    pts = np.empty((N, res, 3))
    for i, x in enumerate(x0):
        pts[i] = odeint(lorenz_ode, x, t,
                        args=(sigma, beta, rho), rtol=rtol, atol=atol)
    
    # Select the colors for the different curves.
    colors = np.zeros((N, 3))
    colors[:,1] = np.linspace(0, 1, N)
    colors = map(tuple, colors.tolist())
    
    # Plot the different trajectories.
    contours = [mlab.plot3d(x[:1,0], x[:1,1], x[:1,2], tube_radius=.15, color=color)
                for x, color in zip(pts, colors)]
    
    # Position the view for the plot.
    mlab.gcf().scene.camera.position = [127.23761585, -108.28736806, 6.35191272]
    mlab.gcf().scene.camera.focal_point = [-1.7792501449584961, -3.6287221908569336, 23.397351264953613]
    mlab.gcf().scene.camera.view_up = [-0.078467260964232038, -0.20339450183237351, 0.97594752194015633]
    mlab.gcf().scene.camera.clipping_range = [128.64624663718814, 328.22549479639167]
    
    # Define the animation.
    @mlab.show
    @mlab.animate(delay=delay)
    def trace_curve():
        for i in xrange(step, res, step):
            for c, x, color in zip(contours, pts, colors):
                c.mlab_source.reset(x=x[:i,0], y=x[:i,1], z=x[:i,2])
            yield
    
    # Run the animation.
    trace_curve()

def make_fractal_movie(seed="X", iter=6):
    """
    """
    f = hilbert.VoxelisedFractal.fromSeed(seed, iter)
    f.center_fractal()
    pos = np.array([vox.pos for vox in f.fractal])
    max_pos = np.max(pos, axis=0)
    mask = lambda x: np.sum((x[0:3]/max_pos[0:3])**2) <= 1
    f = f.to_pretty_plot(refine=5, mayavi=True, mask=mask)
    fig = mlab.gcf()
    fig.scene.set_size([1024, 768])
    fig.scene.render()
    mlab.savefig("example/img/fractal_rotate0000.png")
    step = 5
    for ii in range(1, int(720/step)):
        fig.scene.camera.azimuth(step)
        fig.scene.render()
        mlab.savefig("example/img/fractal_rotate{:04d}.png".format(ii))

    mlab.savefig("example/img/fractal_zoom0000.png")
    nsteps = 200
    for ii in range(nsteps):
        fig.scene.camera.zoom(1.02)
        fig.scene.render()
        mlab.savefig("example/img/fractal_zoom{:04d}.png".format(ii))

    for ii in range(nsteps, nsteps + int(720/step)):
        fig.scene.camera.azimuth(step)
        fig.scene.render()
        mlab.savefig("example/img/fractal_zoom{:04d}.png".format(ii))
    return None

    def action(u, x, xv, y, yv, t, n):
        #print 'action, t=',t,'\nu=',u, '\Nx=',x, '\Ny=', y
        if plot == 1:
            mesh(xv, yv, u, title='t=%g' %t[n])
            time.sleep(0.2) # pause between frames

        elif plot == 2:
            # mayavi plotting
            mlab.clf()
            extent1 = (0, 20, 0, 20,-2, 2)
            s = mlab.surf(x , y, u, colormap='Blues', warp_scale=5,extent=extent1)
            mlab.axes(s, color=(.7, .7, .7), extent=extent1,
            ranges=(0, 10, 0, 10, -1, 1), xlabel='', ylabel='',
            zlabel='',
            x_axis_visibility=False, z_axis_visibility=False)
            mlab.outline(s, color=(0.7, .7, .7), extent=extent1)
            mlab.text(2, -2.5, '', z=-4, width=0.1)
            mlab.colorbar(object=None, title=None, orientation='horizontal', nb_labels=None, nb_colors=None, label_fmt=None)
            mlab.title('test 1D t=%g' % t[n])

            mlab.view(142, -72, 50)
            f = mlab.gcf()
            camera = f.scene.camera
            camera.yaw(0)
        
        
        if plot > 0:
            path = 'Figures_wave2D'
            time.sleep(0) # pause between frames
            if save_plot and plot != 2:
                filename = '%s/%08d.png' % (path, n)
                savefig(filename)  # time consuming!
            elif save_plot and plot == 2:
                filename = '%s/%08d.png' % (path,n)
                mlab.savefig(filename)  # time consuming!

def run_mlab_file(filename, image_file):
    ## XXX: Monkey-patch mlab.show, so that we keep control of the
    ## the mainloop
    old_show = mlab.show
    def my_show(func=None):
        pass
    mlab.show = my_show
    mlab.clf()
    e = mlab.get_engine()
    e.close_scene(mlab.gcf())
    execfile(filename, {'__name__': '__main__'})
    mlab.savefig(image_file)
    size = mlab.gcf().scene.get_size()
    for scene in e.scenes:
        e.close_scene(scene)
    mlab.show = old_show

 def _render_model(self):
   figure = mlab.gcf()
   mlab.clf()
   s = mlab.triangular_mesh(self.coords[:,0], self.coords[:,1],
                            self.coords[:,2], self.tri_index,
                            color=(0.5,0.5,0.5))
   return s.scene

	def _draw(self):
		"""Update the Mayavi objects with new particle information.
		This is called periodically in the GUI thread"""
		if self.data is None:
			return
		
		assert isinstance(threading.current_thread(), threading._MainThread)

		f = mlab.gcf()
		visual.set_viewer(f)

		coords, types, radii, N_changed, bonds, Nbonds_changed, boxl, box_changed = self.data
		self.data = None

		if box_changed or not self.running:
			self.box.set(bounds=(0,boxl[0], 0,boxl[1], 0,boxl[2]))

		if not N_changed:
			self.points.mlab_source.set(x=coords[:,0]%boxl[0], y=coords[:,1]%boxl[1], z=coords[:,2]%boxl[2], u=radii, v=radii, w=radii, scalars=types)
		else:
			self.points.mlab_source.reset(x=coords[:,0]%boxl[0], y=coords[:,1]%boxl[1], z=coords[:,2]%boxl[2], u=radii, v=radii, w=radii, scalars=types)
		if not self.running:
			f.scene.reset_zoom()
			self.running = True

		if not Nbonds_changed:
			if bonds.shape[0] > 0:
				self.arrows.mlab_source.set  (x=bonds[:,0], y=bonds[:,1], z=bonds[:,2], u=bonds[:,3], v=bonds[:,4], w=bonds[:,5])
		else:
			self.arrows.mlab_source.reset(x=bonds[:,0], y=bonds[:,1], z=bonds[:,2], u=bonds[:,3], v=bonds[:,4], w=bonds[:,5], scalars=bonds[:,6])

def example():
    # setup movie
    m = movie(frame_rate=22.5)
    
    # setup example
    from mayavi import mlab
    mlab.options.offscreen = True
    mlab.test_contour3d()
    f = mlab.gcf() 

    # save frames
    for i in range(36*3):
        f.scene.camera.azimuth(10) # rotate
        f.scene.render()
        mlab.savefig(m.next_frame(), figure=f)
    
    # encode frames
    encoded = m.encode()

    # remove tempdir
    del m
    
    if not encoded:
        # report error
        exit(m.output)

def m2screenshot(mayavi_fig=None, mpl_axes=None, autocrop=True):
    """ Capture a screeshot of the Mayavi figure and display it in the
        matplotlib axes.
    """
    import pylab as pl
    # Late import to avoid triggering wx imports before needed.
    try:
        from mayavi import mlab
    except ImportError:
        # Try out old install of Mayavi, with namespace packages
        from enthought.mayavi import mlab

    if mayavi_fig is None:
        mayavi_fig = mlab.gcf()
    else:
        mlab.figure(mayavi_fig)
    if mpl_axes is not None:
        pl.axes(mpl_axes)

    filename = tempfile.mktemp('.png')
    mlab.savefig(filename, figure=mayavi_fig)
    image3d = pl.imread(filename)
    if autocrop:
        bg_color = mayavi_fig.scene.background
        image3d = autocrop_img(image3d, bg_color)
    pl.imshow(image3d)
    pl.axis('off')
    os.unlink(filename)

def main():
    from basic_move import Box
    from enable.api import Container
    container = Container()
    box = Box(bounds=[30,30], position=[20,20], padding=5)
    container.add(box)

    # Create the mlab test mesh and get references to various parts of the
    # VTK pipeline
    m = mlab.test_mesh()
    scene = mlab.gcf().scene
    render_window = scene.render_window
    renderer = scene.renderer
    rwi = scene.interactor

    # Create the Enable Window
    window = EnableVTKWindow(rwi, renderer,
            component=container,
            #istyle_class = tvtk.InteractorStyleSwitch,
            istyle_class = tvtk.InteractorStyle,
            resizable = "v",
            bounds = [100, 100],
            padding_top = 20,
            padding_bottom = 20,
            padding_left = 20,
            )

    #rwi.render()
    #rwi.start()
    mlab.show()
    return window, render_window

 def plotvtk3D(self):
     """
     3D plot using the vtk libraries
     """
     from tvtk.api import tvtk
     from mayavi import mlab
     # Plot the data on a 3D grid
     xy = np.column_stack((self.grd.xp,self.grd.yp))
     dvp = self.F(xy)
     vertexag = 50.0
     
     points = np.column_stack((self.grd.xp,self.grd.yp,dvp*vertexag))
     tri_type = tvtk.Triangle().cell_type
     #tet_type = tvtk.Tetra().cell_type
     ug = tvtk.UnstructuredGrid(points=points)
     ug.set_cells(tri_type, self.grd.cells)
     
     ug.cell_data.scalars = self.grd.dv
     ug.cell_data.scalars.name = 'depths'
     
     f=mlab.gcf()
     f.scene.background = (0.,0.,0.)
     d = mlab.pipeline.add_dataset(ug)
     h=mlab.pipeline.surface(d,colormap='gist_earth')
     mlab.colorbar(object=h,orientation='vertical')
     mlab.view(0,0)
     
     outfile = self.suntanspath+'/depths.png'
     f.scene.save(outfile)      
     print 'Figure saved to %s.'%outfile

def plot_targets(tasks, color=(1., 0., 0.), scale_factor=0.005):
    '''Plot start and end points of tasks.

    Starting points are designated by spheres, and end points by cubes.
    '''
    # draw target positions
    p3d = mlab.pipeline
    fig = mlab.gcf()
    fig.scene.disable_render = True
    
    starts, stops = tasks[:,:3], tasks[:,3:]

    # plot start points
    x,y,z = starts.T
    start_src = p3d.scalar_scatter(x, y, z)
    gl1 = p3d.glyph(start_src)
    gl1.glyph.glyph.scale_factor = scale_factor         # smaller
    gl1.actor.actor.property.color = color # red
    gl1.actor.actor.property.opacity = 0.5 # red

    # plot stop points
    #x,y,z = stops.T
    #stop_src = p3d.scalar_scatter(x, y, z)
    #gl2 = p3d.glyph(stop_src)
    #gl2.glyph.glyph.scale_factor = scale_factor         # smaller
    #gl2.actor.actor.property.color = color # red
    #cube = gl2.glyph.glyph_source.glyph_dict['cube_source']
    #gl2.glyph.glyph_source.glyph_source = cube
    #gl2.actor.actor.property.opacity = 0.5

    fig.scene.disable_render = False
    return gl1

def plot_surface(surf_name, data_orig, hemi='rh', cmin=None, cmax=None, colorbar=True, smooth=5):
    ''' Plots data in a 3D brain surface using the current Mayavi's mlab window. surf_name is a string with the name of the surface, data is a vector of the same length of the number of points in the surface. Function performs smoothing by default, and set the color of the brain to gray by default for points we don't have data. '''

    surf = loadmat('/Users/sudregp/Documents/surfaces/IMAGING_TOOLS/' + surf_name + '.mat')

    nbr = surf['nbr_' + hemi].astype(int)

    # making sure we don't change the original data
    data = data_orig.copy()
    num_voxels = len(data)
    # smoothing data to neighboring voxels with zero value. Algorithm described in MNE-manual-2.7, page 208/356.
    for p in range(smooth):
        print 'Smoothing step ' + str(p + 1) + '/' + str(smooth)
        S = np.zeros([num_voxels, num_voxels])
        for j in range(num_voxels):
            my_neighbors = nbr[j, :] - 1
            # remove entries that are -1
            my_neighbors = np.delete(my_neighbors, np.nonzero(my_neighbors == -1))
            # number of immediate neighbors with non zero values
            Nj = np.sum(data[my_neighbors] != 0)
            if Nj > 0:
                pdb.set_trace()
                S[j, my_neighbors] = 1 / float(Nj)
        data = np.dot(S, data)

    pdb.set_trace()
    # replacing all values that are still 0 by something that is not in the data range.
    data[data == 0] = np.inf

    if cmin is None:
        cmin = np.min(data)
    if cmax is None:
        # check whether we have empty points in the brain. If we do, it has Inf value. If not, the max is the actual maximum of the data
        if len(np.nonzero(np.isinf(data) == True)[0]) == 0:
            cmax = np.max(data)
            add_grey = False
        else:
            cmax = np.unique(np.sort(data))[-2]
            add_grey = True
    else:
        # if cmax was specified, let the person deal with it
        add_grey = False

    surf = mlab.triangular_mesh(surf['coord_' + hemi][0, :], surf['coord_' + hemi][1, :], surf['coord_' + hemi][2, :], surf['tri_' + hemi] - 1, scalars=data, vmax=cmax, vmin=cmin)

    if add_grey:
        # add grey color to scale, and make everything that's infinity that color
        surf.module_manager.scalar_lut_manager.number_of_colors += 1
        lut = surf.module_manager.scalar_lut_manager.lut.table.to_array()
        grey_row = [192, 192, 192, 255]
        # sets the max value in the data range to be gray
        lut = np.vstack((lut, grey_row))
        surf.module_manager.scalar_lut_manager.lut.table = lut

    fig = mlab.gcf()
    # fig.on_mouse_pick(picker_callback)
    mlab.colorbar()

    return surf

def load_nvm_and_register_images(nvm_file, 
                                 images_path, 
                                 max_sleep_seconds,
                                 every_nth,
                                 single_image = None):

    cams, points, cams_mlab, points_mlab = mayaviu.load_and_plot_nvm_cams(nvm_file)    
    fig = mlab.gcf()
    # signal.signal(signal.SIGUSR1, plot_localized_image_after_signal(fig))

    # if single_image is not None:
    #     tmp_dir = 'tmp_im_dir'
    #     if not os.path.isdir(tmp_dir):
    #         os.mkdir(tmp_dir)
    #     shutil.copy(
    
    # signal.signal(signal.SIGUSR1, localized_image_signal)

    model = vsfm_model.create_from_nvm(nvm_file)
    vsfm_interface = vsfmu.vsfm_interface()
    vsfm_interface.sfm_more_less_visualization_data()
    
    image_files = ['{}/{}'.format(images_path, i) \
                      for i in os.listdir(images_path) if i.find('.jpg') > 0]
    image_basenames = sorted([os.path.basename(_) for _ in image_files])
    image_files = sorted([os.path.abspath(_) for _ in image_files])

    vsfm_interface.sfm_load_nview_match(nvm_file)
    vsfm_interface.view_image_thumbnails()

    mlab.show()
    rtst = rts_thread(vsfm_interface, 
                      nvm_file,
                      model, 
                      model_image_dir = '/home/nrhineha/dev/activity/data/kitchen_seq_x1',
                      image_files = image_files,
                      max_sleep_seconds = max_sleep_seconds)
    rtst.start()
    mlab.show()
    
    # seq = register_temporal_sequence(vsfm_interface, 
    #                                  nvm_file,
    #                                  model, 
    #                                  model_image_dir = '/home/nrhineha/dev/activity/data/kitchen_seq_x1',
    #                                  image_files = image_files,
    #                                  max_sleep_seconds = max_sleep_seconds)
    
    # cams, new_cams = register_image(vsfm_interface, image_files[0])

    # # vsfm_interface.sfm_delete_selected_camera()

    # near_cams= model.lookup_nearby_cameras(new_cams[0])
    # mfn = write_specified_match_file('/home/nrhineha/Desktop/test/loc2.jpg',
    #                                  near_cams, 
    #                                  '/home/nrhineha/dev/activity/data/kitchen_seq_x1')

    # register_image(vsfm_interface, images_path, image_files[1], image_basenames[1], match_specified_fn = os.path.abspath(mfn))
    
    ipdb.set_trace()

def plot_lines(lines, scalars=None, style="tube", figure=None,
               name="Lines", tube_radius=0.05, tube_sides=6, **kwargs):
    """Make 3D mayavi plot of lines

    Scalars can be a bunch of single values, or a bunch of rgb data
    to set the color of each line / vertex explicitly. This is
    explained in :py:func:`lines2source`.

    Example:
        A common use case of setting the line color from a topology
        will want to use :py:func:`viscid.topology2color`::

            >>> import viscid
            >>> from viscid.plot import mvi
            >>>
            >>> B = viscid.vlab.get_dipole()
            >>> seeds = viscid.Line([-4, 0, 0], [4, 0, 0])
            >>> lines, topology = viscid.calc_streamlines(B, seeds,
            >>>                                           ibound=0.05)
            >>> scalars = viscid.topology2color(topology)
            >>> mvi.plot_lines(lines, scalars, tube_radius=0.02)
            >>> mvi.mlab.savefig("dipole.x3d")
            >>> viscid.vutil.meshlab_convert("dipole.x3d", "dae")
            >>> mvi.show()

    Parameters:
        lines (list): See :py:func:`lines2source`
        scalar (ndarray): See :py:func:`lines2source`
        style (str): 'strip' or 'tube'
        figure (mayavi.core.scene.Scene): specific figure, or
            :py:func:`mayavi.mlab.gcf`
        name (str): Description
        tube_radius (float): Radius if style == 'tube'
        tube_sides (int): Angular resolution if style == 'tube'
        **kwargs: passed to :meth:`mayavi.mlab.pipeline.surface`. This
            is useful for setting a colormap among other things.

    Returns:
        Mayavi surface module

    Raises:
        ValueError: if style is neither tube nor strip
    """
    style = style.lower()
    if not figure:
        figure = mlab.gcf()

    src = lines2source(lines, scalars=scalars, name=name)

    if style == "tube":
        lines = mlab.pipeline.tube(src, figure=figure, tube_radius=tube_radius,
                                   tube_sides=tube_sides)
    elif style == "strip":
        lines = mlab.pipeline.stripper(src, figure=figure)
    else:
        raise ValueError("Unknown style for lines: {0}".format(style))

    surface = mlab.pipeline.surface(lines, **kwargs)
    return surface

def disp_pts(pts1, pts2, color1=(1,0,0), color2=(0,1,0)):
    figure = mlab.gcf()
    mlab.clf()
    figure.scene.disable_render = True

    pts1_glyphs   = mlab.points3d(pts1[:,0], pts1[:,1], pts1[:,2], color=color1, resolution=20, scale_factor=0.001)
    pts2_glyphs   = mlab.points3d(pts2[:,0], pts2[:,1], pts2[:,2], color=color2, resolution=20, scale_factor=0.001)
    glyph_points1 = pts1_glyphs.glyph.glyph_source.glyph_source.output.points.to_array()
    glyph_points2 = pts2_glyphs.glyph.glyph_source.glyph_source.output.points.to_array()

    dd = 0.001

    outline1 = mlab.outline(pts1_glyphs, line_width=3)
    outline1.outline_mode = 'full'
    p1x, p1y, p1z = pts1[0,:]
    outline1.bounds = (p1x-dd, p1x+dd,
                       p1y-dd, p1y+dd,
                       p1z-dd, p1z+dd)

    pt_id1 = mlab.text(0.8, 0.2, '0 .', width=0.1, color=color1)

    outline2 = mlab.outline(pts2_glyphs, line_width=3)
    outline2.outline_mode = 'full'
    p2x, p2y, p2z = pts2[0,:]
    outline2.bounds = (p2x-dd, p2x+dd,
                       p2y-dd, p2y+dd,
                       p2z-dd, p2z+dd)  
    pt_id2 = mlab.text(0.8, 0.01, '0 .', width=0.1, color=color2)
    
    figure.scene.disable_render = False


    def picker_callback(picker):
        """ Picker callback: this gets called during pick events.
        """
        if picker.actor in pts1_glyphs.actor.actors:
            point_id = picker.point_id/glyph_points1.shape[0]
            if point_id != -1:
                ### show the point id
                pt_id1.text = '%d .'%point_id
                #mlab.title('%d'%point_id)
                x, y, z = pts1[point_id,:]
                outline1.bounds = (x-dd, x+dd,
                                   y-dd, y+dd,
                                   z-dd, z+dd)
        elif picker.actor in pts2_glyphs.actor.actors:
            point_id = picker.point_id/glyph_points2.shape[0]
            if point_id != -1:
                ### show the point id
                pt_id2.text = '%d .'%point_id
                x, y, z = pts2[point_id,:]
                outline2.bounds = (x-dd, x+dd,
                                   y-dd, y+dd,
                                   z-dd, z+dd)


    picker = figure.on_mouse_pick(picker_callback)
    picker.tolerance = dd/2.
    mlab.show()

 def anim():
     f = mlab.gcf()
     for count, i in enumerate(range(2,361,2)):
     #while 1:
         f.scene.camera.azimuth(degree_step)
         f.scene.render()
         mlab.savefig('/Users/jessebrown/Desktop/mayavi_figs/r_hipp_network%03d.png' %count)
         yield

def lorenz_perturbed(N=10, res=10000, step=5, t=50, seed_=120, atol=1E-15,
                     rtol=1E-13, epsilon=2.2e-16, delay=10,
                     sigma=10., beta=8./3, rho=28.):
    """ Animate the trajectories given by the Lorenz equations.
    Plot two trajectories, one with the initial value given by the
    random number generator after you seed it,
    and another that is equal to (1 + epsilon) times the other initial value.
    Choose random x, y, and z values between -15 and 15.
    Seed the random number generator with 'seed_'.
    Use a resolution of 'res' for the points in the plot.
    Plot the time values between 0 ant 't'.
    Pass the tolerances 'atol' and 'rtol' to the ODE solver.
    At each update, add 'step' points to the plot.
    Use a delay of 'delay' at each update in the animation.
    Use different colors for each trajectory.
    Use the values of 'sigma', 'beta', and 'rho' in the Lorenz ODE. """
    # Get initial conditions.
    seed(seed_)
    x1 = -15 + 30 * rand(3)
    x2 = x1 * (1. + epsilon)
    
    # Solve for the trajectories.
    # Plot them.
    t = np.linspace(0, t, res)
    y1 = odeint(lorenz_ode, x1, t, args=(sigma, beta, rho), atol=atol, rtol=rtol)
    c1 = mlab.plot3d(y1[:1,0], y1[:1,1], y1[:1,2], tube_radius=.2, color=(1, 0, 0))
    y2 = odeint(lorenz_ode, x2, t, args=(sigma, beta, rho), atol=atol, rtol=rtol)
    c2 = mlab.plot3d(y2[:1,0], y2[:1,1], y2[:1,2], tube_radius=.2, color=(0, 0, 1))
    
    # Position the view for the plot.
    mlab.gcf().scene.camera.position = [127.23761585, -108.28736806, 6.35191272]
    mlab.gcf().scene.camera.focal_point = [-1.7792501449584961, -3.6287221908569336, 23.397351264953613]
    mlab.gcf().scene.camera.view_up = [-0.078467260964232038, -0.20339450183237351, 0.97594752194015633]
    mlab.gcf().scene.camera.clipping_range = [128.64624663718814, 328.22549479639167]
    
    # Define the animation.
    @mlab.show
    @mlab.animate(delay=delay)
    def trace_curve():
        for i in xrange(2, res, step):
            c1.mlab_source.reset(x=y1[:i,0], y=y1[:i,1], z=y1[:i,2])
            c2.mlab_source.reset(x=y2[:i,0], y=y2[:i,1], z=y2[:i,2])
            yield
    
    # Run the animation.
    trace_curve()

    def plot_u(u, x, xv, y, yv, t, n):
        """User action function for plotting."""
        if t[n] == 0:
            time.sleep(2)
        if plot_method == 1:
            # Works well with Gnuplot backend, not with Matplotlib
            st.mesh(x, y, u, title='t=%g' % t[n], zlim=[-1,1],
                    caxis=[-1,1])
        elif plot_method == 2:
            # Works well with Gnuplot backend, not with Matplotlib
            st.surfc(xv, yv, u, title='t=%g' % t[n], zlim=[-1, 1],
                  colorbar=True, colormap=st.hot(), caxis=[-1,1],
                  shading='flat')
        elif plot_method == 3:
            print 'Experimental 3D matplotlib...under development...'
            # Probably too slow
            #plt.clf()
            ax = fig.add_subplot(111, projection='3d')
            u_surf = ax.plot_surface(xv, yv, u, alpha=0.3)
            #ax.contourf(xv, yv, u, zdir='z', offset=-100, cmap=cm.coolwarm)
            #ax.set_zlim(-1, 1)
            # Remove old surface before drawing
            if u_surf is not None:
                ax.collections.remove(u_surf)
            plt.draw()
            time.sleep(1)
        elif plot_method == 4:
	    # Mayavi visualization
            mlab.clf()
            extent1 = (0, 20, 0, 20,-2, 2)
            s = mlab.surf(x , y, u,
                          colormap='Blues',
                          warp_scale=5,extent=extent1)
            mlab.axes(s, color=(.7, .7, .7), extent=extent1,
                      ranges=(0, 10, 0, 10, -1, 1),
                      xlabel='', ylabel='', zlabel='',
                      x_axis_visibility=False,
                      z_axis_visibility=False)
            mlab.outline(s, color=(0.7, .7, .7), extent=extent1)
            mlab.text(6, -2.5, '', z=-4, width=0.14)
            mlab.colorbar(object=None, title=None,
                          orientation='horizontal',
                          nb_labels=None, nb_colors=None,
                          label_fmt=None)
            mlab.title('Gaussian t=%g' % t[n])
            mlab.view(142, -72, 50)
            f = mlab.gcf()
            camera = f.scene.camera
            camera.yaw(0)

        if plot_method > 0:
            time.sleep(0) # pause between frames
            if save_plot:
                filename = 'tmp_%04d.png' % n
		if plot_method == 4:
                    mlab.savefig(filename)  # time consuming!
		elif plot_method in (1,2):
                    st.savefig(filename)  # time consuming!

def close():
    """Close the scene."""
    f = mlab.gcf()
    e = mlab.get_engine()
    e.window.workbench.prompt_on_exit = False
    e.window.close()
    mlab.options.backend = 'auto'
    # Hack: on Linux the splash screen does not go away so we force it.
    GUI.invoke_after(500, e.window.workbench.application.gui.stop_event_loop)